diff options
Diffstat (limited to 'gsl-1.9/doc/randist.texi')
| -rw-r--r-- | gsl-1.9/doc/randist.texi | 2296 |
1 files changed, 2296 insertions, 0 deletions
diff --git a/gsl-1.9/doc/randist.texi b/gsl-1.9/doc/randist.texi new file mode 100644 index 0000000..e193a09 --- /dev/null +++ b/gsl-1.9/doc/randist.texi @@ -0,0 +1,2296 @@ +@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 |