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