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 STANDARD-DEVIATION is simply the root
  70 ;; (for interesting definitions of "simply"), and VARIANCE is a bit
  71 ;; more general.
  72 (defgeneric variance (x)
  73   (:documentation "compute the variance of the entity X, i.e. if a
  74   scalar, vector, or (covariance) matrix.")
  75   (:method ((x list))
  76     (let ((n (length x))
  77           (r (- x (mean x))))
  78       (sqrt (* (mean (* r r)) (/ n (- n 1))))))
  79   (:method ((x vector-like))
  80     ;; FIXME!!!
  81     (let* ((negone nil) ;; FIXME!
  82            (negresid  (axpy (mean x) negone x))) ; -mu + x = x - mu
  83       (/ (loop for i from 0 to (- (nelts x) 1)
  84             summing (* (vref negresid i)
  85                        (vref negresid i)))
  86          (- (nelts negresid) 1))))
  87   (:method ((x matrix-like))
  88     (error "FIXME: define variance for matrices as covariance (?).")))
  89
  90 (defun covariance-matrix (&rest args)
  91 "Args: (&rest args)
  92 Returns the sample covariance matrix of the data columns in ARGS. ARGS may
  93 consist of lists, vectors or matrices."
  94   (let ((columns (apply #'append
  95                         (mapcar #'(lambda (x)
  96                                     (if (typep x 'matrix-like)
  97                                         (list-of-columns x)
  98                                         (list x)))
  99                                 args))))
 100     (/ (cross-product (reduce #'bind2
 101                               (- columns (mapcar #'mean columns))))
 102        (- (length (car columns)) 1))))
 103
 104 (defgeneric standard-deviation (x)
 105   (:documentation "Compute standard deivation of entity X.")
 106   (:method ((x vector-like)) (sqrt (variance x)))
 107   (:method ((x list))
 108     ;; if elements are not  numeric, error, otherwise
 109     (sqrt (variance x)))
 110   (:method ((x matrix-like))
 111     (error "FIXME: define SD for matrix-like objects")))
 112
 113 (defun standard-deviation-fn (x)
 114   "Function front-end to generic.  Do we really need this?"
 115   (standard-deviation x))
 116
 117 (defun standard-deviation-old (x)
 118 "Args: (x)
 119 Returns the standard deviation of the elements x. Vector reducing.
 120
 121 FIXME AJR: This should be redone as square-root of the variance, and
 122 defer to that structure for computation."
 123   (let ((n (count-elements x))
 124         (r (- x (mean x))))
 125     (sqrt (* (mean (* r r)) (/ n (- n 1))))))
 126
 127 (defgeneric quantile (x p)
 128   (:documentation "Returns the P-th quantile(s) of sequence X as a
 129   native lisp list (which could be fed into a vector-like or similar
 130   for the final format).")
 131   (:method ((x sequence) (p real))
 132     (let ((np (* p (- (length x) 1))))
 133       (mean (list (aref x (floor np))
 134                   (aref x (ceiling np)))))) ;; aref work in general for lists too, or use nth?
 135   (:method ((x vector-like) (p real))  ;; average of sorted elements.  Could store compile
 136     (let ((np (* p (- (nelts x) 1))))
 137       (mean (list (vref (sort x) (floor np))
 138                   (vref (sort x) (ceiling np))))))
 139   (:method ((x sequence) (p sequence))
 140     (error "FIXME: generalize.   Basically, do mapcar or similar across a vector."))
 141   (:method ((x vector-like) (p sequence))
 142     (error "FIXME: generalize."))
 143   (:method ((x vector-like) (p vector-like))
 144     (error "FIXME: generalize.")))
 145
 146 ;;; things to build on top of quantile.
 147 ;;
 148 ;; Do as generics if possible, keeping in mind that quantile returns a
 149 ;; raw lisp list.  Nearly all of these use native lisp structures for speed,
 150 ;; such as raw lists for single and vector values, and arrays
 151 ;; (simple?) for matrix-values  and array-valued results.
 152
 153 (defun median (x)
 154   "Return median of X, using whatever the quantile generic function supports."
 155   (quantile x 0.5))
 156
 157 (defun interquartile-range (x)
 158   "Returns the interquartile range of the elements of X."
 159   (reduce #'- (quantile x (list 0.75 0.25))))
 160
 161 (defun fivnum (x)
 162   "Returns the five number summary (min, 1st quartile, medinan, 3rd quartile,
 163 max) of the elements X."
 164   (quantile x (list 0 .25 .5 .75 1)))
 165
 166 ;;;; Sampling / Resampling
 167
 168 (defgeneric sample (x n &optional replace weights)
 169   (:documentation "Draw a sample of size n from sequence x, with/out
 170   replacement, with weights per element of x as described as a
 171   probability mass distribution vector.  Length of weights should
 172   match length of X.  Unclear if (= X (uniq X))? (possible issue).
 173   Return value shoud match the input values.")
 174   (:method ((x sequence) (n fixnum)
 175             &optional replace weights)
 176     (sample-fn x n replace))
 177   (:method ((x vector-like) (n fixnum)
 178             &optional replace weights)
 179     (error "SAMPLE for vector-likes not implemented yet.")))
 180
 181 (defun sample-fn (x ssize &optional replace)
 182 "Args: (x n &optional (replace nil))
 183 Returns a list of a random sample of size N from sequence X drawn with or
 184 without replacement."
 185   (check-sequence x)
 186   (let ((n (length x))
 187         (x (if (consp x) (coerce x 'vector) (copy-vector x)))
 188         (result nil))
 189     (if (< 0 n)
 190         (dotimes (i ssize result)
 191                  (let ((j (if replace (random n) (+ i (random (- n i))))))
 192                    (setf result (cons (aref x j) result))
 193                    (unless replace     ;; swap elements i and j
 194                            (let ((temp (aref x i)))
 195                              (setf (aref x i) (aref x j))
 196                              (setf (aref x j) temp))))))))
 197
 198
 199