use complete forms to initialize environment
[CommonLispStat.git] / statistics.lsp
blob9e64a74042f671dbeb9ca4eca7f20903d4407be8
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.
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.
13 (in-package :cl-user)
15 (defpackage :lisp-stat-descriptive-statistics
16 (:use :common-lisp
17 :lisp-stat-data
18 :lisp-stat-math)
19 (:shadowing-import-from :lisp-stat-math
20 expt + - * / ** mod rem abs 1+ 1- log exp sqrt sin cos tan
21 asin acos atan sinh cosh tanh asinh acosh atanh float random
22 truncate floor ceiling round minusp zerop plusp evenp oddp
23 < <= = /= >= > complex conjugate realpart imagpart phase
24 min max logand logior logxor lognot ffloor fceiling
25 ftruncate fround signum cis)
26 (:export ;; descriptive stats
27 standard-deviation quantile median interquartile-range
28 fivnum
30 sample
32 ;; the following are more matrix-centric
33 covariance-matrix matrix print-matrix solve
34 backsolve eigenvalues eigenvectors accumulate cumsum combine
35 lowess))
37 (in-package :lisp-stat-descriptive-statistics)
39 ;;;;
40 ;;;; Basic Summary Statistics
41 ;;;;
43 (defun standard-deviation (x)
44 "Args: (x)
45 Returns the standard deviation of the elements x. Vector reducing."
46 (let ((n (count-elements x))
47 (r (- x (mean x))))
48 (sqrt (* (mean (* r r)) (/ n (- n 1))))))
51 ;; FIXME the following assume that we are using the vector based functions
52 (defun quantile (x p)
53 "Args: (x p)
54 Returns the P-th quantile(s) of sequence X. P can be a number or a sequence."
55 (let* ((x (sort-data x))
56 (n (length x))
57 (np (* p (- n 1)))
58 (low (floor np))
59 (high (ceiling np)))
60 (/ (+ (select x low) (select x high)) 2)))
62 (defun median (x)
63 "Args: (x)
64 Returns the median of the elements of X."
65 (quantile x 0.5))
67 (defun interquartile-range (x)
68 "Args: (number-data)
69 Returns the interquartile range of the elements of X."
70 (reduce #'- (quantile x '(0.75 0.25))))
72 (defun fivnum (x)
73 "Args: (number-data)
74 Returns the five number summary (min, 1st quartile, medinan, 3rd quartile,
75 max) of the elements X."
76 (quantile x '(0 .25 .5 .75 1)))
78 (defun covariance-matrix (&rest args)
79 "Args: (&rest args)
80 Returns the sample covariance matrix of the data columns in ARGS. ARGS may
81 consist of lists, vectors or matrices."
82 (let ((columns (apply #'append
83 (mapcar #'(lambda (x)
84 (if (matrixp x) (column-list x) (list x)))
85 args))))
86 (/ (cross-product (apply #'bind-columns
87 (- columns (mapcar #'mean columns))))
88 (- (length (car columns)) 1))))
90 ;;;; Sampling / Resampling
92 (defun sample (x ssize &optional replace)
93 "Args: (x n &optional (replace nil))
94 Returns a list of a random sample of size N from sequence X drawn with or
95 without replacement."
96 (check-sequence x)
97 (let ((n (length x))
98 (x (if (consp x) (coerce x 'vector) (copy-vector x)))
99 (result nil))
100 (if (< 0 n)
101 (dotimes (i ssize result)
102 (let ((j (if replace (random n) (+ i (random (- n i))))))
103 (setf result (cons (aref x j) result))
104 (unless replace ;; swap elements i and j
105 (let ((temp (aref x i)))
106 (setf (aref x i) (aref x j))
107 (setf (aref x j) temp))))))))
111 (defun mean (x)
112 "Args: (x)
113 Returns the mean of the elements x. Vector reducing."
114 (let ((mean 0.0)
115 (count 0.0))
116 (labels ((add-to-mean (x)
117 (let ((count+1 (+ count 1.0)))
118 (setf mean (+ (* (/ count count+1) mean) (* (/ count+1) x)))
119 (setf count count+1)))
120 (find-mean (x)
121 (if (numberp x)
122 (add-to-mean x)
123 (let ((seq (compound-data-seq x)))
124 (if (consp seq)
125 (dolist (x seq)
126 (if (numberp x) (add-to-mean x) (find-mean x)))
127 (let ((n (length seq)))
128 (dotimes (i n)
129 (declare (fixnum i))
130 (let ((x (aref seq i)))
131 (if (numberp x)
132 (add-to-mean x)
133 (find-mean x))))))))))
134 (find-mean x)
135 mean)))
137 ;;;;
138 ;;;; Linear Algebra Functions
139 ;;;;
141 (defun matrix (dim data)
142 "Args: (dim data)
143 returns a matrix of dimensions DIM initialized using sequence DATA
144 in row major order."
145 (let ((dim (coerce dim 'list))
146 (data (coerce data 'list)))
147 (make-array dim :initial-contents (split-list data (nth 1 dim)))))
149 (defun print-matrix (a &optional (stream *standard-output*))
150 "Args: (matrix &optional stream)
151 Prints MATRIX to STREAM in a nice form that is still machine readable"
152 (unless (matrixp a) (error "not a matrix - ~a" a))
153 (let ((size (min 15 (max (map-elements #'flatsize a)))))
154 (format stream "#2a(~%")
155 (dolist (x (row-list a))
156 (format stream " (")
157 (let ((n (length x)))
158 (dotimes (i n)
159 (let ((y (aref x i)))
160 (cond
161 ((integerp y) (format stream "~vd" size y))
162 ((floatp y) (format stream "~vg" size y))
163 (t (format stream "~va" size y))))
164 (if (< i (- n 1)) (format stream " "))))
165 (format stream ")~%"))
166 (format stream " )~%")
167 nil))
169 (defun solve (a b)
170 "Args: (a b)
171 Solves A x = B using LU decomposition and backsolving. B can be a sequence
172 or a matrix."
173 (let ((lu (lu-decomp a)))
174 (if (matrixp b)
175 (apply #'bind-columns
176 (mapcar #'(lambda (x) (lu-solve lu x)) (column-list b)))
177 (lu-solve lu b))))
179 (defun backsolve (a b)
180 "Args: (a b)
181 Solves A x = B by backsolving, assuming A is upper triangular. B must be a
182 sequence. For use with qr-decomp."
183 (let* ((n (length b))
184 (sol (make-array n)))
185 (dotimes (i n)
186 (let* ((k (- n i 1))
187 (val (elt b k)))
188 (dotimes (j i)
189 (let ((l (- n j 1)))
190 (setq val (- val (* (aref sol l) (aref a k l))))))
191 (setf (aref sol k) (/ val (aref a k k)))))
192 (if (listp b) (coerce sol 'list) sol)))
194 (defun eigenvalues (a)
195 "Args: (a)
196 Returns list of eigenvalues of square, symmetric matrix A"
197 (first (eigen a)))
199 (defun eigenvectors (a)
200 "Args: (a)
201 Returns list of eigenvectors of square, symmetric matrix A"
202 (second (eigen a)))
204 (defun accumulate (f s)
205 "Args: (f s)
206 Accumulates elements of sequence S using binary function F.
207 (accumulate #'+ x) returns the cumulative sum of x."
208 (let* ((result (list (elt s 0)))
209 (tail result))
210 (flet ((acc (dummy x)
211 (rplacd tail (list (funcall f (first tail) x)))
212 (setf tail (cdr tail))))
213 (reduce #'acc s))
214 (if (vectorp s) (coerce result 'vector) result)))
216 (defun cumsum (x)
217 "Args: (x)
218 Returns the cumulative sum of X."
219 (accumulate #'+ x))
221 (defun combine (&rest args)
222 "Args (&rest args)
223 Returns sequence of elements of all arguments."
224 (copy-seq (element-seq args)))
226 (defun lowess (x y &key (f .25) (steps 2) (delta -1) sorted)
227 "Args: (x y &key (f .25) (steps 2) delta sorted)
228 Returns (list X YS) with YS the LOWESS fit. F is the fraction of data used for
229 each point, STEPS is the number of robust iterations. Fits for points within
230 DELTA of each other are interpolated linearly. If the X values setting SORTED
231 to T speeds up the computation."
232 (let ((x (if sorted x (sort-data x)))
233 (y (if sorted y (select y (order x))))
234 (delta (if (> delta 0.0) delta (/ (- (max x) (min x)) 50))))
235 (list x)));; (|base-lowess| x y f steps delta))))