Commit 9418a9b8 authored by liam's avatar liam
Browse files

Random distributions Type 1 and Type 2 Gumbel, and Dirichlet.


git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3151 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent 65dbcf2e
......@@ -3,7 +3,7 @@
; description: Definition of GSLL system
; date:
; author: Liam Healy
; modified: Sun Oct 22 2006 - 22:12
; modified: Sun Oct 29 2006 - 22:17
;********************************************************
;;; $Id: $
......@@ -115,4 +115,7 @@
(:file "logistic" :depends-on (rng-types))
(:file "pareto" :depends-on (rng-types))
(:file "spherical-vector" :depends-on (rng-types))
(:file "weibull" :depends-on (rng-types))))))
(:file "weibull" :depends-on (rng-types))
(:file "gumbel1" :depends-on (rng-types))
(:file "gumbel2" :depends-on (rng-types))
(:file "dirichlet" :depends-on (rng-types))))))
;********************************************************
; file: dirichlet.lisp
; description: Dirichlet distribution
; date: Sun Oct 29 2006
; author: Liam M. Healy
; modified: Sun Oct 29 2006 - 22:00
;********************************************************
;;; $Id: $
(in-package :gsl)
(defun-gsl dirichlet (generator alpha theta)
"gsl_ran_dirichlet"
(((generator generator) :pointer)
((1- (dim0 alpha)) :size)
((gsl-array alpha) :pointer)
;; theta had better be at least as long as alpha, or they'll be trouble
((gsl-array theta) :pointer))
:c-return :void
:documentation
"An array of @var{K} random variates from a Dirichlet
distribution of order @var{K}-1. The distribution function is
p(\theta_1,\ldots,\theta_K) \, d\theta_1 \cdots d\theta_K =
{1 \over Z} \prod_{i=1}^{K} \theta_i^{\alpha_i - 1}
\; \delta(1 -\sum_{i=1}^K \theta_i) d\theta_1 \cdots d\theta_K
@math{theta_i >= 0} and @math{alpha_i >= 0}.
The delta function ensures that @math{\sum \theta_i = 1}.
The normalization factor @math{Z} is
Z = {\prod_{i=1}^K \Gamma(\alpha_i) \over \Gamma( \sum_{i=1}^K \alpha_i)}
The random variates are generated by sampling @var{K} values
from gamma distributions with parameters @math{a=alpha_i, b=1},
and renormalizing.
See A.M. Law, W.D. Kelton, @cite{Simulation Modeling and Analysis}
(1991).")
(defun-gsl dirichlet-pdf (alpha theta)
"gsl_ran_dirichlet_pdf"
(((1- (dim0 alpha)) :size)
((gsl-array alpha) :pointer)
;; theta had better be at least as long as alpha, or they'll be trouble
((gsl-array theta) :pointer))
:c-return :double
:documentation
"The probability density @math{p(\theta_1, ... , \theta_K)}
at @var{theta}[@var{K}] for a Dirichlet distribution with parameters
@var{alpha}[@var{K}], using the formula given for #'dirichlet.")
(defun-gsl dirichlet-log-pdf (alpha theta)
"gsl_ran_dirichlet_lnpdf"
(((1- (dim0 alpha)) :size)
((gsl-array alpha) :pointer)
;; theta had better be at least as long as alpha, or they'll be trouble
((gsl-array theta) :pointer))
:c-return :double
:documentation
"The logarithm of the probability density
@math{p(\theta_1, ... , \theta_K)}
for a Dirichlet distribution with parameters
@var{alpha}[@var{K}].")
;;; Examples and unit test
(lisp-unit:define-test dirichlet
(lisp-unit:assert-equal
'("0.485599280805d-04" "0.578734200082d+00"
"0.421217239990d+00" "0.200000000000d+00")
(lisp-unit:fp-sequence
(progn
(rng-set *rng-mt19937* 0)
(with-data (alpha vector-double 4)
(with-data (theta vector-double 4)
(setf (data alpha) #(1.0d0 2.0d0 3.0d0 4.0d0))
(dirichlet *rng-mt19937* alpha theta)
(data theta))))))
(lisp-unit:assert-first-fp-equal
"0.288000000000d+01"
(with-data (alpha vector-double 4)
(with-data (theta vector-double 4)
(setf (data alpha) #(1.0 2.0 3.0 4.0))
(setf (data theta) #(0.1 0.3 0.4 0.2))
(dirichlet-pdf alpha theta))))
(lisp-unit:assert-first-fp-equal
"0.105779029415d+01"
(with-data (alpha vector-double 4)
(with-data (theta vector-double 4)
(setf (data alpha) #(1.0 2.0 3.0 4.0))
(setf (data theta) #(0.1 0.3 0.4 0.2))
(dirichlet-log-pdf alpha theta)))))
;********************************************************
; file: gumbel1.lisp
; description: Beta distribution
; date: Sun Oct 29 2006
; author: Liam M. Healy
; modified: Sun Oct 29 2006 - 21:07
;********************************************************
;;; $Id: $
(in-package :gsl)
(defun-gsl gumbel1 (generator a b)
"gsl_ran_gumbel1"
(((generator generator) :pointer) (a :double) (b :double))
:c-return :double
:documentation
"A random variate from the Type-1 Gumbel
distribution,
p(x) dx = a b \exp(-(b \exp(-ax) + ax)) dx
for @math{-\infty < x < \infty}.")
(defun-gsl gumbel1-pdf (x a b)
"gsl_ran_gumbel1_pdf" ((x :double) (a :double) (b :double))
:c-return :double
:documentation
"The probability density @math{p(x)} at @var{x}
for a Type-1 Gumbel distribution with parameters @var{a} and @var{b},
using the formula given for #'gumbel1.")
(defun-gsl gumbel1-P (x a b)
"gsl_cdf_gumbel1_P" ((x :double) (a :double) (b :double))
:c-return :double
:documentation "The cumulative distribution functions
@math{P(x)} for the Type-1 Gumbel distribution with
parameters @var{a} and @var{b}.")
(defun-gsl gumbel1-Q (x a b)
"gsl_cdf_gumbel1_Q" ((x :double) (a :double) (b :double))
:c-return :double
:documentation "The cumulative distribution functions
@math{Q(x)} for the Type-1 Gumbel distribution with
parameters @var{a} and @var{b}.")
(defun-gsl gumbel1-Pinv (P a b)
"gsl_cdf_gumbel1_Pinv" ((P :double) (a :double) (b :double))
:c-return :double
:documentation "The inverse cumulative distribution functions
@math{P(x)} for the Type-1 Gumbel distribution with
parameters @var{a} and @var{b}.")
(defun-gsl gumbel1-Qinv (Q a b)
"gsl_cdf_gumbel1_Qinv" ((Q :double) (a :double) (b :double))
:c-return :double
:documentation "The inverse cumulative distribution functions
@math{Q(x)} for the Type-1 Gumbel distribution with
parameters @var{a} and @var{b}.")
;;; Examples and unit test
(lisp-unit:define-test gumbel1
(lisp-unit:assert-equal
'("0.895459625749d+01" "0.973051899751d-01" "0.459135062331d+00"
"0.360741242243d+01" "0.313000274682d+00" "0.101657969497d+01"
"0.382920819366d+01" "0.191289739318d+01" "0.117748457895d+01"
"0.189323210797d+01" "0.198591186168d+01")
(lisp-unit:fp-sequence
(progn
(rng-set *rng-mt19937* 0)
(loop for i from 0 to 10
collect
(gumbel1 *rng-mt19937* 1.0d0 2.0d0)))))
(lisp-unit:assert-first-fp-equal
"0.296257089650d+00"
(gumbel1-pdf 0.1 1.0 2.0))
(lisp-unit:assert-first-fp-equal
"0.163707359877d+00"
(gumbel1-P 0.1 1.0 2.0))
(lisp-unit:assert-first-fp-equal
"0.836292640123d+00"
(gumbel1-Q 0.1 1.0 2.0))
(lisp-unit:assert-first-fp-equal
"0.100000000000d+00"
(gumbel1-Pinv 0.1637073598773166 1.0 2.0))
(lisp-unit:assert-first-fp-equal
"0.100000000000d+00"
(gumbel1-Qinv 0.8362926401226833 1.0 2.0)))
;********************************************************
; file: gumbel2.lisp
; description: Beta distribution
; date: Sun Oct 29 2006
; author: Liam M. Healy
; modified: Sun Oct 29 2006 - 21:07
;********************************************************
;;; $Id: $
(in-package :gsl)
(defun-gsl gumbel2 (generator a b)
"gsl_ran_gumbel2"
(((generator generator) :pointer) (a :double) (b :double))
:c-return :double
:documentation
"A random variate from the Type-2 Gumbel
distribution,
p(x) dx = a b x^@{-a-1@} \exp(-b x^@{-a@}) dx
for @math{0 < x < \infty}.")
(defun-gsl gumbel2-pdf (x a b)
"gsl_ran_gumbel2_pdf" ((x :double) (a :double) (b :double))
:c-return :double
:documentation
"The probability density @math{p(x)} at @var{x}
for a Type-2 Gumbel distribution with parameters @var{a} and @var{b},
using the formula given in #'gumbel2.")
(defun-gsl gumbel2-P (x a b)
"gsl_cdf_gumbel2_P" ((x :double) (a :double) (b :double))
:c-return :double
:documentation "The cumulative distribution functions
@math{P(x)} for the Type-2 Gumbel distribution with
parameters @var{a} and @var{b}.")
(defun-gsl gumbel2-Q (x a b)
"gsl_cdf_gumbel2_Q" ((x :double) (a :double) (b :double))
:c-return :double
:documentation "The cumulative distribution functions
@math{Q(x)} for the Type-2 Gumbel distribution with
parameters @var{a} and @var{b}.")
(defun-gsl gumbel2-Pinv (P a b)
"gsl_cdf_gumbel2_Pinv" ((P :double) (a :double) (b :double))
:c-return :double
:documentation "The inverse cumulative distribution functions
@math{P(x)} for the Type-2 Gumbel distribution with
parameters @var{a} and @var{b}.")
(defun-gsl gumbel2-Qinv (Q a b)
"gsl_cdf_gumbel2_Qinv" ((Q :double) (a :double) (b :double))
:c-return :double
:documentation "The inverse cumulative distribution functions
@math{Q(x)} for the Type-2 Gumbel distribution with
parameters @var{a} and @var{b}.")
;;; Examples and unit test
(lisp-unit:define-test gumbel2
(lisp-unit:assert-equal
'("0.774340085852d+04" "0.110219670168d+01" "0.158270445210d+01"
"0.368705239332d+02" "0.136752190666d+01" "0.276372579456d+01"
"0.460260800603d+02" "0.677268352507d+01" "0.324619836866d+01"
"0.664079780729d+01" "0.728568789702d+01")
(lisp-unit:fp-sequence
(progn
(rng-set *rng-mt19937* 0)
(loop for i from 0 to 10
collect
(gumbel2 *rng-mt19937* 1.0d0 2.0d0)))))
(lisp-unit:assert-first-fp-equal
"0.536256036829d-01"
(gumbel2-pdf 5.0 1.0 2.0))
(lisp-unit:assert-first-fp-equal
"0.818730753078d+00"
(gumbel2-P 10.0 1.0 2.0))
(lisp-unit:assert-first-fp-equal
"0.181269246922d+00"
(gumbel2-Q 10.0 1.0 2.0))
(lisp-unit:assert-first-fp-equal
"1.000000000000d+01"
(gumbel2-Pinv 0.8187307530779818 1.0 2.0))
(lisp-unit:assert-first-fp-equal
"0.100000000000d+02"
(gumbel2-Qinv 0.18126924692201815 1.0 2.0)))
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