mean-variance.lisp 9.72 KB
Newer Older
liam's avatar
liam committed
1 2
;; Mean, standard deviation, and variance    
;; Liam Healy, Sat Dec  2 2006 - 22:15
3
;; Time-stamp: <2008-12-26 12:21:30EST mean-variance.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5 6 7

(in-package :gsl)

Liam Healy's avatar
Liam Healy committed
8
;;; To do: stride other than 1 when that information is available from
9 10
;;; the vector.

11 12 13 14
;;;;****************************************************************************
;;;; Mean and weighted mean
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
15
(defmfun mean ((array both))
16
  ("gsl_stats" :type "_mean")
Liam Healy's avatar
Liam Healy committed
17
  (((c-pointer array) :pointer) (1 :int) ((total-size array) sizet))
18 19 20
  :definition :generic
  :element-types :no-complex
  :c-return :double
21
  :inputs (array)
22
  :documentation			; FDL
Liam Healy's avatar
Liam Healy committed
23
   "The arithmetic mean of the array.
24
   The arithmetic mean, or sample mean, is denoted by
25
   \Hat\mu and defined as \Hat\mu = (1/N) \sum x_i.  Returns a double-float.")
26

Liam Healy's avatar
Liam Healy committed
27
(defmfun weighted-mean ((array both) (weights both))
28 29
  ("gsl_stats" :type "_wmean")
  (((c-pointer weights) :pointer) (1 :int)
Liam Healy's avatar
Liam Healy committed
30 31
   ((c-pointer array) :pointer) (1 :int)
   ((total-size array) sizet))
32 33 34
  :definition :generic
  :element-types :float
  :c-return :element-c-type
35
  :inputs (array weights)
36 37
  :documentation			; FDL
  "The weighted mean of the dataset, using the set of weights
38
    The weighted mean is defined as
39
    \Hat\mu = (\sum w_i x_i) / (\sum w_i).")
40 41 42 43 44

;;;;****************************************************************************
;;;; Variance
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
45
(defmfun variance ((array both) &optional mean)
46 47
  (("gsl_stats" :type "_variance")
   ("gsl_stats" :type "_variance_m"))
Liam Healy's avatar
Liam Healy committed
48 49 50
  ((((c-pointer array) :pointer) (1 :int) ((total-size array) sizet))
   (((c-pointer array) :pointer) (1 :int)
   ((total-size array) sizet) (mean :double)))
51 52 53
  :definition :generic
  :element-types :no-complex
  :c-return :double
54
  :inputs (array)
55 56 57 58 59 60 61 62 63 64 65 66
  :documentation			; FDL
  "The estimated, or sample, variance of data.  The
   estimated variance is denoted by \Hat\sigma^2 and is defined by
   \Hat\sigma^2 = (1/(N-1)) \sum (x_i - \Hat\mu)^2
   where x_i are the elements of the dataset data.  Note that
   the normalization factor of 1/(N-1) results from the derivation
   of \Hat\sigma^2 as an unbiased estimator of the population
   variance \sigma^2.  For samples drawn from a gaussian distribution
   the variance of \Hat\sigma^2 itself is 2 \sigma^4 / N.
   If the mean value is known, it may be supplied which will use more
   efficient routines to compute the variance.")

Liam Healy's avatar
Liam Healy committed
67
(defmfun weighted-variance ((array both) (weights both) &optional mean)
Liam Healy's avatar
Liam Healy committed
68 69 70
  (("gsl_stats" :type "_wvariance")
   ("gsl_stats" :type "_wvariance_m"))
  ((((c-pointer weights) :pointer) (1 :int)
Liam Healy's avatar
Liam Healy committed
71 72
    ((c-pointer array) :pointer) (1 :int)
    ((total-size array) sizet))
Liam Healy's avatar
Liam Healy committed
73
   (((c-pointer weights) :pointer) (1 :int)
Liam Healy's avatar
Liam Healy committed
74 75
    ((c-pointer array) :pointer) (1 :int)
    ((total-size array) sizet) (mean :double)))
Liam Healy's avatar
Liam Healy committed
76 77 78
  :definition :generic
  :element-types :float
  :c-return :double
79
  :inputs (array weights)
Liam Healy's avatar
Liam Healy committed
80 81
  :documentation			; FDL
  "The estimated variance of a weighted dataset is defined as
82 83 84
   \Hat\sigma^2 = ((\sum w_i)/((\sum w_i)^2 - \sum (w_i^2))) 
                \sum w_i (x_i - \Hat\mu)^2
   Note that this expression reduces to an unweighted variance with the
Liam Healy's avatar
Liam Healy committed
85 86 87
   familiar 1/(N-1) factor when there are N equal non-zero weights.
   If the mean value is known, it may be supplied which will use more
   efficient routines to compute the variance.")
88 89 90 91 92

;;;;****************************************************************************
;;;; Standard deviation
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
93
(defmfun standard-deviation ((array both) &optional mean)
Liam Healy's avatar
Liam Healy committed
94 95
  (("gsl_stats" :type "_sd")
   ("gsl_stats" :type "_sd_m"))
Liam Healy's avatar
Liam Healy committed
96 97 98
  ((((c-pointer array) :pointer) (1 :int) ((total-size array) sizet))
   (((c-pointer array) :pointer) (1 :int)
   ((total-size array) sizet) (mean :double)))
Liam Healy's avatar
Liam Healy committed
99 100 101
  :definition :generic
  :element-types :no-complex
  :c-return :double
102
  :inputs (array)
Liam Healy's avatar
Liam Healy committed
103 104 105 106
  :documentation			; FDL
  "The standard deviation, square root of the variance.
   If the mean value is known, it may be supplied which will use more
   efficient routines to compute the variance.")
107

Liam Healy's avatar
Liam Healy committed
108
(defmfun weighted-standard-deviation ((array both) (weights both) &optional mean)
Liam Healy's avatar
Liam Healy committed
109 110 111
  (("gsl_stats" :type "_wsd")
   ("gsl_stats" :type "_wsd_m"))
  ((((c-pointer weights) :pointer) (1 :int)
Liam Healy's avatar
Liam Healy committed
112 113
    ((c-pointer array) :pointer) (1 :int)
    ((total-size array) sizet))
Liam Healy's avatar
Liam Healy committed
114
   (((c-pointer weights) :pointer) (1 :int)
Liam Healy's avatar
Liam Healy committed
115 116
    ((c-pointer array) :pointer) (1 :int)
    ((total-size array) sizet) (mean :double)))
Liam Healy's avatar
Liam Healy committed
117 118 119
  :definition :generic
  :element-types :float
  :c-return :double
120
  :inputs (array weights)
Liam Healy's avatar
Liam Healy committed
121 122 123 124
  :documentation			; FDL
  "The weighted standard deviation, square root of the variance.
   If the mean value is known, it may be supplied which will use more
   efficient routines to compute the variance.")
125 126

;;;;****************************************************************************
Liam Healy's avatar
Liam Healy committed
127
;;;; With fixed mean
128 129
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
130
(defmfun variance-with-fixed-mean ((array both) mean)
Liam Healy's avatar
Liam Healy committed
131
  ("gsl_stats" :type "_variance_with_fixed_mean")
Liam Healy's avatar
Liam Healy committed
132 133
  (((c-pointer array) :pointer) (1 :int)
   ((total-size array) sizet) (mean :double))
Liam Healy's avatar
Liam Healy committed
134 135 136
  :definition :generic
  :element-types :no-complex
  :c-return :double
137
  :inputs (array)
Liam Healy's avatar
Liam Healy committed
138 139
  :documentation			; FDL
  "An unbiased estimate of the variance of
Liam Healy's avatar
Liam Healy committed
140
    data when the population mean of the underlying
141 142 143
    distribution is known a priori.  In this case the estimator for
    the variance uses the factor 1/N and the sample mean
    \Hat\mu is replaced by the known population mean \mu,
Liam Healy's avatar
Liam Healy committed
144
    \Hat\sigma^2 = (1/N) \sum (x_i - \mu)^2.")
145

Liam Healy's avatar
Liam Healy committed
146
(defmfun standard-deviation-with-fixed-mean ((array both) mean)
Liam Healy's avatar
Liam Healy committed
147
  ("gsl_stats" :type "_sd_with_fixed_mean")
Liam Healy's avatar
Liam Healy committed
148 149
  (((c-pointer array) :pointer) (1 :int)
   ((total-size array) sizet) (mean :double))
Liam Healy's avatar
Liam Healy committed
150 151 152
  :definition :generic
  :element-types :no-complex
  :c-return :double
153
  :inputs (array)
Liam Healy's avatar
Liam Healy committed
154
  :documentation			; FDL
155 156
   "The standard deviation of data for a fixed population
    mean.  The result is the square root of the
Liam Healy's avatar
Liam Healy committed
157
    corresponding variance function.")
158

159
;;;;****************************************************************************
Liam Healy's avatar
Liam Healy committed
160
;;;; Weighted with fixed mean
161 162
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
163
(defmfun weighted-variance-with-fixed-mean
Liam Healy's avatar
Liam Healy committed
164
    ((array both) (weights both) mean)
Liam Healy's avatar
Liam Healy committed
165 166
  ("gsl_stats" :type "_wvariance_with_fixed_mean")
  (((c-pointer weights) :pointer) (1 :int)
Liam Healy's avatar
Liam Healy committed
167 168
   ((c-pointer array) :pointer) (1 :int)
   ((total-size array) sizet) (mean :double))
Liam Healy's avatar
Liam Healy committed
169 170 171
  :definition :generic
  :element-types :float
  :c-return :double
172
  :inputs (array weights)
Liam Healy's avatar
Liam Healy committed
173 174
  :documentation			; FDL
  "An unbiased estimate of the variance of weighted
175 176 177 178
    dataset when the population mean of the underlying
    distribution is known a priori.  In this case the estimator for
    the variance replaces the sample mean \Hat\mu by the known
    population mean \mu,
Liam Healy's avatar
Liam Healy committed
179
    \Hat\sigma^2 = (\sum w_i (x_i - \mu)^2) / (\sum w_i).")
180

Liam Healy's avatar
Liam Healy committed
181
(defmfun weighted-standard-deviation-with-fixed-mean
Liam Healy's avatar
Liam Healy committed
182
    ((array both) (weights both) mean)
Liam Healy's avatar
Liam Healy committed
183 184
  ("gsl_stats" :type "_wsd_with_fixed_mean")
  (((c-pointer weights) :pointer) (1 :int)
Liam Healy's avatar
Liam Healy committed
185 186
   ((c-pointer array) :pointer) (1 :int)
   ((total-size array) sizet) (mean :double))
Liam Healy's avatar
Liam Healy committed
187 188 189
  :definition :generic
  :element-types :float
  :c-return :double
190
  :inputs (array weights)
Liam Healy's avatar
Liam Healy committed
191 192 193
  :documentation			; FDL
  "The square root of the corresponding variance
   function #'weighted-variance-with-fixed-mean.")
194 195 196 197 198

;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
199
(generate-all-array-tests vector-mean :no-complex
200
 (let ((v1 (array-default 8)))
Liam Healy's avatar
Liam Healy committed
201 202 203
   (mean v1)))

(generate-all-array-tests matrix-mean :no-complex
204
 (let ((m1 (array-default '(3 3))))
Liam Healy's avatar
Liam Healy committed
205 206 207
   (mean m1)))

(generate-all-array-tests vector-variance :no-complex
208
 (let ((v1 (array-default 8)))
Liam Healy's avatar
Liam Healy committed
209 210 211
   (variance v1)))

(generate-all-array-tests matrix-variance :no-complex
212
 (let ((m1 (array-default '(3 3))))
Liam Healy's avatar
Liam Healy committed
213 214 215
   (variance m1)))

(generate-all-array-tests vector-variance-with-mean :no-complex
216
 (let ((v1 (array-default 8)))
Liam Healy's avatar
Liam Healy committed
217 218 219
   (variance v1 (mean v1))))

(generate-all-array-tests matrix-variance-with-mean :no-complex
220
 (let ((m1 (array-default '(3 3))))
Liam Healy's avatar
Liam Healy committed
221 222 223
   (variance m1 (mean m1))))

(generate-all-array-tests vector-standard-deviation :no-complex
224
 (let ((v1 (array-default 8)))
Liam Healy's avatar
Liam Healy committed
225 226 227
   (standard-deviation v1)))

(generate-all-array-tests matrix-standard-deviation :no-complex
228
 (let ((m1 (array-default '(3 3))))
Liam Healy's avatar
Liam Healy committed
229 230 231
   (standard-deviation m1)))

(generate-all-array-tests vector-standard-deviation-with-mean :no-complex
232
 (let ((v1 (array-default 8)))
Liam Healy's avatar
Liam Healy committed
233 234 235
   (standard-deviation v1 (mean v1))))

(generate-all-array-tests matrix-standard-deviation-with-mean :no-complex
236
 (let ((m1 (array-default '(3 3))))
Liam Healy's avatar
Liam Healy committed
237 238 239
   (standard-deviation m1 (mean m1))))

(generate-all-array-tests vector-variance-with-fixed-mean :no-complex
240
 (let ((v1 (array-default 8)))
Liam Healy's avatar
Liam Healy committed
241 242 243
   (variance-with-fixed-mean v1 (mean v1))))

(generate-all-array-tests matrix-variance-with-fixed-mean :no-complex
244
 (let ((m1 (array-default '(3 3))))
Liam Healy's avatar
Liam Healy committed
245 246 247
   (variance-with-fixed-mean m1 (mean m1))))

(generate-all-array-tests vector-standard-deviation-with-fixed-mean :no-complex
248
 (let ((v1 (array-default 8)))
Liam Healy's avatar
Liam Healy committed
249 250 251
   (standard-deviation-with-fixed-mean v1 (mean v1))))

(generate-all-array-tests matrix-standard-deviation-with-fixed-mean :no-complex
252
 (let ((m1 (array-default '(3 3))))
Liam Healy's avatar
Liam Healy committed
253 254 255 256 257 258
   (standard-deviation-with-fixed-mean m1 (mean m1))))

#|
;;; Weighted mean seems to be in error in GSL
;;; Hold off on all weghted tests until this is resolved.
(generate-all-array-tests weighted-mean :float
259
 (let ((v1 (array-default 8))
Liam Healy's avatar
Liam Healy committed
260 261 262 263 264 265
	(v2 (array-default 8)))
   (loop for i below (first (dimensions v2))
      do (setf (maref v2 i) (abs (maref v2 i))))
   (weighted-mean v1 v2)))
|#