Common Lisp Package: RANDOM-DISTRIBUTIONS

README:

FUNCTION

Public

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).

MAKE-MT-RANDOM-STATE (&OPTIONAL STATE)

Analogous to Common Lisp's MAKE-RANDOM-STATE except that this function works on random states for JMT's Mersenne Twister implementation. AJR-comment: this version is better for statistical computations, since it can be used by an externally defined seed. However, it must be optimized, being the core of any pseudo-stochastic algorithm.

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

RANDOM-BINOMIAL

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

RANDOM-CHI-SQUARE

Generate random variable for chi square distribution: p(x) dx = (1/(2*Gamma(nu/2))) (x/2)^(nu/2 - 1) exp(-x/2) dx

RANDOM-EXPONENTIAL

Random values for: p(x) dx = exp(-x/mu) dx/mu

RANDOM-F

Random value for: p(x) dx = (nu1^(nu1/2) nu2^(nu2/2) Gamma((nu1 + nu2)/2) / Gamma(nu1/2) Gamma(nu2/2)) * x^(nu1/2 - 1) (nu2 + nu1 * x)^(-nu1/2 -nu2/2) dx

RANDOM-GAMMA (A &OPTIONAL (B 1.0d0))

[syntax suggar] Generate a random variable with gamma distribution using the MT method (see random-gamma-mt)

RANDOM-GAMMA-INT (A)

Random variable with gamma distribution with integer parameter.

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.

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.

RANDOM-GIG

Random Generalized Inverse Poisson The algorithm is based on that given by Dagpunar (1989)

RANDOM-GIG-IID

Random Generalized Inverse Poisson (vector version)

RANDOM-MULTINOMIAL%

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

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)

RANDOM-NORMAL (&OPTIONAL (MEAN 0.0d0) (SIGMA 1.0d0))

Generate random variable with normal distribution using ziggurat method

RANDOM-NORMAL-BIVARIATE (SIGMA-X SIGMA-Y &KEY (RHO 0.0d0) (ERROR-LIMIT 500))

Return a pair of numbers with specific correlation coefficent rho and with specified variances sigma-x and sigma-y; a direct port of gsl_ran_bivariate_gaussian from the GNU Scientific Library: void gsl_ran_bivariate_gaussian (const gsl_rng * r, double sigma_x, double sigma_y, double rho, double *x, double *y) { double u, v, r2, scale; do { /* choose x,y in uniform square (-1,-1) to (+1,+1) */ u = -1 + 2 * gsl_rng_uniform (r); v = -1 + 2 * gsl_rng_uniform (r); /* see if it is in the unit circle */ r2 = u * u + v * v; } while (r2 > 1.0 || r2 == 0); scale = sqrt (-2.0 * log (r2) / r2); *x = sigma_x * u * scale; *y = sigma_y * (rho * u + sqrt(1 - rho*rho) * v) * scale; } I need a good test for this one.

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.

RANDOM-PARETO

Random value for parato distribution: p(x) dx = (a/b) / (x/b)^(a+1) dx for x >= b

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.

Undocumented

MAKE-DISCRETE-MONOTONE-RANDOM-VAR (P)

MAKE-RANDOM-VARIABLE-GIG

MAKE-RANDOM-VARIABLE-GIG-POISSON

RANDOM-BERNOULLI (P)

RANDOM-MT (X &OPTIONAL (STATE *RANDOM-STATE*))

RANDOM-MULTINOMIAL

Private

MT-MAKE-RANDOM-STATE-INTEGER (N)

Use the single integer to expand into a bunch of integers to use as an MT-RANDOM-STATE. Copied from the 'sgenrand' function in mt19937int.c. This is mostly an internal function. I recommend using MAKE-MT-RANDOM-STATE unless specific circumstances dictate otherwise.

MT-MAKE-RANDOM-STATE-RANDOM

Generate a new random state from a new, hopefully somewhat random, value. This is mostly an internal function. I recommend using MAKE-MT-RANDOM-STATE unless specific circumstances dictate otherwise.

MT-RANDOM-STATE-ARR (INSTANCE)

@arg[extid]{A @class{extid}} @return[sytemid]{puri:uri or nil} Returns the System ID part of this External ID.

MT-RANDOM-STATE-MTI (INSTANCE)

@arg[extid]{A @class{extid}} @return[sytemid]{puri:uri or nil} Returns the System ID part of this External ID.

MT-REFILL

In the C program mt19937int.c, there is a function called 'genrand', & in that function there is a block of code to execute when the mt[] array is exhausted. This function is that block of code. I've removed it from my MT-GENRAND function for clarity.

