randist.texi

This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY without ev

@noindent
for @math{-\infty < x < +\infty}.
@end deftypefun

@deftypefun double gsl_ran_logistic_pdf (double @var{x}, double @var{a})
This function computes the probability density @math{p(x)} at @var{x}
for a logistic distribution with scale parameter @var{a}, using the
formula given above.
@end deftypefun

@sp 1
@tex
\centerline{\input rand-logistic.tex}
@end tex

@deftypefun double gsl_cdf_logistic_P (double @var{x}, double @var{a})
@deftypefunx double gsl_cdf_logistic_Q (double @var{x}, double @var{a})
@deftypefunx double gsl_cdf_logistic_Pinv (double @var{P}, double @var{a})
@deftypefunx double gsl_cdf_logistic_Qinv (double @var{Q}, double @var{a})
These functions compute the cumulative distribution functions
@math{P(x)}, @math{Q(x)} and their inverses for the logistic
distribution with scale parameter @var{a}.
@end deftypefun

@page
@node The Pareto Distribution
@section The Pareto Distribution
@deftypefun double gsl_ran_pareto (const gsl_rng * @var{r}, double @var{a}, double @var{b})
@cindex Pareto distribution
This function returns a random variate from the Pareto distribution of
order @var{a}.  The distribution function is,
@tex
\beforedisplay
$$
p(x) dx = (a/b) / (x/b)^{a+1} dx
$$
\afterdisplay
@end tex
@ifinfo

@example
p(x) dx = (a/b) / (x/b)^@{a+1@} dx
@end example

@end ifinfo
@noindent
for @c{$x \ge b$}
@math{x >= b}.
@end deftypefun

@deftypefun double gsl_ran_pareto_pdf (double @var{x}, double @var{a}, double @var{b})
This function computes the probability density @math{p(x)} at @var{x}
for a Pareto distribution with exponent @var{a} and scale @var{b}, using
the formula given above.
@end deftypefun

@sp 1
@tex
\centerline{\input rand-pareto.tex}
@end tex

