randist.texi

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

@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})