(load "utils.ss")
(load "math.ss")

(define mean 
  (lambda (l)
    (/ (apply + l) (length l))))

;; computes an estimate of the variance of a list of numbers.
(define varest
  (lambda (l)
    (let* ((m (mean l))
           (func (lambda (x) 
                   (let ((norm (- x m)))
                     (* norm norm)))))
      (mean (map func l)))))


(define erf ;;error function
  (lambda (z)
    (let*((t (/ 1.0 (+ 1.0 (* 0.5 (abs z)))))
          (coeflist '(-1.26551223 1.00002368 0.37409196 0.09678418 -0.18628806 0.27886807
                                  -1.13520398 1.48851587 -0.82215223 0.17087277))
          (mantissa (+ (- (* z z))
                       (horneval2 coeflist t)))
          (ans (- 1.0 (* t (exp mantissa)))))
      (if (>= z 0) ans (- ans)))))
      
(define Phi ;; cumulative normal distribution
  (lambda (z)
    (* 0.5 (+ 1.0 (erf (/ z (sqrt 2.0)))))))


(define bicum
  (lambda (a b rho)
    (cond ((and (<= a 0) (<= b 0) (<= rho 0))
           (let ((alist '(0.3253030 0.4211071 0.1334425 0.006374323))
                 (blist '(0.1337764 0.6243247 1.3425378 2.2626645))
                 (auxfunc (lambda (x y) (faux x y a b rho)))
                 (themul (/ (sqrt (- 1 (* rho rho))) pi)))
             (let ((aout (flatten (outer * alist alist)))
                   (bout (flatten (outer auxfunc blist blist))))
               (* themul (ldot aout bout)))))
          ((and (<= a 0) (>= b 0) (>= rho 0))
           (- (Phi a) (bicum a (- b) (- rho))))
          ((and (>= a 0) (<= b 0) (>= rho 0))
           (- (Phi b) (bicum (- a) b (- rho))))
          ((and (>= a 0) (>= b 0) (<= rho 0))
           (+ (Phi a) (Phi b) -1 (bicum (- a) (- b) rho)))
          (else 
           (let ((tmp (/ 1 (sqrt (+ (* a a) (* -2 rho a b) (* b b))))))
             (let ((rho1 (* (sgn a) (- (* rho a) b) tmp))
                   (rho2 (* (sgn b) (- (* rho b) a) tmp))
                   (d (/ (- 1 (* (sgn a) (sgn b)) 4))))
               (+ (bicum a 0 rho1) (bicum b 0 rho2) (- d))))))))