Commit 703c7591 authored by liam's avatar liam
Browse files

Added the rest of the statistics definitions.

Miscellaneous minor fixes.


git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3162 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent 4729f28f
......@@ -13,6 +13,18 @@ Features:
- This library is pure common lisp. No C is used; the only requirement
is that the GSL library must be loaded and accessible.
---- policies:
- All arrays and vectors that need to be passed to and from are
converted to or from GSL vectors or matrices; for GSL functions that
take/return raw C arrays (or pointers), a layer is wrapped
around them that do the conversion.
- Each defun-gsl corresponds to one C function. Where multiple
such functions may be combined into a single regular or generic
function.
======================
Needs work:
See Thu Apr 20 2006
......@@ -45,3 +57,7 @@ This means that the following functions have been put aside:
- #'polynomial-solve (requires vector_complex)
- The index functions in sorting that are not returning permutations,
because they need a vector_int.
It would be nice to use e.g. vector-sap for direct vectors in
e.g. SBCL.
......@@ -3,7 +3,7 @@
; description: Using GSL storage.
; date: Sun Mar 26 2006 - 16:32
; author: Liam M. Healy
; modified: Sun Dec 3 2006 - 10:24
; modified: Sun Dec 31 2006 - 17:34
;********************************************************
;;; $Id: $
......@@ -104,6 +104,7 @@
'((FIXNUM . "_int")
(SINGLE . "_float")
(DOUBLE . "")
;;(LONG-DOUBLE . "_long_double")
(COMPLEX . "_complex")))
(defmacro defdata
......
......@@ -3,7 +3,7 @@
; description: Definition of GSLL system
; date:
; author: Liam Healy
; modified: Sun Dec 3 2006 - 21:58
; modified: Sun Dec 31 2006 - 22:07
;********************************************************
;;; $Id: $
......@@ -133,4 +133,10 @@
(:module statistics
:depends-on (init data)
:components
((:file "mean-variance")))))
((:file "mean-variance")
(:file "absolute-deviation")
(:file "higher-moments")
(:file "autocorrelation")
(:file "covariance")
;; minimum and maximum values provided in vector.lisp
(:file "median-percentile")))))
......@@ -3,7 +3,7 @@
; description: Macros to interface GSL functions.
; date: Mon Mar 6 2006 - 22:35
; author: Liam M. Healy
; modified: Fri Jul 7 2006 - 23:24
; modified: Sun Dec 31 2006 - 22:49
;********************************************************
(in-package :gsl)
......@@ -46,6 +46,9 @@
;; (#:RET3504 (:DOUBLE (1+ KMAX))) (#:RET3505 (:DOUBLE (1+ KMAX)))
;; (#:RET3506 :DOUBLE) (#:RET3507 :DOUBLE))
(defun make-st (symbol type)
(list symbol type))
(defun st-symbol (decl)
(first decl))
......@@ -169,7 +172,8 @@
,@(when invalidate `((cl-invalidate ,@invalidate)))
,@(when (or null-pointer-info (eq c-return :pointer))
`((check-null-pointer ,cret-name
,@(or null-pointer-info '(:ENOMEM "No memory allocated")))))
,@(or null-pointer-info
'(:ENOMEM "No memory allocated")))))
,@after
(values
,@(case c-return
......@@ -189,4 +193,22 @@
(map-name ',(or index name) ,gsl-name)
,@(when export `((export ',name))))))
;;;;****************************************************************************
;;;; Macro defun-gsl
;;;;****************************************************************************
(defmacro defun-optionals
(name arglist no-optional optionals &optional documentation)
"Define a function with and without optional arguments."
(let* ((optpos (position '&optional arglist))
(mandatory-arglist (subseq arglist 0 optpos))
(optional-arglist (subseq arglist (1+ optpos))))
`(defun ,name ,arglist
,documentation
(if ,(first optional-arglist)
(,(intern
(concatenate 'string (string name) (string optionals)))
,@mandatory-arglist ,@optional-arglist)
(,(intern
(concatenate 'string (string name) (string no-optional)))
,@mandatory-arglist)))))
......@@ -3,7 +3,7 @@
; description: Quasi-Random Number Generators
; date: Sun Jul 16 2006 - 15:54
; author: Liam M. Healy
; modified: Mon Jul 17 2006 - 10:10
; modified: Sun Dec 31 2006 - 22:59
;********************************************************
;;; $Id: $
......@@ -13,8 +13,12 @@
((dimension :initarg :dimension :accessor qr-dimension))
(:documentation "A generator of quasi-random numbers."))
(defparameter *default-quasi-random-number-generator* nil)
(defun make-quasi-random-number-generator
(dimension &optional (type *niederreiter2*) (generator t))
(dimension
&optional (type *default-quasi-random-number-generator*)
(generator t))
"Make a random number generator; by default it is allocated on creation."
(let ((instance
(make-instance
......@@ -110,6 +114,9 @@
valid up to 12 dimensions."
"gsl_qrng_niederreiter_2")
(eval-when (:load-toplevel :execute)
(setf *default-quasi-random-number-generator* *niederreiter2*))
(def-rng-type *sobol*
"This generator uses the Sobol sequence described in Antonov, Saleev,
@cite{USSR Comput. Maths. Math. Phys.} 19, 252 (1980). It is valid up to
......
;********************************************************
; file: absolute-deviation.lisp
; description: Absolute deviation
; date: Sun Dec 31 2006 - 13:19
; author: Liam M. Healy
; modified: Sun Dec 31 2006 - 21:47
;********************************************************
;;; $Id: $
(in-package :gsl)
;;; To do: stride other than 1 when that information is availble from
;;; the vector.
(defun-gsl absolute-deviation-nom (data)
"gsl_stats_absdev"
(((gsl-array data) :pointer) (1 :int) ((dim0 data) :size))
:c-return :double
:index absolute-deviation
:export nil
:documentation
"The absolute deviation from the mean of @var{data}. The
absolute deviation from the mean is defined as
absdev = (1/N) \sum |x_i - \Hat\mu|
where @math{x_i} are the elements of the dataset @var{data}. The
absolute deviation from the mean provides a more robust measure of the
width of a distribution than the variance. This function computes the
mean of @var{data} via a call to #'mean.")
(defun-gsl absolute-deviation-m (data mean)
"gsl_stats_absdev_m"
(((gsl-array data) :pointer) (1 :int)
((dim0 data) :size) (mean :double))
:c-return :double
:index absolute-deviation
:export nil
:documentation
"The absolute deviation of the dataset @var{data}
relative to the given value of @var{mean},
absdev = (1/N) \sum |x_i - mean|.
This function is useful if you have already computed the mean of
@var{data} (and want to avoid recomputing it), or wish to calculate the
absolute deviation relative to another value (such as zero, or the
median).")
(export 'absolute-deviation)
(defun-optionals absolute-deviation (data &optional mean)
-nom -m
"The absolute deviation from the mean of @var{data}. The
absolute deviation from the mean is defined as
absdev = (1/N) \sum |x_i - \Hat\mu|
where @math{x_i} are the elements of the dataset @var{data}. The
absolute deviation from the mean provides a more robust measure of the
width of a distribution than the variance. This function computes the
mean of @var{data} via a call to #'mean.")
(defun-gsl weighted-absolute-deviation-nom (data weights)
"gsl_stats_wabsdev"
(((gsl-array weights) :pointer) (1 :int)
((gsl-array data) :pointer) (1 :int) ((dim0 data) :size))
:c-return :double
:index weighted-absolute-deviation
:export nil)
(defun-gsl weighted-absolute-deviation-m (data weights mean)
"gsl_stats_wabsdev_m"
(((gsl-array weights) :pointer) (1 :int)
((gsl-array data) :pointer) (1 :int)
((dim0 data) :size) (mean :double))
:c-return :double
:index weighted-absolute-deviation
:export nil)
(export 'weighted-absolute-deviation)
(defun-optionals weighted-absolute-deviation (data weights &optional mean)
-nom -m
"The weighted absolute deviation from the weighted
mean, defined as
absdev = (\sum w_i |x_i - \Hat\mu|) / (\sum w_i).")
(lisp-unit:define-test absolute-deviation
(lisp-unit:assert-equal
'("0.618000000000d+01" "0.664777777778d+01" "0.618000000000d+01")
(with-data (vec vector-double 3)
(setf (data vec) #(-3.21d0 1.0d0 12.8d0))
(with-data (weights vector-double 3)
(setf (data weights) #(3.0d0 1.0d0 2.0d0))
(let ((mean (mean vec)))
(lisp-unit:fp-sequence
(list
(absolute-deviation vec)
(weighted-absolute-deviation vec weights)
(absolute-deviation vec mean))))))))
;********************************************************
; file: autocorrelation.lisp
; description: Autocorrelation
; date: Sun Dec 31 2006 - 13:19
; author: Liam M. Healy
; modified: Sun Dec 31 2006 - 16:51
;********************************************************
;;; $Id: $
(in-package :gsl)
;;; To do: stride other than 1 when that information is availble from
;;; the vector.
(defun-gsl autocorrelation-nom (data)
"gsl_stats_lag1_autocorrelation"
(((gsl-array data) :pointer) (1 :int) ((dim0 data) :size))
:c-return :double
:index autocorrelation
:export nil)
(defun-gsl autocorrelation-m (data mean)
"gsl_stats_lag1_autocorrelation_m"
(((gsl-array data) :pointer) (1 :int)
((dim0 data) :size) (mean :double))
:c-return :double
:index autocorrelation
:export nil)
(export 'autocorrelation)
(defun-optionals autocorrelation (data &optional mean)
-nom -m
"The lag-1 autocorrelation of the dataset @var{data}.
a_1 = {\sum_{i = 1}^{n} (x_{i} - \Hat\mu) (x_{i-1} - \Hat\mu)
\over
\sum_{i = 1}^{n} (x_{i} - \Hat\mu) (x_{i} - \Hat\mu)}.")
(lisp-unit:define-test autocorrelation
(lisp-unit:assert-equal
'("-0.464636683425d-01" "-0.464636683425d-01")
(with-data (vec vector-double 3)
(setf (data vec) #(-3.21d0 1.0d0 12.8d0))
(let ((mean (mean vec)))
(lisp-unit:fp-sequence
(list
(autocorrelation vec)
(autocorrelation vec mean)))))))
;********************************************************
; file: covariance.lisp
; description: Covariance
; date: Sun Dec 31 2006 - 13:19
; author: Liam M. Healy
; modified: Sun Dec 31 2006 - 16:59
;********************************************************
;;; $Id: $
(in-package :gsl)
;;; To do: stride other than 1 when that information is availble from
;;; the vector.
(defun-gsl covariance-nom (data1 data2)
"gsl_stats_covariance"
(((gsl-array data1) :pointer) (1 :int)
((gsl-array data2) :pointer) (1 :int) ((dim0 data2) :size))
:c-return :double
:index covariance
:export nil)
(defun-gsl covariance-m (data1 data2 mean1 mean2)
"gsl_stats_covariance_m"
(((gsl-array data1) :pointer) (1 :int)
((gsl-array data2) :pointer) (1 :int) ((dim0 data2) :size)
(mean1 :double) (mean2 :double))
:c-return :double
:index covariance
:export nil)
(export 'covariance)
(defun-optionals covariance (data1 data2 &optional mean1 mean2)
-nom -m
"The covariance of the datasets @var{data1} and
@var{data2} which must both be of the same length @var{n}.
covar = {1 \over (n - 1)} \sum_{i = 1}^{n}
(x_{i} - \Hat x) (y_{i} - \Hat y).")
(lisp-unit:define-test covariance
(lisp-unit:assert-equal
'("-0.293000000000d+00" "-0.293000000000d+00")
(with-data (vec1 vector-double 3)
(setf (data vec1) #(-3.21d0 1.0d0 12.8d0))
(with-data (vec2 vector-double 3)
(setf (data vec2) #(1.15d0 -1.0d0 0.5d0))
(let ((mean1 (mean vec1))
(mean2 (mean vec2)))
(lisp-unit:fp-sequence
(list
(covariance vec1 vec2)
(covariance vec1 vec2 mean1 mean2))))))))
;********************************************************
; file: higher-moments.lisp
; description: Skewness and kurtosis.
; date: Sun Dec 31 2006 - 14:20
; author: Liam M. Healy
; modified: Sun Dec 31 2006 - 21:32
;********************************************************
;;; $Id: $
(in-package :gsl)
;;; To do: stride other than 1 when that information is availble from
;;; the vector.
(defun-gsl skewness-nomsd (data)
"gsl_stats_skew"
(((gsl-array data) :pointer) (1 :int) ((dim0 data) :size))
:c-return :double
:index skewness
:export nil
:documentation
"The skewness of @var{data}, defined as
skew = (1/N) \sum ((x_i - \Hat\mu)/\Hat\sigma)^3
where @math{x_i} are the elements of the dataset @var{data}.
The skewness measures the asymmetry of the tails of a distribution.")
(defun-gsl skewness-msd (data mean standard-deviation)
"gsl_stats_skew_m_sd"
(((gsl-array data) :pointer) (1 :int)
((dim0 data) :size) (mean :double) (standard-deviation :double))
:c-return :double
:index skewness
:export nil
:documentation
"The skewness of the dataset @var{data} using the
given values of the mean @var{mean} and standard deviation @var{sd},
skew = (1/N) \sum ((x_i - mean)/sd)^3
These functions are useful if you have already computed the mean and
standard deviation of @var{data} and want to avoid recomputing them.")
(export 'skewness)
(defun-optionals skewness
(data &optional mean standard-deviation)
-nomsd -msd
"The skewness of @var{data}, defined as
skew = (1/N) \sum ((x_i - \Hat\mu)/\Hat\sigma)^3
where @math{x_i} are the elements of the dataset @var{data}.
The skewness measures the asymmetry of the tails of a distribution.")
(defun-gsl kurtosis-nomsd (data)
"gsl_stats_kurtosis"
(((gsl-array data) :pointer) (1 :int) ((dim0 data) :size))
:c-return :double
:index kurtosis
:export nil)
(defun-gsl kurtosis-msd (data mean standard-deviation)
"gsl_stats_kurtosis_m_sd"
(((gsl-array data) :pointer) (1 :int)
((dim0 data) :size) (mean :double) (standard-deviation :double))
:c-return :double
:index kurtosis
:export nil)
(export 'kurtosis)
(defun-optionals kurtosis
(data &optional mean standard-deviation)
-nomsd -msd
"The kurtosis of @var{data}, defined as
kurtosis = ((1/N) \sum ((x_i - \Hat\mu)/\Hat\sigma)^4) - 3
The kurtosis measures how sharply peaked a distribution is,
relative to its width. The kurtosis is normalized to zero
for a gaussian distribution.")
(defun-gsl weighted-skewness-nomsd (data weights)
"gsl_stats_wskew"
(((gsl-array weights) :pointer) (1 :int)
((gsl-array data) :pointer) (1 :int) ((dim0 data) :size))
:c-return :double
:index weighted-skewness
:export nil)
(defun-gsl weighted-skewness-msd (data weights mean standard-deviation)
"gsl_stats_wskew_m_sd"
(((gsl-array data) :pointer) (1 :int)
((dim0 data) :size) (mean :double) (standard-deviation :double))
:c-return :double
:index weighted-skewness
:export nil)
(export 'weighted-skewness)
(defun-optionals weighted-skewness
(data weights &optional mean standard-deviation)
-nomsd -msd
"The weighted skewness of the dataset.
skew = (\sum w_i ((x_i - xbar)/\sigma)^3) / (\sum w_i).")
(defun-gsl weighted-kurtosis-nomsd (data weights)
"gsl_stats_wkurtosis"
(((gsl-array data) :pointer) (1 :int) ((dim0 data) :size))
:c-return :double
:index weighted-kurtosis
:export nil)
(defun-gsl weighted-kurtosis-msd (data weights mean standard-deviation)
"gsl_stats_wkurtosis_m_sd"
(((gsl-array data) :pointer) (1 :int)
((dim0 data) :size) (mean :double) (standard-deviation :double))
:c-return :double
:index weighted-kurtosis
:export nil)
(export 'weighted-kurtosis)
(defun-optionals weighted-kurtosis
(data weights &optional mean standard-deviation)
-nomsd -msd
"The weighted kurtosis of the dataset.
kurtosis = ((\sum w_i ((x_i - xbar)/sigma)^4) / (\sum w_i)) - 3.")
(lisp-unit:define-test higher-moments
(lisp-unit:assert-equal
'("0.276511898399d+00" "0.276511898399d+00"
"-0.233333333333d+01" "-0.233333333333d+01")
(with-data (vec vector-double 3)
(setf (data vec) #(-3.21d0 1.0d0 12.8d0))
(let* ((mean (mean vec))
(sd (standard-deviation vec mean)))
(lisp-unit:fp-sequence
(list
(skewness vec)
(skewness vec mean sd)
(kurtosis vec)
(kurtosis vec mean sd)))))))
......@@ -3,7 +3,7 @@
; description: Mean, standard deviation, and variance
; date: Sat Dec 2 2006 - 22:15
; author: Liam M. Healy
; modified: Sun Dec 3 2006 - 18:18
; modified: Sun Dec 31 2006 - 21:45
;********************************************************
;;; $Id: $
......@@ -15,12 +15,18 @@
(defmacro defun-gsl-stats (&rest args)
"A defun-gsl for stats of double, single, and fixnum."
(defun-gsl-all
;;'(double single fixnum long)
'(double single fixnum)
;;'(:double :float :int :long-double)
'(:double :float :int)
"stats"
'gsl-vector
args))
;;;;****************************************************************************
;;;; Mean and weighted mean
;;;;****************************************************************************
(defgeneric mean (gsl-vector)
(:documentation "The arithmetic mean of the vector.
The arithmetic mean, or sample mean, is denoted by
......@@ -32,7 +38,23 @@
(((gsl-array vector) :pointer) (1 :int) ((dim0 vector) :size))
:c-return :double)
(defgeneric variance (gsl-vector)
(defgeneric weighted-mean (gsl-vector weights)
(:documentation
"The weighted mean of the dataset, using the set of weights
The weighted mean is defined as
\Hat\mu = (\sum w_i x_i) / (\sum w_i)."))
(defun-gsl-stats weighted-mean ((vector gsl-vector) weights)
"gsl_stats_wmean"
(((gsl-array weights) :pointer) (1 :int)
((gsl-array vector) :pointer) (1 :int) ((dim0 vector) :size))
:c-return :double)
;;;;****************************************************************************
;;;; Variance
;;;;****************************************************************************
(defgeneric variance-nom (gsl-vector)
(:documentation "The estimated, or @dfn{sample}, variance of
@var{data}, a dataset of length @var{n} with stride @var{stride}. The
estimated variance is denoted by @math{\Hat\sigma^2} and is defined by
......@@ -45,11 +67,13 @@
This function computes the mean via a call to #'mean. If
you have already computed the mean then you can pass it directly to
#'variance-m"))
#'variance-m."))
(defun-gsl-stats variance ((vector gsl-vector))
(defun-gsl-stats variance-nom ((vector gsl-vector))
"gsl_stats_variance"
(((gsl-array vector) :pointer) (1 :int) ((dim0 vector) :size))
:index variance
:export nil
:c-return :double)
(defgeneric variance-m (gsl-vector mean)
......@@ -60,14 +84,74 @@
"gsl_stats_variance_m"
(((gsl-array vector) :pointer) (1 :int)
((dim0 vector) :size) (mean :double))
:index variance
:export nil
:c-return :double)
(defgeneric standard-deviation (gsl-vector)
(export 'variance)
(defun-optionals variance (data &optional mean)
-nom -m
"The estimated, or @dfn{sample}, variance of
@var{data}, a dataset of length @var{n} with stride @var{stride}. The
estimated variance is denoted by @math{\Hat\sigma^2} and is defined by
\Hat\sigma^2 = (1/(N-1)) \sum (x_i - \Hat\mu)^2
where @math{x_i} are the elements of the dataset @var{data}. Note that
the normalization factor of @math{1/(N-1)} results from the derivation
of @math{\Hat\sigma^2} as an unbiased estimator of the population
variance @math{\sigma^2}. For samples drawn from a gaussian distribution
the variance of @math{\Hat\sigma^2} itself is @math{2 \sigma^4 / N}.")
;;;;****************************************************************************
;;;; Weighted variance
;;;;****************************************************************************
(defgeneric weighted-variance-nom (gsl-vector weights)
(:documentation nil))
(defun-gsl-stats weighted-variance-nom ((vector gsl-vector) weights)
"gsl_stats_wvariance"
(((gsl-array weights) :pointer) (1 :int)
((gsl-array vector) :pointer) (1 :int) ((dim0 vector) :size))
:index weighted-variance
:export nil
:c-return :double)
(defgeneric weighted-variance-m (gsl-vector weights mean)
(:documentation "Compute the weighted variance with the mean known to
avoid its recomputation."))
(defun-gsl-stats weighted-variance-m ((vector gsl-vector) weights mean)
"gsl_stats_wvariance_m"
(((gsl-array weights) :pointer) (1 :int)
((gsl-array vector) :pointer) (1 :int)
((dim0 vector) :size) (mean :double))
:index variance
:export nil
:c-return :double)
(export 'weighted-variance)
(defun-optionals weighted-variance (data weights &optional mean)
-nom -m
"The estimated variance using weights
The estimated variance of a weighted dataset is defined as
\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
familiar @math{1/(N-1)} factor when there are @math{N} equal non-zero
weights.")
;;;;****************************************************************************
;;;; Standard deviation
;;;;****************************************************************************
(defgeneric standard-deviation-nom (gsl-vector)
(:documentation "The standard deviation, square root of the variance."))
(defun-gsl-stats standard-deviation ((vector gsl-vector))
(defun-gsl-stats standard-deviation-nom ((vector gsl-vector))
"gsl_stats_sd"
(((gsl-array vector) :pointer) (1 :int) ((dim0 vector) :size))
:index standard-deviation
:export nil
:c-return :double)
(defgeneric standard-deviation-m (gsl-vector mean)
......@@ -80,6 +164,45 @@
((dim0 vector) :size) (mean :double))
:c-return :double)