@cindex random number distributions
@cindex cumulative distribution functions (CDFs)
@cindex CDFs, cumulative distribution functions
@cindex inverse cumulative distribution functions
@cindex quantile functions
This chapter describes functions for generating random variates and
computing their probability distributions.  Samples from the
distributions described in this chapter can be obtained using any of the
random number generators in the library as an underlying source of
randomness.  

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.  This method uses one call to
the random number generator.  More complicated distributions are created
by the @dfn{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.

The library also provides cumulative distribution functions and inverse
cumulative distribution functions, sometimes referred to as quantile
functions.  The cumulative distribution functions and their inverses are
computed separately for the upper and lower tails of the distribution,
allowing full accuracy to be retained for small results.

The functions for random variates and probability density functions
described in this section are declared in @file{gsl_randist.h}.  The
corresponding cumulative distribution functions are declared in
@file{gsl_cdf.h}.

Note that the discrete random variate functions always
return a value of type @code{unsigned int}, and on most platforms this
has a maximum value of @c{$2^{32}-1 \approx 4.29\times10^9$} 
@math{2^32-1 ~=~ 4.29e9}. They should only be called with
a safe range of parameters (where there is a negligible probability of
a variate exceeding this limit) to prevent incorrect results due to
overflow.

@menu
* Random Number Distribution Introduction::  
* The Gaussian Distribution::   
* The Gaussian Tail Distribution::  
* The Bivariate Gaussian Distribution::  
* The Exponential Distribution::  
* The Laplace Distribution::    
* The Exponential Power Distribution::  
* The Cauchy Distribution::     
* The Rayleigh Distribution::   
* The Rayleigh Tail Distribution::  
* The Landau Distribution::     
* The Levy alpha-Stable Distributions::  
* The Levy skew alpha-Stable Distribution::  
* The Gamma Distribution::      
* The Flat (Uniform) Distribution::  
* The Lognormal Distribution::  
* The Chi-squared Distribution::  
* The F-distribution::          
* The t-distribution::          
* The Beta Distribution::       
* The Logistic Distribution::   
* The Pareto Distribution::     
* Spherical Vector Distributions::  
* The Weibull Distribution::    
* The Type-1 Gumbel Distribution::  
* The Type-2 Gumbel Distribution::  
* The Dirichlet Distribution::  
* General Discrete Distributions::  
* The Poisson Distribution::    
* The Bernoulli Distribution::  
* The Binomial Distribution::   
* The Multinomial Distribution::  
* The Negative Binomial Distribution::  
* The Pascal Distribution::     
* The Geometric Distribution::  
* The Hypergeometric Distribution::  
* The Logarithmic Distribution::  
* Shuffling and Sampling::      
* Random Number Distribution Examples::  
* Random Number Distribution References and Further Reading::  
@end menu

@node Random Number Distribution Introduction
@section Introduction

Continuous random number distributions are defined by a probability
density function, @math{p(x)}, such that the probability of @math{x}
occurring in the infinitesimal range @math{x} to @math{x+dx} is @c{$p\,dx$}
@math{p dx}.

The cumulative distribution function for the lower tail @math{P(x)} is
defined by the integral,
@tex
\beforedisplay
$$
P(x) = \int_{-\infty}^{x} dx' p(x')
$$
\afterdisplay
@end tex
@ifinfo

@example
P(x) = \int_@{-\infty@}^@{x@} dx' p(x')
@end example

@end ifinfo
@noindent
and gives the probability of a variate taking a value less than @math{x}.

The cumulative distribution function for the upper tail @math{Q(x)} is
defined by the integral,
@tex
\beforedisplay
$$
Q(x) = \int_{x}^{+\infty} dx' p(x')
$$
\afterdisplay
@end tex
@ifinfo

@example
Q(x) = \int_@{x@}^@{+\infty@} dx' p(x')
@end example

@end ifinfo
@noindent
and gives the probability of a variate taking a value greater than @math{x}.

The upper and lower cumulative distribution functions are related by
@math{P(x) + Q(x) = 1} and satisfy @c{$0 \le P(x) \le 1$}
@math{0 <= P(x) <= 1}, @c{$0 \le Q(x) \le 1$}
@math{0 <= Q(x) <= 1}.

The inverse cumulative distributions, @c{$x=P^{-1}(P)$}
@math{x=P^@{-1@}(P)} and @c{$x=Q^{-1}(Q)$}
@math{x=Q^@{-1@}(Q)} give the values of @math{x}
which correspond to a specific value of @math{P} or @math{Q}.  
They can be used to find confidence limits from probability values.

For discrete distributions the probability of sampling the integer
value @math{k} is given by @math{p(k)}, where @math{\sum_k p(k) = 1}.
The cumulative distribution for the lower tail @math{P(k)} of a
discrete distribution is defined as,
@tex
\beforedisplay
$$
P(k) = \sum_{i \le k} p(i) 
$$
\afterdisplay
@end tex
@ifinfo

@example
P(k) = \sum_@{i <= k@} p(i)
@end example

@end ifinfo
@noindent
where the sum is over the allowed range of the distribution less than
or equal to @math{k}.  

The cumulative distribution for the upper tail of a discrete
distribution @math{Q(k)} is defined as
@tex
\beforedisplay
$$
Q(k) = \sum_{i > k} p(i) 
$$
\afterdisplay
@end tex
@ifinfo

@example
Q(k) = \sum_@{i > k@} p(i)
@end example

@end ifinfo
@noindent
giving the sum of probabilities for all values greater than @math{k}.
These two definitions satisfy the identity @math{P(k)+Q(k)=1}.

If the range of the distribution is 1 to @math{n} inclusive then
@math{P(n)=1}, @math{Q(n)=0} while @math{P(1) = p(1)},
@math{Q(1)=1-p(1)}.

@page
@node The Gaussian Distribution
@section The Gaussian Distribution
@deftypefun double gsl_ran_gaussian (const gsl_rng * @var{r}, double @var{sigma})
@cindex Gaussian distribution
This function returns a Gaussian random variate, with mean zero and
standard deviation @var{sigma}.  The probability distribution for
Gaussian random variates is,
@tex
\beforedisplay
$$
p(x) dx = {1 \over \sqrt{2 \pi \sigma^2}} \exp (-x^2 / 2\sigma^2) dx
$$
\afterdisplay
@end tex
@ifinfo

@example
p(x) dx = @{1 \over \sqrt@{2 \pi \sigma^2@}@} \exp (-x^2 / 2\sigma^2) dx
@end example

@end ifinfo
@noindent
for @math{x} in the range @math{-\infty} to @math{+\infty}.  Use the
transformation @math{z = \mu + x} on the numbers returned by
@code{gsl_ran_gaussian} to obtain a Gaussian distribution with mean
@math{\mu}.  This function uses the Box-Mueller algorithm which requires two
calls to the random number generator @var{r}.
@end deftypefun

@deftypefun double gsl_ran_gaussian_pdf (double @var{x}, double @var{sigma})
This function computes the probability density @math{p(x)} at @var{x}
for a Gaussian distribution with standard deviation @var{sigma}, using
the formula given above.
@end deftypefun

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

