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* ((negone nil) ;; FIXME!
  80            (negresid  (axpy (mean x) negone x))) ; -mu + x = x - mu
  81       (/ (loop for i from 0 to (- (nelts x) 1)
  82             summing (* (vref negresid i)
  83                        (vref negresid i)))
  84          (- (nelts negresid) 1))))
  85   (:method ((x matrix-like))
  86     (error "FIXME: define variance for matrices as covariance (?).")))
  87
  88 (defun covariance-matrix (&rest args)
  89 "Args: (&rest args)
  90 Returns the sample covariance matrix of the data columns in ARGS. ARGS may
  91 consist of lists, vectors or matrices."
  92   (let ((columns (apply #'append
  93                         (mapcar #'(lambda (x)
  94                                     (if (typep x 'matrix-like)
  95                                         (list-of-columns x)
  96                                         (list x)))
  97                                 args))))
  98     (/ (cross-product (reduce #'bind2
  99                               (- columns (mapcar #'mean columns))))
 100        (- (length (car columns)) 1))))
 101
 102 (defgeneric standard-deviation (x)
 103   (:documentation "Compute standard deivation of entity X.")
 104   (:method ((x vector-like)) (sqrt (variance x)))
 105   (:method ((x list))
 106     ;; if elements are not  numeric, error, otherwise
 107     (sqrt (variance x)))
 108   (:method ((x matrix-like))
 109     (error "FIXME: define SD for matrix-like objects")))
 110
 111 (defun standard-deviation-fn (x)
 112 "Args: (x)
 113 Returns the standard deviation of the elements x. Vector reducing.
 114
 115 FIXME AJR: This should be redone as square-root of the variance, and
 116 defer to that structure for computation."
 117   (let ((n (count-elements x))
 118         (r (- x (mean x))))
 119     (sqrt (* (mean (* r r)) (/ n (- n 1))))))
 120
 121 (defgeneric quantile (x p &optional return-type)
 122   (:documentation "Returns the P-th quantile(s) of sequence X.")
 123   (:method ((x sequence) (p real) &optional (return-type  symbol))
 124     (let ((np (* p (- (length x) 1))))
 125       (mean (list (aref x (floor np))
 126                   (aref x (ceiling np)))))) ;; aref work in general for lists too, or use nth?
 127   (:method ((x vector-like) (p real) &optional (return-type symbol))  ;; average of sorted elements.  Could store compile
 128     (let ((np (* p (- (nelts x) 1))))
 129       (mean (list (vref (sort x) (floor np))
 130                   (vref (sort x) (ceiling np))))))
 131   (:method ((x sequence) (p sequence) &optional return-type)
 132     (error "FIXME: generalize.   Basically, do mapcar or similar across a vector."))
 133   (:method ((x vector-like) (p sequence) &optional return-type)
 134     (error "FIXME: generalize."))
 135   (:method ((x vector-like) (p vector-like) &optional return-type)
 136     (error "FIXME: generalize.")))
 137
 138 ;;; things to build on top of quantiles...!
 139 ;; need to incorporate a return-type if possible.
 140
 141 (defun median (x)
 142   "Return median of X, using whatever the quantile generic function supports."
 143   (quantile x 0.5))
 144
 145 (defun interquartile-range (x)
 146   "Returns the interquartile range of the elements of X."
 147   (reduce #'- (quantile x (list 0.75 0.25))))
 148
 149 (defun fivnum (x)
 150   "Returns the five number summary (min, 1st quartile, medinan, 3rd quartile,
 151 max) of the elements X."
 152   (quantile x (list 0 .25 .5 .75 1)))
 153
 154 ;;;; Sampling / Resampling
 155
 156 (defgeneric sample (x n &optional replace weights)
 157   (:documentation "Draw a sample of size n from sequence x, with/out
 158   replacement, with weights per element of x as described as a
 159   probability mass distribution vector.")
 160   (:method ((x sequence) (n fixnum)
 161             &optional (replace sequence) (weights sequence))
 162     (sample-fn x n replace))
 163   (:method ((x vector-like) (n fixnum)
 164             &optional (replace vector-like) (weights vector-like))
 165     (sample-fn x n replace)))
 166
 167 (defun sample-fn (x ssize &optional replace)
 168 "Args: (x n &optional (replace nil))
 169 Returns a list of a random sample of size N from sequence X drawn with or
 170 without replacement."
 171   (check-sequence x)
 172   (let ((n (length x))
 173         (x (if (consp x) (coerce x 'vector) (copy-vector x)))
 174         (result nil))
 175     (if (< 0 n)
 176         (dotimes (i ssize result)
 177                  (let ((j (if replace (random n) (+ i (random (- n i))))))
 178                    (setf result (cons (aref x j) result))
 179                    (unless replace     ;; swap elements i and j
 180                            (let ((temp (aref x i)))
 181                              (setf (aref x i) (aref x j))
 182                              (setf (aref x j) temp))))))))
 183
 184
 185