statistics.lsp

   1 ;;; -*- mode: lisp -*-
   2 ;;; Copyright (c) 2005--2008, 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 :cl-user)
  14
  15 (defpackage :lisp-stat-descriptive-statistics
  16  (:use :common-lisp
  17        :lisp-stat-data
  18        :lisp-stat-math
  19        :lisp-stat-compound-data
  20        :lisp-stat-matrix
  21        :lisp-stat-linalg-data
  22        :lisp-stat-linalg
  23        :lisp-stat-basics)
  24    (:shadowing-import-from :lisp-stat-math ;; life is a vector!
  25       expt + - * / ** mod rem abs 1+ 1- log exp sqrt sin cos tan
  26       asin acos atan sinh cosh tanh asinh acosh atanh float random
  27       truncate floor ceiling round minusp zerop plusp evenp oddp
  28       < <= = /= >= > ;; complex
  29       conjugate realpart imagpart phase
  30       min max logand logior logxor lognot ffloor fceiling
  31       ftruncate fround signum cis)
  32    (:export standard-deviation quantile median interquartile-range
  33             fivnum sample))
  34
  35 (in-package :lisp-stat-descriptive-statistics)
  36
  37 ;;;
  38 ;;; Basic Summary Statistics
  39 ;;;
  40
  41 (defun standard-deviation (x)
  42 "Args: (x)
  43 Returns the standard deviation of the elements x. Vector reducing."
  44   (let ((n (count-elements x))
  45         (r (- x (mean x))))
  46     (sqrt (* (mean (* r r)) (/ n (- n 1))))))
  47
  48
  49 (defgeneric quantile (x p)
  50   (:documentation "Args: (x p)
  51 Returns the P-th quantile(s) of sequence X. P can be a number or a sequence."))
  52
  53 (defmethod quantile ((x sequence) (p sequence))
  54   (let* ((x (sort-data x))
  55          (n (length x))
  56          (np (* p (- n 1)))
  57          (low (floor np))
  58          (high (ceiling np)))
  59     (/ (+ (select x low) (select x high)) 2)))
  60
  61 (defmethod quantile ((x sequence) (p real))
  62   (let* ((x (sort-data x))
  63          (n (length x))
  64          (np (* p (- n 1)))
  65          (low (floor np))
  66          (high (ceiling np)))
  67     (/ (+ (select x low) (select x high)) 2)))
  68
  69
  70 ;;; things to build on top of quantiles...!
  71
  72 ;; Args: (x)
  73 ;; Returns the median of the elements of X.
  74 (defmacro median (x)
  75   `(quantile ,x 0.5))
  76 ;; (macroexpand '(median (list 1 2 3)))
  77 ;; (median (list 1 2 3))
  78
  79
  80 ;; Args: (number-data)
  81 ;; Returns the interquartile range of the elements of X.
  82 (defmacro interquartile-range (x)
  83   '(reduce #'- (quantile ,x '(0.75 0.25))))
  84
  85 (defun fivnum (x)
  86 "Args: (number-data)
  87 Returns the five number summary (min, 1st quartile, medinan, 3rd quartile,
  88 max) of the elements X."
  89   (quantile x '(0 .25 .5 .75 1)))
  90
  91 (defun covariance-matrix (&rest args)
  92 "Args: (&rest args)
  93 Returns the sample covariance matrix of the data columns in ARGS. ARGS may
  94 consist of lists, vectors or matrices."
  95   (let ((columns (apply #'append
  96                         (mapcar #'(lambda (x)
  97                                     (if (matrixp x) (column-list x) (list x)))
  98                                 args))))
  99     (/ (cross-product (apply #'bind-columns
 100                              (- columns (mapcar #'mean columns))))
 101        (- (length (car columns)) 1))))
 102
 103 ;;;; Sampling / Resampling
 104
 105 (defun sample (x ssize &optional replace)
 106 "Args: (x n &optional (replace nil))
 107 Returns a list of a random sample of size N from sequence X drawn with or
 108 without replacement."
 109   (check-sequence x)
 110   (let ((n (length x))
 111         (x (if (consp x) (coerce x 'vector) (copy-vector x)))
 112         (result nil))
 113     (if (< 0 n)
 114         (dotimes (i ssize result)
 115                  (let ((j (if replace (random n) (+ i (random (- n i))))))
 116                    (setf result (cons (aref x j) result))
 117                    (unless replace     ;; swap elements i and j
 118                            (let ((temp (aref x i)))
 119                              (setf (aref x i) (aref x j))
 120                              (setf (aref x j) temp))))))))
 121
 122
 123
 124 (defun mean (x)
 125 "Args: (x)
 126 Returns the mean of the elements x. Vector reducing."
 127   (let ((mean 0.0)
 128         (count 0.0))
 129     (labels ((add-to-mean (x)
 130               (let ((count+1 (+ count 1.0)))
 131                 (setf mean (+ (* (/ count count+1) mean) (* (/ count+1) x)))
 132                 (setf count count+1)))
 133              (find-mean (x)
 134                (if (numberp x)
 135                  (add-to-mean x)
 136                  (let ((seq (compound-data-seq x)))
 137                    (if (consp seq)
 138                      (dolist (x seq)
 139                        (if (numberp x) (add-to-mean x) (find-mean x)))
 140                      (let ((n (length seq)))
 141                        (dotimes (i n)
 142                          (declare (fixnum i))
 143                          (let ((x (aref seq i)))
 144                            (if (numberp x)
 145                                (add-to-mean x)
 146                                (find-mean x))))))))))
 147       (find-mean x)
 148       mean)))
 149