mean-variance.lisp 10.6 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: <2011-01-10 18:16:17EST mean-variance.lisp>
4
;;
5
;; Copyright 2006, 2007, 2008, 2009, 2011 Liam M. Healy
6 7 8 9 10 11 12 13 14 15 16 17 18 19
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <http://www.gnu.org/licenses/>.
20 21 22

(in-package :gsl)

Liam Healy's avatar
Liam Healy committed
23
;;; To do: stride other than 1 when that information is available from
24 25
;;; the vector.

26 27 28 29
;;;;****************************************************************************
;;;; Mean and weighted mean
;;;;****************************************************************************

30
(defmfun mean ((array both))
31
  ("gsl_stats" :type "_mean")
32
  (((grid:foreign-pointer array) :pointer) (1 :int) ((size array) sizet))
33 34 35
  :definition :generic
  :element-types :no-complex
  :c-return :double
36
  :inputs (array)
37
  :documentation			; FDL
38
   "The arithmetic mean of the array.
39
   The arithmetic mean, or sample mean, is denoted by
40
   \Hat\mu and defined as \Hat\mu = (1/N) \sum x_i.  Returns a double-float.")
41

42
(defmfun weighted-mean ((array both) (weights both))
43
  ("gsl_stats" :type "_wmean")
44 45
  (((grid:foreign-pointer weights) :pointer) (1 :int)
   ((grid:foreign-pointer array) :pointer) (1 :int)
46
   ((size array) sizet))
47 48 49
  :definition :generic
  :element-types :float
  :c-return :element-c-type
50
  :inputs (array weights)
51 52
  :documentation			; FDL
  "The weighted mean of the dataset, using the set of weights
53
    The weighted mean is defined as
54
    \Hat\mu = (\sum w_i x_i) / (\sum w_i).")
55 56 57 58 59

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

60
(defmfun variance ((array both) &optional mean)
61 62
  (("gsl_stats" :type "_variance")
   ("gsl_stats" :type "_variance_m"))
63 64
  ((((grid:foreign-pointer array) :pointer) (1 :int) ((size array) sizet))
   (((grid:foreign-pointer array) :pointer) (1 :int)
65
   ((size array) sizet) (mean :double)))
66 67 68
  :definition :generic
  :element-types :no-complex
  :c-return :double
69
  :inputs (array)
70 71 72 73 74 75 76 77 78 79 80 81
  :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.")

82
(defmfun weighted-variance ((array both) (weights both) &optional mean)
Liam Healy's avatar
Liam Healy committed
83 84
  (("gsl_stats" :type "_wvariance")
   ("gsl_stats" :type "_wvariance_m"))
85 86
  ((((grid:foreign-pointer weights) :pointer) (1 :int)
    ((grid:foreign-pointer array) :pointer) (1 :int)
87
    ((size array) sizet))
88 89
   (((grid:foreign-pointer weights) :pointer) (1 :int)
    ((grid:foreign-pointer array) :pointer) (1 :int)
90
    ((size array) sizet) (mean :double)))
Liam Healy's avatar
Liam Healy committed
91 92 93
  :definition :generic
  :element-types :float
  :c-return :double
94
  :inputs (array weights)
Liam Healy's avatar
Liam Healy committed
95 96
  :documentation			; FDL
  "The estimated variance of a weighted dataset is defined as
97 98 99
   \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
100 101 102
   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.")
103 104 105 106 107

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

108
(defmfun standard-deviation ((array both) &optional mean)
Liam Healy's avatar
Liam Healy committed
109 110
  (("gsl_stats" :type "_sd")
   ("gsl_stats" :type "_sd_m"))
111 112
  ((((grid:foreign-pointer array) :pointer) (1 :int) ((size array) sizet))
   (((grid:foreign-pointer array) :pointer) (1 :int)
113
   ((size array) sizet) (mean :double)))
Liam Healy's avatar
Liam Healy committed
114 115 116
  :definition :generic
  :element-types :no-complex
  :c-return :double
117
  :inputs (array)
Liam Healy's avatar
Liam Healy committed
118 119 120 121
  :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.")
122

123
(defmfun weighted-standard-deviation ((array both) (weights both) &optional mean)
Liam Healy's avatar
Liam Healy committed
124 125
  (("gsl_stats" :type "_wsd")
   ("gsl_stats" :type "_wsd_m"))
126 127
  ((((grid:foreign-pointer weights) :pointer) (1 :int)
    ((grid:foreign-pointer array) :pointer) (1 :int)
128
    ((size array) sizet))
129 130
   (((grid:foreign-pointer weights) :pointer) (1 :int)
    ((grid:foreign-pointer array) :pointer) (1 :int)
131
    ((size array) sizet) (mean :double)))
Liam Healy's avatar
Liam Healy committed
132 133 134
  :definition :generic
  :element-types :float
  :c-return :double
135
  :inputs (array weights)
Liam Healy's avatar
Liam Healy committed
136 137 138 139
  :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.")
140 141

;;;;****************************************************************************
Liam Healy's avatar
Liam Healy committed
142
;;;; With fixed mean
143 144
;;;;****************************************************************************

145
(defmfun variance-with-fixed-mean ((array both) mean)
Liam Healy's avatar
Liam Healy committed
146
  ("gsl_stats" :type "_variance_with_fixed_mean")
147
  (((grid:foreign-pointer array) :pointer) (1 :int)
148
   ((size array) sizet) (mean :double))
Liam Healy's avatar
Liam Healy committed
149 150 151
  :definition :generic
  :element-types :no-complex
  :c-return :double
152
  :inputs (array)
Liam Healy's avatar
Liam Healy committed
153 154
  :documentation			; FDL
  "An unbiased estimate of the variance of
155
    data when the population mean of the underlying
156 157 158
    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
159
    \Hat\sigma^2 = (1/N) \sum (x_i - \mu)^2.")
160

161
(defmfun standard-deviation-with-fixed-mean ((array both) mean)
Liam Healy's avatar
Liam Healy committed
162
  ("gsl_stats" :type "_sd_with_fixed_mean")
163
  (((grid:foreign-pointer array) :pointer) (1 :int)
164
   ((size array) sizet) (mean :double))
Liam Healy's avatar
Liam Healy committed
165 166 167
  :definition :generic
  :element-types :no-complex
  :c-return :double
168
  :inputs (array)
Liam Healy's avatar
Liam Healy committed
169
  :documentation			; FDL
170 171
   "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
172
    corresponding variance function.")
173

174
;;;;****************************************************************************
Liam Healy's avatar
Liam Healy committed
175
;;;; Weighted with fixed mean
176 177
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
178
(defmfun weighted-variance-with-fixed-mean
179
    ((array both) (weights both) mean)
Liam Healy's avatar
Liam Healy committed
180
  ("gsl_stats" :type "_wvariance_with_fixed_mean")
181 182
  (((grid:foreign-pointer weights) :pointer) (1 :int)
   ((grid:foreign-pointer array) :pointer) (1 :int)
183
   ((size array) sizet) (mean :double))
Liam Healy's avatar
Liam Healy committed
184 185 186
  :definition :generic
  :element-types :float
  :c-return :double
187
  :inputs (array weights)
Liam Healy's avatar
Liam Healy committed
188 189
  :documentation			; FDL
  "An unbiased estimate of the variance of weighted
190 191 192 193
    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
194
    \Hat\sigma^2 = (\sum w_i (x_i - \mu)^2) / (\sum w_i).")
195

Liam Healy's avatar
Liam Healy committed
196
(defmfun weighted-standard-deviation-with-fixed-mean
197
    ((array both) (weights both) mean)
Liam Healy's avatar
Liam Healy committed
198
  ("gsl_stats" :type "_wsd_with_fixed_mean")
199 200
  (((grid:foreign-pointer weights) :pointer) (1 :int)
   ((grid:foreign-pointer array) :pointer) (1 :int)
201
   ((size array) sizet) (mean :double))
Liam Healy's avatar
Liam Healy committed
202 203 204
  :definition :generic
  :element-types :float
  :c-return :double
205
  :inputs (array weights)
Liam Healy's avatar
Liam Healy committed
206 207 208
  :documentation			; FDL
  "The square root of the corresponding variance
   function #'weighted-variance-with-fixed-mean.")
209 210 211 212 213

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

214
(generate-all-array-tests vector-mean :no-complex
215
 (let ((v1 (array-default 8)))
216 217 218
   (mean v1)))

(generate-all-array-tests matrix-mean :no-complex
219
 (let ((m1 (array-default '(3 3))))
220 221 222
   (mean m1)))

(generate-all-array-tests vector-variance :no-complex
223
 (let ((v1 (array-default 8)))
224 225 226
   (variance v1)))

(generate-all-array-tests matrix-variance :no-complex
227
 (let ((m1 (array-default '(3 3))))
228 229 230
   (variance m1)))

(generate-all-array-tests vector-variance-with-mean :no-complex
231
 (let ((v1 (array-default 8)))
232 233 234
   (variance v1 (mean v1))))

(generate-all-array-tests matrix-variance-with-mean :no-complex
235
 (let ((m1 (array-default '(3 3))))
236 237 238
   (variance m1 (mean m1))))

(generate-all-array-tests vector-standard-deviation :no-complex
239
 (let ((v1 (array-default 8)))
240 241 242
   (standard-deviation v1)))

(generate-all-array-tests matrix-standard-deviation :no-complex
243
 (let ((m1 (array-default '(3 3))))
244 245 246
   (standard-deviation m1)))

(generate-all-array-tests vector-standard-deviation-with-mean :no-complex
247
 (let ((v1 (array-default 8)))
248 249 250
   (standard-deviation v1 (mean v1))))

(generate-all-array-tests matrix-standard-deviation-with-mean :no-complex
251
 (let ((m1 (array-default '(3 3))))
252 253 254
   (standard-deviation m1 (mean m1))))

(generate-all-array-tests vector-variance-with-fixed-mean :no-complex
255
 (let ((v1 (array-default 8)))
256 257 258
   (variance-with-fixed-mean v1 (mean v1))))

(generate-all-array-tests matrix-variance-with-fixed-mean :no-complex
259
 (let ((m1 (array-default '(3 3))))
260 261 262
   (variance-with-fixed-mean m1 (mean m1))))

(generate-all-array-tests vector-standard-deviation-with-fixed-mean :no-complex
263
 (let ((v1 (array-default 8)))
264 265 266
   (standard-deviation-with-fixed-mean v1 (mean v1))))

(generate-all-array-tests matrix-standard-deviation-with-fixed-mean :no-complex
267
 (let ((m1 (array-default '(3 3))))
268 269 270 271 272 273
   (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
274
 (let ((v1 (array-default 8))
275
	(v2 (array-default 8)))
276
   (loop for i below (first (grid:dimensions v2))
277
      do (setf (grid:gref v2 i) (abs (grid:gref v2 i))))
278 279 280
   (weighted-mean v1 v2)))
|#