@deftypefun double gsl_ran_gaussian_ziggurat (const gsl_rng * @var{r}, double @var{sigma})
@deftypefunx double gsl_ran_gaussian_ratio_method (const gsl_rng * @var{r}, double @var{sigma})
This function computes a Gaussian random variate using the alternative
Marsaglia-Tsang ziggurat and Kinderman-Monahan-Leva ratio methods.  The
Ziggurat algorithm is the fastest available algorithm in most cases.
@end deftypefun

@deftypefun double gsl_ran_ugaussian (const gsl_rng * @var{r})
@deftypefunx double gsl_ran_ugaussian_pdf (double @var{x})
@deftypefunx double gsl_ran_ugaussian_ratio_method (const gsl_rng * @var{r})
These functions compute results for the unit Gaussian distribution.  They
are equivalent to the functions above with a standard deviation of one,
@var{sigma} = 1.
@end deftypefun

@deftypefun double gsl_cdf_gaussian_P (double @var{x}, double @var{sigma})
@deftypefunx double gsl_cdf_gaussian_Q (double @var{x}, double @var{sigma})
@deftypefunx double gsl_cdf_gaussian_Pinv (double @var{P}, double @var{sigma})
@deftypefunx double gsl_cdf_gaussian_Qinv (double @var{Q}, double @var{sigma})
These functions compute the cumulative distribution functions
@math{P(x)}, @math{Q(x)} and their inverses for the Gaussian
distribution with standard deviation @var{sigma}.
@end deftypefun

@deftypefun double gsl_cdf_ugaussian_P (double @var{x})
@deftypefunx double gsl_cdf_ugaussian_Q (double @var{x})
@deftypefunx double gsl_cdf_ugaussian_Pinv (double @var{P})
@deftypefunx double gsl_cdf_ugaussian_Qinv (double @var{Q})
These functions compute the cumulative distribution functions
@math{P(x)}, @math{Q(x)} and their inverses for the unit Gaussian
distribution.
@end deftypefun

