lh-statistics.lisp

;;; Statistical functions in Common Lisp.  Version 1.05 Aug 14, 2022
;;;
;;; This code is copyright (c) 2000, 2001, 2002, 2005, 2022 by Lawrence Hunter
;;; (Larry.Hunter@CUAnschutz.edu) except where otherwise noted.

;;; Thanks to Paul Cohen for the original CLASP package.  Thanks also to bug
;;; reports from Rob St. Amant and Lee Ayres, and several bug fixes from
;;; Robert Goldman.
;;;
;;; This program is free software; you can redistribute it and/or modify it
;;; under the terms of the MIT License (see LICENSE file).

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Functions provided:
;;;
;;;  Descriptive statistics
;;;   Mean, median, mode, geometric mean, range, percentile, variance,
;;;   standard-deviation (sd), coefficient-of-variation,
;;;   standard-error-of-the-mean
;;;
;;;  Distributional functions
;;;   Poisson & Binomial
;;;    binomial-probability, binomial-cumulative-probability,
;;;    binomial-ge-probability, poisson-probability,
;;;    poisson-cumulative-probability, poisson-ge-probability, Normal
;;;    normal-pdf, convert-to-standard-normal, phi, z, t-distribution,
;;;    chi-square, chi-square-cdf
;;;
;;;  Confidence Intervals

;;;   binomial-probability-ci, poisson-mu-ci, normal-mean-ci,
;;;   normal-mean-ci-on-sequences, normal-variance-ci,
;;;   normal-variance-ci-on-sequence, normal-sd-ci
;;;
;;;  Hypothesis tests (parametric)
;;;   z-test, z-test-on-sequence, t-test-one-sample,
;;;   t-test-one-sample-on-sequence, t-test-paired,
;;;   t-test-paired-on-sequences, t-test-two-sample,
;;;   t-test-two-sample-on-sequences, chi-square-test-one-sample, f-test,
;;;   binomial-test-one-sample, binomial-test-two-sample, fisher-exact-test,
;;;   mcnemars-test, poisson-test-one-sample
;;;
;;;  Hypothesis tests (non-parametric)
;;;   sign-test, sign-test-on-sequence, wilcoxon-signed-rank-test,
;;;   chi-square-test-rxc, chi-square-test-for-trend

;;;  Sample size estimates
;;;   t-test-one-sample-sse, t-test-two-sample-sse
;;;   t-test-paired-sse, binomial-test-one-sample-sse,
;;;   binomial-test-two-sample-sse,
;;;   binomial-test-paired-sse, correlation-sse

;;;  Correlation and Regression
;;;   linear-regression, correlation-coefficient,
;;;   correlation-test-two-sample, spearman-rank-correlation

;;;  Significance test functions
;;;   t-significance, f-significance (chi square significance is calculated
;;;   from chi-square-cdf in various ways depending on the problem).

;;;  Utilities
;;;   random-pick, bin-and-count, fishers-z-transform,
;;;   mean-sd-n, square, choose, permutations, round-float


(declaim (optimize (speed 3) (safety 1) (debug 1)))

