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.
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 :lisp-stat-descriptive-statistics
)
16 ;;; Basic Summary Statistics
20 (:documentation
"compute the mean of lists, vectors, various objects")
24 (:method
((x vector-like
))
25 ;; (defparameter *x* (make-vector 5 :initial-contents '((1d0 2d0 3d0 4d0 5d0))))
26 (/ (loop for i from
0 to
(- (nelts x
) 1)
29 #|
;; a more traditional versino of the above...
30 (:method
((x vector-like
))
32 (type (if (col-vector-p x
) :column
:row
)))
34 (ecase type
(:row
1) (:column n
))
35 (ecase type
(:row n
) (:column
1))))
42 Returns the mean of the elements x. Vector reducing.
44 FIXME: Why is this so complex? When figure out, put into generic."
47 (labels ((add-to-mean (x)
48 (let ((count+1 (+ count
1.0)))
49 (setf mean
(+ (* (/ count count
+1) mean
) (* (/ count
+1) x
)))
50 (setf count count
+1)))
54 (let ((seq (compound-data-seq x
)))
57 (if (numberp x
) (add-to-mean x
) (find-mean x
)))
58 (let ((n (length seq
)))
61 (let ((x (aref seq i
)))
64 (find-mean x
))))))))))
69 ;; We do the variance, since the STANDARD-DEVIATION is simply the root
70 ;; (for interesting definitions of "simply"), and VARIANCE is a bit
72 (defgeneric variance
(x)
73 (:documentation
"compute the variance of the entity X, i.e. if a
74 scalar, vector, or (covariance) matrix.")
78 (sqrt (* (mean (* r r
)) (/ n
(- n
1))))))
79 (:method
((x vector-like
))
81 (let* ((negone nil
) ;; FIXME!
82 (negresid (axpy (mean x
) negone x
))) ; -mu + x = x - mu
83 (/ (loop for i from
0 to
(- (nelts x
) 1)
84 summing
(* (vref negresid i
)
86 (- (nelts negresid
) 1))))
87 (:method
((x matrix-like
))
88 (error "FIXME: define variance for matrices as covariance (?).")))
90 (defun covariance-matrix (&rest args
)
92 Returns the sample covariance matrix of the data columns in ARGS. ARGS may
93 consist of lists, vectors or matrices."
94 (let ((columns (apply #'append
96 (if (typep x
'matrix-like
)
100 (/ (cross-product (reduce #'bind2
101 (- columns
(mapcar #'mean columns
))))
102 (- (length (car columns
)) 1))))
104 (defgeneric standard-deviation
(x)
105 (:documentation
"Compute standard deivation of entity X.")
106 (:method
((x vector-like
)) (sqrt (variance x
)))
108 ;; if elements are not numeric, error, otherwise
110 (:method
((x matrix-like
))
111 (error "FIXME: define SD for matrix-like objects")))
113 (defun standard-deviation-fn (x)
114 "Function front-end to generic. Do we really need this?"
115 (standard-deviation x
))
117 (defun standard-deviation-old (x)
119 Returns the standard deviation of the elements x. Vector reducing.
121 FIXME AJR: This should be redone as square-root of the variance, and
122 defer to that structure for computation."
123 (let ((n (count-elements x
))
125 (sqrt (* (mean (* r r
)) (/ n
(- n
1))))))
127 (defgeneric quantile
(x p
)
128 (:documentation
"Returns the P-th quantile(s) of sequence X as a
129 native lisp list (which could be fed into a vector-like or similar
130 for the final format).")
131 (:method
((x sequence
) (p real
))
132 (let ((np (* p
(- (length x
) 1))))
133 (mean (list (aref x
(floor np
))
134 (aref x
(ceiling np
)))))) ;; aref work in general for lists too, or use nth?
135 (:method
((x vector-like
) (p real
)) ;; average of sorted elements. Could store compile
136 (let ((np (* p
(- (nelts x
) 1))))
137 (mean (list (vref (sort x
) (floor np
))
138 (vref (sort x
) (ceiling np
))))))
139 (:method
((x sequence
) (p sequence
))
140 (error "FIXME: generalize. Basically, do mapcar or similar across a vector."))
141 (:method
((x vector-like
) (p sequence
))
142 (error "FIXME: generalize."))
143 (:method
((x vector-like
) (p vector-like
))
144 (error "FIXME: generalize.")))
146 ;;; things to build on top of quantile.
148 ;; Do as generics if possible, keeping in mind that quantile returns a
149 ;; raw lisp list. Nearly all of these use native lisp structures for speed,
150 ;; such as raw lists for single and vector values, and arrays
151 ;; (simple?) for matrix-values and array-valued results.
154 "Return median of X, using whatever the quantile generic function supports."
157 (defun interquartile-range (x)
158 "Returns the interquartile range of the elements of X."
159 (reduce #'-
(quantile x
(list 0.75 0.25))))
162 "Returns the five number summary (min, 1st quartile, medinan, 3rd quartile,
163 max) of the elements X."
164 (quantile x
(list 0 .25 .5 .75 1)))
166 ;;;; Sampling / Resampling
168 (defgeneric sample
(x n
&optional replace weights
)
169 (:documentation
"Draw a sample of size n from sequence x, with/out
170 replacement, with weights per element of x as described as a
171 probability mass distribution vector. Length of weights should
172 match length of X. Unclear if (= X (uniq X))? (possible issue).
173 Return value shoud match the input values.")
174 (:method
((x sequence
) (n fixnum
)
175 &optional replace weights
)
176 (sample-fn x n replace
))
177 (:method
((x vector-like
) (n fixnum
)
178 &optional replace weights
)
179 (error "SAMPLE for vector-likes not implemented yet.")))
181 (defun sample-fn (x ssize
&optional replace
)
182 "Args: (x n &optional (replace nil))
183 Returns a list of a random sample of size N from sequence X drawn with or
184 without replacement."
187 (x (if (consp x
) (coerce x
'vector
) (copy-vector x
)))
190 (dotimes (i ssize result
)
191 (let ((j (if replace
(random n
) (+ i
(random (- n i
))))))
192 (setf result
(cons (aref x j
) result
))
193 (unless replace
;; swap elements i and j
194 (let ((temp (aref x i
)))
195 (setf (aref x i
) (aref x j
))
196 (setf (aref x j
) temp
))))))))