cl-randist.html @@ -15,6 +15,7 @@ ;; It needs only one call to RANDOM for every value produced. ;; ;; For the licence, see at the end of the file. + (defun create-alias-method-vectors (probabilities) (let* ((N (length probabilities)) (threshold (/ 1.0d0 N)) @@ -60,6 +61,12 @@ (values p_keep alternatives))) (defun make-discrete-random-var (probabilities &optional values) +" The function MAKE-DISCRETE-RANDOM-VAR takes an array of +probabilities and an (optional) array of values. Produces a +function which returns each of the values with the specified +probability (or the corresponding integer no values have been +given)." + (when (and values (/= (length values) (length probabilities))) (error "different number of values and probabilities.")) alias_method.lisp @@ -17,6 +17,11 @@ (declaim (ftype (function (double-float double-float) double-float) random-beta)) (defun random-beta (a b) +"The beta distribution has the form + + p(x) dx = (Gamma(a + b)/(Gamma(a) Gamma(b))) x^(a-1) (1-x)^(b-1) dx + +The method used here is the one described in Knuth " (declare (double-float a b)) (let ((x1 (random-gamma a 1d0)) (x2 (random-gamma b 1d0))) beta.lisp @@ -47,8 +47,12 @@ ;; } (declaim (ftype (function (double-float integer) integer) random-binomial)) - (defun random-binomial (p n) +"The binomial distribution has the form, + + prob(k) = n!/(k!(n-k)!) * p^k (1-p)^(n-k) for k = 0, 1, ..., n + +This is the algorithm from Knuth" (let ((a 0) (b 0) (k 0) (X 0d0) (p p) (n n)) binomial.lisp @@ -205,6 +205,23 @@ (inline random-gamma1)) (defun random-gamma1 (a b) +"The Gamma distribution of order a>0 is defined by: + + p(x) dx = {1 / \Gamma(a) b^a } x^{a-1} e^{-x/b} dx + + for x>0. If X and Y are independent gamma-distributed random + variables of order a1 and a2 with the same scale parameter b, then + X+Y has gamma distribution of order a1+a2. + + The algorithms below are from Knuth, vol 2, 2nd ed, p. 129. + + + Works only if a > 1, and is most efficient if a is large + + This algorithm, reported in Knuth, is attributed to Ahrens. A + faster one, we are told, can be found in: J. H. Ahrens and + U. Dieter, Computing 12 (1974) 223-246." + (declare (double-float a b)) (assert (> a 0d0)) (multiple-value-bind (na frac) (truncate a) @@ -267,6 +284,10 @@ (declaim (ftype (function (double-float double-float) double-float) random-gamma-mt)) (defun random-gamma-mt (a b) +"New version based on Marsaglia and Tsang, 'A Simple Method for +generating gamma variables', ACM Transactions on Mathematical +Software, Vol 26, No 3 (2000), p363-372." + (declare (double-float a b)) (if (< a 1d0) (* (random-gamma-mt (+ 1d0 a) b) (expt (random-uniform) (/ a))) @@ -299,4 +320,5 @@ (declaim (inline random-gamma)) (defun random-gamma (a &optional (b 1d0)) + "[syntax suggar] Generate a random variable with gamma distribution using the MT method (see random-gamma-mt)" (random-gamma-mt a b)) gamma.lisp @@ -59,6 +59,7 @@ gsl_ran_multinomial (const gsl_rng * r, const size_t K, |# (defun random-multinomial1 (NN p n) + "Return the genrated values in the n vector" (let ((norm 0d0) (k (1- (array-dimension p 0)))) @@ -81,6 +82,22 @@ gsl_ran_multinomial (const gsl_rng * r, const size_t K, (defun random-multinomial (NN p) +" The multinomial distribution has the form + + N! n_1 n_2 n_K + prob(n_1, n_2, ... n_K) = -------------------- p_1 p_2 ... p_K + (n_1! n_2! ... n_K!) + + where n_1, n_2, ... n_K are nonnegative integers, sum_{k=1,K} n_k = N, + and p = (p_1, p_2, ..., p_K) is a probability distribution. + + Random variates are generated using the conditional binomial method. + This scales well with N and does not require a setup step. + + Ref: + C.S. David, The computer generation of multinomial random variates, + Comp. Stat. Data Anal. 16 (1993) 205-217" + (let ((n (make-array (array-dimension p 0) :element-type 'integer :adjustable nil))) multinomial.lisp @@ -20,6 +20,11 @@ (inline random-negative-binomial)) (defun random-negative-binomial (p n) +" The negative binomial distribution has the form, + prob(k) = Gamma(n + k)/(Gamma(n) Gamma(k + 1)) p^n (1-p)^k + for k = 0, 1, ... . Note that n does not have to be an integer. + This is the Leger's algorithm (given in the answers in Knuth)" + (declare (type double-float p) (integer n)) (let ((X (random-gamma (coerce n 'double-float) 1d0))) nbinomial.lisp @@ -193,6 +193,29 @@ (declaim (ftype (function (double-float double-float) double-float) random-normal-ziggurat)) (defun random-normal-ziggurat (mean sigma) +" This routine is based on the following article, with a couple of + modifications which simplify the implementation. + + George Marsaglia, Wai Wan Tsang + The Ziggurat Method for Generating Random Variables + Journal of Statistical Software, vol. 5 (2000), no. 8 + http://www.jstatsoft.org/v05/i08/ + + The modifications are: + + 1) use 128 steps instead of 256 to decrease the amount of static + data necessary. + + 2) use an acceptance sampling from an exponential wedge + exp(-R*(x-R/2)) for the tail of the base strip to simplify the + implementation. The area of exponential wedge is used in + calculating 'v' and the coefficients in ziggurat table, so the + coefficients differ slightly from those in the Marsaglia and Tsang + paper. + + See also Leong et al, 'A Comment on the Implementation of the + Ziggurat Method', Journal of Statistical Software, vol 5 (2005), no 7." + (let ((i 0) (j 0) (sign 0) (x 0d0) (y 0d0)) (declare (double-float mean sigma)) (declare (double-float x y)) @@ -241,5 +264,6 @@ (declaim (inline random-normal)) (defun random-normal (&optional (mean 0d0) (sigma 1d0)) + "[sintax suggar] Generate random variable with normal distribution using ziggurat method" (random-normal-ziggurat mean sigma)) \ No newline at end of file normal.lisp @@ -51,6 +51,12 @@ (declaim (ftype (function (double-float) integer) random-poisson)) (defun random-poisson (mu) +"The poisson distribution has the form + + p(n) = (mu^n / n!) exp(-mu) + + for n = 0, 1, 2, ... . The method used here is the one from Knuth." + (declare (type double-float mu)) (let ((k 0)) (declare (type integer k)) poisson.lisp @@ -8,6 +8,7 @@ (random x state))) (defmacro random-uniform () + "[syntax suggar] Random variable with uniform distribution in interval [0,1]" `(random-mt 1d0)) (declaim (ftype (function () double-float) random-pos)) randist.lisp bc4ac5c99c62e0d3cd792f55f5ff0f1dbce93e57 Leonardo Varuzza varuzza@gmail.com http://github.com/lvaruzza/cl-randist/commit/bae34807469e77168f1889c962efdfb0f0faca37 bae34807469e77168f1889c962efdfb0f0faca37 2008-03-07T19:49:34-08:00 2008-03-07T19:49:34-08:00 Add documentation. 311d5e00eb09b4670d3cfbe22ab49574515042b0 Leonardo Varuzza varuzza@gmail.com