Commit bdffb15c authored by liam's avatar liam
Browse files

Add random distributions: poisson, bernoulli, binomial, multnomial,

negative-binomial (includes Pascal), geometric, hypergeometric,
and logarithmic, and fix syntax of discrete.


git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3157 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent 8abe8ab1
......@@ -3,7 +3,7 @@
; description: Definition of GSLL system
; date:
; author: Liam Healy
; modified: Sun Nov 12 2006 - 00:10
; modified: Sat Nov 25 2006 - 23:01
;********************************************************
;;; $Id: $
......@@ -119,4 +119,12 @@
(:file "gumbel1" :depends-on (rng-types))
(:file "gumbel2" :depends-on (rng-types))
(:file "dirichlet" :depends-on (rng-types))
(:file "discrete" :depends-on (rng-types))))))
(:file "discrete" :depends-on (rng-types))
(:file "poisson" :depends-on (rng-types))
(:file "bernoulli" :depends-on (rng-types))
(:file "binomial" :depends-on (rng-types))
(:file "multnomial" :depends-on (rng-types))
(:file "negative-binomial" :depends-on (rng-types))
(:file "geometric" :depends-on (rng-types))
(:file "hypergeometric" :depends-on (rng-types))
(:file "logarithmic" :depends-on (rng-types))))))
;********************************************************
; file: bernoulli.lisp
; description: Bernoulli distribution
; date: Sat Nov 25 2006 - 16:59
; author: Liam M. Healy
; modified: Sat Nov 25 2006 - 17:01
;********************************************************
;;; $Id: $
(in-package :gsl)
(defun-gsl bernoulli (generator p)
"gsl_ran_bernoulli"
(((generator generator) :pointer) (p :double))
:c-return :uint
:documentation
"Returns either 0 or 1, the result of a Bernoulli trial
with probability @var{p}. The probability distribution for
a Bernoulli trial is
p(0) = 1 - p
p(1) = p.")
(defun-gsl bernoulli-pdf (k p)
"gsl_ran_bernoulli_pdf" ((k :uint) (p :double))
:c-return :double
:documentation "The probability @math{p(k)} of obtaining
@var{k} from a Bernoulli distribution with probability parameter
@var{p}, using the formula given in #'bernoulli.")
;;; Examples and unit test
(lisp-unit:define-test bernoulli
(lisp-unit:assert-equal
'(0 1 1 0 1 1 0 0 0 0 0)
(progn
(rng-set *rng-mt19937* 0)
(loop for i from 0 to 10
collect
(bernoulli *rng-mt19937* 0.5d0))))
(lisp-unit:assert-first-fp-equal
"0.500000000000d+00"
(bernoulli-pdf 0 0.5d0)))
;********************************************************
; file: binomial.lisp
; description: Binomial distribution
; date: Sat Nov 25 2006 - 16:00
; author: Liam M. Healy
; modified: Sat Nov 25 2006 - 18:08
;********************************************************
;;; $Id: $
(in-package :gsl)
(defun-gsl binomial (generator p n)
"gsl_ran_binomial"
(((generator generator) :pointer) (p :double) (n :uint))
:c-return :uint
:documentation
"A random integer from the binomial distribution,
the number of successes in @var{n} independent trials with probability
@var{p}. The probability distribution for binomial variates is,
p(k) = {n! \over k! (n-k)!} p^k (1-p)^{n-k}
@math{0 <= k <= n}.")
(defun-gsl binomial-pdf (k p n)
"gsl_ran_binomial_pdf" ((k :uint) (p :double) (n :uint))
:c-return :double
:documentation
"The probability @math{p(k)} of obtaining @var{k}
from a binomial distribution with parameters @var{p} and @var{n}, using
the formula given in #'binomial.")
(defun-gsl binomial-P (k p n)
"gsl_cdf_binomial_P" ((k :uint) (p :double) (n :uint))
:c-return :double
:documentation "The cumulative distribution functions
@math{P(k)} for the Binomial distribution
with parameters @var{p} and @var{n}.")
(defun-gsl binomial-Q (k p n)
"gsl_cdf_binomial_Q" ((k :uint) (p :double) (n :uint))
:c-return :double
:documentation "The cumulative distribution functions
@math{Q(k)} for the Binomial distribution
with parameters @var{p} and @var{n}.")
;;; Examples and unit test
(lisp-unit:define-test binomial
(lisp-unit:assert-equal
'(11 3 4 8 4 5 8 6 5 6 6)
(progn
(rng-set *rng-mt19937* 0)
(loop for i from 0 to 10
collect
(binomial *rng-mt19937* 0.4d0 12))))
(lisp-unit:assert-first-fp-equal
"0.227030335488d+00"
(binomial-pdf 5 0.4d0 12))
(lisp-unit:assert-first-fp-equal
"0.665208557568d+00"
(binomial-P 5 0.4d0 12))
(lisp-unit:assert-first-fp-equal
"0.334791442432d+00"
(binomial-Q 5 0.4d0 12)))
......@@ -3,7 +3,7 @@
; description: Discrete random variables
; date: Sat Nov 11 2006 - 21:51
; author: Liam M. Healy
; modified: Sun Nov 12 2006 - 00:09
; modified: Sat Nov 25 2006 - 16:46
;********************************************************
;;; $Id: $
......@@ -66,7 +66,7 @@
'(1 0 1 1 0 1 1 2 1 2 2)
(with-data (probabilities vector-double 3)
(setf (data probabilities) #(0.25d0 0.5d0 0.25d0))
(with-discrete-table probabilities table
(with-discrete-table (probabilities table)
(rng-set *rng-mt19937* 0)
(loop for i from 0 to 10
collect
......@@ -76,8 +76,7 @@
(with-data (probabilities vector-double 3)
(setf (data probabilities)
#(0.25d0 0.5d0 0.25d0))
(with-discrete-table probabilities
table
(with-discrete-table (probabilities table)
(discrete-pdf
1
table))))))
;********************************************************
; file: geometric.lisp
; description: Geometric distribution
; date: Sat Nov 25 2006 - 16:00
; author: Liam M. Healy
; modified: Sat Nov 25 2006 - 21:56
;********************************************************
;;; $Id: $
(in-package :gsl)
(defun-gsl geometric (generator p)
"gsl_ran_geometric"
(((generator generator) :pointer) (p :double))
:c-return :uint
:documentation
"A random integer from the geometric distribution,
the number of independent trials with probability @var{p} until the
first success. The probability distribution for geometric variates
is p(k) = p (1-p)^{k-1} for @math{k >= 1}.
Note that the distribution begins with @math{k=1} with this
definition. There is another convention in which the exponent @math{k-1}
is replaced by @math{k}.")
(defun-gsl geometric-pdf (k p)
"gsl_ran_geometric_pdf" ((k :uint) (p :double))
:c-return :double
:documentation
"The probability @math{p(k)} of obtaining @var{k}
from a geometric distribution with probability parameter @var{p}, using
the formula given in #'geometric.")
(defun-gsl geometric-P (k p)
"gsl_cdf_geometric_P" ((k :uint) (p :double))
:c-return :double
:documentation "The cumulative distribution functions
@math{P(k)} for the Geometric distribution with parameter @var{p}.")
(defun-gsl geometric-Q (k p)
"gsl_cdf_geometric_Q" ((k :uint) (p :double))
:c-return :double
:documentation "The cumulative distribution functions
@math{Q(k)} for the Geometric distribution with parameters @var{p}.")
;;; Examples and unit test
(lisp-unit:define-test geometric
(lisp-unit:assert-equal
'(1 4 3 1 3 2 1 1 2 1 1)
(progn
(rng-set *rng-mt19937* 0)
(loop for i from 0 to 10
collect
(geometric *rng-mt19937* 0.4d0))))
(lisp-unit:assert-first-fp-equal
"0.240000000000d+00"
(geometric-pdf 2 0.4d0))
(lisp-unit:assert-first-fp-equal
"0.640000000000d+00"
(geometric-P 2 0.4d0))
(lisp-unit:assert-first-fp-equal
"0.360000000000d+00"
(geometric-Q 2 0.4d0)))
;********************************************************
; file: hypergeometric.lisp
; description: Hypergeometric distribution
; date: Sat Nov 25 2006 - 16:00
; author: Liam M. Healy
; modified: Sat Nov 25 2006 - 22:40
;********************************************************
;;; $Id: $
(in-package :gsl)
(defun-gsl hypergeometric (generator n1 n2 tt)
"gsl_ran_hypergeometric"
(((generator generator) :pointer) (n1 :uint) (n2 :uint)(tt :uint))
:c-return :uint
:documentation
"A random integer from the hypergeometric
distribution. The probability distribution for hypergeometric
random variates is
p(k) = C(n_1, k) C(n_2, t - k) / C(n_1 + n_2, t)
where @math{C(a,b) = a!/(b!(a-b)!)} and
@math{t <= n_1 + n_2}. The domain of @math{k} is
@math{max(0,t-n_2), ..., min(t,n_1)}.
If a population contains @math{n_1} elements of ``type 1'' and
@math{n_2} elements of ``type 2'' then the hypergeometric
distribution gives the probability of obtaining @math{k} elements of
``type 1'' in @math{t} samples from the population without
replacement.")
(defun-gsl hypergeometric-pdf (k n1 n2 tt)
"gsl_ran_hypergeometric_pdf" ((k :uint) (n1 :uint) (n2 :uint)(tt :uint))
:c-return :double
:documentation
"The probability @math{p(k)} of obtaining @var{k}
from a hypergeometric distribution with parameters @var{n1}, @var{n2},
@var{tt}, using the formula given in #'hypergeometric.")
(defun-gsl hypergeometric-P (k n1 n2 tt)
"gsl_cdf_hypergeometric_P" ((k :uint) (n1 :uint) (n2 :uint)(tt :uint))
:c-return :double
:documentation
"The cumulative distribution functions @math{P(k)} for the
hypergeometric distribution with parameters @var{n1}, @var{n2}
and @var{tt}.")
(defun-gsl hypergeometric-Q (k n1 n2 tt)
"gsl_cdf_hypergeometric_Q" ((k :uint) (n1 :uint) (n2 :uint)(tt :uint))
:c-return :double
:documentation
"The cumulative distribution functions @math{Q(k)} for the
hypergeometric distribution with parameters @var{n1}, @var{n2}
and @var{tt}.")
;;; Examples and unit test
(lisp-unit:define-test hypergeometric-randist
(lisp-unit:assert-equal
'(2 1 0 0 1 1 3 1 0 1 3)
(progn
(rng-set *rng-mt19937* 0)
(loop for i from 0 to 10
collect
(hypergeometric *rng-mt19937* 3 6 3))))
(lisp-unit:assert-first-fp-equal
"0.357142857143d+00"
(hypergeometric-pdf 0 2 6 3))
(lisp-unit:assert-first-fp-equal
"0.892857142857d+00"
(hypergeometric-P 1 2 6 3))
(lisp-unit:assert-first-fp-equal
"0.107142857143d+00"
(hypergeometric-Q 1 2 6 3)))
;********************************************************
; file: logarithmic.lisp
; description: Logarithmic distribution
; date: Sat Nov 25 2006 - 16:00
; author: Liam M. Healy
; modified: Sat Nov 25 2006 - 22:57
;********************************************************
;;; $Id: $
(in-package :gsl)
(defun-gsl logarithmic (generator p)
"gsl_ran_logarithmic"
(((generator generator) :pointer) (p :double))
:c-return :uint
:documentation
"A random integer from the logarithmic distribution.
The probability distribution for logarithmic random variates
is p(k) = {-1 \over \log(1-p)} {\left( p^k \over k \right)}
for @math{k >= 1}.")
(defun-gsl logarithmic-pdf (k p)
"gsl_ran_logarithmic_pdf" ((k :uint) (p :double))
:c-return :double
:documentation
"The probability @math{p(k)} of obtaining @var{k}
from a logarithmic distribution with probability parameter @var{p},
using the formula given in #'logarithmic.")
;;; Examples and unit test
(lisp-unit:define-test logarithmic
(lisp-unit:assert-equal
'(1 3 1 4 1 1 2 1 1 5 2)
(progn
(rng-set *rng-mt19937* 0)
(loop for i from 0 to 10
collect
(logarithmic *rng-mt19937* 0.9d0))))
(lisp-unit:assert-first-fp-equal
"0.156609215118d+00"
(logarithmic-pdf 2 0.4d0)))
;********************************************************
; file: multinomial.lisp
; description: Multinomial distribution
; date: Sat Nov 25 2006 - 16:00
; author: Liam M. Healy
; modified: Sat Nov 25 2006 - 21:16
;********************************************************
;;; $Id: $
(in-package :gsl)
(defun-gsl multinomial (generator sum p n)
"gsl_ran_multinomial"
(((generator generator) :pointer)
((dim0 p) :size)
(sum :size)
((gsl-array p) :pointer)
;; technically, n should be a uint array, but integers work
((gsl-array n) :pointer))
:c-return :void
:documentation
"Returns an array n of @var{K} random variates from a
multinomial distribution. The sum of the array n is specified
by sum=N. The distribution function is
P(n_1, n_2, ..., n_K) =
(N!/(n_1! n_2! ... n_K!)) p_1^n_1 p_2^n_2 ... p_K^n_K
where @math{(n_1, n_2, ..., n_K)}
are nonnegative integers with
@math{sum_@{k=1@}^K n_k = N},
and @math{(p_1, p_2, ..., p_K)}
is a probability distribution with @math{\sum p_i = 1}.
If the array @var{p}[@var{K}] is not normalized then its entries will be
treated as weights and normalized appropriately.
Random variates are generated using the conditional binomial method (see
C.S. David, @cite{The computer generation of multinomial random
variates}, Comp. Stat. Data Anal. 16 (1993) 205--217 for details).")
(defun-gsl multinomial-pdf (p n)
"gsl_ran_multinomial_pdf"
(((dim0 p) :uint) ((gsl-array p) :pointer) ((gsl-array n) :pointer))
:c-return :double
:documentation
"Compute the probability @math{P(n_1, n_2, ..., n_K)}
of sampling @var{n}[@var{K}] from a multinomial distribution
with parameters @var{p}[@var{K}], using the formula given for #'multinomial.")
(defun-gsl multinomial-log-pdf (p n)
"gsl_ran_multinomial_lnpdf"
(((dim0 p) :uint) ((gsl-array p) :pointer) ((gsl-array n) :pointer))
:c-return :double
:documentation
"Compute the natural logarithm of the probability @math{P(n_1, n_2, ..., n_K)}
of sampling @var{n}[@var{K}] from a multinomial distribution
with parameters @var{p}[@var{K}], using the formula given for #'multinomial.")
;;; Examples and unit test
(lisp-unit:define-test multinomial
(lisp-unit:assert-equalp
#(5 0 1 2)
(progn
(rng-set *rng-mt19937* 0)
(with-data (p vector-double 4)
(setf (data p) #(0.1d0 0.2d0 0.3d0 0.4d0))
(with-data (n vector-fixnum 4)
(multinomial *rng-mt19937* 8 p n)
(data n)))))
(lisp-unit:assert-first-fp-equal
"0.806400000000d-04"
(with-data (p vector-double 4)
(setf (data p) #(0.1d0 0.2d0 0.3d0 0.4d0))
(with-data (n vector-fixnum 4)
(setf (data n) #(5 0 1 2))
(multinomial-pdf p N))))
(lisp-unit:assert-first-fp-equal
"-0.942551575364d+01"
(with-data (p vector-double 4)
(setf (data p) #(0.1d0 0.2d0 0.3d0 0.4d0))
(with-data (n vector-fixnum 4)
(setf (data n) #(5 0 1 2))
(multinomial-log-pdf p n)))))
;********************************************************
; file: negative-binomial.lisp
; description: Negative binomial and Pascal distributions
; date: Sat Nov 25 2006 - 16:00
; author: Liam M. Healy
; modified: Sat Nov 25 2006 - 21:42
;********************************************************
;;; $Id: $
(in-package :gsl)
;;;;****************************************************************************
;;;; Negative binomial
;;;;****************************************************************************
(defun-gsl negative-binomial (generator p n)
"gsl_ran_negative_binomial"
(((generator generator) :pointer) (p :double) (n :double))
:c-return :uint
:documentation
"A random integer from the negative binomial
distribution, the number of failures occurring before @var{n} successes
in independent trials with probability @var{p} of success. The
probability distribution for negative binomial variates is,
p(k) = {\Gamma(n + k) \over \Gamma(k+1) \Gamma(n) } p^n (1-p)^k
Note that @math{n} is not required to be an integer.")
(defun-gsl negative-binomial-pdf (k p n)
"gsl_ran_negative_binomial_pdf" ((k :uint) (p :double) (n :double))
:c-return :double
:documentation
"The probability @math{p(k)} of obtaining @var{k}
from a negative binomial distribution with parameters @var{p} and
@var{n}, using the formula given in #'negative-binomial.")
(defun-gsl negative-binomial-P (k p n)
"gsl_cdf_negative_binomial_P" ((k :uint) (p :double) (n :double))
:c-return :double
:documentation "The cumulative distribution functions
@math{P(k)} for the negative binomial distribution
with parameters @var{p} and @var{n}.")
(defun-gsl negative-binomial-Q (k p n)
"gsl_cdf_negative_binomial_Q" ((k :uint) (p :double) (n :double))
:c-return :double
:documentation "The cumulative distribution functions
@math{Q(k)} for the negative binomial distribution
with parameters @var{p} and @var{n}.")
;;;;****************************************************************************
;;;; Pascal
;;;;****************************************************************************
(defun-gsl pascal (generator p n)
"gsl_ran_pascal"
(((generator generator) :pointer) (p :double) (n :uint))
:c-return :uint
:documentation
"A random integer from the Pascal distribution. The
Pascal distribution is simply a negative binomial distribution with an
integer value of @math{n}.
p(k) = {(n + k - 1)! \over k! (n - 1)! } p^n (1-p)^k
@math{k >= 0}.")
(defun-gsl pascal-pdf (k p n)
"gsl_ran_pascal_pdf" ((k :uint) (p :double) (n :uint))
:c-return :double
:documentation
"The probability @math{p(k)} of obtaining @var{k}
from a Pascal distribution with parameters @var{p} and
@var{n}, using the formula given in #'pascal.")
(defun-gsl pascal-P (k p n)
"gsl_cdf_pascal_P" ((k :uint) (p :double) (n :uint))
:c-return :double
:documentation "The cumulative distribution functions
@math{P(k)} for the Pascal distribution
with parameters @var{p} and @var{n}.")
(defun-gsl pascal-Q (k p n)
"gsl_cdf_pascal_Q" ((k :uint) (p :double) (n :uint))
:c-return :double
:documentation "The cumulative distribution functions
@math{Q(k)} for the Pascal distribution
with parameters @var{p} and @var{n}.")
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************
(lisp-unit:define-test negative-binomial
(lisp-unit:assert-equal
'(10 7 12 23 20 24 18 12 4 22 15)
(progn
(rng-set *rng-mt19937* 0)
(loop for i from 0 to 10
collect
(negative-binomial *rng-mt19937* 0.4d0 12.0d0))))
(lisp-unit:assert-first-fp-equal
"0.569847670899d-02"
(negative-binomial-pdf 5 0.4d0 12.0d0))
(lisp-unit:assert-first-fp-equal
"0.105942025555d-01"
(negative-binomial-P 5 0.4d0 12.0d0))
(lisp-unit:assert-first-fp-equal
"0.989405797444d+00"
(negative-binomial-Q 5 0.4d0 12.0d0))
(lisp-unit:assert-equal
'(10 7 12 23 20 24 18 12 4 22 15)
(progn
(rng-set *rng-mt19937* 0)
(loop for i from 0 to 10
collect
(pascal *rng-mt19937* 0.4d0 12))))
(lisp-unit:assert-first-fp-equal
"0.569847670899d-02"
(pascal-pdf 5 0.4d0 12))
(lisp-unit:assert-first-fp-equal
"0.105942025555d-01"
(pascal-P 5 0.4d0 12))
(lisp-unit:assert-first-fp-equal
"0.989405797444d+00"
(pascal-Q 5 0.4d0 12)))
;********************************************************
; file: poisson.lisp
; description: Poisson distribution
; date: Sat Nov 25 2006 - 16:00
; author: Liam M. Healy
; modified: Sat Nov 25 2006 - 17:00
;********************************************************
;;; $Id: $
(in-package :gsl)
(defun-gsl poisson (generator mu)
"gsl_ran_poisson"
(((generator generator) :pointer) (mu :double))
:c-return :uint
:documentation
"A random integer from the Poisson distribution with mean @var{mu}.
The probability distribution for Poisson variates is
p(k) = {\mu^k \over k!} \exp(-\mu)
@math{k >= 0}.")
(defun-gsl poisson-pdf (k mu)
"gsl_ran_poisson_pdf" ((k :uint) (mu :double))
:c-return :double
:documentation
"The probability @math{p(k)} of obtaining @var{k}
from a Poisson distribution with mean @var{mu} using the formula
given in #'poisson.")
(defun-gsl poisson-P (k mu)
"gsl_cdf_poisson_P" ((k :uint) (mu :double))
:c-return :double
:documentation "The cumulative distribution functions
@math{P(k)} for the Poisson distribution with parameter @var{mu}.")
(defun-gsl poisson-Q (k mu)
"gsl_cdf_poisson_Q" ((k :uint) (mu :double))
:c-return :double
:documentation "The cumulative distribution functions
@math{Q(k)} for the Poisson distribution with parameter @var{mu}.")
;;; Examples and unit test
(lisp-unit:define-test poisson
(lisp-unit:assert-equal
'(15 6 9 9 5 8 11 9 11 5 10)
(progn
(rng-set *rng-mt19937* 0)
(loop for i from 0 to 10
collect
(poisson *rng-mt19937* 10.0d0))))
(lisp-unit:assert-first-fp-equal
"0.112599032149d+00"
(poisson-pdf 8 10.0d0))
(lisp-unit:assert-first-fp-equal
"0.332819678751d+00"
(poisson-P 8 10.0d0))
(lisp-unit:assert-first-fp-equal
"0.667180321249d+00"
(poisson-Q 8 10.0d0)))
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment