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 SD is simply the root.
70 (defgeneric variance
(x)
71 (:documentation
"compute the variance of the entity X, i.e. if a
72 scalar, vector, or (covariance) matrix.")
76 (sqrt (* (mean (* r r
)) (/ n
(- n
1))))))
77 (:method
((x vector-like
))
79 (let* ((negone nil
) ;; FIXME!
80 (negresid (axpy (mean x
) negone x
))) ; -mu + x = x - mu
81 (/ (loop for i from
0 to
(- (nelts x
) 1)
82 summing
(* (vref negresid i
)
84 (- (nelts negresid
) 1))))
85 (:method
((x matrix-like
))
86 (error "FIXME: define variance for matrices as covariance (?).")))
88 (defun covariance-matrix (&rest args
)
90 Returns the sample covariance matrix of the data columns in ARGS. ARGS may
91 consist of lists, vectors or matrices."
92 (let ((columns (apply #'append
94 (if (typep x
'matrix-like
)
98 (/ (cross-product (reduce #'bind2
99 (- columns
(mapcar #'mean columns
))))
100 (- (length (car columns
)) 1))))
102 (defgeneric standard-deviation
(x)
103 (:documentation
"Compute standard deivation of entity X.")
104 (:method
((x vector-like
)) (sqrt (variance x
)))
106 ;; if elements are not numeric, error, otherwise
108 (:method
((x matrix-like
))
109 (error "FIXME: define SD for matrix-like objects")))
111 (defun standard-deviation-fn (x)
112 "Function front-end to generic. Do we really need this?"
113 (standard-deviation x
))
115 (defun standard-deviation-old (x)
117 Returns the standard deviation of the elements x. Vector reducing.
119 FIXME AJR: This should be redone as square-root of the variance, and
120 defer to that structure for computation."
121 (let ((n (count-elements x
))
123 (sqrt (* (mean (* r r
)) (/ n
(- n
1))))))
125 (defgeneric quantile
(x p
&optional return-type
)
126 (:documentation
"Returns the P-th quantile(s) of sequence X.")
127 (:method
((x sequence
) (p real
) &optional
(return-type symbol
))
128 (let ((np (* p
(- (length x
) 1))))
129 (mean (list (aref x
(floor np
))
130 (aref x
(ceiling np
)))))) ;; aref work in general for lists too, or use nth?
131 (:method
((x vector-like
) (p real
) &optional
(return-type symbol
)) ;; average of sorted elements. Could store compile
132 (let ((np (* p
(- (nelts x
) 1))))
133 (mean (list (vref (sort x
) (floor np
))
134 (vref (sort x
) (ceiling np
))))))
135 (:method
((x sequence
) (p sequence
) &optional return-type
)
136 (error "FIXME: generalize. Basically, do mapcar or similar across a vector."))
137 (:method
((x vector-like
) (p sequence
) &optional return-type
)
138 (error "FIXME: generalize."))
139 (:method
((x vector-like
) (p vector-like
) &optional return-type
)
140 (error "FIXME: generalize.")))
142 ;;; things to build on top of quantiles...!
143 ;; need to incorporate a return-type if possible.
146 "Return median of X, using whatever the quantile generic function supports."
149 (defun interquartile-range (x)
150 "Returns the interquartile range of the elements of X."
151 (reduce #'-
(quantile x
(list 0.75 0.25))))
154 "Returns the five number summary (min, 1st quartile, medinan, 3rd quartile,
155 max) of the elements X."
156 (quantile x
(list 0 .25 .5 .75 1)))
158 ;;;; Sampling / Resampling
160 (defgeneric sample
(x n
&optional replace weights
)
161 (:documentation
"Draw a sample of size n from sequence x, with/out
162 replacement, with weights per element of x as described as a
163 probability mass distribution vector.")
164 (:method
((x sequence
) (n fixnum
)
165 &optional
(replace sequence
) (weights sequence
))
166 (sample-fn x n replace
))
167 (:method
((x vector-like
) (n fixnum
)
168 &optional
(replace vector-like
) (weights vector-like
))
169 (sample-fn x n replace
)))
171 (defun sample-fn (x ssize
&optional replace
)
172 "Args: (x n &optional (replace nil))
173 Returns a list of a random sample of size N from sequence X drawn with or
174 without replacement."
177 (x (if (consp x
) (coerce x
'vector
) (copy-vector x
)))
180 (dotimes (i ssize result
)
181 (let ((j (if replace
(random n
) (+ i
(random (- n i
))))))
182 (setf result
(cons (aref x j
) result
))
183 (unless replace
;; swap elements i and j
184 (let ((temp (aref x i
)))
185 (setf (aref x i
) (aref x j
))
186 (setf (aref x j
) temp
))))))))