source: project/release/5/statistics/trunk/statistics.scm @ 36303

Last change on this file since 36303 was 36303, checked in by iraikov, 11 months ago

statistics port to C5

File size: 64.0 KB
Line 
1;;; Port of cl-stats.lisp to Chicken Scheme, by Peter Lane, 2010.
2
3;;; This program is free software: you can redistribute it and/or modify
4;;; it under the terms of the GNU General Public License as published by
5;;; the Free Software Foundation, either version 3 of the License, or
6;;; (at your option) any later version.
7
8;;; This program is distributed in the hope that it will be useful,
9;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
10;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
11;;; GNU General Public License for more details.
12
13;;; You should have received a copy of the GNU General Public License
14;;; along with this program.  If not, see <http://www.gnu.org/licenses/>.
15
16;;; ------- Source of methods -------
17;;;
18;;; The formulas and methods used are largely taken from Bernard Rosner,
19;;; "Fundamentals of Biostatistics," 5th edition.  "Rosner x" is a page
20;;; number.
21
22;;;  These abreviations used in function and variable names:
23;;;    ci = confidence interval
24;;;    cdf = cumulative density function
25;;;    ge = greater than or equal to
26;;;    le = less than or equal to
27;;;    pdf = probability density function
28;;;    sd = standard deviation
29;;;    rxc = rows by columns
30;;;    sse = sample size estimate
31
32;;; ------- Original copyright notices from cl-stats.lisp -------
33
34;;; Statistical functions in Common Lisp.  Version 1.04 Feb 24, 2005
35;;;
36;;; This code is copyright (c) 2000, 2001, 2002, 2005 by Larry Hunter
37;;; (Larry.Hunter@uchsc.edu) except where otherwise noted.  (Released
38;;; under GPL version 2.)
39
40;;; CLASP utilities marked as such in comment:
41
42;;; Copyright (c) 1990 - 1994 University of Massachusetts
43;;; Department of Computer Science
44;;; Experimental Knowledge Systems Laboratory
45;;; Professor Paul Cohen, Director.
46;;; All rights reserved.
47
48;;; Permission to use, copy, modify and distribute this software and its
49;;; documentation is hereby granted without fee, provided that the above
50;;; copyright notice of EKSL, this paragraph and the one following appear
51;;; in all copies and in supporting documentation.
52;;; EKSL makes no representation about the suitability of this software for any
53;;; purposes.  It is provided "AS IS", without express or implied warranties
54;;; including (but not limited to) all implied warranties of merchantability
55;;; and fitness for a particular purpose, and notwithstanding any other
56;;; provision contained herein.  In no event shall EKSL be liable for any
57;;; special, indirect or consequential damages whatsoever resulting from loss
58;;; of use, data or profits, whether in an action of contract, negligence or
59;;; other tortuous action, arising out of or in connection with the use or
60;;; performance of this software, even if EKSL is advised of the possiblity of
61;;; such damages.
62
63;;; For more information write to clasp-support@cs.umass.edu
64
65;;; -------------------------------------------------------------------------------
66
67(module
68  statistics
69  (
70    ; utilities
71    average-rank
72    beta-incomplete
73    bin-and-count
74    combinations
75    factorial
76    find-critical-value
77    fisher-z-transform
78    gamma-incomplete
79    gamma-ln
80    permutations
81    random-normal
82    random-pick
83    random-sample
84    random-weighted-sample
85    sign
86    square
87    cumsum
88   
89    ; descriptive statistics
90    mean
91    median
92    mode
93    geometric-mean
94    range
95    percentile
96    variance
97    standard-deviation
98    coefficient-of-variation
99    standard-error-of-the-mean
100    mean-sd-n
101
102    ; distributional functions
103    binomial-probability
104    binomial-cumulative-probability
105    binomial-ge-probability
106    binomial-le-probability
107    poisson-probability
108    poisson-cumulative-probability
109    poisson-ge-probability
110    normal-pdf
111    convert-to-standard-normal
112    phi
113    z
114    t-distribution
115    chi-square
116    chi-square-cdf
117
118    ; Confidence intervals
119    binomial-probability-ci
120    poisson-mu-ci
121    normal-mean-ci
122    normal-mean-ci-on-sequence
123    normal-variance-ci
124    normal-variance-ci-on-sequence
125    normal-sd-ci
126    normal-sd-ci-on-sequence
127
128    ; Hypothesis testing
129    ; (parametric)
130    z-test
131    z-test-on-sequence
132    t-test-one-sample
133    t-test-one-sample-on-sequence
134    t-test-paired
135    t-test-paired-on-sequences
136    t-test-two-sample
137    t-test-two-sample-on-sequences
138    f-test
139    chi-square-test-one-sample
140    binomial-test-one-sample
141    binomial-test-two-sample
142    fisher-exact-test
143    mcnemars-test
144    poisson-test-one-sample
145    ; (non parametric)
146    sign-test
147    sign-test-on-sequence
148    wilcoxon-signed-rank-test
149    wilcoxon-signed-rank-test-on-sequences
150    chi-square-test-rxc
151    chi-square-test-for-trend
152
153    ; Sample size estimates
154    t-test-one-sample-sse
155    t-test-two-sample-sse
156    t-test-paired-sse
157    binomial-test-one-sample-sse
158    binomial-test-two-sample-sse
159    binomial-test-paired-sse
160    correlation-sse
161
162    ; Correlation and regression
163    linear-regression
164    correlation-coefficient
165    correlation-test-two-sample
166    correlation-test-two-sample-on-sequences
167    spearman-rank-correlation
168
169    ; Significance test functions
170    t-significance
171    f-significance 
172    ; (chi-square significance is calculated from chi-square-cdf
173    ;  in ways depending on the problem)
174
175    ;; The Lambert W function, also called the omega function or product logarithm.
176    lambert-W0
177    lambert-Wm1
178    )
179
180  (import scheme (chicken base) (chicken foreign) (chicken format) (chicken keyword) (chicken sort)
181          srfi-1 (only srfi-13 string<) srfi-25 srfi-69 vector-lib)
182
183  ;; ---------------------------------------------------------------------------
184  ;; Include functions from GNU Scientific Library
185
186  #>
187  #include <assert.h>
188  #include <gsl/gsl_rng.h>
189  #include <gsl/gsl_sf_erf.h>
190  #include <gsl/gsl_sf_gamma.h>
191  #include <gsl/gsl_sf_lambert.h>
192 
193  const gsl_rng_type *statistics_rng_T;
194  gsl_rng *statistics_rng_st;
195
196  void statistics_rng_init()
197  {
198    gsl_rng_env_setup();
199                   
200    statistics_rng_T = gsl_rng_default;
201    assert(statistics_rng_T != NULL);
202    statistics_rng_st = gsl_rng_alloc (statistics_rng_T);
203    assert(statistics_rng_st != NULL);
204    }
205
206  gsl_rng *statistics_rng_get_state ()
207  {
208     return statistics_rng_st;
209  }
210  <#
211
212  (define rng-init (foreign-lambda void "statistics_rng_init"))
213  (define rng-state (foreign-lambda nonnull-c-pointer "statistics_rng_get_state"))
214  (define gsl-random-uniform (foreign-lambda double "gsl_rng_uniform" nonnull-c-pointer))
215  (define gsl-random-uniform-int (foreign-lambda unsigned-int "gsl_rng_uniform_int" nonnull-c-pointer unsigned-int))
216
217  ;; The principal branch of the Lambert W function
218  (define lambert-W0 
219    (foreign-lambda double "gsl_sf_lambert_W0" double))
220
221  ;;The secondary real-valued branch of the Lambert W function
222  (define lambert-Wm1
223    (foreign-lambda double "gsl_sf_lambert_Wm1" double))
224
225  (define beta-incomplete
226    (foreign-lambda double "gsl_sf_beta_inc" double double double))
227
228  (define (random-uniform)
229    (gsl-random-uniform (rng-state)))
230
231  (define (random-uniform-int n)
232    (gsl-random-uniform-int (rng-state) n))
233
234  ;; returns complementary normalised incomplete Gamma Function
235  (define (gamma-incomplete a x)
236    (define gamma-inc (foreign-lambda double "gsl_sf_gamma_inc_P" double double))
237    (values (if (= x 0.0) 0.0 (gamma-inc a x))
238            (gamma-ln a)))
239
240  (define (gamma-ln x)
241    (define lngamma (foreign-lambda double "gsl_sf_lngamma" double))
242    (cond ((<= x 0)
243           (error "Argument to gamma-ln must be positive"))
244          ((> x 1.0e302)
245           (error "Argument to gamma-ln too large"))
246          (else
247            (lngamma x))))
248
249  ;; ---------------------------------------------------------------------------
250  (define :1-beta 'scm-stats-1-beta)
251  (define :alpha 'scm-stats-alpha)
252  (define :approximate? 'scm-stats-approximate)
253  (define :epsilon 'scm-stats-epsilon)
254  (define :exact? 'scm-stats-exact)
255  (define :mu 'scm-stats-mu)
256  (define :one-tailed? 'scm-stats-one-tailed)
257  (define :positive 'scm-stats-positive)
258  (define :precision 'scm-stats-precision)
259  (define :rate 'scm-stats-rate)
260  (define :sample-ratio 'scm-stats-sample-ratio)
261  (define :sigma 'scm-stats-sigma)
262  (define :tails 'scm-stats-tails)
263  (define :variances-equal? 'scm-stats-variances-equal)
264  (define :variance-significance-cutoff 'scm-stats-variance-significance-cutoff)
265  (define :x-tolerance 'scm-stats-x-tolerance)
266  (define :y-tolerance 'scm-stats-y-tolerance)
267
268  ;; ---------------------------------------------------------------------------
269  ;; Some utility procedures
270
271  (define PI 3.1415926535897932385)
272
273  (define (array-shape arr)
274    (let ((r (array-rank arr)))
275      (let recur ((d 0) (shape '()))
276        (if (< d r)
277            (let ((s (- (array-end arr d) (array-start arr d))))
278              (recur (+ 1 d) (cons s shape)))
279            (reverse shape)))
280      ))
281 
282  ;; Average rank calculation for non-parametric tests.  Ranks are 1 based,
283  ;; but lisp is 0 based, so add 1!
284  (define (average-rank value sorted-values)
285    (let ((first (position value sorted-values))
286          (last (- (- (length sorted-values) 1)
287                   (position value (reverse sorted-values)))))
288      (+ 1 (if (= first last)
289             first
290             (/ (+ first last) 2)))))
291
292  ;; BIN-AND-COUNT
293  ;;
294  ;; Make N equal width bins and count the number of elements of sequence that
295  ;; belong in each.
296  (define (bin-and-count sequence n)
297    (let ((smallest  (fold min (car sequence) (cdr sequence)))
298          (increment (/ (range sequence)
299                        n))
300          (bins (make-vector n 0)))
301      (let loop ((bin 0))
302        (if (= bin n)
303          (begin ; PCL: Make sure largest element is included in highest bin count!
304            (vector-set! bins (- n 1) (+ (vector-ref bins (- n 1)) 1))
305            bins)
306          (begin
307            (vector-set! bins 
308                         bin
309                         (length
310                           (filter 
311                             (lambda (x) (and (>= x (+ smallest (* bin increment)))
312                                              (< x  (+ smallest (* (+ 1 bin) increment)))))
313                             sequence)))
314            (loop (+ 1 bin)))))))
315
316  ;; COMBINATIONS
317  ;; How many ways to take n things taken k at a time, when order does not matter
318  ;; Rosner 90
319  (define (combinations n k)
320    (/ (factorial n)
321       (* (factorial k)
322          (factorial (- n k)))))
323
324  (define (error-function x)
325    (let-values (((erf ignore) (gamma-incomplete 0.5 (square x))))
326                (if (>= x 0.0)
327                  erf
328                  (- erf))))
329
330  (define (factorial n)
331    (define (fact-n n result)
332      (if (<= n 1)
333        result
334        (fact-n (- n 1) (* n result))))
335    (fact-n n 1))
336
337  ;; assumes p-function is monotonic (for this library)
338  ;; -- #:increasing defaults to #f, may be set to #t
339  (define (find-critical-value p-function p-value . args) ; adopted from CLASP 1.4.3
340    (let ((increasing? (get-keyword #:increasing args (lambda () #f)))
341          (x-tolerance (get-keyword #:x-tolerance args (lambda () 0.00001)))
342          (y-tolerance (get-keyword #:y-tolerance args (lambda () 0.00001))))
343      (let* ((x-low 0)
344             (fx-low (p-function x-low))
345             (x-high 1)
346             (fx-high (p-function x-high)))
347        ; set initial bounding values
348        (let loop ()
349          (unless ((if increasing? <= >=) fx-low p-value fx-high)
350            (set! x-low x-high)
351            (set! fx-low fx-high)
352            (set! x-high (* 2.0 x-high))
353            (set! fx-high (p-function x-high))
354            (loop)))
355        ; binary search
356        (let loop ()
357          (let* ((x-mid (/ (+ x-low x-high) 2.0))
358                 (fx-mid (p-function x-mid))
359                 (y-diff (abs (- fx-mid p-value)))
360                 (x-diff (- x-high x-low)))
361            (if (or (< x-diff x-tolerance)
362                    (< y-diff y-tolerance))
363              x-mid
364              ;; depending on whether p-function is monotonically increasing with x,
365              ;; if the function is above the desired p-value ...
366              (begin
367                (if ((if increasing? > <) p-value fx-mid)
368                  ;; then critical value is in the upper half
369                  (begin (set! x-low x-mid)
370                         (set! fx-low fx-mid))
371                  ;; otherwise it's in the lower half
372                  (begin (set! x-high x-mid)
373                         (set! fx-high fx-mid)))
374                (loop))))))))
375
376  ;; FISHER-Z-TRANSFORM
377  ;; Rosner 458
378  ;; Transforms the correlation coefficient to an approximately
379  ;; normal distribution.
380  (define (fisher-z-transform r)
381    (* 0.5 (log (/ (+ 1 r) (- 1 r)))))
382
383  ;; taken from TSPL 3
384  (define (list-sort comparison items)
385    (define (dosort pred? ls n)
386      (if (= n 1)
387        (list (car ls))
388        (let ((i (quotient n 2)))
389          (domerge pred?
390                   (dosort pred? ls i)
391                   (dosort pred? (list-tail ls i) (- n i))))))
392    (define (domerge pred? l1 l2)
393      (cond ((null? l1) l2)
394            ((null? l2) l1)
395            ((pred? (car l2) (car l1))
396             (cons (car l2) (domerge pred? l1 (cdr l2))))
397            (else
398              (cons (car l1) (domerge pred? (cdr l1) l2)))))
399    (if (null? items)
400      items
401      (dosort comparison items (length items))))
402
403  ;; PERMUTATIONS
404  ;; How many ways to take n things taken k at a time, when order matters
405  ;; Rosner 88
406  (define (permutations n k)
407    (fold * 1 (iota k n -1)))
408
409  (define (position item items)
410    (let ((index (list-index (lambda (x) (equal? item x))
411                             items)))
412      (if (number? index)
413        index
414        (error "Position: item not in list"))))
415
416  (define (cumsum sequence)
417    (if (null? sequence)
418        0
419        (reverse (fold (lambda (x ax) (cons (+ x (car ax)) ax)) (list (car sequence)) (cdr sequence)))))
420
421  (define (sign x)
422    (cond ((negative? x) -1)
423          ((positive? x) 1)
424          ((zero? x) 0)
425          (else '())))
426
427  (define (square x) 
428    (* x x))
429
430  ;; TODO
431  (define (underflow-goes-to-zero x)
432    x)
433
434  (define (vector-sort comparison items)
435    (list->vector (list-sort comparison (vector->list items))))
436
437  ;; RANDOM-NORMAL
438  ;; returns a random number with mean and standard-distribution as specified.
439  (define (random-normal mean sd)
440    (+ mean (* sd (/ (random-uniform-int 1000000) 1000000))))
441
442  ;; RANDOM-PICK
443  ;; Random selection from list
444  (define (random-pick items)
445    (if (and (list? items) (not (null? items)))
446      (list-ref items (random-uniform-int (length items)))
447      #f))
448
449  ;; RANDOM-SAMPLE
450  ;; Return a random sample of size N from sequence, without replacement. 
451  ;; If N is equal to or greater than the length of the sequence, return
452  ;; the entire sequence.
453  (define (random-sample n items)
454    (cond ((<= n 0) 
455           '())
456          ((>= n (length items)) 
457           items)
458          (else
459            (let loop ((remaining items)
460                       (kept '()))
461              (if (= (length kept) n)
462                kept
463                (let ((one (random-pick remaining)))
464                  (loop (delete one remaining)
465                        (cons one kept))))))))
466
467  ;; RANDOM-WEIGHTED-SAMPLE Return a random sample of size M from
468  ;; sequence of length N, without replacement, where each element has
469  ;; a defined probability of selection (weight) W.  If M is equal to
470  ;; or greater to N, return the entire sequence.
471  (define (random-weighted-sample m items weights)
472    (let ((n (length items)))
473      (cond ((<= m 0) '())
474            ((>= m n) items)
475            (else
476             (let* ((keys (map (lambda (w) (expt (random-uniform) (/ 1 w))) weights))
477                    (sorted-items (list-sort (lambda (x y) (> (car x) (car y))) (zip keys items))))
478               (map cadr (take sorted-items m)))
479             ))
480      ))
481
482               
483  ;; ---------------------------------------------------------------------------
484  ;; Descriptive statistics
485
486  ;; MEAN
487  ;; Rosner 10
488  (define (mean sequence)
489    (if (null? sequence)
490      0
491      (/ (fold + 0 sequence) (length sequence))))
492
493  ;; MEDIAN
494  ;; Rosner 12 (and 19)
495  (define (median sequence)
496    (percentile sequence 50))
497
498  ;; MODE
499  ;; Rosner 14
500  ;; Returns two values: a list of the modes and the number of times they occur
501  (define (mode sequence)
502    (if (null? sequence)
503      (error "Mode: Sequence must not be null")
504      (let ((count-table (make-hash-table eqv?))
505            (modes '())
506            (mode-count 0))
507        (for-each (lambda (item) 
508                    (hash-table-set! count-table
509                                     item
510                                     (+ 1 (hash-table-ref count-table item (lambda () 0)))))
511                  sequence)
512        (for-each (lambda (key) 
513                    (let ((val (hash-table-ref count-table key (lambda () #f))))
514                      (cond ((> val mode-count) ; keep mode
515                             (set! modes (list key))
516                             (set! mode-count val))
517                            ((= val mode-count) ; store multiple modes
518                             (set! modes (cons key modes))))))
519                  (hash-table-keys count-table))
520        (cond ((every number? modes) (set! modes (list-sort < modes)))
521              ((every string? modes) (set! modes (list-sort string< modes)))
522              )
523        (values modes mode-count))))
524
525  ;; GEOMETRIC-MEAN
526  ;; Rosner 16
527  (define (geometric-mean sequence . args)
528    (let ((base (if (null? args) 10 (car args))))
529      (expt base (mean (map (lambda (x) (/ (log x)
530                                           (log base))) 
531                            sequence)))))
532
533  ;; RANGE
534  ;; Rosner 18
535  (define (range sequence)
536    (if (null? sequence)
537      0
538      (- (fold max (car sequence) (cdr sequence))
539         (fold min (car sequence) (cdr sequence)))))
540
541  ;; PERCENTILE
542  ;; Rosner 19
543  (define (percentile sequence percent)
544    (if (null? sequence)
545      (error "Percentile: Sequence must not be null")
546      (let* ((sorted-vec (vector-sort < (apply vector sequence)))
547             (n (vector-length sorted-vec))
548             (k (* n (/ percent 100)))
549             (floor-k (floor k)))
550        (if (= k floor-k)
551          (/ (+ (vector-ref sorted-vec k)
552                (vector-ref sorted-vec (- k 1)))
553             2)
554          (vector-ref sorted-vec floor-k)))))
555
556  ;; VARIANCE
557  ;; Rosner 21
558  (define (variance sequence)
559    (if (< (length sequence) 2)
560      (error "variance: sequence must contain at least two elements")
561      (let ((mean1 (mean sequence)))
562        (/ (fold + 0 (map (lambda (x) (square (- mean1 x))) sequence))
563           (- (length sequence) 1)))))
564
565  ;; STANDARD-DEVIATION
566  ;; Rosner 21
567  (define (standard-deviation sequence)
568    (sqrt (variance sequence)))
569
570  ;; COEFFICIENT-OF-VARIATION
571  ;; Rosner 24
572  (define (coefficient-of-variation sequence)
573    (* 100 
574       (/ (standard-deviation sequence)
575          (mean sequence))))
576
577  ;; STANDARD-ERROR-OF-THE-MEAN
578  ;; Rosner 172
579  (define (standard-error-of-the-mean sequence)
580    (/ (standard-deviation sequence)
581       (sqrt (length sequence))))
582
583  ;; MEAN-SD-N
584  ;; Combined calculation, returns mean, standard deviation and length of
585  ;; sequence of numbers
586  (define (mean-sd-n sequence)
587    (values (mean sequence)
588            (standard-deviation sequence)
589            (length sequence)))
590
591  ;; ---------------------------------------------------------------------------
592  ;; Distributional functions
593
594  ;; BINOMIAL-PROBABILITY
595  ;; exact: Rosner 93, approximate 105
596  ;;
597  ;; P(X=k) for X a binomial random variable with parameters n & p.
598  ;; Binomial expectations for seeing k events in N trials, each having
599  ;; probability p.  Use the Poisson approximation if N>100 and P<0.01.
600  (define (binomial-probability n k p)
601    (if (and (> n 100)
602             (< p 0.01))
603      (poisson-probability (* n p) k)
604      (* (combinations n k)
605         (expt p k)
606         (expt (- 1 p) (- n k)))))
607
608  ;; BINOMIAL-CUMULATIVE-PROBABILITY
609  ;; Rosner 94
610  ;;
611  ;; P(X<k) for X a binomial random variable with parameters n & p.
612  ;; Binomial expectations for fewer than k evens in N trials, each
613  ;; having probability p.
614  (define (binomial-cumulative-probability n k p)
615    (let loop ((sum-upto-k-1 0.0)
616               (i 0))
617      (if (= i k)
618        sum-upto-k-1
619        (loop (+ sum-upto-k-1 (binomial-probability n i p))
620              (+ 1 i)))))
621
622  ;; BINOMIAL-GE-PROBABILITY
623  ;; Rosner 94
624  ;;
625  ;; The probability of k or more occurrences in N events, each with
626  ;; probability p.
627  (define (binomial-ge-probability n k p)
628    (- 1 (binomial-cumulative-probability n k p)))
629
630  ;; for convenience later
631  (define (binomial-le-probability n k p)
632    (let ((sum-up-to-k 0))
633      (let loop ((i 0))
634        (if (= i (+ 1 k))
635          sum-up-to-k
636          (set! sum-up-to-k (+ sum-up-to-k (binomial-probability n i p)))))))
637
638  ;; POISSON-PROBABILITY
639  ;; Rosner 100
640  ;;
641  ;; Probability of seeing k events over a time period when the expected
642  ;; number of events over that time is mu.
643  (define (poisson-probability mu k)
644    (/ (* (exp (- mu)) (expt mu k))
645       (factorial k)))
646
647  ;; POISSON-CUMULATIVE-PROBABILITY
648  ;; Rosner 197
649  ;;
650  ;; Probability of seeing fewer than K events over a time period when the
651  ;; expected number of events over that time is mu.
652  (define (poisson-cumulative-probability mu k)
653    (if (< k 170)
654      (let loop ((sum 0)
655                 (x 0))
656        (if (= x k)
657          sum
658          (loop (+ sum (poisson-probability mu x))
659                (+ x 1))))
660      (- 1 (let-values (((gi ignore) (gamma-incomplete k mu))) gi))))
661
662  ;; POISSON-GE-PROBABILITY
663  ;; Probability of X or more events when expected number of events is mu.
664  (define (poisson-ge-probability mu x)
665    (- 1 (poisson-cumulative-probability mu x)))
666
667  ;; NORMAL-PDF
668  ;; Rosner 115
669  ;; The probability density function (PDF) for a normal distribution
670  ;; with mean mu and variance sigma-squared at point x.
671  (define (normal-pdf x mu sigma-squared)
672    (* (/ (sqrt (* 2 PI sigma-squared)))
673       (exp (- (/ (square (- x mu))
674                  (* 2 sigma-squared))))))
675
676  ;; CONVERT-TO-STANDARD-NORMAL
677  ;; Rosner 130
678  ;; Convert X from a Normal distribution with mean mu and variance sigma
679  ;; to standard normal
680  (define (convert-to-standard-normal x mu sigma)
681    (/ (- x mu) sigma))
682
683  ;; PHI
684  ;; the CDF of standard normal distribution
685  ;; Rosner 125
686  (define (phi x)
687    (* 0.5 (+ 1.0 (error-function (/ x (sqrt 2.0))))))
688
689  ;; Z
690  ;; The inverse normal function, P(X<Zu) = u where X is distributed as
691  ;; the standard normal.  Uses binary search.
692  ;; Rosner 128
693  (define (z percentile . args)
694    (let ((epsilon (get-keyword #:epsilon args (lambda () 1e-15))))
695      (let loop ((min-value -9e0)
696                 (max-value 9e0)
697                 (guess 0.0))
698        (if (< (- max-value min-value) epsilon)
699          guess
700          (if (< (phi guess) percentile)
701            (loop guess max-value (mean (list guess max-value)))
702            (loop min-value guess (mean (list min-value guess))))))))
703
704  ;; T-DISTRIBUTION
705  ;; Rosner 178
706  ;; Returns the point which is the indicated percentile in the T distribution
707  ;; with dof degrees of freedom
708  (define (t-distribution dof percentile)
709    (find-critical-value
710      (lambda (x) (t-significance x dof #:tails ':positive))
711      (- 1 percentile)))
712
713  ;; CHI-SQUARE
714  ;; Rosner 187
715  ;; Returns the point which is the indicated percentile in the Chi Square
716  ;; distribution with dof degrees of freedom.
717  (define (chi-square dof percentile)
718    (find-critical-value (lambda (x) (chi-square-cdf x dof))
719                         (- 1 percentile)
720                         #:increasing #t))
721
722  ;; CHI-SQUARE-CDF
723  ;; Computes the left hand tail area under the chi square
724  ;; distribution under dof degrees of freedom up to X.
725  (define (chi-square-cdf x dof)
726    (let-values (((cdf ignore) (gamma-incomplete (* 0.5 dof) (* 0.5 x))))
727                cdf))
728
729  ;; ---------------------------------------------------------------------------
730  ;; Confidence intervals
731  ;;
732
733  ;; BINOMIAL-PROBABILITY-CI
734  ;; Rosner 193 (approximate) & 194 (exact)
735
736  ;; Confidence intervals on a binomial probability.  If a binomial
737  ;; probability of p has been observed in N trials, what is the 1-alpha
738  ;; confidence interval around p?  Approximate (using normal theory
739  ;; approximation) when npq >= 10 unless told otherwise
740
741  (define (binomial-probability-ci n p alpha . args)
742    (let ((exact (get-keyword #:exact? args (lambda () #f))))
743      (if (and (> (* n p (- 1 p)) 10) 
744               (not exact))
745        (let ((difference (* (z (- 1 (/ alpha 2)))
746                             (sqrt (/ (* p (- 1 p)) n)))))
747          (values (- p difference) (+ p difference)))
748        (values (find-critical-value
749                  (lambda (p1) (binomial-cumulative-probability n (floor (* p n)) p1))
750                  (- 1 (/ alpha 2)))
751                (find-critical-value
752                  (lambda (p2) (binomial-cumulative-probability n (+ 1 (floor (* p n))) p2))
753                  (/ alpha 2))))))
754
755  ;; POISSON-MU-CI
756  ;; Confidence interval for the Poisson parameter mu
757  ;; Rosner 197
758  ;;
759  ;; Given x observations in a unit of time, what is the 1-alpha confidence
760  ;; interval on the Poisson parameter mu (= lambda*T)?
761  ;;
762  ;; Since find-critical-value assumes that the function is monotonic
763  ;; increasing, adjust the value we are looking for taking advantage of
764  ;; reflectiveness
765  (define (poisson-mu-ci x alpha)
766    (values
767      (find-critical-value
768        (lambda (mu) (poisson-cumulative-probability mu (- x 1)))
769        (- 1 (/ alpha 2)))
770      (find-critical-value
771        (lambda (mu) (poisson-cumulative-probability mu x))
772        (/ alpha 2))))
773
774  ;; NORMAL-MEAN-CI
775  ;; Confidence interval for the mean of a normal distribution
776  ;; Rosner 180
777  ;; The 1-alpha percent confidence interval on the mean of a normal
778  ;; distribution with parameters mean, sd & n.
779  (define (normal-mean-ci mean sd n alpha)
780    (let ((t-value (t-distribution (- n 1) (- 1 (/ alpha 2)))))
781      (values (- mean (* t-value (/ sd (sqrt n))))
782              (+ mean (* t-value (/ sd (sqrt n)))))))
783
784  ;; NORMAL-MEAN-CI-ON-SEQUENCE
785  ;;
786  ;; The 1-alpha confidence interval on the mean of a sequence of numbers
787  ;; drawn from a Normal distribution.
788  (define (normal-mean-ci-on-sequence sequence alpha)
789    (let-values (((mean sd n) (mean-sd-n sequence)))
790                (normal-mean-ci mean sd n alpha)))
791
792  ;; NORMAL-VARIANCE-CI
793  ;; Rosner 190
794  ;; The 1-alpha confidence interval on the variance of a sequence of numbers
795  ;; drawn from a Normal distribution.
796  (define (normal-variance-ci variance n alpha)
797    (let ((chi-square-low (chi-square (- n 1) (- 1 (/ alpha 2))))
798          (chi-square-high (chi-square (- n 1) (/ alpha 2)))
799          (the-numerator (* (- n 1) variance)))
800      (values (/ the-numerator chi-square-low) 
801              (/ the-numerator chi-square-high))))
802
803  ;; NORMAL-VARIANCE-CI-ON-SEQUENCE
804  ;;
805  (define (normal-variance-ci-on-sequence sequence alpha)
806    (let ((variance (variance sequence))
807          (n (length sequence)))
808      (normal-variance-ci variance n alpha)))
809
810  ;; NORMAL-SD-CI
811  ;; Rosner 190
812  ;; As above, but a confidence inverval for the standard deviation.
813  (define (normal-sd-ci sd n alpha)
814    (let-values (((low high) (normal-variance-ci (square sd) n alpha)))
815                (values (sqrt low) (sqrt high))))
816
817  ;; NORMAL-SD-CI-ON-SEQUENCE
818  (define (normal-sd-ci-on-sequence sequence alpha)
819    (let ((sd (standard-deviation sequence))
820          (n (length sequence)))
821      (normal-sd-ci sd n alpha)))
822
823  ;; ---------------------------------------------------------------------------
824  ;; Hypothesis testing
825  ;;
826
827  ;; Z-TEST
828  ;; Rosner 228
829
830  ;; The significance of a one sample Z test for the mean of a normal
831  ;; distribution with known variance.  mu is the null hypothesis mean, x-bar
832  ;; is the observed mean, sigma is the standard deviation and N is the number
833  ;; of observations.  If tails is :both, the significance of a difference
834  ;; between x-bar and mu.  If tails is :positive, the significance of x-bar
835  ;; is greater than mu, and if tails is :negative, the significance of x-bar
836  ;; being less than mu.  Returns a p value.
837  (define (z-test x-bar n . args)
838    (let ((mu (get-keyword #:mu args (lambda () 0)))
839          (sigma (get-keyword #:sigma args (lambda () 1)))
840          (tails (get-keyword #:tails args (lambda () ':both))))
841      (let ((z (/ (- x-bar mu) (/ sigma (sqrt n)))))
842        (case tails
843          ((:both) (if (< z 0)
844                     (* 2 (phi z))
845                     (* 2 (- 1 (phi z)))))
846          ((:negative) (phi z))
847          ((:positive) (- 1 (phi z)))))))
848
849  ;; Z-TEST-ON-SEQUENCE
850  ;; Returns a p value for significance of a sample compared with
851  ;; given mean and standard deviation for the null hypothesis.
852  (define (z-test-on-sequence sequence . args)
853    (let ((mu (get-keyword #:mu args (lambda () 0)))
854          (sigma (get-keyword #:sigma args (lambda () 1)))
855          (tails (get-keyword #:tails args (lambda () ':both))))
856      (let ((x-bar (mean sequence))
857            (n (length sequence)))
858        (z-test x-bar n #:mu mu #:sigma sigma #:tails tails))))
859
860  ;; T-TEST-ONE-SAMPLE
861  ;; T-TEST-ONE-SAMPLE-ON-SEQUENCE
862  ;; Rosner 216
863
864  ;; The significance of a one sample T test for the mean of a normal
865  ;; distribution with unknown variance.  X-bar is the observed mean, sd is
866  ;; the observed standard deviation, N is the number of observations and mu
867  ;; is the test mean.  -ON-SAMPLE is the same, but calculates the observed
868  ;; values from a sequence of numbers.
869  (define (t-test-one-sample x-bar sd n mu . args)
870    (let ((tails (get-keyword #:tails args (lambda () ':both))))
871      (t-significance  (/ (- x-bar mu) (/ sd (sqrt n))) (- n 1) #:tails tails)))
872
873  (define (t-test-one-sample-on-sequence sequence mu . args)
874    (let ((tails (get-keyword #:tails args (lambda () ':both))))
875      (let-values (((mean sd n) (mean-sd-n sequence)))
876                  (t-test-one-sample mean sd n mu #:tails tails))))
877
878  ;; T-TEST-PAIRED
879  ;; Rosner 276
880
881  ;; The significance of a paired t test for the means of two normal
882  ;; distributions in a longitudinal study.  D-bar is the mean difference, sd
883  ;; is the standard deviation of the differences, N is the number of pairs.
884  (define (t-test-paired d-bar sd n . args)
885    (let ((tails (get-keyword #:tails args (lambda () ':both))))
886      (t-significance (/ d-bar (/ sd (sqrt n))) (- n 1) #:tails tails)))
887
888  ;; T-TEST-PAIRED-ON-SEQUENCES
889  ;; Rosner 276
890
891  ;; The significance of a paired t test for means of two normal distributions
892  ;; in a longitudinal study.  Before is a sequence of before values, after is
893  ;; the sequence of paired after values (which must be the same length as the
894  ;; before sequence).
895  (define (t-test-paired-on-sequences before after . args)
896    (let ((tails (get-keyword #:tails args (lambda () ':both))))
897      (let-values (((mean sd n) (mean-sd-n (map - before after))))
898                  (t-test-paired mean sd n #:tailed tails))))
899
900  ;; T-TEST-TWO-SAMPLE
901  ;; Rosner  282, 294, 297
902
903  ;; The significance of the difference of two means (x-bar1 and x-bar2) with
904  ;; standard deviations sd1 and sd2, and sample sizes n1 and n2 respectively.
905  ;; The form of the two sample t test depends on whether the sample variances
906  ;; are equal or not.   If the variable variances-equal? is :test, then we
907  ;; use an F test and the variance-significance-cutoff to determine if they
908  ;; are equal.  If the variances are equal, then we use the two sample t test
909  ;; for equal variances.  If they are not equal, we use the Satterthwaite
910  ;; method, which has good type I error properties (at the loss of some
911  ;; power). 
912  (define (t-test-two-sample x-bar1 sd1 n1 x-bar2 sd2 n2 . args)
913    (let ((variances-equal? (get-keyword #:variances-equal? args (lambda () ':test)))
914          (variances-significance-cutoff (get-keyword #:variance-significance-cutoff args (lambda () 0.05)))
915          (tails (get-keyword #:tails args (lambda () ':both))))
916      (let ((t-statistic #f)
917            (dof #f))
918        (if (case variances-equal?
919              ((:test)
920               (> (f-test (square sd1) n1 (square sd2) n2 #:tails tails)
921                  variances-significance-cutoff))
922              ((#t :yes :equal)
923               #t)
924              ((#f :no :unequal)
925               #f))
926          (let ((s (sqrt (/ (+ (* (- n1 1) (square sd1))
927                               (* (- n2 1) (square sd2)))
928                            (+ n1 n2 -2)))))
929            (set! t-statistic (/ (- x-bar1 x-bar2)
930                                 (* s (sqrt (+ (/ n1) (/ n2))))))
931            (set! dof (+ n1 n2 -2)))
932          (let ((variance-ratio-1 (/ (square sd1) n1))
933                (variance-ratio-2 (/ (square sd2) n2)))
934            (set! t-statistic (/ (- x-bar1 x-bar2)
935                                 (sqrt (+ variance-ratio-1 variance-ratio-2))))
936            (set! dof (round (/ (square (+ variance-ratio-1 variance-ratio-2))
937                                (+ (/ (square variance-ratio-1) (- n1 1))
938                                   (/ (square variance-ratio-2) (- n2 1))))))))
939        (t-significance t-statistic dof #:tails tails))))
940
941  ;; T-TEST-TWO-SAMPLE-ON-SEQUENCES
942  ;; Same as above, but providing the sequences rather than the summaries
943  (define (t-test-two-sample-on-sequences sequence1 sequence2 . args)
944    (let ((variance-significance-cutoff (get-keyword #:variance-significance-cutoff args (lambda () 0.05)))
945          (tails (get-keyword #:tails args (lambda () ':both))))
946      (let-values (((x-bar-1 sd1 n1) (mean-sd-n sequence1))
947                   ((x-bar-2 sd2 n2) (mean-sd-n sequence2)))
948                  (t-test-two-sample x-bar-1 sd1 n1 x-bar-2 sd2 n2
949                                     #:variances-equal? ':test 
950                                     #:variance-significance-cutoff variance-significance-cutoff 
951                                     #:tails tails))))
952
953  ;; F-TEST
954  ;; Rosner 290
955  ;; F test for the equality of two variances
956  (define (f-test variance1 n1 variance2 n2 . args)
957    (let ((tails (get-keyword #:tails args (lambda () ':both))))
958      (let ((significance (f-significance (/ variance1 variance2) (- n1 1) (- n2 1)
959                                          (not (eq? tails ':both)))))
960        (case tails
961          ((:both) significance)
962          ((:positive) (if (> variance1 variance2) significance (- 1 significance)))
963          ((:negative) (if (< variance1 variance2) significance (- 1 significance)))))))
964
965  ;; CHI-SQUARE-TEST-ONE-SAMPLE
966  ;; Rosner 246
967  ;; The significance of a one sample Chi square test for the variance of a
968  ;; normal distribution.  Variance is the observed variance, N is the number
969  ;; of observations, and sigma-squared is the test variance.
970  (define (chi-square-test-one-sample variance n sigma-squared . args)
971    (let ((tails (get-keyword #:tails args (lambda () ':both))))
972      (let ((cdf (chi-square-cdf (/ (* (- n 1) variance)
973                                    sigma-squared)
974                                 (- n 1))))
975        (case tails
976          ((:negative)
977           cdf)
978          ((:positive)
979           (- 1 cdf))
980          ((:both)
981           (if (<= variance sigma-squared)
982             (* 2 cdf)
983             (* 2 (- 1 cdf))))))))
984
985  ;; BINOMIAL-TEST-ONE-SAMPLE
986  ;; Rosner 249
987  ;; The significance of a one sample test for the equality of an observed
988  ;; probability p-hat to an expected probability p under a binomial
989  ;; distribution with N observations.  Use the normal theory approximation if
990  ;; n*p(1-p) > 10 (unless the exact flag is true).
991  (define (binomial-test-one-sample p-hat n p . args)
992    (let ((tails (get-keyword #:tails args (lambda () ':both)))
993          (exact (get-keyword #:exact? args (lambda () #f))))
994      (let ((q (- 1 p)))
995        (if (and (> (* n p q) 10) (not exact))
996          (let ((z (/ (- p-hat p) (sqrt (/ (* p q) n)))))
997            (case tails
998              ((:negative) (phi z))
999              ((:positive) (- 1 (phi z)))
1000              ((:both) (* 2 (if (<= p-hat p) (phi z) (- 1 (phi z)))))))
1001          (let* ((observed (round (* p-hat n)))
1002                 (probability-more-extreme
1003                   (if (<= p-hat p)
1004                     (binomial-cumulative-probability n observed p)
1005                     (binomial-ge-probability n observed p))))
1006            (case tails
1007              ((:negative :positive) probability-more-extreme)
1008              ((:both) (min (* 2 probability-more-extreme) 1.0))))))))
1009
1010  ;; BINOMIAL-TEST-TWO-SAMPLE
1011  ;; Rosner 357
1012  ;; Are the observed probabilities of an event (p-hat1 and p-hat2) in N1/N2
1013  ;; trials different?  The normal theory method implemented here.  The exact
1014  ;; test is Fisher's contingency table method, below.
1015  (define (binomial-test-two-sample p-hat1 n1 p-hat2 n2 . args)
1016    (let ((tails (get-keyword #:tails args (lambda () ':both)))
1017          (exact (get-keyword #:exact? args (lambda () #f))))
1018      (let* ((p-hat (/ (+ (* p-hat1 n1) (* p-hat2 n2)) (+ n1 n2)))
1019             (q-hat (- 1 p-hat))
1020             (z (/ (- (abs (- p-hat1 p-hat2)) (+ (/ (* 2 n1)) (/ (* 2 n2))))
1021                   (sqrt (* p-hat q-hat (+ (/ n1) (/ n2)))))))
1022        (if (and (> (* n1 p-hat q-hat) 5)
1023                 (> (* n2 p-hat q-hat) 5)
1024                 (not exact))
1025          (* (- 1 (phi z)) (if (eq? tails ':both) 2 1))
1026          (fisher-exact-test (inexact->exact (round (* p-hat1 n1)))
1027                             (inexact->exact (round (* (- 1 p-hat1) n1)))
1028                             (inexact->exact (round (* p-hat2 n2)))
1029                             (inexact->exact (round (* (- 1 p-hat2) n2)))
1030                             tails)))))
1031
1032  ;; FISHER-EXACT-TEST
1033  ;; Rosner 371
1034  ;; Fisher's exact test.  Gives a p value for a particular 2x2 contingency table
1035  (define (fisher-exact-test a b c d . args)
1036    (let ((tails (get-keyword #:tails args (lambda () ':both))))
1037      (define (table-probability ta tb tc td)
1038        (let ((tn (+ ta tb tc td)))
1039          (/ (* 1.0 (factorial (+ ta tb)) (factorial (+ tc td))
1040                (factorial (+ ta tc)) (factorial (+ tb td)))
1041             (* 1.0 (factorial tn) (factorial ta) (factorial tb)
1042                (factorial tc) (factorial td)))))
1043      (let* ((row-margin1 (+ a b))
1044             (row-margin2 (+ c d))
1045             (column-margin1 (+ a c))
1046             (column-margin2 (+ b d))
1047             (n (+ a b c d))
1048             (table-probabilities (make-vector (+ 1 (min row-margin1
1049                                                         row-margin2
1050                                                         column-margin1
1051                                                         column-margin2))
1052                                               0)))
1053        ;; rearrange so that the first row and column marginals are
1054        ;; smallest.  Only need to change first margins and a.
1055        (cond ((and (< row-margin2 row-margin1)
1056                    (< column-margin2 column-margin1))
1057               (set! a d)
1058               (set! row-margin1 row-margin2)
1059               (set! column-margin1 column-margin2))
1060              ((< row-margin2 row-margin1)
1061               (set! a c)
1062               (set! row-margin1 row-margin2))
1063              ((< column-margin2 column-margin1)
1064               (set! a b)
1065               (set! column-margin1 column-margin2)))
1066        (let loop ((i 0))
1067          (unless (= i (vector-length table-probabilities))
1068            (let* ((test-a i)
1069                   (test-b (- row-margin1 i))
1070                   (test-c (- column-margin1 i))
1071                   (test-d (- n (+ test-a test-b test-c))))
1072              (vector-set! table-probabilities
1073                           i 
1074                           (table-probability test-a test-b test-c test-d)))
1075            (loop (+ 1 i))))
1076        (let ((above 0.0)
1077              (below 0.0))
1078          (let loop ((i 0))
1079            (unless (= i (vector-length table-probabilities))
1080              (if (< i (+ 1 a))
1081                (set! above (+ above (vector-ref table-probabilities i)))
1082                (set! below (+ below (vector-ref table-probabilities i))))
1083              (loop (+ 1 i))))
1084          (case tails
1085            ((:both) (* 2 (min above below)))
1086            ((:positive) below)
1087            ((:negative) above))))))
1088
1089  ;; MCNEMARS-TEST
1090  ;; Rosner 379 and 381
1091
1092  ;; McNemar's test for correlated proportions, used for longitudinal
1093  ;; studies.  Look only at the number of discordant pairs (one treatment
1094  ;; is effective and the other is not).  If the two treatments are A
1095  ;; and B, a-discordant-count is the number where A worked and B did not,
1096  ;; and b-discordant-count is the number where B worked and A did not.
1097  (define (mcnemars-test a-discordant-count b-discordant-count . args)
1098    (let ((exact (get-keyword #:exact? args (lambda () #f))))
1099      (let ((n (+ a-discordant-count b-discordant-count)))
1100        (if (and (> n 20) (not exact))
1101          (let ((x2 (/ (square (- (abs (- a-discordant-count b-discordant-count)) 1))
1102                       n)))
1103            (- 1 (chi-square-cdf x2 1)))
1104          (cond ((= a-discordant-count b-discordant-count) 
1105                 1.0)
1106                ((< a-discordant-count b-discordant-count)
1107                 (* 2 (binomial-le-probability n a-discordant-count 0.5)))
1108                (else
1109                  (* 2 (binomial-ge-probability n a-discordant-count 0.5))))))))
1110
1111  ;; POISSON-TEST-ONE-SAMPLE
1112  ;; Rosner 256 (approximation on 259)
1113  ;; The significance of a one sample test for the equality of an observed
1114  ;; number of events (observed) and an expected number mu under the
1115  ;; poisson distribution.  Normal theory approximation is not that great,
1116  ;; so don't use it unless told.
1117  (define (poisson-test-one-sample observed mu . args)
1118    (let ((tails (get-keyword #:tails args (lambda () ':both)))
1119          (approximate (get-keyword #:approximate? args (lambda () #f))))
1120      (if approximate
1121        (let ((x-square (/ (square (- observed mu)) mu)))
1122          (- 1 (chi-square-cdf x-square 1)))
1123        (let ((probability-more-extreme 
1124                (if (< observed mu)
1125                  (poisson-cumulative-probability mu observed)
1126                  (poisson-ge-probability mu observed))))
1127          (case tails
1128            ((:negative :positive) probability-more-extreme)
1129            ((:both) (min (* 2 probability-more-extreme) 1.0)))))))
1130
1131  ;; SIGN-TEST
1132  ;; Rosner 335-7
1133  ;; Really, just a special case of the binomial one sample test with p = 1/2.
1134  ;; The normal theory version has a correction factor to make it a better
1135  ;; approximation.
1136  (define (sign-test plus-count minus-count . args)
1137    (let ((exact (get-keyword #:exact? args (lambda () #f)))
1138          (tails (get-keyword #:tails args (lambda () ':both))))
1139      (let* ((n (+ plus-count minus-count))
1140             (p-hat (/ plus-count n)))
1141        (if (or (< n 20) exact)
1142          (binomial-test-one-sample p-hat n 0.5 :tails tails :exact? #t)
1143          (let ((area (- 1 (phi (/ (- (abs (- plus-count minus-count)) 1)
1144                                   (sqrt n))))))
1145            (if (eq? tails ':both)
1146              (* 2 area)
1147              area))))))
1148
1149  ;; SIGN-TEST-ON-SEQUENCE
1150  ;; Same as above, but takes two sequences and tests whether the entries
1151  ;; in one are different (greater or less) than the other.
1152  (define (sign-test-on-sequence sequence1 sequence2 . args)
1153    (let ((exact (get-keyword #:exact? args (lambda () #f)))
1154          (tails (get-keyword #:tails args (lambda () ':both))))
1155      (let* ((differences (map - sequence1 sequence2))
1156             (plus-count (length (filter positive? differences)))
1157             (minus-count (length (filter negative? differences))))
1158        (sign-test plus-count minus-count :exact? exact :tails tails))))
1159
1160  ;; WILCOXON-SIGNED-RANK-TEST
1161  ;; Rosner 341
1162  ;; A test on the ranking of positive and negative differences (are the
1163  ;; positive differences significantly larger/smaller than the negative
1164  ;; ones). Assumes a continuous and symmetric distribution of differences,
1165  ;; although not a normal one.  This is the normal theory approximation,
1166  ;; which is only valid when N > 15.
1167
1168  ;; This test is completely equivalent to the Mann-Whitney test.
1169  (define (wilcoxon-signed-rank-test differences . args)
1170    (let ((tails (get-keyword #:tails args (lambda () ':both))))
1171      (let* ((nonzero-differences (filter (lambda (n) (not (zero? n))) differences))
1172             (sorted-list (list-sort (lambda (x y) (< (car x) (car y)))
1173                                     (map (lambda (dif)
1174                                            (list (abs dif)
1175                                                  (sign dif)))
1176                                          nonzero-differences)))
1177             (distinct-values (delete-duplicates (map car sorted-list)))
1178             (ties '()))
1179        (when (< (length nonzero-differences) 16)
1180          (error 
1181            "This Wilcoxon Signed-Rank Test (normal approximation method) requires nonzero N > 15"))
1182
1183        ; add avg-rank to the sorted values
1184        (for-each (lambda (value)
1185                    (let ((first (position value (map car sorted-list)))
1186                          (last (position value (reverse (map car sorted-list)))))
1187                      (if (= first last)
1188                        (append (find (lambda (item) (= (car item) value))
1189                                      sorted-list)
1190                                (list (+ 1 first)))
1191                        (let ((number-tied (+ 1 (- last first)))
1192                              (avg-rank (+ 1 (/ (+ first last) 2)))) ; + 1 since 0 based
1193                          (set! ties (cons number-tied ties))
1194                          (let loop ((i 0)
1195                                     (result '()))
1196                            (if (= i number-tied)
1197                              (reverse result)
1198                              (loop (+ 1 i)
1199                                    (cons (cons (list-ref sorted-list (+ first i))
1200                                                (list avg-rank))
1201                                          result))))))))
1202                  distinct-values)
1203        (set! ties (reverse ties))
1204        (let* ((direction (if (eq? tails ':negative) -1 1))
1205               (r1 (fold + 0
1206                         (map (lambda (entry) 
1207                                (if (= (cadr entry) direction)
1208                                  (caddr entry)
1209                                  0))
1210                              sorted-list)))
1211               (n (length nonzero-differences))
1212               (expected-r1 (/ (* n (+ 1 n)) 4))
1213               (ties-factor (if ties
1214                              (/ (fold + 0 
1215                                       (map (lambda (ti) (- (* ti ti ti) ti))
1216                                            ties))
1217                                 48)
1218                              0))
1219               (var-r1 (- (/ (* n (+ 1 n) (+ 1 (* 2 n))) 24) ties-factor))
1220               (T-score (/ (- (abs (- r1 expected-r1)) 0.5) (sqrt var-r1))))
1221          (* (if (eq? tails ':both) 2 1)
1222             (- 1 (phi T-score)))))))
1223
1224  (define (wilcoxon-signed-rank-test-on-sequences sequence1 sequence2 . args)
1225    (let ((tails (get-keyword #:tails args (lambda () ':both))))
1226      (wilcoxon-signed-rank-test (map - sequence1 sequence2) tails)))
1227
1228  ;; CHI-SQUARE-TEST-RXC
1229  ;; Rosner 395
1230  ;; Takes contingency-table, an RxC array, and returns the significance of
1231  ;; the relationship between the row variable and the column variable.  Any
1232  ;; difference in proportion will cause this test to be significant --
1233  ;; consider using the test for trend instead if you are looking for a
1234  ;; consistent change.
1235  ;; contingency-table is an array, as defined using srfi-25
1236  (define (chi-square-test-rxc contingency-table)
1237    (let* ((shape (array-shape contingency-table))
1238           (rows (car shape))
1239           (columns (cadr shape))
1240           (row-marginals (make-vector rows 0.0))
1241           (column-marginals (make-vector columns 0.0))
1242           (total 0.0)
1243           (expected-lt-5 0)
1244           (expected-lt-1 0)
1245           (expected-values (make-array '#(0) rows columns))
1246           (x2 0.0))
1247      (let loop ((i 0)
1248                 (j 0))
1249        (cond ((= j columns)
1250               (loop (+ 1 i) 0))
1251              ((= i rows)
1252               '())
1253              (else
1254                (let ((cell (array-ref contingency-table i j)))
1255                  (vector-set! row-marginals i cell)
1256                  (vector-set! column-marginals j cell)
1257                  (set! total (+ total cell))))))
1258      (let loop ((i 0)
1259                 (j 0))
1260        (cond ((= j columns)
1261               (loop (+ 1 i) 0))
1262              ((= i rows)
1263               '())
1264              (else
1265                (let ((expected (/ (* (vector-ref row-marginals i) (vector-ref column-marginals j))
1266                                   total)))
1267                  (when (< expected 1) (set! expected-lt-1 (+ 1 expected-lt-1)))
1268                  (when (< expected 5) (set! expected-lt-5 (+ 1 expected-lt-5)))
1269                  (array-set! expected-values i j expected)))))
1270      (when (positive? expected-lt-1)
1271        (error "This test cannot be used when an expected value is less than one"))
1272      (when (> expected-lt-5 (/ (* rows columns) 5))
1273        (error "This test cannot be used when more than 1/5 of the expected values are less than 5."))
1274      (let loop ((i 0)
1275                 (j 0))
1276        (cond ((= j columns)
1277               (loop (+ 1 i) 0))
1278              ((= i rows) 
1279               '())
1280              (else
1281                (set! x2 (/ (square (- (array-ref contingency-table i j)
1282                                       (array-ref expected-values i j)))
1283                            (array-ref expected-values i j))))))
1284      (- 1 (chi-square-cdf x2 (* (- rows 1) (- columns 1))))))
1285
1286  ;; CHI-SQUARE-TEST-FOR-TREND
1287  ;; Rosner 398
1288  ;; This test works on a 2xk table and assesses if there is an increasing or
1289  ;; decreasing trend.  Arguments are equal sized lists counts.  Optionally,
1290  ;; provide a list of scores, which represent some numeric attribute of the
1291  ;; group.  If not provided, scores are assumed to be 1 to k.
1292  (define (chi-square-test-for-trend row1-counts row2-counts . args)
1293    (let ((scores (if (and (not (null? args))
1294                           (= 1 (length args)))
1295                    (car args)
1296                    (let loop ((i 0)
1297                               (res '()))
1298                      (if (= i (length row1-counts))
1299                        (reverse res)
1300                        (loop (+ 1 i)
1301                              (cons (+ 1 i) res)))))))
1302      (let* ((ns (map + row1-counts row2-counts))
1303             (p-hats (map / row1-counts ns))
1304             (n (fold + 0 ns))
1305             (p-bar (/ (fold + 0 row1-counts) n))
1306             (q-bar (- 1 p-bar))
1307             (s-bar (mean scores))
1308             (a (fold + 0.0
1309                      (map (lambda (p-hat ni s)
1310                             (* ni (- p-hat p-bar) (- s s-bar)))
1311                           p-hats ns scores)))
1312             (b (* 1.0 p-bar q-bar (- (fold + 0 (map (lambda (ni s) (* ni (square s)))
1313                                                 ns scores))
1314                                  (/ (square (fold + 0 (map (lambda (ni s) (* ni s))
1315                                                            ns scores)))
1316                                     n))))
1317             (x2 (/ (square a) b))
1318             (significance (- 1 (chi-square-cdf x2 1))))
1319        (when (< (* p-bar q-bar n) 5)
1320          (error 
1321            "This test is only applicable when N * p-bar * q-bar >= 5"))
1322        (format #t "~%The trend is ~a, p=~a~&" (if (< a 0) "decreasing" "increasing") significance)
1323        significance)))
1324
1325  ;; ---------------------------------------------------------------------------
1326  ;; Sample size estimates
1327  ;;
1328
1329  ;; T-TEST-ONE-SAMPLE-SSE
1330  ;; Rosner 238
1331  ;; Returns the number of subjects needed to test whether the mean of a
1332  ;; normally distributed sample mu is different from a null hypothesis mean
1333  ;; mu-null and variance variance, with alpha, 1-beta and tails as specified.
1334  (define (t-test-one-sample-sse mu mu-null variance . args)
1335    (let ((alpha (get-keyword #:alpha args (lambda () 0.05)))
1336          (one-beta (get-keyword #:1-beta args (lambda () 0.95)))
1337          (tails (get-keyword #:tails args (lambda () ':both))))
1338      (let ((z-beta (z one-beta))
1339            (z-alpha (z (- 1 (if (eq? tails ':both) (/ alpha 2) alpha)))))
1340        (inexact->exact
1341          (ceiling (/ (* variance (square (+ z-beta z-alpha)))
1342                      (square (- mu-null mu))))))))
1343
1344  ;; T-TEST-TWO-SAMPLE-SSE
1345  ;; Rosner 308
1346  ;; Returns the number of subjects needed to test whether the mean mu1 of a
1347  ;; normally distributed sample (with variance variance1) is different from a
1348  ;; second sample with mean mu2 and variance variance2, with alpha, 1-beta
1349  ;; and tails as specified.  It is also possible to set a sample size ratio
1350  ;; of sample 1 to sample 2.
1351  (define (t-test-two-sample-sse mu1 variance1 mu2 variance2 . args)
1352    (let ((sample-ratio (get-keyword #:sample-ratio args (lambda () 1)))
1353          (alpha (get-keyword #:alpha args (lambda () 0.05)))
1354          (one-beta (get-keyword #:1-beta args (lambda () 0.95)))
1355          (tails (get-keyword #:tails args (lambda () ':both))))
1356      (let* ((delta2 (square (- mu1 mu2)))
1357             (z-term (square (+ (z one-beta)
1358                                (z (- 1 (if (eq? tails ':both)
1359                                          (/ alpha 2)
1360                                          alpha))))))
1361             (n1 (ceiling (/ (* (+ variance1 (/ variance2 sample-ratio)) z-term)
1362                             delta2))))
1363        (values (inexact->exact n1)
1364                (inexact->exact (ceiling (* sample-ratio n1)))))))
1365
1366  ;; T-TEST-PAIRED-SSE
1367  ;; Rosner 311
1368  ;; Returns the number of subjects needed to test whether the differences
1369  ;; with mean difference-mu and variance difference-variance, with alpha,
1370  ;; 1-beta and tails as specified.
1371  (define (t-test-paired-sse difference-mu difference-variance . args)
1372    (let ((alpha (get-keyword #:alpha args (lambda () 0.05)))
1373          (one-beta (get-keyword #:1-beta args (lambda () 0.95)))
1374          (tails (get-keyword #:tails args (lambda () ':both))))
1375      (ceiling (/ (* 2 difference-variance
1376                     (square (+ (z one-beta)
1377                                (z (- 1 (if (eq? tails ':both)
1378                                          (/ alpha 2)
1379                                          alpha))))))
1380                  (square difference-mu)))))
1381
1382  ;; BINOMIAL-TEST-ONE-SAMPLE-SSE
1383  ;; Rosner 254
1384  ;; Returns the number of subjects needed to test whether an observed
1385  ;; probability is significantly different from a particular binomial
1386  ;; null hypothesis with a significance alpha and a power 1-beta.
1387  (define (binomial-test-one-sample-sse p-estimated p-null . args)
1388    (let ((alpha (get-keyword #:alpha args (lambda () 0.05)))
1389          (one-beta (get-keyword #:1-beta args (lambda () 0.95)))
1390          (tails (get-keyword #:tails args (lambda () ':both))))
1391      (let ((q-null (- 1 p-null))
1392            (q-estimated (- 1 p-estimated)))
1393        (inexact->exact
1394          (ceiling
1395            (/ (* p-null
1396                  q-null
1397                  (square (+ (z (- 1 (if (eq? tails ':both)
1398                                       (/ alpha 2)
1399                                       alpha)))
1400                             (* (z one-beta)
1401                                (sqrt (/ (* p-estimated q-estimated)
1402                                         (* p-null q-null)))))))
1403               (square (- p-estimated p-null))))))))
1404
1405  ;; BINOMIAL-TEST-TWO-SAMPLE-SSE
1406  ;; Rosner 384
1407  ;; The number of subjects needed to test if two binomial probabilities
1408  ;; are different at a given significance alpha and power 1-beta.  The
1409  ;; sample sizes can be unequal; the p2 sample is sample-sse-ratio * the
1410  ;; size of the p1 sample.  It can be a one tailed or two tailed test.
1411  (define (binomial-test-two-sample-sse p1 p2 . args)
1412    (let ((alpha (get-keyword #:alpha args (lambda () 0.05)))
1413          (sample-ratio (get-keyword #:sample-ratio args (lambda () 1)))
1414          (one-beta (get-keyword #:1-beta args (lambda () 0.95)))
1415          (tails (get-keyword #:tails args (lambda () ':both))))
1416      (let* ((q1 (- 1 p1))
1417             (q2 (- 1 p2))
1418             (delta (abs (- p1 p2)))
1419             (p-bar (/ (+ p1 (* sample-ratio p2)) (+ 1 sample-ratio)))
1420             (q-bar (- 1 p-bar))
1421             (z-alpha (z (- 1 (if (eq? tails ':both) (/ alpha 2) alpha))))
1422             (z-beta (z one-beta))
1423             (n1 (ceiling
1424                   (/ (square (+ (* (sqrt (* p-bar q-bar (+ 1 (/ sample-ratio))))
1425                                    z-alpha)
1426                                 (* (sqrt (+ (* p1 q1) (/ (* p2 q2) sample-ratio)))
1427                                    z-beta)))
1428                      (square delta)))))
1429        (values n1 (ceiling (* sample-ratio n1))))))
1430
1431  ;; BINOMIAL-TEST-PAIRED-SSE
1432  ;; Rosner 387
1433
1434  ;; Sample size estimate for the McNemar (discordant pairs) test.  Pd is the
1435  ;; projected proportion of discordant pairs among all pairs, and Pa is the
1436  ;; projected proportion of type A pairs among discordant pairs.  alpha,
1437  ;; 1-beta and tails are as above.  Returns the number of individuals
1438  ;; necessary; that is twice the number of matched pairs necessary.
1439  (define (binomial-test-paired-sse pd pa . args)
1440    (let ((alpha (get-keyword #:alpha args (lambda () 0.05)))
1441          (one-beta (get-keyword #:1-beta args (lambda () 0.95)))
1442          (tails (get-keyword #:tails args (lambda () ':both))))
1443      (let ((qa (- 1 pa))
1444            (z-alpha (z (- 1 (if (eq? tails ':both) (/ alpha 2) alpha))))
1445            (z-beta (z one-beta)))
1446        (ceiling (/ (square (+ z-alpha (* 2 z-beta (sqrt (* pa qa)))))
1447                    (* 2 (square (- pa 0.5)) pd))))))
1448
1449  ;; CORRELATION-SSE
1450  ;; Rosner 463
1451  ;;
1452  ;; Returns the size of a sample necessary to find a correlation of
1453  ;; expected value rho with significance alpha and power 1-beta.
1454  (define (correlation-sse rho . args)
1455    (let ((alpha (get-keyword #:alpha args (lambda () 0.05)))
1456          (one-beta (get-keyword #:1-beta args (lambda () 0.95))))
1457      (inexact->exact
1458        (ceiling (+ 3 (/ (square (+ (z (- 1 alpha)) (z one-beta)))
1459                         (square (fisher-z-transform rho))))))))
1460
1461  ;; ---------------------------------------------------------------------------
1462  ;; Correlation and regression
1463  ;;
1464
1465  ;; LINEAR-REGRESSION
1466  ;; Rosner 31, 441 for t-test
1467  ;; Computes the regression equation for a least squares fit of a line to a
1468  ;; sequence of points (each a list of two numbers, e.g. '((1.0 0.1) (2.0 0.2)))
1469  ;; and reports the intercept, slope, correlation coefficient r, R^2, and
1470  ;; the significance of the difference of the slope from 0.
1471  (define (linear-regression points)
1472    (unless (> (length points) 2)
1473      (error "Requires at least three points"))
1474    (let ((xs (map car points))
1475          (ys (map cadr points)))
1476      (let* ((x-bar (mean xs))
1477             (y-bar (mean ys))
1478             (n (length points))
1479             (Lxx (fold + 0
1480                        (map (lambda (xi) (square (- xi x-bar))) xs)))
1481             (Lyy (fold + 0
1482                        (map (lambda (yi) (square (- yi y-bar))) ys)))
1483             (Lxy (fold + 0
1484                        (map (lambda (xi yi) (* (- xi x-bar) (- yi y-bar)))
1485                             xs ys)))
1486             (b (if (zero? Lxx) 0 (/ Lxy Lxx)))
1487             (a (- y-bar (* b x-bar)))
1488             (reg-ss (* b Lxy))
1489             (res-ms (/ (- Lyy reg-ss) (- n 2)))
1490             (r (/ Lxy (sqrt (* Lxx Lyy))))
1491             (r2 (/ reg-ss Lyy))
1492             (t-test (/ b (sqrt (/ res-ms Lxx))))
1493             (t-sign (t-significance t-test (- n 2) #:tails ':both)))
1494        (format #t "~%Intercept = ~A, slope = ~A, r = ~A, R^2 = ~A, p = ~A~%"
1495                a b r r2 t-sign)
1496        (values a b r r2 t-sign))))
1497
1498  ;; CORRELATION-COEFFICIENT
1499  ;; Also called Pearson Correlation
1500  (define (correlation-coefficient points)
1501    (let* ((xs (map car points))
1502           (ys (map cadr points))
1503           (x-bar (mean xs))
1504           (y-bar (mean ys)))
1505      (/ (fold + 0 (map (lambda (xi yi) (* (- xi x-bar)
1506                                           (- yi y-bar)))
1507                        xs
1508                        ys))
1509         (sqrt (* (fold + 0 (map (lambda (xi) (square (- xi x-bar)))
1510                                 xs))
1511                  (fold + 0 (map (lambda (yi) (square (- yi y-bar)))
1512                                 ys)))))))
1513
1514
1515  ;; CORRELATION-TEST-TWO-SAMPLE
1516  ;; Rosner 464
1517  ;; Test if two correlation coefficients are different.  Uses Fisher's Z test.
1518  (define (correlation-test-two-sample r1 n1 r2 n2 . args)
1519    (let ((tails (get-keyword #:tails args (lambda () ':both))))
1520      (let* ((z1 (fisher-z-transform r1))
1521             (z2 (fisher-z-transform r2))
1522             (my-lambda (/ (- z1 z2)
1523                           (sqrt (+ (/ (- n1 3))
1524                                    (/ (- n2 3)))))))
1525        (case tails
1526          ((:both) (* 2 (if (<= my-lambda 0) (phi my-lambda) (- 1 (phi my-lambda)))))
1527          ((:positive) (- 1 (phi my-lambda)))
1528          ((:negative) (phi my-lambda))))))
1529
1530  (define (correlation-test-two-sample-on-sequences points1 points2 . args)
1531    (let ((tails (get-keyword #:tails args (lambda () ':both))))
1532      (let ((r1 (correlation-coefficient points1))
1533            (n1 (length points1))
1534            (r2 (correlation-coefficient points2))
1535            (n2 (length points2)))
1536        (correlation-test-two-sample r1 n1 r2 n2 tails))))
1537
1538  ;; SPEARMAN-RANK-CORRELATION
1539  ;; Rosner 498
1540
1541  ;; Spearman rank correlation computes the relationship between a pair of
1542  ;; variables when one or both are either ordinal or have a distribution that
1543  ;; is far from normal.   It takes a list of points (same format as
1544  ;; linear-regression) and returns the spearman rank correlation coefficient
1545  ;; and its significance.
1546  (define (spearman-rank-correlation points)
1547    (let ((xis (map car points))
1548          (yis (map cadr points)))
1549      (let* ((n (length points))
1550             (sorted-xis (list-sort < xis))
1551             (sorted-yis (list-sort < yis))
1552             (average-x-ranks (map (lambda (x) (average-rank x sorted-xis)) xis))
1553             (average-y-ranks (map (lambda (y) (average-rank y sorted-yis)) yis))
1554             (mean-x-rank (mean average-x-ranks))
1555             (mean-y-rank (mean average-y-ranks))
1556             (Lxx (fold + 0 
1557                        (map (lambda (xi-rank) (square (- xi-rank mean-x-rank)))
1558                             average-x-ranks)))
1559             (Lyy (fold + 0
1560                        (map (lambda (yi-rank) (square (- yi-rank mean-y-rank)))
1561                             average-y-ranks)))
1562             (Lxy (fold + 0
1563                        (map (lambda (xi-rank yi-rank)
1564                               (* (- xi-rank mean-x-rank)
1565                                  (- yi-rank mean-y-rank)))
1566                             average-x-ranks average-y-ranks)))
1567             (rs (if (> (* Lxx Lyy) 0) (/ Lxy (sqrt (* Lxx Lyy))) 0)) ; TODO: think about rs = 1
1568             (ts (if (= 1 rs) 1 (/ (* rs (sqrt (- n 2))) (sqrt (- 1 (square rs))))))
1569             (p (t-significance ts (- n 2) #:tails ':both)))
1570        (format #t "~%Spearman correlation coefficient ~A, p = ~A~%" rs p)
1571        (values rs p))))
1572
1573  ;; ---------------------------------------------------------------------------
1574  ;; Significance test functions
1575  ;;
1576
1577  (define (t-significance t-statistic dof . args)
1578    (set! t-statistic (exact->inexact t-statistic))
1579    (set! dof (exact->inexact dof))
1580    (let ((tails (get-keyword #:tails args (lambda () ':both))))
1581      (let ((a (beta-incomplete (* 0.5 dof) 0.5 (/ dof (+ dof (square t-statistic))))))
1582        (case tails
1583          ((:both) 
1584           a)
1585          ((:positive)
1586           (if (positive? t-statistic)
1587             (* 0.5 a)
1588             (- 1.0 (* 0.5 a))))
1589          ((:negative)
1590           (if (positive? t-statistic)
1591             (- 1.0 (* 0.5 a))
1592             (* 0.5 a)))))))
1593
1594  (define (f-significance f-statistic numerator-dof denominator-dof . args)
1595    (let ((one-tailed? (get-keyword #:one-tailed? args (lambda () #f))))
1596      (let ((tail-area (beta-incomplete (* 0.5 denominator-dof)
1597                                        (* 0.5 numerator-dof)
1598                                        (/ denominator-dof
1599                                           (+ denominator-dof 
1600                                              (* numerator-dof f-statistic))))))
1601        (if one-tailed? 
1602          (if (< f-statistic 1)
1603            (- 1 tail-area)
1604            tail-area)
1605          (begin (set! tail-area (* 2.0 tail-area))
1606                 (if (> tail-area 1.0)
1607                   (- 2.0 tail-area)
1608                   tail-area))))))
1609
1610  (rng-init)
1611 
1612  ) ; end of library
1613
Note: See TracBrowser for help on using the repository browser.