RANDOM-MULTINOMIAL1

Return the genrated values in the n vector

RANDOM-POS

Create the sign, i.e. random positive or negative, similar to a binary result.

RANDOM-VECTOR-IID (N VARIABLE)

Return a vector with n IID instances of variable

ZEROIN

zero of the function f(x) is computed in the interval ax,bx . input.. ax left endpoint of initial interval bx right endpoint of initial interval f function subprogram which evaluates f(x) for any x in the interval ax,bx tol desired length of the interval of uncertainty of the final result ( .ge. 0.0d0) output.. zeroin abcissa approximating a zero of f in the interval ax,bx it is assumed that f(ax) and f(bx) have opposite signs without a check. zeroin returns a zero x in the given interval ax,bx to within a tolerance 4*macheps*abs(x) + tol, where macheps is the relative machine precision. this function subprogram is a slightly modified translation of the algol 60 procedure zero given in richard brent, algorithms for minimization without derivatives, prentice - hall, inc. (1973).

Undocumented

%RANDOM-GIG

COEFFICIENT-OF-VARIATION (SEQUENCE)

CONVERT-TO-DOUBLE-FLOAT-VECTOR

COPY-MT-RANDOM-STATE (INSTANCE)

CREATE-ALIAS-METHOD-VECTORS (PROBABILITIES)

GAMMA-FRAC (A)

GAMMA-LARGE (A)

GEN-UNIFORM (&OPTIONAL (N 100) (MAX 1.0))

GENLIST (F N)

GEOMETRIC-MEAN (SEQUENCE &OPTIONAL (BASE 10))

GIG-SETUP

MEAN (SEQUENCE)

MEDIAN (SEQUENCE)

MODE (SEQUENCE)

MT-GENRAND

MT-INTERNAL-MAKE-RANDOM-STATE (&KEY ((MTI DUM0) NIL) ((ARR DUM1) NIL))

SETFMT-RANDOM-STATE-ARR (NEW-VALUE INSTANCE)

SETFMT-RANDOM-STATE-MTI (NEW-VALUE INSTANCE)

MT-RANDOM-STATE-P (OBJECT)

MT-TEMPERING-SHIFT-L (N)

MT-TEMPERING-SHIFT-S (N)

MT-TEMPERING-SHIFT-T (N)

MT-TEMPERING-SHIFT-U (N)

PERCENTILE (SEQUENCE PERCENT)

PROFILE-GAMMA (&OPTIONAL (N (* 100 1000)))

RANDOM-DIRICHLET1 (ALPHA THETA)

RANDOM-GENERALIZED-INVERSE-POISSON

RANDOM-T (NU)

RANDOM-T-NU<=2 (NU)

RANDOM-T-NU>2 (NU)

RANGE (SEQUENCE)

REPORT-UNIFORM-COUNT (C N)

SETUP-CUT-POINT-RANDIST (P)

SQUARE (X)

STANDARD-DEVIATION (SEQUENCE)

STANDARD-ERROR-OF-THE-MEAN (SEQUENCE)

TEST-ALIAS-METHOD (N RUNS)

TEST-ALIAS-METHOD-UNIFORM (N)

TEST-ALIAS-METHOD-ZEROS

TEST-BINOMIAL (&OPTIONAL (N *N*))

TEST-DIST (F &OPTIONAL (N 1000))

TEST-GAMMA-MT-SPEED (&OPTIONAL (N *N*))

TEST-GAMMA-SPEED (&OPTIONAL (N *N*))

TEST-GIG-RANGE

TEST-GIG-SPEED

TEST-MULTINOMIAL

TEST-MULTINOMIAL1

TEST-ZIGG-SPEED (&OPTIONAL (N *N*))

TRANSFER-SIGN

UNIFORM-REFERENCE (N)

VARIANCE (SEQUENCE)

MACRO

Public

RANDOM-UNIFORM

[syntax suggar] Random variable with uniform distribution in interval [0,1]

Private

Undocumented

SD (&REST ARGS)

VAR (&REST ARGS)

VARIABLE

Public

*MT-RANDOM-STATE*

Unlike the reference implementation, we'll initialize the random state to a hopefully somewhat random & unique value.

Private

Undocumented

*KTAB*

*N*

*WTAB*

*YTAB*

CLASS

Private

Undocumented

MT-RANDOM-STATE

CONSTANT

Private

+MT-K-INVERSE-2^32F+

1/(2^32), as a floating-point number

+MT-LOWER-MASK+

least significant r bits

+MT-UPPER-MASK+

most significant w-r bits

Undocumented

+E+

+MT-K2^32+

+MT-M+

+MT-N+

+R+