src/describe/statistics.lsp

   1 ;;; -*- mode: lisp -*-
   2 ;;; Copyright (c) 2005--2009, by A.J. Rossini <blindglobe@gmail.com>
   3 ;;; See COPYRIGHT file for any additional restrictions (BSD license).
   4 ;;; Since 1991, ANSI was finally finished.  Edited for ANSI Common Lisp.
   5
   6 ;;; XLisp-ism's removed to focus on Common Lisp.  Original source from:
   7 ;;;; statistics.lsp XLISP-STAT statistics functions
   8 ;;;; XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney
   9 ;;;; Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz
  10 ;;;; You may give out copies of this software; for conditions see the file
  11 ;;;; COPYING included with this distribution.
  12
  13 (in-package :lisp-stat-descriptive-statistics)
  14
  15 ;;;
  16 ;;; Basic Summary Statistics
  17 ;;;
  18
  19 (defgeneric mean (x)
  20   (:documentation "compute the mean of lists, vectors, various objects")
  21   (:method ((x list))
  22     (/ (reduce #'+ x)
  23        (length x)))
  24   (:method ((x vector-like))
  25       ;; (defparameter *x* (make-vector 5 :initial-contents '((1d0 2d0 3d0 4d0 5d0))))
  26     (/ (loop for i from 0 to (- (nelts x) 1)
  27           summing (vref x i))
  28        (nelts x)))
  29 #| ;; a more traditional versino of the above...
  30   (:method ((x vector-like))
  31     (let ((n (nelts x))
  32           (type (if (col-vector-p x) :column :row)))
  33       (/ (gemm x (ones
  34                   (ecase type (:row 1) (:column n))
  35                   (ecase type (:row n) (:column 1))))
  36          n)))
  37 |#
  38   )
  39
  40 (defun mean-fn (x)
  41 "Args: (x)
  42 Returns the mean of the elements x. Vector reducing.
  43
  44 FIXME: Why is this so complex?  When figure out, put into generic."
  45   (let ((mean 0.0)
  46         (count 0.0))
  47     (labels ((add-to-mean (x)
  48               (let ((count+1 (+ count 1.0)))
  49                 (setf mean (+ (* (/ count count+1) mean) (* (/ count+1) x)))
  50                 (setf count count+1)))
  51              (find-mean (x)
  52                (if (numberp x)
  53                  (add-to-mean x)
  54                  (let ((seq (compound-data-seq x)))
  55                    (if (consp seq)
  56                      (dolist (x seq)
  57                        (if (numberp x) (add-to-mean x) (find-mean x)))
  58                      (let ((n (length seq)))
  59                        (dotimes (i n)
  60                          (declare (fixnum i))
  61                          (let ((x (aref seq i)))
  62                            (if (numberp x)
  63                                (add-to-mean x)
  64                                (find-mean x))))))))))
  65       (find-mean x)
  66       mean)))
  67
  68
  69 ;; We do the variance, since the SD is simply the root.
  70 (defgeneric variance (x)
  71   (:documentation "compute the variance of the entity X, i.e. if a
  72   scalar, vector, or (covariance) matrix.")
  73   (:method ((x list))
  74     (let ((n (length x))
  75           (r (- x (mean x))))
  76       (sqrt (* (mean (* r r)) (/ n (- n 1))))))
  77   (:method ((x vector-like))
  78     ;; FIXME!!!
  79     (let ((negresid  (axpy (mean x) negone x)))  ; -mu + x = x - mu
  80       (/ (loop for i from 0 to (- (nelts x) 1)
  81             summing (* (vref negresid i)
  82                        (vref negresid i)))
  83          (- (nelts negresid) 1))))
  84   (:method ((x matrix-like))
  85     (error "FIXME: define variance for matrices as covariance (?).")))
  86
  87 (defun covariance-matrix (&rest args)
  88 "Args: (&rest args)
  89 Returns the sample covariance matrix of the data columns in ARGS. ARGS may
  90 consist of lists, vectors or matrices."
  91   (let ((columns (apply #'append
  92                         (mapcar #'(lambda (x)
  93                                     (if (typep x 'matrix-like)
  94                                         (list-of-columns x)
  95                                         (list x)))
  96                                 args))))
  97     (/ (cross-product (reduce #'bind2
  98                               (- columns (mapcar #'mean columns))))
  99        (- (length (car columns)) 1))))
 100
 101 (defgeneric standard-deviation (x)
 102   (:documentation "Compute standard deivation of entity X.")
 103   (:method ((x vector-like)) (sqrt (variance x)))
 104   (:method ((x list))
 105     ;; if elements are not  numeric, error, otherwise
 106     (sqrt (variance x)))
 107   (:method ((x matrix-like))
 108     (error "FIXME: define SD for matrix-like objects")))
 109
 110 (defun standard-deviation-fn (x)
 111 "Args: (x)
 112 Returns the standard deviation of the elements x. Vector reducing.
 113
 114 FIXME AJR: This should be redone as square-root of the variance, and
 115 defer to that structure for computation."
 116   (let ((n (count-elements x))
 117         (r (- x (mean x))))
 118     (sqrt (* (mean (* r r)) (/ n (- n 1))))))
 119
 120 (defgeneric quantile (x p)
 121   (:documentation "Returns the P-th quantile(s) of sequence X.")
 122   (:method ((x sequence) (p real))
 123     (let ((np (* p (- n 1))))
 124       (mean (aref x (floor np))
 125             (aref x (ceiling np))))) ;; aref work in general for lists too, or use nth?
 126   (:method ((x vector-like) (p real))  ;; average of sorted elements.  Could store compile
 127     (let ((np (* p (- (nelts x) 1))))
 128       (mean (list (vref (sort x) (floor np))
 129                   (vref (sort x) (ceiling np))))))
 130   (:method ((x sequence) (p sequence))
 131     (error "FIXME: generalize.   Basically, do mapcar or similar across a vector."))
 132   (:method ((x vector-like) (p sequence))
 133     (error "FIXME: generalize."))
 134   (:method ((x vector-like) (p vector-like))
 135     (error "FIXME: generalize.")))
 136
 137
 138 ;;; things to build on top of quantiles...!
 139
 140 ;; Args: (x)
 141 ;; Returns the median of the elements of X.
 142 (defmacro median (x)
 143   `(quantile ,x 0.5))
 144 ;; (macroexpand '(median (list 1 2 3)))
 145 ;; (median (list 1 2 3))
 146
 147
 148 ;; Args: (number-data)
 149 ;; Returns the interquartile range of the elements of X.
 150 (defmacro interquartile-range (x)
 151   `(reduce #'- (quantile ,x '(0.75 0.25))))
 152
 153 (defun fivnum (x)
 154 "Args: (number-data)
 155 Returns the five number summary (min, 1st quartile, medinan, 3rd quartile,
 156 max) of the elements X."
 157   (quantile x '(0 .25 .5 .75 1)))
 158
 159 ;;;; Sampling / Resampling
 160
 161 (defgeneric sample (x n &optional replace weights)
 162   (:documentation "Draw a sample of size n from sequence x, with/out
 163   replacement, with weights per element of x as described as a
 164   probability mass distribution vector.")
 165   (:method ((x list) (n fixnum)
 166             &optional (replace list) (weights list))
 167     )
 168   (:method ((x vector-like) (n fixnum)
 169             &optional (replace vector-like) (weights vector-like))
 170     )
 171   )
 172
 173 (defun sample-fn (x ssize &optional replace)
 174 "Args: (x n &optional (replace nil))
 175 Returns a list of a random sample of size N from sequence X drawn with or
 176 without replacement."
 177   (check-sequence x)
 178   (let ((n (length x))
 179         (x (if (consp x) (coerce x 'vector) (copy-vector x)))
 180         (result nil))
 181     (if (< 0 n)
 182         (dotimes (i ssize result)
 183                  (let ((j (if replace (random n) (+ i (random (- n i))))))
 184                    (setf result (cons (aref x j) result))
 185                    (unless replace     ;; swap elements i and j
 186                            (let ((temp (aref x i)))
 187                              (setf (aref x i) (aref x j))
 188                              (setf (aref x j) temp))))))))
 189
 190
 191