source: project/release/4/statistics/trunk/statistics.scm @ 33028

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

statistics: added gsl bindings for lambert's w function; added maintainer to meta file

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