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.
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.
15 (defpackage :lisp-stat-descriptive-statistics
19 :lisp-stat-compound-data
21 :lisp-stat-linalg-data
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
35 (in-package :lisp-stat-descriptive-statistics
)
38 ;;; Basic Summary Statistics
41 (defun standard-deviation (x)
43 Returns the standard deviation of the elements x. Vector reducing."
44 (let ((n (count-elements x
))
46 (sqrt (* (mean (* r r
)) (/ n
(- n
1))))))
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."))
53 (defmethod quantile ((x sequence
) (p sequence
))
54 (let* ((x (sort-data x
))
59 (/ (+ (select x low
) (select x high
)) 2)))
61 (defmethod quantile ((x sequence
) (p real
))
62 (let* ((x (sort-data x
))
67 (/ (+ (select x low
) (select x high
)) 2)))
70 ;;; things to build on top of quantiles...!
73 ;; Returns the median of the elements of X.
76 ;; (macroexpand '(median (list 1 2 3)))
77 ;; (median (list 1 2 3))
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))))
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)))
91 (defun covariance-matrix (&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
97 (if (matrixp x
) (column-list x
) (list x
)))
99 (/ (cross-product (apply #'bind-columns
100 (- columns
(mapcar #'mean columns
))))
101 (- (length (car columns
)) 1))))
103 ;;;; Sampling / Resampling
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."
111 (x (if (consp x
) (coerce x
'vector
) (copy-vector x
)))
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
))))))))
126 Returns the mean of the elements x. Vector reducing."
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)))
136 (let ((seq (compound-data-seq x
)))
139 (if (numberp x
) (add-to-mean x
) (find-mean x
)))
140 (let ((n (length seq
)))
143 (let ((x (aref seq i
)))
146 (find-mean x
))))))))))