statistics.lsp

   1 ;;; -*- mode: lisp -*-
   2 ;;; Copyright (c) 2005--2007, 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 (defpackage :lisp-stat
  14  (:use :common-lisp
  15        :lisp-stat-data)
  16  (:export standard-deviation quantile median interquartile-range
  17           fivnum covariance-matrix difference rseq matrix print-matrix solve
  18           backsolve eigenvalues eigenvectors accumulate cumsum combine
  19           lowess))
  20
  21 (:in-package :lisp-stat)
  22
  23 ;;;;
  24 ;;;; Basic Summary Statistics
  25 ;;;;
  26
  27 (defun standard-deviation (x)
  28 "Args: (x)
  29 Returns the standard deviation of the elements x. Vector reducing."
  30   (let ((n (count-elements x))
  31         (r (- x (mean x))))
  32     (sqrt (* (mean (* r r)) (/ n (- n 1))))))
  33
  34 (defun quantile (x p)
  35 "Args: (x p)
  36 Returns the P-th quantile(s) of sequence X. P can be a number or a sequence."
  37   (let* ((x (sort-data x))
  38          (n (length x))
  39          (np (* p (- n 1)))
  40          (low (floor np))
  41          (high (ceiling np)))
  42     (/ (+ (select x low) (select x high)) 2)))
  43
  44 (defun median (x)
  45 "Args: (x)
  46 Returns the median of the elements of X."
  47   (quantile x 0.5))
  48
  49 (defun interquartile-range (x)
  50 "Args: (number-data)
  51 Returns the interquartile range of the elements of X."
  52   (apply #'- (quantile x '(0.75 0.25))))
  53
  54 (defun fivnum (x)
  55 "Args: (number-data)
  56 Returns the five number summary (min, 1st quartile, medinan, 3rd quartile,
  57 max) of the elements X."
  58   (quantile x '(0 .25 .5 .75 1)))
  59
  60 (defun covariance-matrix (&rest args)
  61 "Args: (&rest args)
  62 Returns the sample covariance matrix of the data columns in ARGS. ARGS may
  63 consist of lists, vectors or matrices."
  64   (let ((columns (apply #'append
  65                         (mapcar #'(lambda (x)
  66                                     (if (matrixp x) (column-list x) (list x)))
  67                                 args))))
  68     (/ (cross-product (apply #'bind-columns
  69                              (- columns (mapcar #'mean columns))))
  70        (- (length (car columns)) 1))))
  71
  72 ;;;;
  73 ;;;; Basic Sequence Operations
  74 ;;;;
  75
  76 (defun difference (x)
  77 "Args: (x)
  78 Returns differences for a sequence X."
  79   (let ((n (length x)))
  80     (- (select x (iseq 1 (1- n))) (select x (iseq 0 (- n 2))))))
  81
  82 (defun rseq (a b num)
  83 "Args: (a b num)
  84 Returns a list of NUM equally spaced points starting at A and ending at B."
  85   (+ a (* (iseq 0 (1- num)) (/ (float (- b a)) (1- num)))))
  86
  87 ;;;;
  88 ;;;; Linear Algebra Functions
  89 ;;;;
  90
  91 (defun matrix (dim data)
  92 "Args: (dim data)
  93 returns a matrix of dimensions DIM initialized using sequence DATA
  94 in row major order."
  95   (let ((dim (coerce dim 'list))
  96         (data (coerce data 'list)))
  97     (make-array dim :initial-contents (split-list data (nth 1 dim)))))
  98
  99 (defun print-matrix (a &optional (stream *standard-output*))
 100 "Args: (matrix &optional stream)
 101 Prints MATRIX to STREAM in a nice form that is still machine readable"
 102   (unless (matrixp a) (error "not a matrix - ~a" a))
 103   (let ((size (min 15 (max (map-elements #'flatsize a)))))
 104     (format stream "#2a(~%")
 105     (dolist (x (row-list a))
 106             (format stream "    (")
 107             (let ((n (length x)))
 108               (dotimes (i n)
 109                        (let ((y (aref x i)))
 110                          (cond
 111                            ((integerp y) (format stream "~vd" size y))
 112                            ((floatp y) (format stream "~vg" size y))
 113                            (t (format stream "~va" size y))))
 114                        (if (< i (- n 1)) (format stream " "))))
 115             (format stream ")~%"))
 116     (format stream "   )~%")
 117     nil))
 118
 119 (defun solve (a b)
 120 "Args: (a b)
 121 Solves A x = B using LU decomposition and backsolving. B can be a sequence
 122 or a matrix."
 123   (let ((lu (lu-decomp a)))
 124     (if (matrixp b)
 125         (apply #'bind-columns
 126                (mapcar #'(lambda (x) (lu-solve lu x)) (column-list b)))
 127         (lu-solve lu b))))
 128
 129 (defun backsolve (a b)
 130 "Args: (a b)
 131 Solves A x = B by backsolving, assuming A is upper triangular. B must be a
 132 sequence. For use with qr-decomp."
 133   (let* ((n (length b))
 134          (sol (make-array n)))
 135     (dotimes (i n)
 136              (let* ((k (- n i 1))
 137                     (val (elt b k)))
 138                (dotimes (j i)
 139                         (let ((l (- n j 1)))
 140                           (setq val (- val (* (aref sol l) (aref a k l))))))
 141                (setf (aref sol k) (/ val (aref a k k)))))
 142     (if (listp b) (coerce sol 'list) sol)))
 143
 144 (defun eigenvalues (a)
 145 "Args: (a)
 146 Returns list of eigenvalues of square, symmetric matrix A"
 147   (first (eigen a)))
 148
 149 (defun eigenvectors (a)
 150 "Args: (a)
 151 Returns list of eigenvectors of square, symmetric matrix A"
 152   (second (eigen a)))
 153
 154 (defun accumulate (f s)
 155 "Args: (f s)
 156 Accumulates elements of sequence S using binary function F.
 157 (accumulate #'+ x) returns the cumulative sum of x."
 158   (let* ((result (list (elt s 0)))
 159          (tail result))
 160     (flet ((acc (dummy x)
 161                 (rplacd tail (list (funcall f (first tail) x)))
 162                 (setf tail (cdr tail))))
 163       (reduce #'acc s))
 164     (if (vectorp s) (coerce result 'vector) result)))
 165
 166 (defun cumsum (x)
 167   "Args: (x)
 168 Returns the cumulative sum of X."
 169   (accumulate #'+ x))
 170
 171 (defun combine (&rest args)
 172   "Args (&rest args)
 173 Returns sequence of elements of all arguments."
 174   (copy-seq (element-seq args)))
 175
 176 (defun lowess (x y &key (f .25) (steps 2) (delta -1) sorted)
 177 "Args: (x y &key (f .25) (steps 2) delta sorted)
 178 Returns (list X YS) with YS the LOWESS fit. F is the fraction of data used for
 179 each point, STEPS is the number of robust iterations. Fits for points within
 180 DELTA of each other are interpolated linearly. If the X values setting SORTED
 181 to T speeds up the computation."
 182   (let ((x (if sorted x (sort-data x)))
 183         (y (if sorted y (select y (order x))))
 184         (delta (if (> delta 0.0) delta (/ (- (max x) (min x)) 50))))
 185     (list x)));; (|base-lowess| x y f steps delta))))