higher-moments.lisp 3.35 KB
Newer Older
1 2
;; Skewness and kurtosis.
;; Liam Healy, Sun Dec 31 2006 - 14:20
3
;; Time-stamp: <2008-12-26 12:20:08EST higher-moments.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5 6 7 8 9

(in-package :gsl)
;;; To do: stride other than 1 when that information is availble from
;;; the vector.

Liam Healy's avatar
Liam Healy committed
10 11 12 13 14 15 16 17
(defmfun skewness ((data vector) &optional mean standard-deviation)
  (("gsl_stats" :type "_skew")
   ("gsl_stats" :type "_skew_m_sd"))
  ((((c-pointer data) :pointer) (1 :int) ((dim0 data) sizet))
   (((c-pointer data) :pointer) (1 :int) ((dim0 data) sizet)
    (mean :double) (standard-deviation :double)))
  :definition :generic
  :element-types :no-complex
18
  :c-return :double
19
  :inputs (data)
20
  :documentation			; FDL
Liam Healy's avatar
Liam Healy committed
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
  "The skewness of data, defined as skew = (1/N) \sum ((x_i -
  \Hat\mu)/\Hat\sigma)^3 where x_i are the elements of the dataset
  data.  The skewness measures the asymmetry of the tails of a
  distribution.  If mean and standard deviation are supplied, compute
  skewness of the dataset data using the given values skew = (1/N)
  \sum ((x_i - mean)/sd)^3.  This is useful if you have
  already computed the mean and standard deviation of data and want to
  avoid recomputing them.")

(defmfun kurtosis ((data vector) &optional mean standard-deviation)
  (("gsl_stats" :type "_kurtosis")
   ("gsl_stats" :type "_kurtosis_m_sd"))
  ((((c-pointer data) :pointer) (1 :int) ((dim0 data) sizet))
   (((c-pointer data) :pointer) (1 :int) ((dim0 data) sizet)
    (mean :double) (standard-deviation :double)))
  :definition :generic
  :element-types :no-complex
38
  :c-return :double
39
  :inputs (data)
40 41
  :documentation			; FDL
  "The kurtosis of data defined as
42 43 44 45 46
   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.")

Liam Healy's avatar
Liam Healy committed
47 48 49 50 51 52 53 54 55 56 57 58 59
(defmfun weighted-skewness
    ((data vector) (weights vector) &optional mean standard-deviation)
  (("gsl_stats" :type "_wskew")
   ("gsl_stats" :type "_wskew_m_sd"))
  ((((c-pointer weights) :pointer) (1 :int)
    ((c-pointer data) :pointer) (1 :int)
    ((dim0 data) sizet))
   (((c-pointer weights) :pointer) (1 :int)
    ((c-pointer data) :pointer) (1 :int)
    ((dim0 data) sizet)
    (mean :double) (standard-deviation :double)))
  :definition :generic
  :element-types :float
60
  :c-return :double
61
  :inputs (data weights)
Liam Healy's avatar
Liam Healy committed
62
  :documentation			; FDL
63 64 65
  "The weighted skewness of the dataset.
   skew = (\sum w_i ((x_i - xbar)/\sigma)^3) / (\sum w_i).")

Liam Healy's avatar
Liam Healy committed
66 67 68 69 70 71 72 73 74 75 76 77 78
(defmfun weighted-kurtosis
    ((data vector) (weights vector) &optional mean standard-deviation)
  (("gsl_stats" :type "_wkurtosis")
   ("gsl_stats" :type "_wkurtosis_m_sd"))
  ((((c-pointer weights) :pointer) (1 :int)
    ((c-pointer data) :pointer) (1 :int)
    ((dim0 data) sizet))
   (((c-pointer weights) :pointer) (1 :int)
    ((c-pointer data) :pointer) (1 :int)
    ((dim0 data) sizet)
    (mean :double) (standard-deviation :double)))
  :definition :generic
  :element-types :float
79
  :c-return :double
80
  :inputs (data weights)
Liam Healy's avatar
Liam Healy committed
81
  :documentation			; FDL
82 83 84
  "The weighted kurtosis of the dataset.
   kurtosis = ((\sum w_i ((x_i - xbar)/sigma)^4) / (\sum w_i)) - 3.")

liam's avatar
liam committed
85 86
;;; Examples and unit test

Liam Healy's avatar
Liam Healy committed
87
(save-test higher-moments
88
  (let ((vec #m(-3.21d0 1.0d0 12.8d0)))
89 90
      (let* ((mean (mean vec))
	     (sd (standard-deviation vec mean)))
91 92 93 94
	(list
	 (skewness vec)
	 (skewness vec mean sd)
	 (kurtosis vec)
liam's avatar
liam committed
95
	 (kurtosis vec mean sd)))))