@deftypefun double gsl_cdf_pareto_P (double @var{x}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_pareto_Q (double @var{x}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_pareto_Pinv (double @var{P}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_pareto_Qinv (double @var{Q}, double @var{a}, double @var{b})
These functions compute the cumulative distribution functions
@math{P(x)}, @math{Q(x)} and their inverses for the Pareto
distribution with exponent @var{a} and scale @var{b}.
@end deftypefun

@page
@node Spherical Vector Distributions
@section Spherical Vector Distributions

The spherical distributions generate random vectors, located on a
spherical surface.  They can be used as random directions, for example in
the steps of a random walk.

@deftypefun void gsl_ran_dir_2d (const gsl_rng * @var{r}, double * @var{x}, double * @var{y})
@deftypefunx void gsl_ran_dir_2d_trig_method (const gsl_rng * @var{r}, double * @var{x}, double * @var{y})
@cindex 2D random direction vector
@cindex direction vector, random 2D
@cindex spherical random variates, 2D
This function returns a random direction vector @math{v} =
(@var{x},@var{y}) in two dimensions.  The vector is normalized such that
@math{|v|^2 = x^2 + y^2 = 1}.  The obvious way to do this is to take a
uniform random number between 0 and @math{2\pi} and let @var{x} and
@var{y} be the sine and cosine respectively.  Two trig functions would
have been expensive in the old days, but with modern hardware
implementations, this is sometimes the fastest way to go.  This is the
case for the Pentium (but not the case for the Sun Sparcstation).
One can avoid the trig evaluations by choosing @var{x} and
@var{y} in the interior of a unit circle (choose them at random from the
interior of the enclosing square, and then reject those that are outside
the unit circle), and then dividing by @c{$\sqrt{x^2 + y^2}$}
@math{\sqrt@{x^2 + y^2@}}.
A much cleverer approach, attributed to von Neumann (See Knuth, v2, 3rd
ed, p140, exercise 23), requires neither trig nor a square root.  In
this approach, @var{u} and @var{v} are chosen at random from the
interior of a unit circle, and then @math{x=(u^2-v^2)/(u^2+v^2)} and
@math{y=2uv/(u^2+v^2)}.
@end deftypefun

@deftypefun void gsl_ran_dir_3d (const gsl_rng * @var{r}, double * @var{x}, double * @var{y}, double * @var{z})
@cindex 3D random direction vector
@cindex direction vector, random 3D
@cindex spherical random variates, 3D
This function returns a random direction vector @math{v} =
(@var{x},@var{y},@var{z}) in three dimensions.  The vector is normalized
such that @math{|v|^2 = x^2 + y^2 + z^2 = 1}.  The method employed is
due to Robert E. Knop (CACM 13, 326 (1970)), and explained in Knuth, v2,
3rd ed, p136.  It uses the surprising fact that the distribution
projected along any axis is actually uniform (this is only true for 3
dimensions).
@end deftypefun

@deftypefun void gsl_ran_dir_nd (const gsl_rng * @var{r}, size_t @var{n}, double * @var{x})
@cindex N-dimensional random direction vector
@cindex direction vector, random N-dimensional
@cindex spherical random variates, N-dimensional

This function returns a random direction vector
@c{$v = (x_1,x_2,\ldots,x_n)$}
@math{v = (x_1,x_2,...,x_n)} in @var{n} dimensions.  The vector is normalized
such that 
@c{$|v|^2 = x_1^2 + x_2^2 + \cdots + x_n^2 = 1$}
@math{|v|^2 = x_1^2 + x_2^2 + ... + x_n^2 = 1}.  The method
uses the fact that a multivariate gaussian distribution is spherically
symmetric.  Each component is generated to have a gaussian distribution,
and then the components are normalized.  The method is described by
Knuth, v2, 3rd ed, p135--136, and attributed to G. W. Brown, Modern
Mathematics for the Engineer (1956).
@end deftypefun

@page
@node The Weibull Distribution
@section The Weibull Distribution
@deftypefun double gsl_ran_weibull (const gsl_rng * @var{r}, double @var{a}, double @var{b})
@cindex Weibull distribution
This function returns a random variate from the Weibull distribution.  The
distribution function is,
@tex
\beforedisplay
$$
p(x) dx = {b \over a^b} x^{b-1}  \exp(-(x/a)^b) dx
$$
\afterdisplay
@end tex
@ifinfo

@example
p(x) dx = @{b \over a^b@} x^@{b-1@}  \exp(-(x/a)^b) dx
@end example

@end ifinfo
@noindent
for @c{$x \ge 0$}
@math{x >= 0}.
@end deftypefun

@deftypefun double gsl_ran_weibull_pdf (double @var{x}, double @var{a}, double @var{b})
This function computes the probability density @math{p(x)} at @var{x}
for a Weibull distribution with scale @var{a} and exponent @var{b},
using the formula given above.
@end deftypefun

@sp 1
@tex
\centerline{\input rand-weibull.tex}
@end tex

@deftypefun double gsl_cdf_weibull_P (double @var{x}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_weibull_Q (double @var{x}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_weibull_Pinv (double @var{P}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_weibull_Qinv (double @var{Q}, double @var{a}, double @var{b})
These functions compute the cumulative distribution functions
@math{P(x)}, @math{Q(x)} and their inverses for the Weibull
distribution with scale @var{a} and exponent @var{b}.
@end deftypefun


@page
@node The Type-1 Gumbel Distribution
@section  The Type-1 Gumbel Distribution
@deftypefun double gsl_ran_gumbel1 (const gsl_rng * @var{r}, double @var{a}, double @var{b})
@cindex Gumbel distribution (Type 1)
@cindex Type 1 Gumbel distribution, random variates
This function returns  a random variate from the Type-1 Gumbel
distribution.  The Type-1 Gumbel distribution function is,
@tex
\beforedisplay
$$
p(x) dx = a b \exp(-(b \exp(-ax) + ax)) dx
$$
\afterdisplay
@end tex
@ifinfo

@example
p(x) dx = a b \exp(-(b \exp(-ax) + ax)) dx
@end example

@end ifinfo
@noindent
for @math{-\infty < x < \infty}. 
@end deftypefun

@deftypefun double gsl_ran_gumbel1_pdf (double @var{x}, double @var{a}, double @var{b})
This function computes 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 above.
@end deftypefun

@sp 1
@tex
\centerline{\input rand-gumbel1.tex}
@end tex

@deftypefun double gsl_cdf_gumbel1_P (double @var{x}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_gumbel1_Q (double @var{x}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_gumbel1_Pinv (double @var{P}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_gumbel1_Qinv (double @var{Q}, double @var{a}, double @var{b})
These functions compute the cumulative distribution functions
@math{P(x)}, @math{Q(x)} and their inverses for the Type-1 Gumbel
distribution with parameters @var{a} and @var{b}.
@end deftypefun


@page
@node The Type-2 Gumbel Distribution
@section  The Type-2 Gumbel Distribution
@deftypefun double gsl_ran_gumbel2 (const gsl_rng * @var{r}, double @var{a}, double @var{b})
@cindex Gumbel distribution (Type 2)
@cindex Type 2 Gumbel distribution
This function returns a random variate from the Type-2 Gumbel
distribution.  The Type-2 Gumbel distribution function is,
@tex
\beforedisplay
$$
p(x) dx = a b x^{-a-1} \exp(-b x^{-a}) dx
$$
\afterdisplay
@end tex
@ifinfo

@example
p(x) dx = a b x^@{-a-1@} \exp(-b x^@{-a@}) dx
@end example

@end ifinfo
@noindent
for @math{0 < x < \infty}.
@end deftypefun

@deftypefun double gsl_ran_gumbel2_pdf (double @var{x}, double @var{a}, double @var{b})
This function computes 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 above.
@end deftypefun

@sp 1
@tex
\centerline{\input rand-gumbel2.tex}
@end tex

@deftypefun double gsl_cdf_gumbel2_P (double @var{x}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_gumbel2_Q (double @var{x}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_gumbel2_Pinv (double @var{P}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_gumbel2_Qinv (double @var{Q}, double @var{a}, double @var{b})
These functions compute the cumulative distribution functions
@math{P(x)}, @math{Q(x)} and their inverses for the Type-2 Gumbel
distribution with parameters @var{a} and @var{b}.
@end deftypefun


@page
@node The Dirichlet Distribution
@section The Dirichlet Distribution
@deftypefun void gsl_ran_dirichlet (const gsl_rng * @var{r}, size_t @var{K}, const double @var{alpha}[], double @var{theta}[])
@cindex Dirichlet distribution 
This function returns an array of @var{K} random variates from a Dirichlet
distribution of order @var{K}-1. The distribution function is
@tex
\beforedisplay
$$
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
$$
\afterdisplay
@end tex
@ifinfo

@example
p(\theta_1, ..., \theta_K) d\theta_1 ... d\theta_K = 
  (1/Z) \prod_@{i=1@}^K \theta_i^@{\alpha_i - 1@} \delta(1 -\sum_@{i=1@}^K \theta_i) d\theta_1 ... d\theta_K
@end example

@end ifinfo
@noindent
for @c{$\theta_i \ge 0$} 
@math{theta_i >= 0}
and @c{$\alpha_i \ge 0$} 
@math{alpha_i >= 0}.  The delta function ensures that @math{\sum \theta_i = 1}.
The normalization factor @math{Z} is
@tex
\beforedisplay
$$
Z = {\prod_{i=1}^K \Gamma(\alpha_i) \over \Gamma( \sum_{i=1}^K \alpha_i)}
$$
\afterdisplay
@end tex
@ifinfo

@example
Z = @{\prod_@{i=1@}^K \Gamma(\alpha_i)@} / @{\Gamma( \sum_@{i=1@}^K \alpha_i)@}
@end example
@end ifinfo

The random variates are generated by sampling @var{K} values 
from gamma distributions with parameters 
@c{$a=\alpha_i$, $b=1$} 
@math{a=alpha_i, b=1}, 
and renormalizing. 
See A.M. Law, W.D. Kelton, @cite{Simulation Modeling and Analysis} (1991).
@end deftypefun

@deftypefun double gsl_ran_dirichlet_pdf (size_t @var{K}, const double @var{alpha}[], const double @var{theta}[]) 
This function computes the probability density 
@c{$p(\theta_1, \ldots , \theta_K)$}
@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 above.
@end deftypefun

@deftypefun double gsl_ran_dirichlet_lnpdf (size_t @var{K}, const double @var{alpha}[], const double @var{theta}[]) 
This function computes the logarithm of the probability density 
@c{$p(\theta_1, \ldots , \theta_K)$}
@math{p(\theta_1, ... , \theta_K)}
for a Dirichlet distribution with parameters 
@var{alpha}[@var{K}].
@end deftypefun

@page
@node General Discrete Distributions
@section General Discrete Distributions

Given @math{K} discrete events with different probabilities @math{P[k]},
produce a random value @math{k} consistent with its probability.

The obvious way to do this is to preprocess the probability list by
generating a cumulative probability array with @math{K+1} elements:
@tex
\beforedisplay
$$
\eqalign{
C[0] & = 0 \cr
C[k+1] &= C[k]+P[k].
}
$$
\afterdisplay
@end tex
@ifinfo

@example
  C[0] = 0 
C[k+1] = C[k]+P[k].
@end example

@end ifinfo
@noindent
Note that this construction produces @math{C[K]=1}.  Now choose a
uniform deviate @math{u} between 0 and 1, and find the value of @math{k}
such that @c{$C[k] \le u < C[k+1]$}
@math{C[k] <= u < C[k+1]}.  
Although this in principle requires of order @math{\log K} steps per
random number generation, they are fast steps, and if you use something
like @math{\lfloor uK \rfloor} as a starting point, you can often do
pretty well.

But faster methods have been devised.  Again, the idea is to preprocess
the probability list, and save the result in some form of lookup table;
then the individual calls for a random discrete event can go rapidly.
An approach invented by G. Marsaglia (Generating discrete random numbers
in a computer, Comm ACM 6, 37--38 (1963)) is very clever, and readers
interested in examples of good algorithm design are directed to this
short and well-written paper.  Unfortunately, for large @math{K},
Marsaglia's lookup table can be quite large.  

A much better approach is due to Alastair J. Walker (An efficient method
for generating discrete random variables with general distributions, ACM
Trans on Mathematical Software 3, 253--256 (1977); see also Knuth, v2,
3rd ed, p120--121,139).  This requires two lookup tables, one floating
point and one integer, but both only of size @math{K}.  After
preprocessing, the random numbers are generated in O(1) time, even for
large @math{K}.  The preprocessing suggested by Walker requires
@math{O(K^2)} effort, but that is not actually necessary, and the
implementation provided here only takes @math{O(K)} effort.  In general,
more preprocessing leads to faster generation of the individual random
numbers, but a diminishing return is reached pretty early.  Knuth points
out that the optimal preprocessing is combinatorially difficult for
large @math{K}.

This method can be used to speed up some of the discrete random number
generators below, such as the binomial distribution.  To use it for