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     (let ((n (nelts x))
  26           (type (if (column-vector-p x) :column :row)))
  27       (/ (gemm x (ones
  28                   (ecase type (:row 1) (:column n))
  29                   (ecase type (:row n) (:column 1))))
  30          n))))
  31
  32 (defun mean-fn (x)
  33 "Args: (x)
  34 Returns the mean of the elements x. Vector reducing.
  35
  36 FIXME: Why is this so complex?  When figure out, put into generic."
  37   (let ((mean 0.0)
  38         (count 0.0))
  39     (labels ((add-to-mean (x)
  40               (let ((count+1 (+ count 1.0)))
  41                 (setf mean (+ (* (/ count count+1) mean) (* (/ count+1) x)))
  42                 (setf count count+1)))
  43              (find-mean (x)
  44                (if (numberp x)
  45                  (add-to-mean x)
  46                  (let ((seq (compound-data-seq x)))
  47                    (if (consp seq)
  48                      (dolist (x seq)
  49                        (if (numberp x) (add-to-mean x) (find-mean x)))
  50                      (let ((n (length seq)))
  51                        (dotimes (i n)
  52                          (declare (fixnum i))
  53                          (let ((x (aref seq i)))
  54                            (if (numberp x)
  55                                (add-to-mean x)
  56                                (find-mean x))))))))))
  57       (find-mean x)
  58       mean)))
  59
  60
  61 ;; We do the variance, since the SD is simply the root.
  62 (defgeneric variance (x)
  63   (:documentation "compute the variance of the entity X, i.e. if a
  64   scalar, vector, or matrix.")
  65   (:method ((x list))
  66     (let ((n (length x))
  67           (r (- x (mean x))))
  68       (sqrt (* (mean (* r r)) (/ n (- n 1))))))
  69   (:method ((x vector-like))))
  70
  71 (defgeneric standard-deviation (x)
  72   (:documentation "Compute standard deivation of entity X.")
  73   (:method ((x vector-like)) (sqrt (variance x)))
  74   (:method ((x list))
  75     ;; if elements are not  numeric, error, otherwise
  76     (sqrt (variance x)))
  77   (:method ((x matrix-like))
  78     (error "FIXME: define SD for matrix-like objects")))
  79
  80 (defun standard-deviation-fn (x)
  81 "Args: (x)
  82 Returns the standard deviation of the elements x. Vector reducing.
  83
  84 FIXME AJR: This should be redone as square-root of the variance, and
  85 defer to that structure for computation."
  86   (let ((n (count-elements x))
  87         (r (- x (mean x))))
  88     (sqrt (* (mean (* r r)) (/ n (- n 1))))))
  89
  90 (defgeneric quantile (x p)
  91   (:documentation "Returns the P-th quantile(s) of sequence X.")
  92   (:method ((x sequence) (p sequence))
  93     (let* ((x (sort-data x))
  94            (n (length x))
  95            (np (* p (- n 1)))
  96            (low (floor np))
  97            (high (ceiling np)))
  98       (/ (+ (select x low) (select x high)) 2)))
  99   (:method ((x sequence) (p real))
 100     (error "FIXME: implement.")))
 101
 102
 103 ;;; things to build on top of quantiles...!
 104
 105 ;; Args: (x)
 106 ;; Returns the median of the elements of X.
 107 (defmacro median (x)
 108   `(quantile ,x 0.5))
 109 ;; (macroexpand '(median (list 1 2 3)))
 110 ;; (median (list 1 2 3))
 111
 112
 113 ;; Args: (number-data)
 114 ;; Returns the interquartile range of the elements of X.
 115 (defmacro interquartile-range (x)
 116   `(reduce #'- (quantile ,x '(0.75 0.25))))
 117
 118 (defun fivnum (x)
 119 "Args: (number-data)
 120 Returns the five number summary (min, 1st quartile, medinan, 3rd quartile,
 121 max) of the elements X."
 122   (quantile x '(0 .25 .5 .75 1)))
 123
 124 (defun covariance-matrix (&rest args)
 125 "Args: (&rest args)
 126 Returns the sample covariance matrix of the data columns in ARGS. ARGS may
 127 consist of lists, vectors or matrices."
 128   (let ((columns (apply #'append
 129                         (mapcar #'(lambda (x)
 130                                     (if (typep x 'matrix-like)
 131                                         (list-of-columns x)
 132                                         (list x)))
 133                                 args))))
 134     (/ (cross-product (reduce #'bind-columns
 135                              (- columns (mapcar #'mean columns))))
 136        (- (length (car columns)) 1))))
 137
 138 ;;;; Sampling / Resampling
 139
 140 (defun sample (x ssize &optional replace)
 141 "Args: (x n &optional (replace nil))
 142 Returns a list of a random sample of size N from sequence X drawn with or
 143 without replacement."
 144   (check-sequence x)
 145   (let ((n (length x))
 146         (x (if (consp x) (coerce x 'vector) (copy-vector x)))
 147         (result nil))
 148     (if (< 0 n)
 149         (dotimes (i ssize result)
 150                  (let ((j (if replace (random n) (+ i (random (- n i))))))
 151                    (setf result (cons (aref x j) result))
 152                    (unless replace     ;; swap elements i and j
 153                            (let ((temp (aref x i)))
 154                              (setf (aref x i) (aref x j))
 155                              (setf (aref x j) temp))))))))
 156
 157
 158