(uiop:define-package :lh.statistics
  (:use :common-lisp)
  (:export #:mean #:median #:mode #:geometric-mean #:sample-range #:percentile
           #:variance #:standard-deviation #:sd #:coefficient-of-variation
           #:standard-error-of-the-mean #:permutations #:choose
           #:binomial-probability #:binomial-cumulative-probability
           #:binomial-ge-probability #:poisson-probability
           #:poisson-cumulative-probability #:poisson-ge-probability
           #:poisson-cdf #:normal-pdf #:convert-to-standard-normal #:phi #:z
           #:t-distribution #:chi-square #:chi-square-cdf
           #:binomial-probability-ci #:poisson-mu-ci #:normal-mean-ci
           #:normal-mean-ci-on-sequence #:normal-variance-ci
           #:normal-variance-ci-on-sequence #:normal-sd-ci
           #:normal-sd-ci-on-sequence #:z-test #:z-test-on-sequence
           #:t-test-one-sample #:t-test-one-sample-on-sequence
           #:t-test-paired #:t-test-paired-on-sequences #:t-test-two-sample
           #:t-test-two-sample-on-sequences #:chi-square-test-one-sample
           #:f-test #:binomial-test-one-sample #:binomial-test-two-sample
           #:fisher-exact-test #:mcnemars-test #:poisson-test-one-sample
           #:sign-test #:sign-test-on-sequences #:wilcoxon-signed-rank-test
           #:wilcoxon-signed-rank-test-on-sequences
           #:chi-square-test-rxc #:chi-square-test-for-trend
           #:t-test-one-sample-sse #:t-test-two-sample-sse
           #:t-test-paired-sse #:binomial-test-one-sample-sse
           #:binomial-test-two-sample-sse #:binomial-test-paired-sse
           #:correlation-sse #:linear-regression #:correlation-coefficient
           #:correlation-test-two-sample
           #:correlation-test-two-sample-on-sequences #:spearman-rank-correlation
           #:t-significance #:f-significance #:random-pick #:test-variables
           #:bin-and-count #:fisher-z-transform #:mean-sd-n #:square
           #:round-float #:false-discovery-correction
           #:random-normal)
  (:documentation "The formulas and methods used are largely taken from Bernard Rosner, *Fundamentals of Biostatistics* 5th Edition.  'Rosner x' is a page number.  Some numeric functions were taken from CLASP, a 1994 common lisp package that implemented some of the statistical functions from *Numeric recipes in C* For CLASP functions, see copyright notice below.

These abreviations used in function and variable names:
   ci = confidence interval
   cdf = cumulative density function
   ge = greater than or equal to
   le = less than or equal to
   pdf = probability density function
   sd = standard deviation
   rxc = rows by columns
   sse = sample size estimate"))


(in-package :lh.statistics)

;;;;; Macros

;; This macro makes assertions more readable.  There are several special
;; types defined: :probability (:prob), :positive-integer (:posint),
;; :positive-number (:posnum), :number-sequence (:numseq),
;; :positive-integer-sequence (:posintseq), :probability-sequence
;; (:probseq), :nonzero-number-sequence (:nonzero-numseq) and :percentage
;; (:percent).  Other assertions are assumed to be internal types.  The
;; arguments to test-variables are lists.  The first element of the list is
;; a variable name, and the second element is either a special or built-in
;; type.  If the variable binding is not of the type specified, and error is
;; signalled indicating the problem.  One variable may have multiple type
;; requirements, which are conjunctive.

(eval-when (:compile-toplevel :load-toplevel :execute)

(defmacro test-variables (&rest args)
  (let ((assertions nil))
    (dolist (arg args (append `(or ,@(nreverse assertions))))
      (let* ((name (first arg))
             (type (second arg))
             (test (case type
                     ((:probability :prob)
                      `(and (numberp ,name) (not (minusp ,name)) (<= ,name 1)))
                     ((:positive-integer :posint)
                      `(and (integerp ,name) (plusp ,name)))
                     ((:positive-number :posnum)
                      `(and (numberp ,name) (plusp ,name)))
                     ((:number-sequence :numseq)
                      `(and (typep ,name 'sequence) (every #'numberp ,name)
                            (not (null ,name))))
                     ((:nonzero-number-sequence :nonzero-numseq)
                      `(and (typep ,name 'sequence) (not (null ,name))
                        (every #'(lambda (x) (and (numberp x) (not (= 0 x))))
                         ,name)))
                     ((:probability-sequence :probseq)
                      `(and (typep ,name 'sequence) (not (null ,name))
                        (every #'(lambda (x) (and (numberp x) (not (minusp x))
                                                  (<= x 1.0))) ,name)))
                     ((:positive-integer-sequence :posintseq)
                      `(and (typep ,name 'sequence) (not (null ,name))
                        (every #'(lambda (x) (and (typep x 'integer) (plusp
                                                                      x)))
                         ,name)))
                     (:percentage
                      `(and (numberp ,name) (plusp ,name) (<= ,name 100)))
                     (:test (third arg))
                     (t `(typep ,name ',type))))
             (message `(error
                        ,(if (eql type :test)
                             name
                             (format nil "~a = ~~a is not a ~a" name
                                     (case type
                                       ((:positive-integer :posint)
                                        "positive integer")
                                       ((:positive-number :posnum)
                                        "positive number")
                                       ((:probability :prob) "probability")
                                       ((:number-sequence :numseq)
                                        "sequence of numbers")
                                       ((:nonzero-number-sequence
                                         :nonzero-numseq)
                                        "sequence of non-zero numbers")
                                       ((:positive-integer-sequence :posintseq)
                                        "sequence of positive integers")
                                       ((:probability-sequence :probseq)
                                        "sequence of probabilities")
                                       ((:percent :percentile) "percent")
                                       (t type))))
                        ,@(unless (eql type :test) `(,name)))))
        (push `(unless ,test ,message) assertions)))))

(defmacro square (x)
  `(* ,x ,x))

(defmacro underflow-goes-to-zero (&body body)
  "Protects against floating point underflow errors and sets the value to 0.0 instead."
  `(handler-case
       (progn ,@body)
       (floating-point-underflow (condition)
	(declare (ignore condition))
	(values 0.0d0))))


) ;end eval-when



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;;  Descriptive statistics
;;;

;;; Measures of location

;; Rosner 10
(defun mean (sequence)
  "Returns the mean of SEQUENCE"
  (test-variables (sequence :numseq))
  (/ (reduce #'+ sequence) (length sequence)))

;; Rosner 12 (and 19)
(defun median (sequence)
  "Returns the median of SEQUENCE"
  (test-variables (sequence :numseq))
  (percentile sequence 50))

;; Rosner 14
(defun mode (sequence)
  "Returns two values: a list of the modes and the number of times they occur"

  ;; Rob St. Amant <stamant@csc.ncsu.edu> suggested using a hash
  ;; table instead of an alist, and Lee Ayres <ayres@acm.org> noted that
  ;; my revision failed to handle multiple modes properly.

  (test-variables (sequence :numseq))
  (let ((count-table (make-hash-table :test #'eql))
	(modes nil)
	(mode-count 0))
    (map nil (lambda (elt) (incf (gethash elt count-table 0))) sequence)
    (maphash (lambda (key value)
               (cond ((> value mode-count)
		      (setf modes (list key)
			    mode-count value))
		     ((= value mode-count)
		      (push key modes))))
             count-table)
    (values modes mode-count)))

;; Rosner 16
(defun geometric-mean (sequence &optional (base 10))
  "Returns the geometric mean of SEQUENCE
The geometric mean is a mean or average, which indicates the central tendency or typical value of a set of numbers by using the product of their values (as opposed to the arithmetic mean which uses their sum)"
  (test-variables (sequence :nonzero-numseq) (base :posnum))
  (expt base (mean (map 'list #'(lambda (x) (log x base)) sequence))))



;;; Measures of spread

;; Rosner 18
(defun sample-range (sequence)
  "Return the difference between the largest and smallest values in SEQUENCE"
  (test-variables (sequence :numseq))
  (- (reduce #'max sequence) (reduce #'min sequence)))

;; Rosner 19
(defun percentile (sequence percent)
  "Return an element from SEQUENCE at percentile PERCENT
This function is also known as quantile."
  (test-variables (sequence :numseq) (percent :percentage))
  (let* ((sorted-vect (coerce (sort (copy-seq sequence) #'<) 'simple-vector))
         (n (length sorted-vect))
         (k (* n (/ percent 100)))
         (floor-k (floor k)))
    (if (= k floor-k)
        (/ (+ (aref sorted-vect k)
              (aref sorted-vect (1- k)))
           2)
        (aref sorted-vect floor-k))))

;; Rosner 21
(defun variance (sequence)
  "Return variance of SEQUENCE"
  (test-variables (sequence :numseq))
  (let ((mean (mean sequence))
        (n (length sequence)))
    (/ (reduce #'+ (map 'list #'(lambda (x) (square (- mean x))) sequence))
       (1- n))))

;; Rosner 21
(defun standard-deviation (sequence)
  "Return the standard deviation of SEQUENCE"
  (test-variables (sequence :numseq))
  (let ((mean (mean sequence))
        (n (length sequence)))
    (sqrt (/ (reduce #'+ (map 'list #'(lambda (x) (square (- mean x))) sequence))
             (1- n)))))

(defun sd (sequence)
  (test-variables (sequence :numseq))
  (let ((mean (mean sequence))
        (n (length sequence)))
    (sqrt (/ (reduce #'+ (map 'list #'(lambda (x) (square (- mean x))) sequence))
             (1- n)))))


;;; The coefficient of variation (CV) is a standardized measure of
;;; dispersion of a probability distribution or frequency
;;; distribution.

;; Rosner 24
(defun coefficient-of-variation (sequence)
  "Return coefficient of variation"
  (* 100 (/ (standard-deviation sequence) (mean sequence))))

;; Rosner 172
(defun standard-error-of-the-mean (sequence)
  "Return the estimated standard deviation obtained from a set of sample means from repeated samples"
  (/ (standard-deviation sequence) (sqrt (length sequence))))




;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;;  Distributional functions
;;;

;;;
;;; Binomial and Poisson distributions
;;;

;; exact: Rosner 93, approximate 105
(defun binomial-probability (n k p)
  "Return P(X=k) for X a binomial random variable with parameters n & p.  Binomial expectations for seeing k events in N trials, each having probability p.  Use the Poisson approximation if N>100 and P<0.01."
  (test-variables (n :posint) (p :prob)
                  ("K must be between 0 and N (inclusive)" :test (and (>= k 0) (<= k n))))
  (if (and (> n 100) (< p 0.01))
      (poisson-probability (* n p) k)
      (let ((p (coerce p 'double-float)))
        (* (choose n k)
           (expt p k)
           (expt (- 1 p) (- n k))))))

;; Rosner 94
(defun binomial-cumulative-probability (n k p)
  "Return P(X<k) for X a binomial random variable with parameters n & p. Bionomial expecations for fewer than k events in N trials, each having probability p.  This is also known as probability mass function (PMF), the probability of getting exactly k successes in n independent Bernoulli trials."
  (test-variables (n :posint) (k :posint) (p :prob)
             ("K must be less than or equal to N" :test (<= k n)))
  (let ((sum-up-to-k-1 0d0))
    (dotimes (i k sum-up-to-k-1)
      (incf sum-up-to-k-1 (binomial-probability n i p)))))

;; Rosner 94
(defun binomial-ge-probability (n k p)
  "The probability of k or more occurances in N events, each with probability p."
  (- 1 (binomial-cumulative-probability n k p)))

;; Just for convenience later, binomial-le-probability
(defun binomial-le-probability (n k p)
  (test-variables (n :posint) (k :posint) (p :prob)
             ("K must be less than or equal to N" :test (<= k n)))
  (let ((sum-up-to-k 0d0))
    (dotimes (i (1+ k) sum-up-to-k)
      (incf sum-up-to-k (binomial-probability n i p)))))

;; Rosner 100
(defun poisson-probability (mu k)
  "Probability of seeing k events over a time period when the expected number of events over that time is mu."
  (test-variables (mu :posnum) (k :posint))
  (let ((mu (coerce mu 'double-float)))
    (/ (* (exp (- mu)) (expt mu k))
       (factorial k))))

;; Rosner 197
(defun poisson-cumulative-probability (mu k)
  "Probability of seeing fewer than K events over a time period when the expected number events over that time is mu."
  (test-variables (mu :posnum) (k :posint))
  (if (< k 170)
      (let ((sum 0d0))
        (dotimes (x k sum)
          (incf sum (poisson-probability mu x))))
      (let ((mu (coerce mu 'double-float))
            (k (coerce k 'double-float)))
        (- 1d0 (gamma-incomplete k mu)))))

(defun poisson-ge-probability (mu x)
  "Probability of X or more events when expected is mu."
  (- 1 (poisson-cumulative-probability mu x)))



;;;  Normal distribution functions

;; Rosner 115
(defun normal-pdf (x mu sigma)
  "The probability density function (PDF) for a normal distribution with mean mu and variance sigma at point x."
  (test-variables (x number) (mu number) (sigma :posnum))
  (* (/ (* (sqrt (* 2 pi)) sigma))
     (exp (* (- (/ (* 2 (square sigma)))) (square (- x mu))))))

;; Rosner 130
(defun convert-to-standard-normal (x mu sigma)
  "Convert X from a Normal distribution with mean mu and variance sigma to standard normal"
  (test-variables (x number) (mu number) (sigma :posnum))
  (/ (- x mu) sigma))

;; Rosner 125
;; Adopted from CLASP 1.4.3, see copyright notice at http://eksl-www.cs.umass.edu/clasp.html"
(defun phi (x)
  "the CDF of standard normal distribution"
  (test-variables (x number))
  (setf x (coerce x 'double-float))
  (locally (declare (type double-float x))
    (* 0.5d0 (+ 1.0d0 (error-function (/ x (sqrt 2.0d0)))))))

;; Rosner 128.
(defun z (percentile &key (epsilon 1d-15))
  "The inverse normal function, P(X<Zu) = u where X is distributed as the standard normal.  Uses binary search.("
  (test-variables (percentile :prob))
  (let ((target (coerce percentile 'double-float)))
    (do ((min -9d0 min)
         (max 9d0 max)
         (guess 0d0 (+ min (/ (- max min) 2d0))))
        ((< (- max min) epsilon) guess)
      (let ((result (coerce (phi guess) 'double-float)))
        (if (< result target)
            (setq min guess)
            (setq max guess))))))

;; Rosner 178
;;Adopted from CLASP 1.4.3, http://eksl-www.cs.umass.edu/clasp.html"
(defun t-distribution (dof percentile)
  "Returns the point which is the indicated percentile in the T distribution with dof degrees of freedom"
  (test-variables (dof :posint) (percentile :prob))
  (find-critical-value
   #'(lambda (x) (t-significance x dof :tails :positive))
   (- 1 percentile)))

;; Rosner 187
(defun chi-square (dof percentile)
  ";; Returns the point which is the indicated percentile in the Chi Square distribution with dof degrees of freedom."
  (test-variables (dof :posint) (percentile :prob))
  (find-critical-value #'(lambda (x) (chi-square-cdf x dof))
                       (- 1 percentile)))

;; Adopted from CLASP 1.4.3, http://eksl-www.cs.umass.edu/clasp.html"
(defun chi-square-cdf (x dof)
  " Chi-square-cdf computes the left hand tail area under the chi square distribution under dof degrees of freedom up to X."
  (test-variables (x :posnum) (dof :posint))
  (multiple-value-bind (cdf ignore)
      (gamma-incomplete (* 0.5 dof) (* 0.5 x))
    (declare (ignore ignore))
    cdf))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;;  Confidence intervals
;;;

;; Rosner 193 (approximate) & 194 (exact)
(defun binomial-probability-ci (n p alpha &key exact?)
  "Confidence intervals on a binomial probability.  If a binomial probability of p has been observed in N trials, what is the 1-alpha confidence interval around p?  Approximate (using normal theory approximation) when npq >= 10 unless told otherwise"
  (test-variables (n :posint) (p :prob) (alpha :prob))
  (if (and (> (* n p (- 1 p)) 10) (not exact?))
      (let ((difference (* (z (- 1 (/ alpha 2)))
                           (sqrt (/ (* p (- 1 p)) n)))))
        (values (- p difference) (+ p difference)))
      (values (find-critical-value
               (lambda (p1) (binomial-cumulative-probability n (floor (* p n)) p1))
               (- 1 (/ alpha 2)))
              (find-critical-value
               (lambda (p2) (binomial-cumulative-probability n (1+ (floor (* p n))) p2))
               (/ alpha 2)))))

;; Rosner 197
(defun poisson-mu-ci (x alpha)
  "Confidence interval for the Poisson parameter mu

Given x observations in a unit of time, what is the 1-alpha confidence interval on the Poisson parameter mu (= lambda*T)?

Since find-critical-value assumes that the function is monotonic increasing, adjust the value we are looking for taking advantage of reflectiveness"
  (test-variables (x :posnum) (alpha :prob))
  (values
   (find-critical-value
    #'(lambda (mu) (poisson-cumulative-probability mu (1- x)))
    (- 1 (/ alpha 2)))
   (find-critical-value
    #'(lambda (mu) (poisson-cumulative-probability mu x))
    (/ alpha 2))))

;; Rosner 180
(defun normal-mean-ci (mean sd n alpha)
  "Confidence interval for the mean of a normal distribution.

The 1-alpha percent confidence interval on the mean of a normal distribution with parameters mean, sd & n."
  (test-variables (mean number) (sd :posnum) (n :posint) (alpha :prob))
  (let ((t-value (t-distribution (1- n) (- 1 (/ alpha 2)))))
    (values (- mean (* t-value (/ sd (sqrt n))))
            (+ mean (* t-value (/ sd (sqrt n)))))))

(defun normal-mean-ci-on-sequence (sequence alpha)
  "The 1-alpha confidence interval on the mean of a sequence of numbers drawn from a Normal distribution."
  (test-variables (alpha :prob)) ; sequence tested by mean-sd-n
  (multiple-value-call #'normal-mean-ci (mean-sd-n sequence) alpha))

;; Rosner 190
(defun normal-variance-ci (variance n alpha)
  "The 1-alpha confidence interval on the variance of a sequence of numbers drawn from a Normal distribution."
  (test-variables (variance :posnum) (n :posint) (alpha :prob))
  (let ((chi-square-low (chi-square (1- n) (- 1 (/ alpha 2))))
        (chi-square-high (chi-square (1- n) (/ alpha 2)))
        (numerator (* (1- n) variance)))
    (values (/ numerator chi-square-low) (/ numerator chi-square-high))))

(defun normal-variance-ci-on-sequence (sequence alpha)
  (test-variables (sequence :numseq) (alpha :prob))
  (let ((variance (variance sequence))
        (n (length sequence)))
    (normal-variance-ci variance n alpha)))

;; Rosner 190
(defun normal-sd-ci (sd n alpha)
  "The 1-alpha confidence interval on the standard deviation of a sequence of numbers drawn from a Normal distribution."
  (multiple-value-bind (low high)
      (normal-variance-ci (square sd) n alpha)
    (values (sqrt low) (sqrt high))))

(defun normal-sd-ci-on-sequence (sequence alpha)
  (test-variables (sequence :numseq) (alpha :prob))
  (let ((sd (standard-deviation sequence))
        (n (length sequence)))
    (normal-sd-ci sd n alpha)))



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Hypothesis testing
;;;

;; Rosner 228
(defun z-test (x-bar n &key (mu 0) (sigma 1) (tails :both))
  "The significance of a one sample Z test for the mean of a normal distribution with known variance.  mu is the null hypothesis mean, x-bar is the observed mean, sigma is the standard deviation and N is the number of observations.  If tails is :both, the significance of a difference between x-bar and mu.  If tails is :positive, the significance of x-bar is greater than mu, and if tails is :negative, the significance of x-bar being less than mu.  Returns a p value."
  (test-variables (x-bar number) (n :posint) (mu number) (sigma :posnum))
  (let ((z (/ (- x-bar mu) (/ sigma (sqrt n)))))
    (ecase tails
      (:both (if (< z 0)
                 (* 2 (phi z))
                 (* 2 (- 1 (phi z)))))
      (:negative (phi z))
      (:positive (- 1 (phi z))))))

(defun z-test-on-sequence (sequence &key (mu 0) (sigma 1) (tails :both))
  (test-variables (sequence :numseq)) ; the rest handled by z-test
  (let ((x-bar (mean sequence))
        (n (length sequence)))
    (z-test x-bar n :mu mu :sigma sigma :tails tails)))

;; Rosner 216
(defun t-test-one-sample (x-bar sd n mu &key (tails :both))
  "The significance of a one sample T test for the mean of a normal distribution with unknown variance.  X-bar is the observed mean, sd is the observed standard deviation, N is the number of observations and mu is the test mean."
  (test-variables (x-bar number) (sd :posnum) (n :posint) (mu number))
  (t-significance  (/ (- x-bar mu) (/ sd (sqrt n))) (1- n) :tails tails))

;; Rosner 216
(defun t-test-one-sample-on-sequence (sequence mu &key (tails :both))
  "The significance of a one sample T test for the mean of a normal sequence of numbers with unknown variance.  X-bar is the observed mean, sd is the observed standard deviation, N is the number of observations and mu is the test mean."
  (multiple-value-call #'t-test-one-sample
    (mean-sd-n sequence) mu :tails tails))

;; Rosner 276
(defun t-test-paired (d-bar sd n &key (tails :both))
  "The significance of a paired t test for the means of two normal distributions in a longitudinal study.  D-bar is the mean difference, sd is the standard deviation of the differences, N is the number of pairs."
  (test-variables (d-bar number) (sd :posnum) (n :posint))
  (t-significance (/ d-bar (/ sd (sqrt n))) (1- n) :tails tails))

;; Rosner 276
(defun t-test-paired-on-sequences (before after &key (tails :both))
  "The significance of a paired t test for means of two normal distributions in a longitudinal study.  Before is a sequence of before values, after is the sequence of paired after values (which must be the same length as the before sequence)."
  (test-variables (before :numseq) (after :numseq)
             ("Before and after sequences must be of equal length"
              :test (= (length before) (length after))))
  (multiple-value-call #'t-test-paired
     (mean-sd-n (map 'list #'- before after)) :tails tails))

;; Rosner  282, 294, 297
(defun t-test-two-sample (x-bar1 sd1 n1 x-bar2 sd2 n2 &key
                          (variances-equal? :test)
                          (variance-significance-cutoff 0.05)
							(tails :both))
  "The significance of the difference of two means (x-bar1 and x-bar2) with standard deviations sd1 and sd2, and sample sizes n1 and n2 respectively. The form of the two sample t test depends on whether the sample variances are equal or not.   If the variable variances-equal? is :test, then we use an F test and the variance-significance-cutoff to determine if they are equal.  If the variances are equal, then we use the two sample t test for equal variances.  If they are not equal, we use the Satterthwaite method, which has good type I error properties (at the loss of some power)."
  (test-variables (x-bar1 number) (sd1 :posnum) (n1 :posint)
             (x-bar2 number) (sd2 :posnum) (n2 :posint))
  (let (t-statistic dof)
    (if (ecase variances-equal?
          (:test (> (f-test (square sd1) n1 (square sd2) n2 :tails tails)
                    variance-significance-cutoff))
          ((t :yes :equal) t)
          ((nil :no :unequal) nil))
        (let ((s (sqrt (/ (+ (* (1- n1) (square sd1))
                             (* (1- n2) (square sd2)))
                          (+ n1 n2 -2)))))
          (setq t-statistic  (/ (- x-bar1 x-bar2)
                                (* s (sqrt (+ (/ n1) (/ n2)))))
                dof (+ n1 n2 -2)))
        (let* ((variance-ratio1 (/ (square sd1) n1))
               (variance-ratio2 (/ (square sd2) n2)))
          (setq t-statistic (/ (- x-bar1 x-bar2)
                               (sqrt (+ variance-ratio1 variance-ratio2)))
                dof (round (/ (square (+ variance-ratio1 variance-ratio2))
                              (+ (/ (square variance-ratio1) (1- n1))
                                 (/ (square variance-ratio2) (1- n2))))))))
    (t-significance t-statistic dof :tails tails)))


;; Same as above, but providing the sequences rather than the summaries.
(defun t-test-two-sample-on-sequences (sequence1 sequence2 &key
                                       (variance-significance-cutoff 0.05)
							     (tails :both))
  "The significance of the difference of two means of SEQUENCE1 and SEQUENCE2 with standard deviations sd1 and sd2, and sample sizes n1 and n2 respectively. The form of the two sample t test depends on whether the sample variances are equal or not.   If the variable variances-equal? is :test, then we use an F test and the variance-significance-cutoff to determine if they are equal.  If the variances are equal, then we use the two sample t test for equal variances.  If they are not equal, we use the Satterthwaite method."
  (multiple-value-call #'t-test-two-sample
    (mean-sd-n sequence1) (mean-sd-n sequence2) :tails tails
    :variance-significance-cutoff variance-significance-cutoff))

;; Rosner 290
(defun f-test (variance1 n1 variance2 n2 &key (tails :both))
  "F test for the equality of two variances"
  (test-variables (variance1 :posnum) (n1 :posint) (variance2 :posnum) (n2 :posint))
  (let ((significance (f-significance (/ variance1 variance2) (1- n1) (1- n2)
                                      (not (eql tails :both)))))
    (ecase tails
      (:both significance)
      (:positive (if (> variance1 variance2) significance (- 1 significance)))
      (:negative (if (< variance1 variance2) significance (- 1 significance))))))

;; Rosner 246
(defun chi-square-test-one-sample (variance n sigma-squared &key (tails :both))
  "The significance of a one sample Chi square test for the variance of a normal distribution.  Variance is the observed variance, N is the number of observations, and sigma-squared is the test variance."
  (test-variables (variance :posnum) (n :posint) (sigma-squared :posnum))
  (let ((cdf (chi-square-cdf (/ (* (1- n) variance) sigma-squared) (1- n))))
    (ecase tails
      (:negative cdf)
      (:positive (- 1 cdf))
      (:both (if (<= variance sigma-squared)
                 (* 2 cdf)
                 (* 2 (- 1 cdf)))))))

;; Rosner 249
(defun binomial-test-one-sample (p-hat n p &key (tails :both) (exact? nil))
  "The significance of a one sample test for the equality of an observed probability p-hat to an expected probability p under a binomial distribution with N observations.  Use the normal theory approximation if n*p*(1-p) > 10 (unless the exact flag is true)."
  (test-variables (p-hat :prob) (n :posint) (p :prob))
  (let ((q (- 1 p)))
    (if (and (> (* n p q) 10) (not exact?))
        (let ((z (/ (- p-hat p) (sqrt (/ (* p q) n)))))
          (ecase tails
            (:negative (phi z))
            (:positive (- 1 (phi z)))
            (:both (* 2 (if (<= p-hat p) (phi z) (- 1 (phi z)))))))
        (let* ((observed (round (* p-hat n)))
               (probability-more-extreme
                (if (<= p-hat p)
                    (binomial-cumulative-probability n observed p)
                    (binomial-ge-probability n observed p))))
          (ecase tails
            ((:negative :positive) probability-more-extreme)
            (:both (min (* 2 probability-more-extreme) 1.0)))))))

;; Rosner 357
(defun binomial-test-two-sample (p-hat1 n1 p-hat2 n2 &key (tails :both)
						       (exact? nil))
  "Are the observed probabilities of an event (p-hat1 and p-hat2) in N1/N2 trials different? The normal theory method implemented here.  The exact test is Fisher's contingency table method, below."
  (test-variables (p-hat1 :prob) (n1 :posint) (p-hat2 :prob) (n2 :posint))
  (let* ((p-hat (/ (+ (* p-hat1 n1) (* p-hat2 n2)) (+ n1 n2)))
         (q-hat (- 1 p-hat))
         (z (/ (- (abs (- p-hat1 p-hat2)) (+ (/ (* 2 n1)) (/ (* 2 n2))))
               (sqrt (* p-hat q-hat (+ (/ n1) (/ n2)))))))
    (if (and (> (* n1 p-hat q-hat) 5)
             (> (* n2 p-hat q-hat) 5)
             (not exact?))
        (* (- 1 (phi z)) (if (eql tails :both) 2 1))
        (let ((contingency-table (make-array '(2 2))))
          (setf (aref contingency-table 0 0) (* p-hat1 n1)
                (aref contingency-table 0 1) (- 1 (* p-hat1 n1))
                (aref contingency-table 1 0) (* p-hat2 n2)
                (aref contingency-table 1 1) (- 1 (* p-hat2 n2)))
          (fisher-exact-test contingency-table :tails tails)))))

;; Rosner 371
(defun fisher-exact-test (contingency-table &key (tails :both))
  "Fisher's exact test.  Gives a p value for a particular 2x2 contingency table"
  (flet ((table-probability (a b c d)
           (let ((n (+ a b c d)))
             (/ (* (factorial (+ a b)) (factorial (+ c d))
                   (factorial (+ a c)) (factorial (+ b d)))
                (* (factorial n) (factorial a) (factorial b)
                   (factorial c) (factorial d))))))

    (let ((a (aref contingency-table 0 0))
          (b (aref contingency-table 0 1))
          (c (aref contingency-table 1 0))
          (d (aref contingency-table 1 1)))
      (test-variables (a number) (b number) (c number) (d number))
      (let* ((row-margin1 (+ a b))
             (row-margin2 (+ c d))
             (column-margin1 (+ a c))
             (column-margin2 (+ b d))
             (n (+ a b c d))
             (table-probabilities
              (make-array (1+ (min row-margin1 row-margin2 column-margin1
                                   column-margin2)))))

        ;; rearrange so that the first row and column marginals are
        ;; smallest.  Only need to change first margins and a.

        (cond ((and (< row-margin2 row-margin1) (< column-margin2 column-margin1))
               (psetq a d
                      row-margin1 row-margin2
                      column-margin1 column-margin2))
              ((< row-margin2 row-margin1)
               (psetq a c
                      row-margin1 row-margin2))
              ((< column-margin2 column-margin1)
               (psetq a b
                      column-margin1 column-margin2)))
        (dotimes (i (length table-probabilities))
          (let* ((test-a i)
                 (test-b (- row-margin1 i))
                 (test-c (- column-margin1 i))
                 (test-d (- n (+ test-a test-b test-c))))
            (setf (aref table-probabilities i)
                  (table-probability test-a test-b test-c test-d))))
        (let ((above (reduce #'+ (subseq table-probabilities 0 (1+ a))))
              (below (reduce #'+ (subseq table-probabilities a))))
          (float
           (ecase tails
             ((:both) (* 2 (min above below)))
             ((:positive) below)
             ((:negative) above))
           1d0))))))

;; Rosner 379 and 381
(defun mcnemars-test (a-discordant-count b-discordant-count &key (exact? nil))
  "McNemar's test for correlated proportions, used for longitudinal studies. Look only at the number of discordant pairs (one treatment is effective and the other is not).  If the two treatments are A and B, a-discordant-count is the number where A worked and B did not, and b-discordant-count is the number where B worked and A did not."
  (test-variables (a-discordant-count :posint) (b-discordant-count :posint))
  (let ((n (+ a-discordant-count b-discordant-count)))
    (if (and (> n 20) (not exact?))
        (let ((x2 (/ (square
                      (- (abs (- a-discordant-count b-discordant-count)) 1))
                     n)))
          (- 1 (chi-square-cdf x2 1)))
        (cond ((= a-discordant-count b-discordant-count) 1.0)
              ((< a-discordant-count b-discordant-count)
               (* 2 (binomial-le-probability n a-discordant-count 1/2)))
              (t (* 2 (binomial-ge-probability n a-discordant-count 1/2)))))))

;; Rosner 256 (approximation on 259)
(defun poisson-test-one-sample (observed mu &key (tails :both) (approximate? nil))
  "The significance of a one sample test for the equality of an observed number of events (observed) and an expected number mu under the poisson distribution.  Normal theory approximation is not that great, so don't use it unless told."
  (test-variables (observed :posnum) (mu :posnum))
  (if approximate?
      (let ((x-square (/ (square (- observed mu)) mu)))
        (- 1 (chi-square-cdf x-square 1)))
      (let ((probability-more-extreme
             (if (< observed mu)
                 (poisson-cumulative-probability mu observed)
                 (poisson-ge-probability mu observed))))
        (ecase tails
          ((:negative :positive) probability-more-extreme)
          (:both (min (* 2 probability-more-extreme) 1.0))))))


;;; Non-parametric hypothesis testing

;; Rosner 335-7.
(defun sign-test (plus-count minus-count &key (exact? nil) (tails :both))
  "Really just a special case of the binomial one sample test with p = 1/2. The normal theory version has a correction factor to make it a better approximation."
  (test-variables (plus-count :posint) (minus-count :posint))
  (let* ((n (+ plus-count minus-count))
         (p-hat (/ plus-count n)))
    (if (or (< n 20) exact?)
        (binomial-test-one-sample p-hat n 0.5 :tails tails :exact? t)
        (let ((area (- 1 (phi (/ (1- (abs (- plus-count minus-count)))
                                 (sqrt n))))))
          (if (eql tails :both)
              (* 2 area)
              area)))))

(defun sign-test-on-sequences (sequence1 sequence2 &key (exact? nil) (tails :both))
  "Same as SIGN-TEST, but takes two sequences and tests whether the entries in one are different (greater or less) than the other."
  (test-variables (sequence1 :numseq) (sequence2 :numseq)
              ("Sequences must be of equal length"
               :test (= (length sequence1) (length sequence2))))
  (let* ((differences (map 'list #'- sequence1 sequence2))
         (plus-count (count #'plusp differences))
         (minus-count (count #'minusp differences)))
    (sign-test plus-count minus-count :exact? exact? :tails tails)))

;; Rosner 341
(defun wilcoxon-signed-rank-test (differences &optional (tails :both))
  "A test on the ranking of positive and negative differences (are the positive differences significantly larger/smaller than the negative ones). Assumes a continuous and symmetric distribution of differences, although not a normal one.  This is the normal theory approximation, which is only valid when N > 15. This test is equivalent to the Mann-Whitney test."
  (let* ((nonzero-differences (remove 0 differences :test #'=))
         (sorted-list (sort (mapcar #'(lambda (dif)
                                        (list (abs dif) (sign dif)))
                                    nonzero-differences)
                            #'<
                            :key #'first))
         (distinct-values (delete-duplicates (mapcar #'first sorted-list)))
         (ties nil))

    (when (< (length nonzero-differences) 16)
      (error "This Wilcoxon Signed-Rank Test (normal approximation method) requires nonzero N > 15"))

    (unless (member tails '(:positive :negative :both))
      (error "tails must be one of :positive, :negative or :both, not ~a" tails))

    ; add avg-rank to the sorted values

    (dolist (value distinct-values)
      (let ((first (position value sorted-list :key #'first))
            (last (position value sorted-list :key #'first :from-end t)))
        (if (= first last)
            (nconc (find value sorted-list :key #'first) (list (1+ first)))
            (let ((number-tied (1+ (- last first)))
                  (avg-rank (1+ (/ (+ first last) 2)))) ; +1 since 0 based
              (push number-tied ties)
              (dotimes (i number-tied)
                (nconc (nth (+ first i) sorted-list) (list avg-rank)))))))
    (setq ties (nreverse ties))
    (let* ((direction (if (eq tails :negative) -1 1))
           (r1 (reduce #'+ (mapcar #'(lambda (entry)
                                       (if (= (second entry) direction)
                                           (third entry)
                                           0))
                                   sorted-list)))
           (n (length nonzero-differences))
           (expected-r1 (/ (* n (1+ n)) 4))
           (ties-factor (if ties
                            (/ (reduce #'+ (mapcar #'(lambda (ti)
                                                       (- (* ti ti ti) ti))
                                                   ties))
                               48)
                            0))
           (var-r1 (- (/ (* n (1+ n) (1+ (* 2 n))) 24) ties-factor))
           (T-score (/ (- (abs (- r1 expected-r1)) 1/2) (sqrt var-r1))))
      (* (if (eq tails :both) 2 1) (- 1 (phi T-score))))))

(defun wilcoxon-signed-rank-test-on-sequences (sequence1 sequence2
                                               &optional (tails :both))
  (test-variables (sequence1 :numseq) (sequence2 :numseq)
              ("Sequences must be of equal length"
               :test (= (length sequence1) (length sequence2))))
  (wilcoxon-signed-rank-test (map 'list #'- sequence1 sequence2) tails))

;; Rosner 395
(defun chi-square-test-rxc (contingency-table)
  "Takes contingency-table, an RxC array, and returns the significance of the relationship between the row variable and the column variable.  Any difference in proportion will cause this test to be significant -- consider using the test for trend instead if you are looking for a consistent change."
  (let* ((rows (array-dimension contingency-table 0))
         (columns (array-dimension contingency-table 1))
         (row-marginals (make-array rows :initial-element 0.0))
         (column-marginals (make-array columns :initial-element 0.0))
         (total 0.0)
         (expected-lt-5 0)
         (expected-lt-1 0)
         (expected-values (make-array (list rows columns)
                                      :element-type 'single-float))
         (x2 0.0))
    (dotimes (i rows)
      (dotimes (j columns)
        (let ((cell (aref contingency-table i j)))
          (incf (svref row-marginals i) cell)
          (incf (svref column-marginals j) cell)
          (incf total cell))))
    (dotimes (i rows)
      (dotimes (j columns)
        (let ((expected (/ (* (aref row-marginals i) (aref column-marginals j))
                           total)))
          (when (< expected 1) (incf expected-lt-1))
          (when (< expected 5) (incf expected-lt-5))
          (setf (aref expected-values i j) expected))))
    (when (plusp expected-lt-1)
      (error "This test cannot be used when an expected value is less than one"))
    (when (> expected-lt-5 (/ (* rows columns) 5))
      (error "This test cannot be used when more than 1/5 of the expected values are less than 5."))
    (dotimes (i rows)
      (dotimes (j columns)
        (incf x2 (/ (square (- (aref contingency-table i j)
                               (aref expected-values i j)))
                    (aref expected-values i j)))))
    (- 1 (chi-square-cdf x2 (* (1- rows) (1- columns))))))

;; Rosner 398
(defun chi-square-test-for-trend (row1-counts row2-counts &optional scores)
  "This test works on a 2xk table and assesses if there is an increasing or decreasing trend.  Arguments are equal sized lists counts.  Optionally, provide a list of scores, which represent some numeric attribute of the group.  If not provided, scores are assumed to be 1 to k."
  (unless scores (setq scores (dotimes (i (length row1-counts) (nreverse scores))
                                (push (1+ i) scores))))
  (test-variables (row1-counts :posintseq) (row2-counts :posintseq) (scores :numseq)
              ("Sequences must be of equal length"
               :test (= (length row1-counts) (length row2-counts))))
  (let* ((ns (map 'list #'+ row1-counts row2-counts))
         (p-hats (map 'list #'/ row1-counts ns))
         (n (reduce #'+ ns))
         (p-bar (/ (reduce #'+ row1-counts) n))
         (q-bar (- 1 p-bar))
         (s-bar (mean scores))
         (a (reduce #'+ (mapcar (lambda (p-hat ni s)
                                  (* ni (- p-hat p-bar) (- s s-bar)))
                                p-hats ns scores)))
         (b (* p-bar q-bar (- (reduce #'+ (mapcar (lambda (ni s) (* ni (square s)))
                                                  ns scores))
                              (/ (square (reduce #'+ (mapcar (lambda (ni s) (* ni s))
                                                             ns scores)))
                                 n))))
         (x2 (/ (square a) b))
         (significance (- 1 (chi-square-cdf (float x2) 1))))
    (when (< (* p-bar q-bar n) 5)
      (error "This test is only applicable when N * p-bar * q-bar >= 5"))
    (format t "~%The trend is ~a, p = ~f"
            (if (< a 0) "decreasing" "increasing")
            significance)
    significance))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Sample size estimates
;;;

;; Rosner 238
(defun t-test-one-sample-sse (mu mu-null variance &key
						    (alpha 0.05) (1-beta .95) (tails :both))
  "Returns the number of subjects needed to test whether the mean of a normally distributed sample mu is different from a null hypothesis mean mu-null and variance variance, with alpha, 1-beta and tails as specified."
  (test-variables (mu number) (mu-null number) (variance :posnum)
              (alpha :prob) (1-beta :prob))
  (let ((z-beta (z 1-beta))
        (z-alpha (z (- 1 (if (eql tails :both) (/ alpha 2) alpha)))))
    (round-up (/ (* variance (square (+ z-beta z-alpha)))
                (square (- mu-null mu))))))

;; Rosner 308
(defun t-test-two-sample-sse (mu1 variance1 mu2 variance2 &key
                                        (sample-ratio 1) (alpha 0.05)
							    (1-beta .95) (tails :both))
  "Returns the number of subjects needed to test whether the mean mu1 of a normally distributed sample (with variance variance1) is different from a second sample with mean mu2 and variance variance2, with alpha, 1-beta and tails as specified.  It is also possible to set a sample size ratio of sample 1 to sample 2."
  (test-variables (mu1 number) (variance1 :posnum) (mu2 number)
              (variance2 :posnum) (sample-ratio :posnum) (alpha :prob)
              (1-beta :prob))
  (let* ((delta2 (square (- mu1 mu2)))
         (z-term (square (+ (z 1-beta)
                           (z (- 1 (if (eql tails :both)
                                       (/ alpha 2)
                                       alpha))))))
         (n1 (round-up (/ (* (+ variance1 (/ variance2 sample-ratio)) z-term)
                          delta2))))
    (values n1 (round-up (* sample-ratio n1)))))

;; Rosner 311
(defun t-test-paired-sse (difference-mu difference-variance
                                    &key (alpha 0.05) (1-beta 0.95)
                                      (tails :both))
  "Returns the number of subjects needed to test whether the differences with mean difference-mu and variance difference-variance, with alpha, 1-beta and tails as specified."
  (test-variables (difference-mu number) (difference-variance :posnum)
              (alpha :prob) (1-beta :prob))
  (round-up (/ (* 2 difference-variance
                 (square (+ (z 1-beta)
                            (z (- 1 (if (eql tails :both)
                                        (/ alpha 2)
                                        alpha))))))
              (square difference-mu))))

;; Rosner 254
(defun binomial-test-one-sample-sse (p-estimated p-null &key
                                               (alpha 0.05) (1-beta 0.95)
							  (tails :both))
  "Returns the number of subjects needed to test whether an observed probability is significantly different from a particular binomial null hypothesis with a significance alpha and a power 1-beta."
  (test-variables (p-estimated :prob) (p-null :prob) (alpha :prob) (1-beta :prob))
  (let ((q-null (- 1 p-null))
        (q-estimated (- 1 p-estimated)))
    (round-up
     (/ (* p-null q-null
           (square (+ (z (- 1 (if (eql tails :both) (/ alpha 2) alpha)))
                      (* (z 1-beta)
                         (sqrt (/ (* p-estimated q-estimated)
                                  (* p-null q-null)))))))
        (square (- p-estimated p-null))))))

;; Rosner 384
(defun binomial-test-two-sample-sse (p1 p2 &key (alpha 0.05)
                                               (sample-ratio 1)
                                             (1-beta .95) (tails :both))
  "The number of subjects needed to test if two binomial probabilities are different at a given significance alpha and power 1-beta.  The sample sizes can be unequal; the p2 sample is sample-sse-ratio * the size of the p1 sample.  It can be a one tailed or two tailed test."
  (test-variables (p1 :prob) (p2 :prob) (alpha :prob) (1-beta :prob)
              (sample-ratio :posnum))
  (let* ((q1 (- 1 p1))
         (q2 (- 1 p2))
         (delta (abs (- p1 p2)))
         (p-bar (/ (+ p1 (* sample-ratio p2)) (1+ sample-ratio)))
         (q-bar (- 1 p-bar))
         (z-alpha (z (- 1 (if (eql tails :both) (/ alpha 2) alpha))))
         (z-beta (z 1-beta))
         (n1 (round-up
              (/ (square (+ (* (sqrt (* p-bar q-bar (1+ (/ sample-ratio))))
                               z-alpha)
                            (* (sqrt (+ (* p1 q1) (/ (* p2 q2) sample-ratio)))
                               z-beta)))
                 (square delta)))))
    (values n1 (round-up (* sample-ratio n1)))))

;; Rosner 387
(defun binomial-test-paired-sse (pd pa &key (alpha 0.05)
					 (1-beta 0.95) (tails :both))
  "Sample size estimate for the McNemar (discordant pairs) test.  Pd is the projected proportion of discordant pairs among all pairs, and Pa is the projected proportion of type A pairs among discordant pairs.  alpha, 1-beta and tails are as above.  Returns the number of individuals necessary; that is twice the number of matched pairs necessary."
  (test-variables (pd :prob) (pa :posnum) (alpha :prob) (1-beta :prob))
  (let ((qa (- 1 pa))
        (z-alpha (z (- 1 (if (eql tails :both) (/ alpha 2) alpha))))
        (z-beta (z 1-beta)))
    (round-up (/ (square (+ z-alpha (* 2 z-beta (sqrt (* pa qa)))))
                 (* 2 (square (- pa 1/2)) pd)))))

;; Rosner 463
(defun correlation-sse (rho &key (alpha 0.05) (1-beta 0.95))
  "Returns the size of a sample necessary to find a correlation of expected value rho with significance alpha and power 1-beta."
  (test-variables (rho :prob) (alpha :prob) (1-beta :prob))
  (round-up (+ 3 (/ (square (+ (z (- 1 alpha)) (z 1-beta)))
                    (square (fisher-z-transform rho))))))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Correlation and Regression
;;;

;; Rosner 431, 441 for t-test
(defun linear-regression (points)
  "Computes the regression equation for a least squares fit of a line to a sequence of points (each a list of two numbers, e.g. '((1.0 0.1) (2.0 0.2))) and report the intercept, slope, correlation coefficient r, R^2, and the significance of the difference of the slope from 0."
  (test-variables (points sequence))
  (let  ((xs (map 'list #'first points))
         (ys (map 'list #'second points)))
    (test-variables (xs :numseq) (ys :numseq))
    (let* ((x-bar (mean xs))
           (y-bar (mean ys))
           (n (length points))
           (Lxx (reduce #'+ (mapcar (lambda (xi) (square (- xi x-bar))) xs)))
           (Lyy (reduce #'+ (mapcar (lambda (yi) (square (- yi y-bar))) ys)))
           (Lxy (reduce #'+ (mapcar (lambda (xi yi) (* (- xi x-bar) (- yi y-bar)))
                                    xs ys)))
           (b (unless (zerop Lxx) (/ Lxy Lxx)))
           (a (when b (- y-bar (* b x-bar))))
           (reg-ss (* b Lxy))
           (res-ms (when (> n 2) (/ (- Lyy reg-ss) (- n 2))))
           (r (unless (or (zerop Lxx) (zerop Lxy)) (/ Lxy (sqrt (* Lxx Lyy)))))
           (r2 (unless (zerop Lyy) (/ reg-ss Lyy)))
           (t-test (unless (or (zerop Lxx) (zerop res-ms)) (/ b (sqrt (/ res-ms Lxx)))))
           (t-significance (when t-test (t-significance t-test (- n 2) :tails :both))))
      (values a b r r2 t-significance))))

(defun correlation-coefficient (points)
  "Pearson correlation"
  (test-variables (points sequence))
  (let ((xs (map 'list #'first points))
        (ys (map 'list #'second points)))
    (test-variables (xs :numseq) (ys :numseq))
    (let ((x-bar (mean xs))
          (y-bar (mean ys)))
      (/ (reduce #'+ (mapcar #'(lambda (xi yi) (* (- xi x-bar) (- yi y-bar)))
                             xs ys))
         (sqrt (* (reduce #'+ (mapcar #'(lambda (xi) (square (- xi x-bar)))
                                      xs))
                  (reduce #'+ (mapcar #'(lambda (yi) (square (- yi y-bar)))
                                      ys))))))))

;; Rosner 464
(defun correlation-test-two-sample (r1 n1 r2 n2 &key (tails :both))
  "Test if two correlation coefficients are different.  Users Fisher's Z test."
  (test-variables (r1 :prob) (n1 :posint) (r2 :prob) (n2 :posint))
  (let* ((z1 (fisher-z-transform r1))
         (z2 (fisher-z-transform r2))
         (lambda (/ (- z1 z2) (sqrt (+ (/ (- n1 3)) (/ (- n2 3)))))))
    (ecase tails
      (:both (* 2 (if (<= lambda 0) (phi lambda) (- 1 (phi lambda)))))
      (:positive (- 1 (phi lambda)))
      (:negative (phi lambda)))))

(defun correlation-test-two-sample-on-sequences (points1 points2 &key (tails :both))
  (test-variables (points1 sequence) (points2 sequence))
  (let ((r1 (correlation-coefficient points1))
        (n1 (length points1))
        (r2 (correlation-coefficient points2))
        (n2 (length points2)))
    (correlation-test-two-sample r1 n1 r2 n2 :tails tails)))

;; Rosner 498
(defun spearman-rank-correlation (points)
  "Spearman rank correlation computes the relationship between a pair of variables when one or both are either ordinal or have a distribution that is far from normal.  It takes a list of points (same format as linear-regression) and returns the spearman rank correlation coefficient and its significance."
  (test-variables (points sequence))
  (let ((xis (mapcar #'first points))
        (yis (mapcar #'second points)))
    (test-variables (xis :numseq) (yis :numseq))
    (let* ((n (length points))
           (sorted-xis (sort (copy-seq xis) #'<))
           (sorted-yis (sort (copy-seq yis) #'<))
           (average-x-ranks (mapcar (lambda (x) (average-rank x sorted-xis)) xis))
           (average-y-ranks (mapcar (lambda (y) (average-rank y sorted-yis)) yis))
           (mean-x-rank (mean average-x-ranks))
           (mean-y-rank (mean average-y-ranks))
           (Lxx (reduce #'+ (mapcar (lambda (xi-rank) (square (- xi-rank mean-x-rank)))
                                    average-x-ranks)))
           (Lyy (reduce #'+ (mapcar (lambda (yi-rank) (square (- yi-rank mean-y-rank)))
                                    average-y-ranks)))
           (Lxy (reduce #'+ (mapcar (lambda (xi-rank yi-rank)
                                      (* (- xi-rank mean-x-rank)
                                         (- yi-rank mean-y-rank)))
                                    average-x-ranks average-y-ranks)))
           (rs (/ Lxy (sqrt (* Lxx Lyy))))
           (ts (/ (* rs (sqrt (- n 2))) (sqrt (- 1 (square rs)))))
           (p (t-significance ts (- n 2) :tails :both)))
      (format t "~%Spearman correlation coefficient ~f, p = ~f" rs p)
      (values rs p))))



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Significance functions
;;;

;;  Lookup table in Rosner; this is adopted from CLASP/Numeric Recipes
(defun t-significance (t-statistic dof &key (tails :both))
  "Adopted from CLASP 1.4.3, http://eksl-www.cs.umass.edu/clasp.html"
  (test-variables (t-statistic number) (dof :posint))
  (setf dof (float dof t-statistic))
  (let ((a (beta-incomplete (* 0.5 dof) 0.5 (/ dof (+ dof (square t-statistic))))))
    ;; A is 2*Integral from (abs t-statistic) to Infinity of t-distribution
    (ecase tails
      (:both a)
      (:positive (if (plusp t-statistic)
		     (* .5 a)
		     (- 1.0 (* .5 a))))
      (:negative (if (plusp t-statistic)
		     (- 1.0 (* .5 a))
		     (* .5 a))))))


;; This implementation follows Numerical Recipes in C, section 6.3 and
;; the `ftest' function in section 13.4."
(defun f-significance
       (f-statistic numerator-dof denominator-dof &optional one-tailed-p)
  "Adopted from CLASP, but changed to handle F < 1 correctly in the
one-tailed case.  The `f-statistic' must be a positive number.  The degrees
of freedom arguments must be positive integers.  The `one-tailed-p' argument
is treated as a boolean."
  (setq f-statistic (float f-statistic))
  (test-variables (f-statistic :posnum) (numerator-dof :posint) (denominator-dof :posint))
  (let ((tail-area (beta-incomplete
		     (* 0.5d0 denominator-dof)
		     (* 0.5d0 numerator-dof)
		     (float (/ denominator-dof
			       (+ denominator-dof
				  (* numerator-dof f-statistic)))
			    1d0))))
    (if one-tailed-p
        (if (< f-statistic 1) (- 1 tail-area) tail-area)
	(progn (setf tail-area (* 2.0 tail-area))
	       (if (> tail-area 1.0)
		   (- 2.0 tail-area)
		   tail-area)))))

;; CHI-SQUARE and NORMAL (Gaussian) significance are calculated from their
;; respective CDF functions.



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Utilities of potential external use
;;;
#|
(defun random-sample (n sequence)
  "Return a random sample of size N from sequence, without replacement.  If N is equal to or greater than the length of the sequence, return the entire sequence.
DERECATED in favour of ls-statistics:random-sample"
  (test-variables (n integer) (sequence sequence))
  (cond ((null sequence) nil)
        ((<= n 0) nil)
        ((< (length sequence) n) sequence)
        (t (let ((one (random-pick sequence)))
             (cons one (random-sample (1- n) (remove one sequence :count 1)))))))
|#
(defun random-pick (sequence)
  "Random selection from sequence"
  (test-variables (sequence sequence))
  (when sequence (elt sequence (random (length sequence)))))

(defun random-normal (&key (mean 0) (sd 1))
  "Returns a random number with mean and standard-distribution as specified."
  (test-variables (mean number) (sd :posnum))
  (let ((random-standard (z (random 1d0))))
    (+ mean (* sd random-standard))))

(defun bin-and-count (sequence n)
  "Make N equal width bins and count the number of elements of sequence that belong in each."
  (let* ((min (reduce #'min sequence))
         (increment (/ (- (reduce #'max sequence) min) n))
         (bins (make-array n :initial-element 0)))
    (dotimes (bin n bins)
      (setf (aref bins bin)
            (count-if #'(lambda (x) (and (>= x (+ min (* bin increment)))
                                         (< x (+ min (* (1+ bin) increment)))))
                      sequence)))))

;; Rosner 458
(defun fisher-z-transform (r)
  "Transforms the correlation coefficient to an approximately normal distribution."
  (test-variables (r :prob))
  (* 1/2 (log (/ (1+ r) (- 1 r)))))

;; Rosner 88
(defun permutations (n k)
  "How many ways to take n things taken k at a time, when order matters"
  (test-variables (n :posint) (k :posint)
             ("K must be less than or equal to N" :test (<= k n)))
  (let ((p 1))
    (dotimes (i (1+ k) p)
      (setq p (* p (- n i))))))

;; Rosner 90
(defun choose (n k)
  "How may ways to take n things taken k at a time, when order doesn't matter"
  (test-variables (n :posint)
             ("K must be between 0 and N (inclusive)" :test (and (>= k 0) (<= k n))))
  (/ (factorial n) (* (factorial k) (factorial (- n k)))))

(defun mean-sd-n (sequence)
  "A combined calculation that is often useful.  Takes a sequence and returns three values: mean, standard deviation and N."
  (test-variables (sequence :numseq))
  (values (mean sequence) (standard-deviation sequence) (length sequence)))

(defun round-float (x &key (precision 5))
  "Rounds a floating point number to a specified number of digits precision."
  (test-variables (x number) (precision :posint))
  (/ (fround x (expt 10 (- precision))) (expt 10 precision)))


;;; Benjamini Y and Hochberg Y. "Controlling the false discovery
;;; rate: a practical and powerful approach to multiple testing." J
;;; R Stat Soc Ser B 57: 289 300, 1995.
(defun false-discovery-correction (p-values &key (rate 0.05))
  "A multiple testing correction that is less conservative than Bonferroni. Takes a list of p-values and a false discovery rate, and returns the number of p-values that are likely to be good enough to reject the null at that rate.  Returns a second value which is the p-value cutoff."
  (let ((number-of-tests (length p-values))
        (sorted-p-values (sort p-values #'>)))
    (do ((p-value (pop sorted-p-values) (pop sorted-p-values))
         (tests-to-go number-of-tests (1- tests-to-go)))
        ((or (null p-value)
             (<= p-value (* rate (/ tests-to-go number-of-tests))))
         (values tests-to-go (* rate (/ tests-to-go number-of-tests)))))))



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Internal functions
;;;

(defun round-up (x)
  (multiple-value-bind (rounded ignore) (ceiling x)
    (declare (ignore ignore))
    rounded))

(defun sign (x)
  (cond ((minusp x) -1)
        ((plusp x) 1)
        ((zerop x) 0)
        (t nil)))

(defun factorial (number)
  (if (not (and (integerp number) (>= number 0)))
      (error "factorial: ~a is not a positive integer" number)
    (labels ((fact (num) (if (= 0 num) 1 (* num (fact (1- num))))))
       (fact number))))

(defun average-rank (value sorted-values)
  "Average rank calculation for non-parametric tests.  Ranks are 1 based, but lisp is 0 based, so add 1!"
  (let ((first (position value sorted-values))
        (last (position value sorted-values :from-end t)))
    (1+ (if (= first last)
        first
        (/ (+ first last) 2)))))



;;; CLASP utilities:

;;; Copyright (c) 1990 - 1994 University of Massachusetts
;;; Department of Computer Science
;;; Experimental Knowledge Systems Laboratory
;;; Professor Paul Cohen, Director.
;;; All rights reserved.

;;; Permission to use, copy, modify and distribute this software and its
;;; documentation is hereby granted without fee, provided that the above
;;; copyright notice of EKSL, this paragraph and the one following appear
;;; in all copies and in supporting documentation.
;;; EKSL makes no representation about the suitability of this software for any
;;; purposes.  It is provided "AS IS", without express or implied warranties
;;; including (but not limited to) all implied warranties of merchantability
;;; and fitness for a particular purpose, and notwithstanding any other
;;; provision contained herein.  In no event shall EKSL be liable for any
;;; special, indirect or consequential damages whatsoever resulting from loss
;;; of use, data or profits, whether in an action of contract, negligence or
;;; other tortuous action, arising out of or in connection with the use or
;;; performance of this software, even if EKSL is advised of the possiblity of
;;; such damages.

;;; For more information write to clasp-support@cs.umass.edu

(defun error-function (x)
  "Adopted from CLASP 1.4.3, http://eksl-www.cs.umass.edu/clasp.html"
  ;; per CMUCL type-checker, the results of this function may not be
  ;; correct, if the input isn't a double-float.  For now, I think
  ;; it's easiest to coerce, but later it would be better to ensure
  ;; that the callers do the right thing.
  (let ((erf (gamma-incomplete 0.5d0 (square (coerce x 'double-float)))))
    (if (>= x 0.0d0) erf (- erf))))

;; Bug fix from rpgoldman@siftech.com.  Have to coerce an to a double-float, not a float.

(defun gamma-incomplete (a x)
  "Adopted from CLASP 1.4.3, http://eksl-www.cs.umass.edu/clasp.html"
  (declare (optimize (safety 3)))
  (setq a (coerce a 'double-float))
  (let ((gln (the double-float (gamma-ln a))))
    (when (= x 0.0)
      (return-from gamma-incomplete (values 0.0d0 gln)))
    (if (< x (+ a 1.0d0))
        ;; Use series representation.  The following is the code of what
        ;; Numerical Recipes in C calls ``GSER'
        (let* ((itmax 1000)
               (eps   3.0d-7)
               (ap    a)
               (sum   (/ 1d0 a))
               (del sum))
          (declare (type double-float ap sum del)
		   (type fixnum itmax))
          (dotimes (i itmax)
            (incf ap 1.0d0)
            (setf del (* del (/ x ap)))
            (incf sum del)
            (if (< (abs del) (* eps (abs sum)))
                (let ((result (underflow-goes-to-zero
                               (* sum (safe-exp (- (* a (log x)) x gln))))))
                  (return-from gamma-incomplete (values result gln)))))
          (error "Series didn't converge:~%~
                  Either a=~s is too large, or ITMAX=~d is too small." a itmax))
        ;; Use the continued fraction representation.  The following is the
        ;; code of what Numerical Recipes in C calls ``GCF.'' Their code
        ;; computes the complement of the desired result, so we subtract from
        ;; 1.0 at the end.
        (let ((itmax 1000)
              (eps   3.0e-7)
              (gold 0d0) (g 0d0) (fac 1d0) (b1 1d0) (b0 0d0)
              (anf 0d0) (ana 0d0) (an 0d0) (a1 x) (a0 1d0))
          (declare (type double-float gold g fac b1 b0 anf ana an a1 a0))
          (dotimes (i itmax)
            (setf an  (coerce (1+ i) 'double-float)
                  ana (- an a)
                  a0  (* fac (+ a1 (* a0 ana)))
                  b0  (* fac (+ b1 (* b0 ana)))
                  anf (* fac an)
                  a1  (+ (* x a0) (* anf a1))
                  b1  (+ (* x b0) (* anf b1)))
            (unless (zerop a1)
              (setf fac (/ 1.0d0 a1)
                    g   (* b1 fac))
              (if (< (abs (/ (- g gold) g)) eps)
                  (let ((result (underflow-goes-to-zero
                                 (* (safe-exp (- (* a (log x)) x gln)) g))))
                    (return-from
                     gamma-incomplete (values (- 1.0d0 result) gln)))
                  (setf gold g))))
          (error "Continued Fraction didn't converge:~%~
                  Either a=~s is too large, or ITMAX=~d is too small." a
                  itmax)))))


(defun gamma-ln (x)
  "Adopted from CLASP 1.4.3, http://eksl-www.cs.umass.edu/clasp.html"
  (cond ((<= x 0) (error "arg to gamma-ln must be positive:  ~s" x))
	((> x 1.0d302)
	 (error "Argument too large:  ~e" x))
	((= x 0.5d0)
	 ;; special case for this arg, since it is used by the error-function
	 (log (sqrt pi)))
	((< x 1)
	 ;; Use reflection formula:  Gamma(1-z) = z*pi/(Gamma(1+z)sin(pi*z))
	 (let ((z (- 1.0d0 x)))
	   (- (+ (log z) (log pi)) (+ (gamma-ln (+ 1.0 z)) (log (sin (* pi z)))))))
	(t (let* ((xx  (- x 1.0d0))
		  (tmp (+ xx 5.5d0))
		  (ser 1.0d0))
	     (declare (type double-float xx tmp ser))
	     (decf tmp (* (+ xx 0.5d0) (log tmp)))
	     (dolist (coef '(76.18009173d0 -86.50532033d0 24.01409822d0
					   -1.231739516d0 0.120858003d-2 -0.536382d-5))
	       (declare (type double-float coef))
	       (incf xx 1.0d0)
	       (incf ser (/ coef xx)))
	     (- (log (* 2.50662827465d0 ser)) tmp)))))

(defun error-function-complement (x)
  "Adopted from CLASP 1.4.3, http://eksl-www.cs.umass.edu/clasp.html"
  (let* ((z   (abs x))
	 (y   (/ 1d0 (+ 1d0 (* 0.5 z))))   ; instead of t
	 (ans
	      (* y (exp (+ (* (- z) z)
			   -1.26551223
			   (* y
			      (+ 1.00002368
				 (* y
				    (+ 0.37409196
				       (* y
					  (+ 0.09678418
					     (* y
						(+ -0.18628806
						   (* y
						      (+ 0.27886807
							 (* y
							    (+ -1.13520398
							       (* y
								  (+ 1.48851587
								     (* y
									(+ -0.82215223
									   (* y 0.17087277))))))))))))))))))))))
    (declare (type double-float z y ans))
    (if (>= x 0.0)
	ans
	(- 2.0 ans))))

(defun find-critical-value
       (p-function p-value &optional (x-tolerance .00001) (y-tolerance .00001))
  "Adopted from CLASP 1.4.3, http://eksl-www.cs.umass.edu/clasp.html"
  (let* ((x-low 0d0)
	 (fx-low 1d0)
	 (x-high 1d0)
	 (fx-high (coerce (funcall p-function x-high) 'double-float)))
    ;; double up
    (declare (type double-float x-low fx-low x-high fx-high))
    (do () (nil)
      ;; for general functions, we'd have to try the other way of bracketing,
      ;; and probably have another way to terminate if, say, y is not in the
      ;; range of f.
      (when (>= fx-low p-value fx-high)
	(return))
      (setf x-low x-high
	    fx-low fx-high
	    x-high (* 2.0 x-high)
	    fx-high (funcall p-function x-high)))
    ;; binary search
    (do () (nil)
      (let* ((x-mid  (/ (+ x-low x-high) 2.0))
	     (fx-mid (funcall p-function x-mid))
	     (y-diff (abs (- fx-mid p-value)))
	     (x-diff (- x-high x-low)))
	(when (or (< x-diff x-tolerance)
		  (< y-diff y-tolerance))
	  (return-from find-critical-value x-mid))
	;; Because significance is monotonically decreasing with x, if the
	;; function is above the desired p-value...
	(if (< p-value fx-mid)
	    ;; then the critical x is in the upper half
	    (setf x-low x-mid
		  fx-low fx-mid)
	    ;; otherwise, it's in the lower half
	    (setf x-high x-mid
		  fx-high fx-mid))))))

(defun beta-incomplete (a b x)
  "Adopted from CLASP 1.4.3, http://eksl-www.cs.umass.edu/clasp.html"
   (flet ((betacf (a b x)
	    ;; straight from Numerical Recipes in C, section 6.3
	    (declare (type double-float a b x))
	    (let ((ITMAX 1000)
		  (EPS   3.0d-7)
		  (qap 0d0) (qam 0d0) (qab 0d0) (em  0d0) (tem 0d0) (d 0d0)
		  (bz  0d0) (bm  1d0) (bp  0d0) (bpp 0d0)
		  (az  1d0) (am  1d0) (ap  0d0) (app 0d0) (aold 0d0))
	      (declare (type double-float qap qam qab tem d
			     bz bm bp bpp az am ap app aold em))
	      (setf qab (+ a b)
		    qap (+ a 1d0)
		    qam (- a 1d0)
		    bz  (- 1d0 (/ (* qab x) qap)))
	      (dotimes (m ITMAX)
		(setf em   (coerce (float (1+ m)) 'double-float))
                (setf
		      tem  (+ em em)
		      d    (/ (* em (- b em) x)
			      (* (+ qam tem) (+ a tem)))
		      ap   (+ az (* d am))
		      bp   (+ bz (* d bm))
		      d    (/ (* (- (+ a em)) (+ qab em) x)
			      (* (+ qap tem) (+ a tem)))
		      app  (+ ap (* d az))
		      bpp  (+ bp (* d bz))
		      aold az
		      am   (/ ap bpp)
		      bm   (/ bp bpp)
		      az   (/ app bpp)
		      bz   1d0)
		(if (< (abs (- az aold)) (* EPS (abs az)))
		    (return-from betacf az)))
	      (error "a=~s or b=~s too big, or ITMAX too small in BETACF"
		     a b))))
      (declare (notinline betacf))
      (setq a (coerce a 'double-float) b (coerce b 'double-float)
            x (coerce x 'double-float))
      (when (or (< x 0d0) (> x 1d0))
	 (error "x must be between 0d0 and 1d0:  ~f" x))
      ;; bt is the factors in front of the continued fraction
      (let ((bt (if (or (= x 0d0) (= x 1d0))
		    0d0
		    (exp (+ (gamma-ln (+ a b))
			    (- (gamma-ln a))
			    (- (gamma-ln b))
			    (* a (log x))
			    (* b (log (- 1d0 x))))))))
	 (if (< x (/ (+ a 1d0) (+ a b 2.0)))
	     ;; use continued fraction directly
	     (/ (* bt (betacf a b x)) a)
	     ;; use continued fraction after making the symmetry transformation
	     (- 1d0 (/ (* bt (betacf b a (- 1d0 x))) b))))))

(defun safe-exp (x)
  "Eliminates floating point underflow for the exponential function.
Instead, it just returns 0.0d0"
  (setf x (coerce x 'double-float))
  (if (< x (log least-positive-double-float))
      0.0d0
      (exp x)))

;;; </code>