Random Number Distributions

This section describes the rnd module, which contains the functions for generating random variates from a specified distribution. Samples from the distributions described in this section can be obtained using any of the random number generators from the module rng as an underlying source of randomness (see Random Number Generator Algorithms).

In the simplest cases a non-uniform distribution can be obtained analytically from the uniform distribution of a random number generator by applying an appropriate transformation, which requires the inverse cumulative distribution function. This method uses one call to the random number generator. More complicated distributions are created by the “acceptance-rejection” method, which compares the desired distribution against a distribution which is similar and known analytically. This usually requires several samples from the generator.

LGSL also provides cumulative distribution functions and inverse cumulative distribution functions, sometimes referred to as quantile functions; these are in the module randist.

rnd.gaussian(r, sigma)

This function returns a Gaussian random variate, with mean zero and standard deviation sigma. The probability distribution for Gaussian random variates is,

\[p(x) dx = {1 \over \sqrt{2 \pi \sigma^2}} \exp (-x^2 / 2\sigma^2) dx\]

for \(x\) in the range \(-\infty\) to \(+\infty\). Use the transformation \(z = \mu + x\) on the numbers returned by gaussian() to obtain a Gaussian distribution with mean \(\mu\). This function uses the Box-Mueller algorithm which requires two calls to the random number generator r.

rnd.exponential(r, mu)

This function returns a random variate from the exponential distribution with mean mu. The distribution is,

\[p(x) dx = {1 \over \mu} \exp(-x/\mu) dx\]

for \(x \ge 0\).

rnd.chisq(r, nu)

The chi-squared (\(\chi^2\)) distribution arises in statistics. If \(Y_i\) are \(n\) independent gaussian random variates with unit variance then the sum-of-squares,

\[X_i = \sum_i Y_i^2\]

has a chi-squared distribution with \(n\) degrees of freedom.

This function returns a random variate from the chi-squared distribution with \(\nu\) degrees of freedom. The distribution function is,

\[p(x) dx = {1 \over 2 \Gamma(\nu/2) } (x/2)^{\nu/2 - 1} \exp(-x/2) dx\]

for \(x \ge 0\).

rnd.laplace(r, a)

This function returns a random variate from the Laplace distribution with width a. The distribution is,

\[p(x) dx = {1 \over 2 a} \exp(-|x/a|) dx\]

for \(-\infty < x < +\infty\).

rnd.tdist(r, nu)

The t-distribution arises in statistics. If \(Y_1\) has a normal distribution and \(Y_2\) has a chi-squared distribution with \(\nu\) degrees of freedom then the ratio,

\[X = { Y_1 \over \sqrt{Y_2 / \nu} }\]

has a t-distribution \(t(x;\nu)\) with \(\nu\) degrees of freedom.

This function returns a random variate from the t-distribution. The distribution function is,

\[p(x) dx = {\Gamma((\nu + 1)/2) \over \sqrt{\pi \nu} \Gamma(\nu/2)} (1 + x^2/\nu)^{-(\nu + 1)/2} dx\]

for \(-\infty < x < +\infty\).

rnd.cauchy(r, a)

This function returns a random variate from the Cauchy distribution with scale parameter a. The probability distribution for Cauchy random variates is,

\[p(x) dx = {1 \over a\pi (1 + (x/a)^2) } dx\]

for \(-\infty < x < +\infty\). The Cauchy distribution is also known as the Lorentz distribution.

rnd.rayleigh(r, sigma)

This function returns a random variate from the Rayleigh distribution with scale parameter sigma. The distribution is,

\[p(x) dx = {x \over \sigma^2} \exp(- x^2/(2 \sigma^2)) dx\]

for \(x > 0\).

rnd.fdist(r, nu1, nu2)

The F-distribution arises in statistics. If \(Y_1\) and \(Y_2\) are chi-squared deviates with \(\nu_1\) and \(\nu_2\) degrees of freedom then the ratio,

\[X = { (Y_1 / \nu_1) \over (Y_2 / \nu_2) }\]

has an F-distribution \(F(x; \nu_1, \nu_2)\).

This function returns a random variate from the F-distribution with degrees of freedom \(\nu_1\) and \(\nu_2\). The distribution function is,

\[p(x) dx = { \Gamma((\nu_1 + \nu_2)/2) \over \Gamma(\nu_1/2) \Gamma(\nu_2/2) } \nu_1^{\nu_1/2} \nu_2^{\nu_2/2} x^{\nu_1/2 - 1} (\nu_2 + \nu_1 x)^{-\nu_1/2 -\nu_2/2}\]

for \(x \ge 0\).

rnd.gamma(r, a, b)

This function returns a random variate from the gamma distribution. The distribution function is,

\[p(x) dx = {1 \over \Gamma(a) b^a} x^{a-1} e^{-x/b} dx\]

for \(x > 0\).

The gamma distribution with an integer parameter a is known as the Erlang distribution. The variates are computed using the Marsaglia-Tsang fast gamma method.

rnd.beta(r, a, b)

This function returns a random variate from the beta distribution. The distribution function is,

\[p(x) dx = {\Gamma(a+b) \over \Gamma(a) \Gamma(b)} x^{a-1} (1-x)^{b-1} dx\]

for \(0 \le x \le 1\).

rnd.gaussian_tail(r, a, sigma)

This function provides random variates from the upper tail of a Gaussian distribution with standard deviation sigma. The values returned are larger than the lower limit a, which must be positive. The method is based on Marsaglia’s famous rectangle-wedge-tail algorithm (Ann. Math. Stat. 32, 894–899 (1961)), with this aspect explained in Knuth, v2, 3rd ed, p139,586 (exercise 11).

The probability distribution for Gaussian tail random variates is,

\[p(x) dx = {1 \over N(a;\sigma) \sqrt{2 \pi \sigma^2}} \exp \left(- \frac{x^2}{2 \sigma^2}\right) dx\]

for \(x > a\) where \(N(a; \sigma)\) is the normalization constant,

\[N(a; \sigma) = (1/2) \textrm{erfc}(a / \sqrt{2 \sigma^2}).\]

rnd.exppow(r, a, b)

This function returns a random variate from the exponential power distribution with scale parameter a and exponent b. The distribution is,

\[p(x) dx = {1 \over 2 a \Gamma(1+1/b)} \exp(-|x/a|^b) dx\]

for \(x \ge 0\). For b = 1 this reduces to the Laplace distribution. For b = 2 it has the same form as a gaussian distribution, but with \(a = \sqrt{2} \sigma\).

rnd.lognormal(r, zeta, sigma)

This function returns a random variate from the lognormal distribution. The distribution function is,

\[p(x) dx = {1 \over x \sqrt{2 \pi \sigma^2} } \exp(-(\ln(x) - \zeta)^2/2 \sigma^2) dx\]

for \(x > 0\).

rnd.binomial(r, p, n)

This function returns a random integer from the binomial distribution, the number of successes in n independent trials with probability p. The probability distribution for binomial variates is,

\[p(k) = {n! \over k! (n-k)! } p^k (1-p)^{n-k}\]

for \(0 \le k \le n\).

rnd.poisson(r, mu)

This function returns a random integer from the Poisson distribution with mean mu. The probability distribution for Poisson variates is,

\[p(k) = {\mu^k \over k!} \exp(-\mu)\]

for k >= 0.