;;
;; Basic test of normal-pdf.
;;
;; This test creates a one-dimensional normal distribution with a
;; given mean and variance. It then samples from this distribution
;; and calculates the sample mean and variance for comparison.
;;

(use srfi-4 normal-pdf random-mtzig test)

;; Number of samples 
(define N 100000)
;(define N 10)

;;
;;  Computes the arithmetic mean of a f64vector using the recurrence relation 
;;     mean_(n) = mean(n-1) + (v[n] - mean(n-1))/(n+1)
;;
(define (mean v)
  (let ((n (f64vector-length v)))
    (let loop ((i 0) (mean 0))
      (if (< i n)
	  (loop (fx+ 1 i) (+ mean (/ (- (f64vector-ref v i) mean) (+ 1 i))))
	  mean))))

(define (variance v mean)
  (let ((n (f64vector-length v)))
    (let loop ((i 0) (sqsum 0))
      (if (< i n)
	  (let ((delta (- (f64vector-ref v i) mean)))
	    (loop (fx+ 1 i) (+ sqsum (* delta delta))))
	  (* (/ 1 (- n 1)) sqsum)))))


;;  distribution parameters
(define mu     1.58798)
(define sigma  4.55175)

(define gpdf (make-normal-pdf 1 mu sigma)) 

(define rng      (random-mtzig:init))
(define input    (random-mtzig:f64vector-randn! N rng))
(define samples  (make-f64vector N 0))
(define density  (make-f64vector N 0)) 
      
(test-group "univariate normal PDF test"

   (let loop ((i 0))
     (if (< i N)
	 (begin (f64vector-set! samples i (normal-pdf:sample gpdf (f64vector-ref input i)))
		(f64vector-set! density i (normal-pdf:density gpdf (f64vector-ref input i)))
		(loop (fx+ 1 i)))))

   ;; test for mean and variance equality to one decimal place
   (test "sample-mean"
	 (truncate (* 10 mu))
	 (truncate (* 10 (mean samples))))

   (test "sample-variance"
	 (truncate (* 10 sigma))
	 (truncate (* 10 (variance samples (mean samples)))))

)