;; some tests for statistics package

(require-extension srfi-42)
(require-extension srfi-78)
(require-extension statistics)

;; -- to help the tests

(define (to-4-dp f)
  (/ (round (* f 10000)) 10000))

(define (to-5-dp f)
  (/ (round (* f 100000)) 100000))

(define (=4 n1 n2)
  (= (to-4-dp n1) (to-4-dp n2)))

(define (=5 n1 n2)
  (= (to-5-dp n1) (to-5-dp n2)))

;; -- tests of utilities

(check (average-rank 3 '(1 2 3 4 5))
       => 3)

(check (average-rank 3 '(1 2 3 3 4 5))
       => 7/2)

(check (beta-incomplete 0.2 0.1 0.3)
       (=> =5) 0.28367)

(check (bin-and-count '(1 2 3 4 5 6 7 8 9 10) 2)
       => (vector 5 5))

(check (bin-and-count '(1 2 3 4 5 6 7 8 9 10) 5)
       => (vector 2 2 2 2 2))

(check (bin-and-count '(1 2 2 2 3 9) 2)
       => (vector 5 1))

(check (combinations 10 3) 
       => 120)

(check (factorial 10)
       => 3628800)

; find-critical-value

(check (fisher-z-transform 0.1)
       (=> =5) 0.10034)

(check (fisher-z-transform 0.5)
       (=> =5) 0.54931)

(check (let-values (((a b) (gamma-incomplete 2.0 1.5))) (list (to-5-dp a) b))
       => (list 0.44217 0.0))

(check (gamma-ln 0.5)
       (=> =5) 0.57236)

(check (permutations 10 3)
       => 720)

; random-normal: unsure how to test

(check (random-pick '(1 2 3 4 5))
       (=> member) '(1 2 3 4 5))

(check (random-sample 3 '(1 2 3 4 5))
       (=> (lambda (result expected) (= (length result) expected))) 3)

(check (sign 3.5)
       => 1)

(check (sign -3.5)
       => -1)

(check (sign 0)
       => 0)

(check (square 10)
       => 100)

;; -- tests of descriptive statistics

(check (mean '())
       => 0)

(check (mean '(1))
       => 1)

(check (mean '(1 2 3 4 5))
       => 3)

(check (median '(1))
       => 1)

(check (median '(1 2 3 4 5))
       => 3)

(check (median '(1 1 2 3 4 5 6 7 8))
       => 4)

(check (let-values (((modes counts) (mode '(1 1 1 2 3 4 5)))) 
                   (list modes counts))
       => '((1) 3))

(check (let-values (((modes counts) (mode '(1 1 1 2 3 4 4 4 5)))) 
                   (list modes counts))
       => '((1 4) 3))

(check (geometric-mean '(1 100))
       => 10.0)

(check (geometric-mean '(1 2 3 4 5))
       (=> =5) 2.60517)

(check (range '(1 2 3 4 5))
       => 4)

(check (range '(1 1 1 2 3 4 10))
       => 9)

(check (percentile '(1 2 3 4 5 6) 50)
       => 7/2)

(check (percentile '(1 2 3 4 5 6) 30)
       => 2)

(check (variance '(1 2 3 4 5))
       => 5/2)

(check (standard-deviation '(1 2 3 4 5))
       (=> =5) 1.58113883)

(check (coefficient-of-variation '(1 2 3 4 5))
       (=> =5) 52.704627)

(check (standard-error-of-the-mean '(1 2 3 4 5))
       (=> =5) 0.707106781186548)

(let-values (((mean sd n) (mean-sd-n '(1 2 3 4 5))))
            (check mean => 3)
            (check sd (=> =5) 1.58113883008419)
            (check n => 5))

;; -- tests of distributional functions

(check-ec (: pair '((0 0.0009765625) (1 0.009765625) (2 0.0439453125) (3 0.1171875)
                                     (4 0.205078125) (5 0.24609375) (6 0.205078125)
                                     (7 0.1171875) (8 0.0439453125) (9 0.009765625)
                                     (10 0.0009765625)))
 (binomial-probability 10 (car pair) 0.5)
 (=> =5) (cadr pair))

(check-ec (: pair '((0 0.0) (1 0.0009765625) (2 0.0107421875) (3 0.0546875)
                            (4 0.171875) (5 0.376953125) (6 0.623046875)
                            (7 0.828125) (8 0.9453125) (9 0.9892578125)
                            (10 0.9990234375)))
          (binomial-cumulative-probability 10 (car pair) 0.5)
          (=> =5) (cadr pair))

(check-ec (: pair '((0 0.0000) (1 0.0005) (2 0.0023) (3 0.0076) (4 0.0189) (5 0.0378)
                               (6 0.0631) (7 0.0901) (8 0.1126) (9 0.1251) (10 0.1251)
                               (11 0.1137) (12 0.0948) (13 0.0729) (14 0.0521)
                               (15 0.0347) (16 0.0217) (17 0.0128) (18 0.0071)
                               (19 0.0037)))
          (poisson-probability 10 (car pair))
          (=> =4) (cadr pair))

(check-ec (: pair '((0 0.0000) (1 0.0000) (2 0.0005) (3 0.0028) (4 0.0103) (5 0.0293)
             (6 0.0671) (7 0.1301) (8 0.2202) (9 0.3328) (10 0.4579) 
             (11 0.5830) (12 0.6968) (13 0.7916) (14 0.8645)
             (15 0.9165) (16 0.9513) (17 0.9730) (18 0.9857) (19 0.9928)))
          (poisson-cumulative-probability 10 (car pair))
          (=> =4) (cadr pair))

(check-ec (: pair '((0 0.0088) (1 0.0270) (2 0.0648) (3 0.1210) (4 0.1760)
                               (5 0.1995) (6 0.1760) (7 0.1210) (8 0.0648)
                               (9 0.0270)(10 0.0088)))
          (normal-pdf (car pair) 5 4)
          (=> =4) (cadr pair))

(check (convert-to-standard-normal 5 6 2)
       => -1/2)

(check-ec (: pair '((-2.0 0.0228) (-1.6 0.0548) (-1.2 0.1151) (-0.8 0.2119)
                                  (-0.4 0.3446) ( 0.0 0.5000) ( 0.4 0.6554)
                                  ( 0.8 0.7881) ( 1.2 0.8849) ( 1.6 0.9452)))
          (phi (car pair))
          (=> =4) (cadr pair))

(check-ec (: pair '((0.1 -1.2815515713806909) (0.2 -0.8416212245351449)))
          (z (car pair))
          (=> =5) (cadr pair))

(check (t-distribution 1 0.9)
       (=> =5) 3.07763671875)

(check (t-distribution 10 0.9)
       (=> =5) 1.3721923828125)

(check (chi-square 10 0.4405)
       (=> =5) 10)

(check (chi-square 3 0.1718)
       (=> =5) 5)

(check (chi-square-cdf 10 10)
       (=> =5) 0.559506714934786)

(check (chi-square-cdf 5 3)
       (=> =5) 0.828202855703266)

;; -- tests of confidence intervals

(let-values (((upper-bound lower-bound) (binomial-probability-ci 10 0.8 0.9)))
            (check upper-bound (=> =5) 0.724273681640625)
            (check lower-bound (=> =5) 0.851547241210938))

(let-values (((upper-bound lower-bound) (poisson-mu-ci 10 0.9)))
            (check upper-bound (=> =5) 8.305419921875)
            (check lower-bound (=> =5) 10.0635986328125))

(let-values (((upper-bound lower-bound) (normal-mean-ci 0.5 0.1 10 0.8)))
            (check upper-bound (=> =5) 0.491747852700165)
            (check lower-bound (=> =5) 0.508252147299835))

(let-values (((upper-bound lower-bound) (normal-mean-ci-on-sequence '(1 2 3 4 5) 0.9)))
            (check upper-bound (=> =5) 2.90535368828478)
            (check lower-bound (=> =5) 3.09464631171522))

;; -- tests of hypothesis testing

(check (z-test 40 1 #:mu 50 #:sigma 10 #:tails ':negative)
       (=> =5) 0.15865525)
(check (z-test 40 10 #:mu 50 #:sigma 10 #:tails ':negative)
       (=> =5) 0.000783)
(check (z-test 40 5 #:mu 50 #:sigma 10)
       (=> =5) 0.025347)

(check (z-test 11/5 5 #:mu 2 #:sigma 1)
       (=> =5) 0.65472085)

(check (z-test-on-sequence '(1 1 2 3 4) #:mu 2 #:sigma 1)
       (=> =5) 0.65472085)

(check (t-test-one-sample 5 0.8 3 6)
       (=> =5) 0.162781641721079)

(check (t-test-one-sample-on-sequence '(5 6 7) 5.9)
       (=> =5) 0.878433865229034)

(check (correlation-test-two-sample 0.9 100 0.85 50)
       (=> =5) 0.224083300908794)

(check (t-test-paired 7/3 0.57735 3)
       (=> =5) 0.0198039411803931)

(check (t-test-paired-on-sequences '(4 3 5) '(1 1 3))
       (=> =5) 0.0198039411803931)

; correlation-test-two-sample-on-sequences

(let-values (((s p) (spearman-rank-correlation '((4 5) (10 8) (3 6) (1 2) (9 10) (2 3) (6 9) (7 4) (8 7) (5 1)))))
            (check s (=> =5) 0.684848484848485)
            (check p (=> =5) 0.0288827975067328))

;; -- tests of sample size estimates

(check (t-test-one-sample-sse 5.0 5.2 0.5)
       => 163)

(let-values (((n1 n2) (t-test-two-sample-sse 5.1 0.5 5.2 0.3)))
            (check n1 => 1040)
            (check n2 => 1040))

(check (correlation-sse 0.80 #:alpha 0.05 #:1-beta 0.9)
       => 11) ;; NOTE: Not same as some on-line sites describe ...

;; -- tests of correlation and regression

(let-values (((a b r r2 p) (linear-regression '((1.0 0.1) (2.0 0.3) (3.0 0.8)))))
            (check a (=> =5) -0.3)
            (check b (=> =5) 0.35)
            (check r (=> =5) 0.970725343394151)
            (check r2 (=> =5) 0.942307692307692)
            (check p (=> =5) 0.154420958311267))

(check (correlation-coefficient '((1.0 0.1) (2.0 0.3) (3.0 0.8)))
       (=> =5) 0.970725343394151)

;; -- tests of significance tests

(check (t-significance 0.2 5)
       (=> =5) 0.849360513995829)

(check (t-significance 0.2 5 #:tails ':both)
       (=> =5) 0.849360513995829)

(check (t-significance 0.2 5 #:tails ':positive)
       (=> =5) 0.424680256997915)

(check (t-significance 0.2 5 #:tails ':negative)
       (=> =5) 0.575319743002086)

(check (f-significance 1.5 8 2)
       (=> =5) 0.920449812578091)

(check (f-significance 1.5 8 2 #:one-tailed? #t)
       (=> =5) 0.460224906289046)

(check (binomial-test-two-sample 0.2 100 0.3 50)
       (=> =5) 0.245930683028145)

(check (fisher-exact-test 10 20 30 40)
       (=> =5) 0.5066621427235114)

;; -- summarise results of tests

(check-report)

;; uncomment for released version - replace N with number of tests
;(if (check-passed? N)
;  'ok
;  (error "failed test suite"))