source: project/release/5/statistics/branches/collections/statistics.scm @ 36371

Last change on this file since 36371 was 36371, checked in by Ivan Raikov, 13 months ago

statistics: yasos collections branch

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