@page
@node The Gaussian Tail Distribution
@section The Gaussian Tail Distribution
@deftypefun double gsl_ran_gaussian_tail (const gsl_rng * @var{r}, double @var{a}, double @var{sigma})
@cindex Gaussian Tail distribution
This function provides random variates from the upper tail of a Gaussian
distribution with standard deviation @var{sigma}.  The values returned
are larger than the lower limit @var{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,
@tex
\beforedisplay
$$
p(x) dx = {1 \over N(a;\sigma) \sqrt{2 \pi \sigma^2}} \exp (- x^2 / 2\sigma^2) dx
$$
\afterdisplay
@end tex
@ifinfo

@example
p(x) dx = @{1 \over N(a;\sigma) \sqrt@{2 \pi \sigma^2@}@} \exp (- x^2/(2 \sigma^2)) dx
@end example

@end ifinfo
@noindent
for @math{x > a} where @math{N(a;\sigma)} is the normalization constant,
@tex
\beforedisplay
$$
N(a;\sigma) = {1 \over 2} \hbox{erfc}\left({a \over \sqrt{2 \sigma^2}}\right).
$$
\afterdisplay
@end tex
@ifinfo

@example
N(a;\sigma) = (1/2) erfc(a / sqrt(2 sigma^2)).
@end example
@end ifinfo

@end deftypefun

@deftypefun double gsl_ran_gaussian_tail_pdf (double @var{x}, double @var{a}, double @var{sigma})
This function computes the probability density @math{p(x)} at @var{x}
for a Gaussian tail distribution with standard deviation @var{sigma} and
lower limit @var{a}, using the formula given above.
@end deftypefun

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

@deftypefun double gsl_ran_ugaussian_tail (const gsl_rng * @var{r}, double @var{a})
@deftypefunx double gsl_ran_ugaussian_tail_pdf (double @var{x}, double @var{a})
These functions compute results for the tail of a unit Gaussian
distribution.  They are equivalent to the functions above with a standard
deviation of one, @var{sigma} = 1.
@end deftypefun


@page
@node The Bivariate Gaussian Distribution
@section The Bivariate Gaussian Distribution

@deftypefun void gsl_ran_bivariate_gaussian (const gsl_rng * @var{r}, double @var{sigma_x}, double @var{sigma_y}, double @var{rho}, double * @var{x}, double * @var{y})
@cindex Bivariate Gaussian distribution
@cindex two dimensional Gaussian distribution
@cindex Gaussian distribution, bivariate
This function generates a pair of correlated Gaussian variates, with
mean zero, correlation coefficient @var{rho} and standard deviations
@var{sigma_x} and @var{sigma_y} in the @math{x} and @math{y} directions.
The probability distribution for bivariate Gaussian random variates is,
@tex
\beforedisplay
$$
p(x,y) dx dy = {1 \over 2 \pi \sigma_x \sigma_y \sqrt{1-\rho^2}} \exp \left(-{(x^2/\sigma_x^2 + y^2/\sigma_y^2 - 2 \rho x y/(\sigma_x\sigma_y)) \over 2(1-\rho^2)}\right) dx dy
$$
\afterdisplay
@end tex
@ifinfo

@example
p(x,y) dx dy = @{1 \over 2 \pi \sigma_x \sigma_y \sqrt@{1-\rho^2@}@} \exp (-(x^2/\sigma_x^2 + y^2/\sigma_y^2 - 2 \rho x y/(\sigma_x\sigma_y))/2(1-\rho^2)) dx dy
@end example

@end ifinfo
@noindent
for @math{x,y} in the range @math{-\infty} to @math{+\infty}.  The
correlation coefficient @var{rho} should lie between @math{1} and
@math{-1}.
@end deftypefun

@deftypefun double gsl_ran_bivariate_gaussian_pdf (double @var{x}, double @var{y}, double @var{sigma_x}, double @var{sigma_y}, double @var{rho})
This function computes the probability density @math{p(x,y)} at
(@var{x},@var{y}) for a bivariate Gaussian distribution with standard
deviations @var{sigma_x}, @var{sigma_y} and correlation coefficient
@var{rho}, using the formula given above.
@end deftypefun

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

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

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

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

@deftypefun double gsl_ran_exponential_pdf (double @var{x}, double @var{mu})
This function computes the probability density @math{p(x)} at @var{x}
for an exponential distribution with mean @var{mu}, using the formula
given above.
@end deftypefun

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

@deftypefun double gsl_cdf_exponential_P (double @var{x}, double @var{mu})
@deftypefunx double gsl_cdf_exponential_Q (double @var{x}, double @var{mu})
@deftypefunx double gsl_cdf_exponential_Pinv (double @var{P}, double @var{mu})
@deftypefunx double gsl_cdf_exponential_Qinv (double @var{Q}, double @var{mu})
These functions compute the cumulative distribution functions
@math{P(x)}, @math{Q(x)} and their inverses for the exponential
distribution with mean @var{mu}.
@end deftypefun

@page
@node The Laplace Distribution
@section The Laplace Distribution
@deftypefun double gsl_ran_laplace (const gsl_rng * @var{r}, double @var{a})
@cindex two-sided exponential distribution
@cindex Laplace distribution
This function returns a random variate from the Laplace distribution
with width @var{a}.  The distribution is,
@tex
\beforedisplay
$$
p(x) dx = {1 \over 2 a}  \exp(-|x/a|) dx
$$
\afterdisplay
@end tex
@ifinfo

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

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

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

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

@deftypefun double gsl_cdf_laplace_P (double @var{x}, double @var{a})
@deftypefunx double gsl_cdf_laplace_Q (double @var{x}, double @var{a})
@deftypefunx double gsl_cdf_laplace_Pinv (double @var{P}, double @var{a})
@deftypefunx double gsl_cdf_laplace_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 Laplace
distribution with width @var{a}.
@end deftypefun


@page
@node The Exponential Power Distribution
@section The Exponential Power Distribution
@deftypefun double gsl_ran_exppow (const gsl_rng * @var{r}, double @var{a}, double @var{b})
@cindex Exponential power distribution
This function returns a random variate from the exponential power distribution
with scale parameter @var{a} and exponent @var{b}.  The distribution is,
@tex
\beforedisplay
$$
p(x) dx = {1 \over 2 a \Gamma(1+1/b)} \exp(-|x/a|^b) dx
$$
\afterdisplay
@end tex
@ifinfo

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

@end ifinfo
@noindent
for @c{$x \ge 0$}
@math{x >= 0}.  For @math{b = 1} this reduces to the Laplace
distribution.  For @math{b = 2} it has the same form as a gaussian
distribution, but with @c{$a = \sqrt{2} \sigma$}
@math{a = \sqrt@{2@} \sigma}.
@end deftypefun

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

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

@deftypefun double gsl_cdf_exppow_P (double @var{x}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_exppow_Q (double @var{x}, double @var{a}, double @var{b})
These functions compute the cumulative distribution functions
@math{P(x)}, @math{Q(x)} for the exponential power distribution with
parameters @var{a} and @var{b}.
@end deftypefun


@page
@node The Cauchy Distribution
@section The Cauchy Distribution
@deftypefun double gsl_ran_cauchy (const gsl_rng * @var{r}, double @var{a})
@cindex Cauchy distribution
This function returns a random variate from the Cauchy distribution with
scale parameter @var{a}.  The probability distribution for Cauchy
random variates is,
@tex
\beforedisplay
$$
p(x) dx = {1 \over a\pi (1 + (x/a)^2) } dx
$$
\afterdisplay
@end tex
@ifinfo

@example
p(x) dx = @{1 \over a\pi (1 + (x/a)^2) @} dx
@end example

@end ifinfo
@noindent
for @math{x} in the range @math{-\infty} to @math{+\infty}.  The Cauchy
distribution is also known as the Lorentz distribution.
@end deftypefun

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

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

@deftypefun double gsl_cdf_cauchy_P (double @var{x}, double @var{a})
@deftypefunx double gsl_cdf_cauchy_Q (double @var{x}, double @var{a})
@deftypefunx double gsl_cdf_cauchy_Pinv (double @var{P}, double @var{a})
@deftypefunx double gsl_cdf_cauchy_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 Cauchy
distribution with scale parameter @var{a}.
@end deftypefun


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

@example
p(x) dx = @{x \over \sigma^2@} \exp(- x^2/(2 \sigma^2)) dx
@end example

@end ifinfo
@noindent
for @math{x > 0}.
@end deftypefun

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

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

@deftypefun double gsl_cdf_rayleigh_P (double @var{x}, double @var{sigma})
@deftypefunx double gsl_cdf_rayleigh_Q (double @var{x}, double @var{sigma})
@deftypefunx double gsl_cdf_rayleigh_Pinv (double @var{P}, double @var{sigma})
@deftypefunx double gsl_cdf_rayleigh_Qinv (double @var{Q}, double @var{sigma})
These functions compute the cumulative distribution functions
@math{P(x)}, @math{Q(x)} and their inverses for the Rayleigh
distribution with scale parameter @var{sigma}.
@end deftypefun


@page
@node The Rayleigh Tail Distribution
@section The Rayleigh Tail Distribution
@deftypefun double gsl_ran_rayleigh_tail (const gsl_rng * @var{r}, double @var{a}, double @var{sigma})
@cindex Rayleigh Tail distribution
This function returns a random variate from the tail of the Rayleigh
distribution with scale parameter @var{sigma} and a lower limit of
@var{a}.  The distribution is,
@tex
\beforedisplay
$$
p(x) dx = {x \over \sigma^2} \exp ((a^2 - x^2) /(2 \sigma^2)) dx
$$
\afterdisplay
@end tex
@ifinfo

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

@end ifinfo
@noindent
for @math{x > a}.
@end deftypefun

@deftypefun double gsl_ran_rayleigh_tail_pdf (double @var{x}, double @var{a}, double @var{sigma})
This function computes the probability density @math{p(x)} at @var{x}
for a Rayleigh tail distribution with scale parameter @var{sigma} and
lower limit @var{a}, using the formula given above.
@end deftypefun

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

@page
@node The Landau Distribution
@section The Landau Distribution
@deftypefun double gsl_ran_landau (const gsl_rng * @var{r})
@cindex Landau distribution
This function returns a random variate from the Landau distribution.  The
probability distribution for Landau random variates is defined
analytically by the complex integral,
@tex
\beforedisplay
$$
p(x) = 
{1 \over {2 \pi i}} \int_{c-i\infty}^{c+i\infty} ds\, \exp(s \log(s) + x s) 
$$
\afterdisplay
@end tex
@ifinfo

@example
p(x) = (1/(2 \pi i)) \int_@{c-i\infty@}^@{c+i\infty@} ds exp(s log(s) + x s) 
@end example
@end ifinfo
For numerical purposes it is more convenient to use the following
equivalent form of the integral,
@tex
\beforedisplay
$$
p(x) = (1/\pi) \int_0^\infty dt \exp(-t \log(t) - x t) \sin(\pi t).
$$
\afterdisplay
@end tex
@ifinfo

@example
p(x) = (1/\pi) \int_0^\infty dt \exp(-t \log(t) - x t) \sin(\pi t).
@end example
@end ifinfo
@end deftypefun

@deftypefun double gsl_ran_landau_pdf (double @var{x})
This function computes the probability density @math{p(x)} at @var{x}
for the Landau distribution using an approximation to the formula given
above.
@end deftypefun

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

@page
@node The Levy alpha-Stable Distributions
@section The Levy alpha-Stable Distributions
@deftypefun double gsl_ran_levy (const gsl_rng * @var{r}, double @var{c}, double @var{alpha})
@cindex Levy distribution
This function returns a random variate from the Levy symmetric stable
distribution with scale @var{c} and exponent @var{alpha}.  The symmetric
stable probability distribution is defined by a fourier transform,
@tex
\beforedisplay
$$
p(x) = {1 \over 2 \pi} \int_{-\infty}^{+\infty} dt \exp(-it x - |c t|^\alpha)
$$
\afterdisplay
@end tex
@ifinfo

@example
p(x) = @{1 \over 2 \pi@} \int_@{-\infty@}^@{+\infty@} dt \exp(-it x - |c t|^alpha)
@end example

@end ifinfo
@noindent
There is no explicit solution for the form of @math{p(x)} and the
library does not define a corresponding @code{pdf} function.  For
@math{\alpha = 1} the distribution reduces to the Cauchy distribution.  For
@math{\alpha = 2} it is a Gaussian distribution with @c{$\sigma = \sqrt{2} c$} 
@math{\sigma = \sqrt@{2@} c}.  For @math{\alpha < 1} the tails of the
distribution become extremely wide.

The algorithm only works for @c{$0 < \alpha \le 2$}
@math{0 < alpha <= 2}.
@end deftypefun

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

@page
@node The Levy skew alpha-Stable Distribution
@section The Levy skew alpha-Stable Distribution

@deftypefun double gsl_ran_levy_skew (const gsl_rng * @var{r}, double @var{c}, double @var{alpha}, double @var{beta})
@cindex Levy distribution, skew
@cindex Skew Levy distribution
This function returns a random variate from the Levy skew stable
distribution with scale @var{c}, exponent @var{alpha} and skewness
parameter @var{beta}.  The skewness parameter must lie in the range
@math{[-1,1]}.  The Levy skew stable probability distribution is defined
by a fourier transform,
@tex
\beforedisplay
$$
p(x) = {1 \over 2 \pi} \int_{-\infty}^{+\infty} dt \exp(-it x - |c t|^\alpha (1-i \beta \sign(t) \tan(\pi\alpha/2)))
$$
\afterdisplay
@end tex
@ifinfo

@example
p(x) = @{1 \over 2 \pi@} \int_@{-\infty@}^@{+\infty@} dt \exp(-it x - |c t|^alpha (1-i beta sign(t) tan(pi alpha/2)))
@end example

@end ifinfo
@noindent
When @math{\alpha = 1} the term @math{\tan(\pi \alpha/2)} is replaced by
@math{-(2/\pi)\log|t|}.  There is no explicit solution for the form of
@math{p(x)} and the library does not define a corresponding @code{pdf}
function.  For @math{\alpha = 2} the distribution reduces to a Gaussian
distribution with @c{$\sigma = \sqrt{2} c$} 
@math{\sigma = \sqrt@{2@} c} and the skewness parameter has no effect.  
For @math{\alpha < 1} the tails of the distribution become extremely
wide.  The symmetric distribution corresponds to @math{\beta =
0}.

The algorithm only works for @c{$0 < \alpha \le 2$}
@math{0 < alpha <= 2}.
@end deftypefun

The Levy alpha-stable distributions have the property that if @math{N}
alpha-stable variates are drawn from the distribution @math{p(c, \alpha,
\beta)} then the sum @math{Y = X_1 + X_2 + \dots + X_N} will also be
distributed as an alpha-stable variate,
@c{$p(N^{1/\alpha} c, \alpha, \beta)$}  
@math{p(N^(1/\alpha) c, \alpha, \beta)}.

@comment PDF not available because there is no analytic expression for it
@comment
@comment @deftypefun double gsl_ran_levy_pdf (double @var{x}, double @var{mu})
@comment This function computes the probability density @math{p(x)} at @var{x}
@comment for a symmetric Levy distribution with scale parameter @var{mu} and
@comment exponent @var{a}, using the formula given above.
@comment @end deftypefun

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

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

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

@end ifinfo
@noindent
for @math{x > 0}.
@comment If @xmath{X} and @xmath{Y} are independent gamma-distributed random
@comment variables of order @xmath{a} and @xmath{b}, then @xmath{X+Y} has a gamma
@comment distribution of order @xmath{a+b}.

@cindex Erlang distribution
The gamma distribution with an integer parameter @var{a} is known as the Erlang distribution.

The variates are computed using the Marsaglia-Tsang fast gamma method.
This function for this method was previously called
@code{gsl_ran_gamma_mt} and can still be accessed using this name.
@end deftypefun

@deftypefun double gsl_ran_gamma_knuth (const gsl_rng * @var{r}, double @var{a}, double @var{b})
This function returns a gamma variate using the algorithms from Knuth (vol 2).
@end deftypefun

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

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

@deftypefun double gsl_cdf_gamma_P (double @var{x}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_gamma_Q (double @var{x}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_gamma_Pinv (double @var{P}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_gamma_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 gamma
distribution with parameters @var{a} and @var{b}.
@end deftypefun

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

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

@end ifinfo
@noindent
if @c{$a \le x < b$}
@math{a <= x < b} and 0 otherwise.
@end deftypefun

@deftypefun double gsl_ran_flat_pdf (double @var{x}, double @var{a}, double @var{b})
This function computes the probability density @math{p(x)} at @var{x}
for a uniform distribution from @var{a} to @var{b}, using the formula
given above.
@end deftypefun

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

@deftypefun double gsl_cdf_flat_P (double @var{x}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_flat_Q (double @var{x}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_flat_Pinv (double @var{P}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_flat_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 a uniform distribution
from @var{a} to @var{b}.
@end deftypefun


@page
@node The Lognormal Distribution
@section The Lognormal Distribution
@deftypefun double gsl_ran_lognormal (const gsl_rng * @var{r}, double @var{zeta}, double @var{sigma})
@cindex Lognormal distribution
This function returns a random variate from the lognormal
distribution.  The distribution function is,
@tex
\beforedisplay
$$
p(x) dx = {1 \over x \sqrt{2 \pi \sigma^2}} \exp(-(\ln(x) - \zeta)^2/2 \sigma^2) dx
$$
\afterdisplay
@end tex
@ifinfo

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

@end ifinfo
@noindent
for @math{x > 0}.
@end deftypefun

@deftypefun double gsl_ran_lognormal_pdf (double @var{x}, double @var{zeta}, double @var{sigma})
This function computes the probability density @math{p(x)} at @var{x}
for a lognormal distribution with parameters @var{zeta} and @var{sigma},
using the formula given above.
@end deftypefun

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

@deftypefun double gsl_cdf_lognormal_P (double @var{x}, double @var{zeta}, double @var{sigma})
@deftypefunx double gsl_cdf_lognormal_Q (double @var{x}, double @var{zeta}, double @var{sigma})
@deftypefunx double gsl_cdf_lognormal_Pinv (double @var{P}, double @var{zeta}, double @var{sigma})
@deftypefunx double gsl_cdf_lognormal_Qinv (double @var{Q}, double @var{zeta}, double @var{sigma})
These functions compute the cumulative distribution functions
@math{P(x)}, @math{Q(x)} and their inverses for the lognormal
distribution with parameters @var{zeta} and @var{sigma}.
@end deftypefun


@page
@node The Chi-squared Distribution
@section The Chi-squared Distribution
The chi-squared distribution arises in statistics.  If @math{Y_i} are
@math{n} independent gaussian random variates with unit variance then the
sum-of-squares,
@tex
\beforedisplay
$$
X_i = \sum_i Y_i^2
$$
\afterdisplay
@end tex
@ifinfo

@example
X_i = \sum_i Y_i^2
@end example

@end ifinfo
@noindent
has a chi-squared distribution with @math{n} degrees of freedom.

@deftypefun double gsl_ran_chisq (const gsl_rng * @var{r}, double @var{nu})
@cindex Chi-squared distribution
This function returns a random variate from the chi-squared distribution
with @var{nu} degrees of freedom. The distribution function is,
@tex
\beforedisplay
$$
p(x) dx = {1 \over 2 \Gamma(\nu/2) } (x/2)^{\nu/2 - 1} \exp(-x/2) dx
$$
\afterdisplay
@end tex
@ifinfo

@example
p(x) dx = @{1 \over 2 \Gamma(\nu/2) @} (x/2)^@{\nu/2 - 1@} \exp(-x/2) dx
@end example

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

@deftypefun double gsl_ran_chisq_pdf (double @var{x}, double @var{nu})
This function computes the probability density @math{p(x)} at @var{x}
for a chi-squared distribution with @var{nu} degrees of freedom, using
the formula given above.
@end deftypefun

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

@deftypefun double gsl_cdf_chisq_P (double @var{x}, double @var{nu})
@deftypefunx double gsl_cdf_chisq_Q (double @var{x}, double @var{nu})
@deftypefunx double gsl_cdf_chisq_Pinv (double @var{P}, double @var{nu})
@deftypefunx double gsl_cdf_chisq_Qinv (double @var{Q}, double @var{nu})
These functions compute the cumulative distribution functions
@math{P(x)}, @math{Q(x)} and their inverses for the chi-squared
distribution with @var{nu} degrees of freedom.
@end deftypefun


@page
@node The F-distribution
@section The F-distribution
The F-distribution arises in statistics.  If @math{Y_1} and @math{Y_2}
are chi-squared deviates with @math{\nu_1} and @math{\nu_2} degrees of
freedom then the ratio,
@tex
\beforedisplay
$$
X = { (Y_1 / \nu_1) \over (Y_2 / \nu_2) }
$$
\afterdisplay
@end tex
@ifinfo

@example
X = @{ (Y_1 / \nu_1) \over (Y_2 / \nu_2) @}
@end example

@end ifinfo
@noindent
has an F-distribution @math{F(x;\nu_1,\nu_2)}.

@deftypefun double gsl_ran_fdist (const gsl_rng * @var{r}, double @var{nu1}, double @var{nu2})
@cindex F-distribution
This function returns a random variate from the F-distribution with degrees of freedom @var{nu1} and @var{nu2}. The distribution function is,
@tex
\beforedisplay
$$
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}
$$
\afterdisplay
@end tex
@ifinfo

@example
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@}
@end example

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

@deftypefun double gsl_ran_fdist_pdf (double @var{x}, double @var{nu1}, double @var{nu2})
This function computes the probability density @math{p(x)} at @var{x}
for an F-distribution with @var{nu1} and @var{nu2} degrees of freedom,
using the formula given above.
@end deftypefun

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

@deftypefun double gsl_cdf_fdist_P (double @var{x}, double @var{nu1}, double @var{nu2})
@deftypefunx double gsl_cdf_fdist_Q (double @var{x}, double @var{nu1}, double @var{nu2})
@deftypefunx double gsl_cdf_fdist_Pinv (double @var{P}, double @var{nu1}, double @var{nu2})
@deftypefunx double gsl_cdf_fdist_Qinv (double @var{Q}, double @var{nu1}, double @var{nu2})
These functions compute the cumulative distribution functions
@math{P(x)}, @math{Q(x)} and their inverses for the F-distribution
with @var{nu1} and @var{nu2} degrees of freedom.
@end deftypefun

@page
@node The t-distribution
@section The t-distribution
The t-distribution arises in statistics.  If @math{Y_1} has a normal
distribution and @math{Y_2} has a chi-squared distribution with
@math{\nu} degrees of freedom then the ratio,
@tex
\beforedisplay
$$
X = { Y_1 \over \sqrt{Y_2 / \nu} }
$$
\afterdisplay
@end tex
@ifinfo

@example
X = @{ Y_1 \over \sqrt@{Y_2 / \nu@} @}
@end example

@end ifinfo
@noindent
has a t-distribution @math{t(x;\nu)} with @math{\nu} degrees of freedom.

@deftypefun double gsl_ran_tdist (const gsl_rng * @var{r}, double @var{nu})
@cindex t-distribution
@cindex Student t-distribution
This function returns a random variate from the t-distribution.  The
distribution function is,
@tex
\beforedisplay
$$
p(x) dx = {\Gamma((\nu + 1)/2) \over \sqrt{\pi \nu} \Gamma(\nu/2)}
   (1 + x^2/\nu)^{-(\nu + 1)/2} dx
$$
\afterdisplay
@end tex
@ifinfo

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

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

@deftypefun double gsl_ran_tdist_pdf (double @var{x}, double @var{nu})
This function computes the probability density @math{p(x)} at @var{x}
for a t-distribution with @var{nu} degrees of freedom, using the formula
given above.
@end deftypefun

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

@deftypefun double gsl_cdf_tdist_P (double @var{x}, double @var{nu})
@deftypefunx double gsl_cdf_tdist_Q (double @var{x}, double @var{nu})
@deftypefunx double gsl_cdf_tdist_Pinv (double @var{P}, double @var{nu})
@deftypefunx double gsl_cdf_tdist_Qinv (double @var{Q}, double @var{nu})
These functions compute the cumulative distribution functions
@math{P(x)}, @math{Q(x)} and their inverses for the t-distribution
with @var{nu} degrees of freedom.
@end deftypefun

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

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

@end ifinfo
@noindent
for @c{$0 \le x \le 1$}
@math{0 <= x <= 1}.
@end deftypefun

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

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

@deftypefun double gsl_cdf_beta_P (double @var{x}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_beta_Q (double @var{x}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_beta_Pinv (double @var{P}, double @var{a}, double @var{b})
@deftypefunx double gsl_cdf_beta_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 beta
distribution with parameters @var{a} and @var{b}.
@end deftypefun

@page
@node The Logistic Distribution
@section The Logistic Distribution

@deftypefun double gsl_ran_logistic (const gsl_rng * @var{r}, double @var{a})
@cindex Logistic distribution
This function returns a random variate from the logistic
distribution.  The distribution function is,
@tex
\beforedisplay
$$
p(x) dx = { \exp(-x/a) \over a (1 + \exp(-x/a))^2 } dx
$$
\afterdisplay
@end tex
@ifinfo

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

@end ifinfo
@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
something like the Poisson Distribution, a modification would have to
be made, since it only takes a finite set of @math{K} outcomes.

@deftypefun {gsl_ran_discrete_t *} gsl_ran_discrete_preproc (size_t @var{K}, const double * @var{P})
@cindex Discrete random numbers
@cindex Discrete random numbers, preprocessing
This function returns a pointer to a structure that contains the lookup
table for the discrete random number generator.  The array @var{P}[] contains
the probabilities of the discrete events; these array elements must all be 
positive, but they needn't add up to one (so you can think of them more
generally as ``weights'')---the preprocessor will normalize appropriately.
This return value is used
as an argument for the @code{gsl_ran_discrete} function below.
@end deftypefun

@deftypefun {size_t} gsl_ran_discrete (const gsl_rng * @var{r}, const gsl_ran_discrete_t * @var{g})
@cindex Discrete random numbers
After the preprocessor, above, has been called, you use this function to
get the discrete random numbers.
@end deftypefun

@deftypefun {double} gsl_ran_discrete_pdf (size_t @var{k}, const gsl_ran_discrete_t * @var{g})
@cindex Discrete random numbers
Returns the probability @math{P[k]} of observing the variable @var{k}.
Since @math{P[k]} is not stored as part of the lookup table, it must be
recomputed; this computation takes @math{O(K)}, so if @var{K} is large
and you care about the original array @math{P[k]} used to create the
lookup table, then you should just keep this original array @math{P[k]}
around.
@end deftypefun

@deftypefun {void} gsl_ran_discrete_free (gsl_ran_discrete_t * @var{g})
@cindex Discrete random numbers
De-allocates the lookup table pointed to by @var{g}.
@end deftypefun

@page
@node The Poisson Distribution
@section The Poisson Distribution
@deftypefun {unsigned int} gsl_ran_poisson (const gsl_rng * @var{r}, double @var{mu})
@cindex Poisson random numbers
This function returns a random integer from the Poisson distribution
with mean @var{mu}.  The probability distribution for Poisson variates is,
@tex
\beforedisplay
$$
p(k) = {\mu^k \over k!} \exp(-\mu)
$$
\afterdisplay
@end tex
@ifinfo

@example
p(k) = @{\mu^k \over k!@} \exp(-\mu)
@end example

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

@deftypefun double gsl_ran_poisson_pdf (unsigned int @var{k}, double @var{mu})
This function computes the probability @math{p(k)} of obtaining  @var{k}
from a Poisson distribution with mean @var{mu}, using the formula
given above.
@end deftypefun

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

@deftypefun double gsl_cdf_poisson_P (unsigned int @var{k}, double @var{mu})
@deftypefunx double gsl_cdf_poisson_Q (unsigned int @var{k}, double @var{mu})
These functions compute the cumulative distribution functions
@math{P(k)}, @math{Q(k)} for the Poisson distribution with parameter
@var{mu}.
@end deftypefun


@page
@node The Bernoulli Distribution
@section The Bernoulli Distribution
@deftypefun {unsigned int} gsl_ran_bernoulli (const gsl_rng * @var{r}, double @var{p})
@cindex Bernoulli trial, random variates
This function returns either 0 or 1, the result of a Bernoulli trial
with probability @var{p}.  The probability distribution for a Bernoulli
trial is,
@tex
\beforedisplay
$$
\eqalign{
p(0) & = 1 - p \cr
p(1) & = p
}
$$
\afterdisplay
@end tex
@ifinfo

@example
p(0) = 1 - p
p(1) = p
@end example
@end ifinfo

@end deftypefun

@deftypefun double gsl_ran_bernoulli_pdf (unsigned int @var{k}, double @var{p})
This function computes the probability @math{p(k)} of obtaining
@var{k} from a Bernoulli distribution with probability parameter
@var{p}, using the formula given above.
@end deftypefun

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

@page
@node The Binomial Distribution
@section The Binomial Distribution
@deftypefun {unsigned int} gsl_ran_binomial (const gsl_rng * @var{r}, double @var{p}, unsigned int @var{n})
@cindex Binomial random variates
This function returns a random integer from the binomial distribution,
the number of successes in @var{n} independent trials with probability
@var{p}.  The probability distribution for binomial variates is,
@tex
\beforedisplay
$$
p(k) = {n! \over k! (n-k)!} p^k (1-p)^{n-k}
$$
\afterdisplay
@end tex
@ifinfo

@example
p(k) = @{n! \over k! (n-k)! @} p^k (1-p)^@{n-k@}
@end example

@end ifinfo
@noindent
for @c{$0 \le k \le n$}
@math{0 <= k <= n}.
@end deftypefun

@deftypefun double gsl_ran_binomial_pdf (unsigned int @var{k}, double @var{p}, unsigned int @var{n})
This function computes the probability @math{p(k)} of obtaining @var{k}
from a binomial distribution with parameters @var{p} and @var{n}, using
the formula given above.
@end deftypefun

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

@deftypefun double gsl_cdf_binomial_P (unsigned int @var{k}, double @var{p}, unsigned int @var{n})
@deftypefunx double gsl_cdf_binomial_Q (unsigned int @var{k}, double @var{p}, unsigned int @var{n})
These functions compute the cumulative distribution functions
@math{P(k)}, @math{Q(k)}  for the binomial
distribution with parameters @var{p} and @var{n}.
@end deftypefun


@page
@node The Multinomial Distribution
@section The Multinomial Distribution
@deftypefun void gsl_ran_multinomial (const gsl_rng * @var{r}, size_t @var{K}, unsigned int @var{N}, const double @var{p}[], unsigned int @var{n}[])
@cindex Multinomial distribution 

This function computes a random sample @var{n}[] from the multinomial
distribution formed by @var{N} trials from an underlying distribution
@var{p}[@var{K}]. The distribution function for @var{n}[] is,
@tex
\beforedisplay
$$
P(n_1, n_2,\cdots, n_K) = {{ N!}\over{n_1 ! n_2 ! \cdots n_K !}} \,
  p_1^{n_1} p_2^{n_2} \cdots p_K^{n_K}
$$
\afterdisplay
@end tex
@ifinfo

@example
P(n_1, n_2, ..., n_K) = 
  (N!/(n_1! n_2! ... n_K!)) p_1^n_1 p_2^n_2 ... p_K^n_K
@end example

@end ifinfo
@noindent
where @c{($n_1$, $n_2$, $\ldots$, $n_K$)}
@math{(n_1, n_2, ..., n_K)} 
are nonnegative integers with 
@c{$\sum_{k=1}^{K} n_k =N$} 
@math{sum_@{k=1@}^K n_k = N},
and
@c{$(p_1, p_2, \ldots, p_K)$} 
@math{(p_1, p_2, ..., p_K)}
is a probability distribution with @math{\sum p_i = 1}.  
If the array @var{p}[@var{K}] is not normalized then its entries will be
treated as weights and normalized appropriately.  The arrays @var{n}[]
and @var{p}[] must both be of length @var{K}.

Random variates are generated using the conditional binomial method (see
C.S. David, @cite{The computer generation of multinomial random
variates}, Comp. Stat. Data Anal. 16 (1993) 205--217 for details).
@end deftypefun

@deftypefun double gsl_ran_multinomial_pdf (size_t @var{K}, const double @var{p}[], const unsigned int @var{n}[]) 
This function computes the probability 
@c{$P(n_1, n_2, \ldots, n_K)$}
@math{P(n_1, n_2, ..., n_K)}
of sampling @var{n}[@var{K}] from a multinomial distribution 
with parameters @var{p}[@var{K}], using the formula given above.
@end deftypefun

@deftypefun double gsl_ran_multinomial_lnpdf (size_t @var{K}, const double @var{p}[], const unsigned int @var{n}[]) 
This function returns the logarithm of the probability for the
multinomial distribution @c{$P(n_1, n_2, \ldots, n_K)$}
@math{P(n_1, n_2, ..., n_K)} with parameters @var{p}[@var{K}].
@end deftypefun

@page
@node The Negative Binomial Distribution
@section The Negative Binomial Distribution
@deftypefun {unsigned int} gsl_ran_negative_binomial (const gsl_rng * @var{r}, double @var{p}, double @var{n})
@cindex Negative Binomial distribution, random variates
This function returns a random integer from the negative binomial
distribution, the number of failures occurring before @var{n} successes
in independent trials with probability @var{p} of success.  The
probability distribution for negative binomial variates is,
@tex
\beforedisplay
$$
p(k) = {\Gamma(n + k) \over \Gamma(k+1) \Gamma(n) } p^n (1-p)^k
$$
\afterdisplay
@end tex
@ifinfo

@example
p(k) = @{\Gamma(n + k) \over \Gamma(k+1) \Gamma(n) @} p^n (1-p)^k
@end example

@end ifinfo
@noindent
Note that @math{n} is not required to be an integer.
@end deftypefun

@deftypefun double gsl_ran_negative_binomial_pdf (unsigned int @var{k}, double @var{p}, double  @var{n})
This function computes the probability @math{p(k)} of obtaining @var{k}
from a negative binomial distribution with parameters @var{p} and
@var{n}, using the formula given above.
@end deftypefun

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

@deftypefun double gsl_cdf_negative_binomial_P (unsigned int @var{k}, double @var{p}, double @var{n})
@deftypefunx double gsl_cdf_negative_binomial_Q (unsigned int @var{k}, double @var{p}, double @var{n})
These functions compute the cumulative distribution functions
@math{P(k)}, @math{Q(k)} for the negative binomial distribution with
parameters @var{p} and @var{n}.
@end deftypefun

@page
@node The Pascal Distribution
@section The Pascal Distribution

@deftypefun {unsigned int} gsl_ran_pascal (const gsl_rng * @var{r}, double @var{p}, unsigned int @var{n})
This function returns a random integer from the Pascal distribution.  The
Pascal distribution is simply a negative binomial distribution with an
integer value of @math{n}.
@tex
\beforedisplay
$$
p(k) = {(n + k - 1)! \over k! (n - 1)! } p^n (1-p)^k
$$
\afterdisplay
@end tex
@ifinfo

@example
p(k) = @{(n + k - 1)! \over k! (n - 1)! @} p^n (1-p)^k
@end example

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

@deftypefun double gsl_ran_pascal_pdf (unsigned int @var{k}, double @var{p}, unsigned int @var{n})
This function computes the probability @math{p(k)} of obtaining @var{k}
from a Pascal distribution with parameters @var{p} and
@var{n}, using the formula given above.
@end deftypefun

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

@deftypefun double gsl_cdf_pascal_P (unsigned int @var{k}, double @var{p}, unsigned int @var{n})
@deftypefunx double gsl_cdf_pascal_Q (unsigned int @var{k}, double @var{p}, unsigned int @var{n})
These functions compute the cumulative distribution functions
@math{P(k)}, @math{Q(k)} for the Pascal distribution with
parameters @var{p} and @var{n}.
@end deftypefun

@page
@node The Geometric Distribution
@section The Geometric Distribution
@deftypefun {unsigned int} gsl_ran_geometric (const gsl_rng * @var{r}, double @var{p})
@cindex Geometric random variates
This function returns a random integer from the geometric distribution,
the number of independent trials with probability @var{p} until the
first success.  The probability distribution for geometric variates
is,
@tex
\beforedisplay
$$
p(k) =  p (1-p)^{k-1}
$$
\afterdisplay
@end tex
@ifinfo

@example
p(k) =  p (1-p)^(k-1)
@end example

@end ifinfo
@noindent
for @c{$k \ge 1$}
@math{k >= 1}.  Note that the distribution begins with @math{k=1} with this
definition.  There is another convention in which the exponent @math{k-1} 
is replaced by @math{k}.
@end deftypefun

@deftypefun double gsl_ran_geometric_pdf (unsigned int @var{k}, double @var{p})
This function computes the probability @math{p(k)} of obtaining @var{k}
from a geometric distribution with probability parameter @var{p}, using
the formula given above.
@end deftypefun

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

@deftypefun double gsl_cdf_geometric_P (unsigned int @var{k}, double @var{p})
@deftypefunx double gsl_cdf_geometric_Q (unsigned int @var{k}, double @var{p})
These functions compute the cumulative distribution functions
@math{P(k)}, @math{Q(k)} for the geometric distribution with parameter
@var{p}.
@end deftypefun


@page
@node The Hypergeometric Distribution
@section The Hypergeometric Distribution
@cindex hypergeometric random variates
@deftypefun {unsigned int} gsl_ran_hypergeometric (const gsl_rng * @var{r}, unsigned int @var{n1}, unsigned int @var{n2}, unsigned int @var{t})
@cindex Geometric random variates
This function returns a random integer from the hypergeometric
distribution.  The probability distribution for hypergeometric
random variates is,
@tex
\beforedisplay
$$
p(k) =  C(n_1, k) C(n_2, t - k) / C(n_1 + n_2, t)
$$
\afterdisplay
@end tex
@ifinfo

@example
p(k) =  C(n_1, k) C(n_2, t - k) / C(n_1 + n_2, t)
@end example

@end ifinfo
@noindent
where @math{C(a,b) = a!/(b!(a-b)!)} and 
@c{$t \leq n_1 + n_2$}
@math{t <= n_1 + n_2}.  The domain of @math{k} is 
@c{$\hbox{max}(0,t-n_2), \ldots, \hbox{min}(t,n_1)$} 
@math{max(0,t-n_2), ..., min(t,n_1)}.

If a population contains @math{n_1} elements of ``type 1'' and
@math{n_2} elements of ``type 2'' then the hypergeometric
distribution gives the probability of obtaining @math{k} elements of
``type 1'' in @math{t} samples from the population without
replacement.
@end deftypefun

@deftypefun double gsl_ran_hypergeometric_pdf (unsigned int @var{k}, unsigned int @var{n1}, unsigned int @var{n2}, unsigned int @var{t})
This function computes the probability @math{p(k)} of obtaining @var{k}
from a hypergeometric distribution with parameters @var{n1}, @var{n2},
@var{t}, using the formula given above.
@end deftypefun

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

@deftypefun double gsl_cdf_hypergeometric_P (unsigned int @var{k}, unsigned int @var{n1}, unsigned int @var{n2}, unsigned int @var{t})
@deftypefunx double gsl_cdf_hypergeometric_Q (unsigned int @var{k}, unsigned int @var{n1}, unsigned int @var{n2}, unsigned int @var{t})
These functions compute the cumulative distribution functions
@math{P(k)}, @math{Q(k)} for the hypergeometric distribution with
parameters @var{n1}, @var{n2} and @var{t}.
@end deftypefun


@page
@node The Logarithmic Distribution
@section The Logarithmic Distribution
@deftypefun {unsigned int} gsl_ran_logarithmic (const gsl_rng * @var{r}, double @var{p})
@cindex Logarithmic random variates
This function returns a random integer from the logarithmic
distribution.  The probability distribution for logarithmic random variates
is,
@tex
\beforedisplay
$$
p(k) = {-1 \over \log(1-p)} {\left( p^k \over k \right)}
$$
\afterdisplay
@end tex
@ifinfo

@example
p(k) = @{-1 \over \log(1-p)@} @{(p^k \over k)@}
@end example

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

@deftypefun double gsl_ran_logarithmic_pdf (unsigned int @var{k}, double @var{p})
This function computes the probability @math{p(k)} of obtaining @var{k}
from a logarithmic distribution with probability parameter @var{p},
using the formula given above.
@end deftypefun

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

@page
@node Shuffling and Sampling
@section Shuffling and Sampling

The following functions allow the shuffling and sampling of a set of
objects.  The algorithms rely on a random number generator as a source of
randomness and a poor quality generator can lead to correlations in the
output.  In particular it is important to avoid generators with a short
period.  For more information see Knuth, v2, 3rd ed, Section 3.4.2,
``Random Sampling and Shuffling''.

@deftypefun void gsl_ran_shuffle (const gsl_rng * @var{r}, void * @var{base}, size_t @var{n}, size_t @var{size})

This function randomly shuffles the order of @var{n} objects, each of
size @var{size}, stored in the array @var{base}[0..@var{n}-1].  The
output of the random number generator @var{r} is used to produce the
permutation.  The algorithm generates all possible @math{n!}
permutations with equal probability, assuming a perfect source of random
numbers.

The following code shows how to shuffle the numbers from 0 to 51,

@example
int a[52];

for (i = 0; i < 52; i++)
  @{
    a[i] = i;
  @}

gsl_ran_shuffle (r, a, 52, sizeof (int));
@end example

@end deftypefun

@deftypefun int gsl_ran_choose (const gsl_rng * @var{r}, void * @var{dest}, size_t @var{k}, void * @var{src}, size_t @var{n}, size_t @var{size})
This function fills the array @var{dest}[k] with @var{k} objects taken
randomly from the @var{n} elements of the array
@var{src}[0..@var{n}-1].  The objects are each of size @var{size}.  The
output of the random number generator @var{r} is used to make the
selection.  The algorithm ensures all possible samples are equally
likely, assuming a perfect source of randomness.

The objects are sampled @emph{without} replacement, thus each object can
only appear once in @var{dest}[k].  It is required that @var{k} be less
than or equal to @code{n}.  The objects in @var{dest} will be in the
same relative order as those in @var{src}.  You will need to call
@code{gsl_ran_shuffle(r, dest, n, size)} if you want to randomize the
order.

The following code shows how to select a random sample of three unique
numbers from the set 0 to 99,

@example
double a[3], b[100];

for (i = 0; i < 100; i++)
  @{
    b[i] = (double) i;
  @}

gsl_ran_choose (r, a, 3, b, 100, sizeof (double));
@end example

@end deftypefun

@deftypefun void gsl_ran_sample (const gsl_rng * @var{r}, void * @var{dest}, size_t @var{k}, void * @var{src}, size_t @var{n}, size_t @var{size})
This function is like @code{gsl_ran_choose} but samples @var{k} items
from the original array of @var{n} items @var{src} with replacement, so
the same object can appear more than once in the output sequence
@var{dest}.  There is no requirement that @var{k} be less than @var{n}
in this case.
@end deftypefun


@node Random Number Distribution Examples
@section Examples

The following program demonstrates the use of a random number generator
to produce variates from a distribution.  It prints 10 samples from the
Poisson distribution with a mean of 3.

@example
@verbatiminclude examples/randpoisson.c
@end example

@noindent
If the library and header files are installed under @file{/usr/local}
(the default location) then the program can be compiled with these
options,

@example
$ gcc -Wall demo.c -lgsl -lgslcblas -lm
@end example

@noindent
Here is the output of the program,

@example
$ ./a.out 
@verbatiminclude examples/randpoisson.out
@end example

@noindent
The variates depend on the seed used by the generator.  The seed for the
default generator type @code{gsl_rng_default} can be changed with the
@code{GSL_RNG_SEED} environment variable to produce a different stream
of variates,

@example
$ GSL_RNG_SEED=123 ./a.out 
@verbatiminclude examples/randpoisson.2.out
@end example

@noindent
The following program generates a random walk in two dimensions.

@example
@verbatiminclude examples/randwalk.c
@end example

@noindent
Here is the output from the program, three 10-step random walks from the origin,

@tex
\centerline{\input random-walk.tex}
@end tex

The following program computes the upper and lower cumulative
distribution functions for the standard normal distribution at
@math{x=2}.

@example
@verbatiminclude examples/cdf.c
@end example

@noindent
Here is the output of the program,

@example
@verbatiminclude examples/cdf.out
@end example

@node Random Number Distribution References and Further Reading
@section References and Further Reading

For an encyclopaedic coverage of the subject readers are advised to
consult the book @cite{Non-Uniform Random Variate Generation} by Luc
Devroye.  It covers every imaginable distribution and provides hundreds
of algorithms.

@itemize @asis
@item
Luc Devroye, @cite{Non-Uniform Random Variate Generation},
Springer-Verlag, ISBN 0-387-96305-7.
@end itemize

@noindent
The subject of random variate generation is also reviewed by Knuth, who
describes algorithms for all the major distributions.

@itemize @asis
@item
Donald E. Knuth, @cite{The Art of Computer Programming: Seminumerical
Algorithms} (Vol 2, 3rd Ed, 1997), Addison-Wesley, ISBN 0201896842.
@end itemize

@noindent
The Particle Data Group provides a short review of techniques for
generating distributions of random numbers in the ``Monte Carlo''
section of its Annual Review of Particle Physics.

@itemize @asis
@item
@cite{Review of Particle Properties}
R.M. Barnett et al., Physical Review D54, 1 (1996)
@uref{http://pdg.lbl.gov/}.
@end itemize

@noindent
The Review of Particle Physics is available online in postscript and pdf
format.

@noindent
An overview of methods used to compute cumulative distribution functions
can be found in @cite{Statistical Computing} by W.J. Kennedy and
J.E. Gentle. Another general reference is @cite{Elements of Statistical
Computing} by R.A. Thisted.

@itemize @asis
@item
William E. Kennedy and James E. Gentle, @cite{Statistical Computing} (1980),
Marcel Dekker, ISBN 0-8247-6898-1.
@end itemize

@itemize @asis
@item
Ronald A. Thisted, @cite{Elements of Statistical Computing} (1988), 
Chapman & Hall, ISBN 0-412-01371-1.
@end itemize

@noindent
The cumulative distribution functions for the Gaussian distribution
are based on the following papers,

@itemize @asis
@item
@cite{Rational Chebyshev Approximations Using Linear Equations},
W.J. Cody, W. Fraser, J.F. Hart. Numerische Mathematik 12, 242--251 (1968).
@end itemize

@itemize @asis
@item
@cite{Rational Chebyshev Approximations for the Error Function},
W.J. Cody. Mathematics of Computation 23, n107, 631--637 (July 1969).
@end itemize