source: project/release/4/statistics/tags/0.6/statistics.scm @ 32986

Last change on this file since 32986 was 32986, checked in by Ivan Raikov, 6 years ago

statistics v 0.6 release

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