1 @cindex random number distributions
2 @cindex cumulative distribution functions (CDFs)
3 @cindex CDFs, cumulative distribution functions
4 @cindex inverse cumulative distribution functions
5 @cindex quantile functions
6 This chapter describes functions for generating random variates and
7 computing their probability distributions. Samples from the
8 distributions described in this chapter can be obtained using any of the
9 random number generators in the library as an underlying source of
12 In the simplest cases a non-uniform distribution can be obtained
13 analytically from the uniform distribution of a random number generator
14 by applying an appropriate transformation. This method uses one call to
15 the random number generator. More complicated distributions are created
16 by the @dfn{acceptance-rejection} method, which compares the desired
17 distribution against a distribution which is similar and known
18 analytically. This usually requires several samples from the generator.
20 The library also provides cumulative distribution functions and inverse
21 cumulative distribution functions, sometimes referred to as quantile
22 functions. The cumulative distribution functions and their inverses are
23 computed separately for the upper and lower tails of the distribution,
24 allowing full accuracy to be retained for small results.
26 The functions for random variates and probability density functions
27 described in this section are declared in @file{gsl_randist.h}. The
28 corresponding cumulative distribution functions are declared in
31 Note that the discrete random variate functions always
32 return a value of type @code{unsigned int}, and on most platforms this
33 has a maximum value of @c{$2^{32}-1 \approx 4.29\times10^9$}
34 @math{2^32-1 ~=~ 4.29e9}. They should only be called with
35 a safe range of parameters (where there is a negligible probability of
36 a variate exceeding this limit) to prevent incorrect results due to
40 * Random Number Distribution Introduction::
41 * The Gaussian Distribution::
42 * The Gaussian Tail Distribution::
43 * The Bivariate Gaussian Distribution::
44 * The Exponential Distribution::
45 * The Laplace Distribution::
46 * The Exponential Power Distribution::
47 * The Cauchy Distribution::
48 * The Rayleigh Distribution::
49 * The Rayleigh Tail Distribution::
50 * The Landau Distribution::
51 * The Levy alpha-Stable Distributions::
52 * The Levy skew alpha-Stable Distribution::
53 * The Gamma Distribution::
54 * The Flat (Uniform) Distribution::
55 * The Lognormal Distribution::
56 * The Chi-squared Distribution::
57 * The F-distribution::
58 * The t-distribution::
59 * The Beta Distribution::
60 * The Logistic Distribution::
61 * The Pareto Distribution::
62 * Spherical Vector Distributions::
63 * The Weibull Distribution::
64 * The Type-1 Gumbel Distribution::
65 * The Type-2 Gumbel Distribution::
66 * The Dirichlet Distribution::
67 * General Discrete Distributions::
68 * The Poisson Distribution::
69 * The Bernoulli Distribution::
70 * The Binomial Distribution::
71 * The Multinomial Distribution::
72 * The Negative Binomial Distribution::
73 * The Pascal Distribution::
74 * The Geometric Distribution::
75 * The Hypergeometric Distribution::
76 * The Logarithmic Distribution::
77 * Shuffling and Sampling::
78 * Random Number Distribution Examples::
79 * Random Number Distribution References and Further Reading::
82 @node Random Number Distribution Introduction
85 Continuous random number distributions are defined by a probability
86 density function, @math{p(x)}, such that the probability of @math{x}
87 occurring in the infinitesimal range @math{x} to @math{x+dx} is @c{$p\,dx$}
90 The cumulative distribution function for the lower tail @math{P(x)} is
91 defined by the integral,
95 P(x) = \int_{-\infty}^{x} dx' p(x')
102 P(x) = \int_@{-\infty@}^@{x@} dx' p(x')
107 and gives the probability of a variate taking a value less than @math{x}.
109 The cumulative distribution function for the upper tail @math{Q(x)} is
110 defined by the integral,
114 Q(x) = \int_{x}^{+\infty} dx' p(x')
121 Q(x) = \int_@{x@}^@{+\infty@} dx' p(x')
126 and gives the probability of a variate taking a value greater than @math{x}.
128 The upper and lower cumulative distribution functions are related by
129 @math{P(x) + Q(x) = 1} and satisfy @c{$0 \le P(x) \le 1$}
130 @math{0 <= P(x) <= 1}, @c{$0 \le Q(x) \le 1$}
131 @math{0 <= Q(x) <= 1}.
133 The inverse cumulative distributions, @c{$x=P^{-1}(P)$}
134 @math{x=P^@{-1@}(P)} and @c{$x=Q^{-1}(Q)$}
135 @math{x=Q^@{-1@}(Q)} give the values of @math{x}
136 which correspond to a specific value of @math{P} or @math{Q}.
137 They can be used to find confidence limits from probability values.
139 For discrete distributions the probability of sampling the integer
140 value @math{k} is given by @math{p(k)}, where @math{\sum_k p(k) = 1}.
141 The cumulative distribution for the lower tail @math{P(k)} of a
142 discrete distribution is defined as,
146 P(k) = \sum_{i \le k} p(i)
153 P(k) = \sum_@{i <= k@} p(i)
158 where the sum is over the allowed range of the distribution less than
159 or equal to @math{k}.
161 The cumulative distribution for the upper tail of a discrete
162 distribution @math{Q(k)} is defined as
166 Q(k) = \sum_{i > k} p(i)
173 Q(k) = \sum_@{i > k@} p(i)
178 giving the sum of probabilities for all values greater than @math{k}.
179 These two definitions satisfy the identity @math{P(k)+Q(k)=1}.
181 If the range of the distribution is 1 to @math{n} inclusive then
182 @math{P(n)=1}, @math{Q(n)=0} while @math{P(1) = p(1)},
186 @node The Gaussian Distribution
187 @section The Gaussian Distribution
188 @deftypefun double gsl_ran_gaussian (const gsl_rng * @var{r}, double @var{sigma})
189 @cindex Gaussian distribution
190 This function returns a Gaussian random variate, with mean zero and
191 standard deviation @var{sigma}. The probability distribution for
192 Gaussian random variates is,
196 p(x) dx = {1 \over \sqrt{2 \pi \sigma^2}} \exp (-x^2 / 2\sigma^2) dx
203 p(x) dx = @{1 \over \sqrt@{2 \pi \sigma^2@}@} \exp (-x^2 / 2\sigma^2) dx
208 for @math{x} in the range @math{-\infty} to @math{+\infty}. Use the
209 transformation @math{z = \mu + x} on the numbers returned by
210 @code{gsl_ran_gaussian} to obtain a Gaussian distribution with mean
211 @math{\mu}. This function uses the Box-Mueller algorithm which requires two
212 calls to the random number generator @var{r}.
215 @deftypefun double gsl_ran_gaussian_pdf (double @var{x}, double @var{sigma})
216 This function computes the probability density @math{p(x)} at @var{x}
217 for a Gaussian distribution with standard deviation @var{sigma}, using
218 the formula given above.
223 \centerline{\input rand-gaussian.tex}
226 @deftypefun double gsl_ran_gaussian_ziggurat (const gsl_rng * @var{r}, double @var{sigma})
227 @deftypefunx double gsl_ran_gaussian_ratio_method (const gsl_rng * @var{r}, double @var{sigma})
228 This function computes a Gaussian random variate using the alternative
229 Marsaglia-Tsang ziggurat and Kinderman-Monahan-Leva ratio methods. The
230 Ziggurat algorithm is the fastest available algorithm in most cases.
233 @deftypefun double gsl_ran_ugaussian (const gsl_rng * @var{r})
234 @deftypefunx double gsl_ran_ugaussian_pdf (double @var{x})
235 @deftypefunx double gsl_ran_ugaussian_ratio_method (const gsl_rng * @var{r})
236 These functions compute results for the unit Gaussian distribution. They
237 are equivalent to the functions above with a standard deviation of one,
241 @deftypefun double gsl_cdf_gaussian_P (double @var{x}, double @var{sigma})
242 @deftypefunx double gsl_cdf_gaussian_Q (double @var{x}, double @var{sigma})
243 @deftypefunx double gsl_cdf_gaussian_Pinv (double @var{P}, double @var{sigma})
244 @deftypefunx double gsl_cdf_gaussian_Qinv (double @var{Q}, double @var{sigma})
245 These functions compute the cumulative distribution functions
246 @math{P(x)}, @math{Q(x)} and their inverses for the Gaussian
247 distribution with standard deviation @var{sigma}.
250 @deftypefun double gsl_cdf_ugaussian_P (double @var{x})
251 @deftypefunx double gsl_cdf_ugaussian_Q (double @var{x})
252 @deftypefunx double gsl_cdf_ugaussian_Pinv (double @var{P})
253 @deftypefunx double gsl_cdf_ugaussian_Qinv (double @var{Q})
254 These functions compute the cumulative distribution functions
255 @math{P(x)}, @math{Q(x)} and their inverses for the unit Gaussian
260 @node The Gaussian Tail Distribution
261 @section The Gaussian Tail Distribution
262 @deftypefun double gsl_ran_gaussian_tail (const gsl_rng * @var{r}, double @var{a}, double @var{sigma})
263 @cindex Gaussian Tail distribution
264 This function provides random variates from the upper tail of a Gaussian
265 distribution with standard deviation @var{sigma}. The values returned
266 are larger than the lower limit @var{a}, which must be positive. The
267 method is based on Marsaglia's famous rectangle-wedge-tail algorithm (Ann.
268 Math. Stat. 32, 894--899 (1961)), with this aspect explained in Knuth, v2,
269 3rd ed, p139,586 (exercise 11).
271 The probability distribution for Gaussian tail random variates is,
275 p(x) dx = {1 \over N(a;\sigma) \sqrt{2 \pi \sigma^2}} \exp (- x^2 / 2\sigma^2) dx
282 p(x) dx = @{1 \over N(a;\sigma) \sqrt@{2 \pi \sigma^2@}@} \exp (- x^2/(2 \sigma^2)) dx
287 for @math{x > a} where @math{N(a;\sigma)} is the normalization constant,
291 N(a;\sigma) = {1 \over 2} \hbox{erfc}\left({a \over \sqrt{2 \sigma^2}}\right).
298 N(a;\sigma) = (1/2) erfc(a / sqrt(2 sigma^2)).
304 @deftypefun double gsl_ran_gaussian_tail_pdf (double @var{x}, double @var{a}, double @var{sigma})
305 This function computes the probability density @math{p(x)} at @var{x}
306 for a Gaussian tail distribution with standard deviation @var{sigma} and
307 lower limit @var{a}, using the formula given above.
312 \centerline{\input rand-gaussian-tail.tex}
315 @deftypefun double gsl_ran_ugaussian_tail (const gsl_rng * @var{r}, double @var{a})
316 @deftypefunx double gsl_ran_ugaussian_tail_pdf (double @var{x}, double @var{a})
317 These functions compute results for the tail of a unit Gaussian
318 distribution. They are equivalent to the functions above with a standard
319 deviation of one, @var{sigma} = 1.
324 @node The Bivariate Gaussian Distribution
325 @section The Bivariate Gaussian Distribution
327 @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})
328 @cindex Bivariate Gaussian distribution
329 @cindex two dimensional Gaussian distribution
330 @cindex Gaussian distribution, bivariate
331 This function generates a pair of correlated Gaussian variates, with
332 mean zero, correlation coefficient @var{rho} and standard deviations
333 @var{sigma_x} and @var{sigma_y} in the @math{x} and @math{y} directions.
334 The probability distribution for bivariate Gaussian random variates is,
338 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
345 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
350 for @math{x,y} in the range @math{-\infty} to @math{+\infty}. The
351 correlation coefficient @var{rho} should lie between @math{1} and
355 @deftypefun double gsl_ran_bivariate_gaussian_pdf (double @var{x}, double @var{y}, double @var{sigma_x}, double @var{sigma_y}, double @var{rho})
356 This function computes the probability density @math{p(x,y)} at
357 (@var{x},@var{y}) for a bivariate Gaussian distribution with standard
358 deviations @var{sigma_x}, @var{sigma_y} and correlation coefficient
359 @var{rho}, using the formula given above.
364 \centerline{\input rand-bivariate-gaussian.tex}
368 @node The Exponential Distribution
369 @section The Exponential Distribution
370 @deftypefun double gsl_ran_exponential (const gsl_rng * @var{r}, double @var{mu})
371 @cindex Exponential distribution
372 This function returns a random variate from the exponential distribution
373 with mean @var{mu}. The distribution is,
377 p(x) dx = {1 \over \mu} \exp(-x/\mu) dx
384 p(x) dx = @{1 \over \mu@} \exp(-x/\mu) dx
393 @deftypefun double gsl_ran_exponential_pdf (double @var{x}, double @var{mu})
394 This function computes the probability density @math{p(x)} at @var{x}
395 for an exponential distribution with mean @var{mu}, using the formula
401 \centerline{\input rand-exponential.tex}
404 @deftypefun double gsl_cdf_exponential_P (double @var{x}, double @var{mu})
405 @deftypefunx double gsl_cdf_exponential_Q (double @var{x}, double @var{mu})
406 @deftypefunx double gsl_cdf_exponential_Pinv (double @var{P}, double @var{mu})
407 @deftypefunx double gsl_cdf_exponential_Qinv (double @var{Q}, double @var{mu})
408 These functions compute the cumulative distribution functions
409 @math{P(x)}, @math{Q(x)} and their inverses for the exponential
410 distribution with mean @var{mu}.
414 @node The Laplace Distribution
415 @section The Laplace Distribution
416 @deftypefun double gsl_ran_laplace (const gsl_rng * @var{r}, double @var{a})
417 @cindex two-sided exponential distribution
418 @cindex Laplace distribution
419 This function returns a random variate from the Laplace distribution
420 with width @var{a}. The distribution is,
424 p(x) dx = {1 \over 2 a} \exp(-|x/a|) dx
431 p(x) dx = @{1 \over 2 a@} \exp(-|x/a|) dx
436 for @math{-\infty < x < \infty}.
439 @deftypefun double gsl_ran_laplace_pdf (double @var{x}, double @var{a})
440 This function computes the probability density @math{p(x)} at @var{x}
441 for a Laplace distribution with width @var{a}, using the formula
447 \centerline{\input rand-laplace.tex}
450 @deftypefun double gsl_cdf_laplace_P (double @var{x}, double @var{a})
451 @deftypefunx double gsl_cdf_laplace_Q (double @var{x}, double @var{a})
452 @deftypefunx double gsl_cdf_laplace_Pinv (double @var{P}, double @var{a})
453 @deftypefunx double gsl_cdf_laplace_Qinv (double @var{Q}, double @var{a})
454 These functions compute the cumulative distribution functions
455 @math{P(x)}, @math{Q(x)} and their inverses for the Laplace
456 distribution with width @var{a}.
461 @node The Exponential Power Distribution
462 @section The Exponential Power Distribution
463 @deftypefun double gsl_ran_exppow (const gsl_rng * @var{r}, double @var{a}, double @var{b})
464 @cindex Exponential power distribution
465 This function returns a random variate from the exponential power distribution
466 with scale parameter @var{a} and exponent @var{b}. The distribution is,
470 p(x) dx = {1 \over 2 a \Gamma(1+1/b)} \exp(-|x/a|^b) dx
477 p(x) dx = @{1 \over 2 a \Gamma(1+1/b)@} \exp(-|x/a|^b) dx
483 @math{x >= 0}. For @math{b = 1} this reduces to the Laplace
484 distribution. For @math{b = 2} it has the same form as a gaussian
485 distribution, but with @c{$a = \sqrt{2} \sigma$}
486 @math{a = \sqrt@{2@} \sigma}.
489 @deftypefun double gsl_ran_exppow_pdf (double @var{x}, double @var{a}, double @var{b})
490 This function computes the probability density @math{p(x)} at @var{x}
491 for an exponential power distribution with scale parameter @var{a}
492 and exponent @var{b}, using the formula given above.
497 \centerline{\input rand-exppow.tex}
500 @deftypefun double gsl_cdf_exppow_P (double @var{x}, double @var{a}, double @var{b})
501 @deftypefunx double gsl_cdf_exppow_Q (double @var{x}, double @var{a}, double @var{b})
502 These functions compute the cumulative distribution functions
503 @math{P(x)}, @math{Q(x)} for the exponential power distribution with
504 parameters @var{a} and @var{b}.
509 @node The Cauchy Distribution
510 @section The Cauchy Distribution
511 @deftypefun double gsl_ran_cauchy (const gsl_rng * @var{r}, double @var{a})
512 @cindex Cauchy distribution
513 This function returns a random variate from the Cauchy distribution with
514 scale parameter @var{a}. The probability distribution for Cauchy
519 p(x) dx = {1 \over a\pi (1 + (x/a)^2) } dx
526 p(x) dx = @{1 \over a\pi (1 + (x/a)^2) @} dx
531 for @math{x} in the range @math{-\infty} to @math{+\infty}. The Cauchy
532 distribution is also known as the Lorentz distribution.
535 @deftypefun double gsl_ran_cauchy_pdf (double @var{x}, double @var{a})
536 This function computes the probability density @math{p(x)} at @var{x}
537 for a Cauchy distribution with scale parameter @var{a}, using the formula
543 \centerline{\input rand-cauchy.tex}
546 @deftypefun double gsl_cdf_cauchy_P (double @var{x}, double @var{a})
547 @deftypefunx double gsl_cdf_cauchy_Q (double @var{x}, double @var{a})
548 @deftypefunx double gsl_cdf_cauchy_Pinv (double @var{P}, double @var{a})
549 @deftypefunx double gsl_cdf_cauchy_Qinv (double @var{Q}, double @var{a})
550 These functions compute the cumulative distribution functions
551 @math{P(x)}, @math{Q(x)} and their inverses for the Cauchy
552 distribution with scale parameter @var{a}.
557 @node The Rayleigh Distribution
558 @section The Rayleigh Distribution
559 @deftypefun double gsl_ran_rayleigh (const gsl_rng * @var{r}, double @var{sigma})
560 @cindex Rayleigh distribution
561 This function returns a random variate from the Rayleigh distribution with
562 scale parameter @var{sigma}. The distribution is,
566 p(x) dx = {x \over \sigma^2} \exp(- x^2/(2 \sigma^2)) dx
573 p(x) dx = @{x \over \sigma^2@} \exp(- x^2/(2 \sigma^2)) dx
581 @deftypefun double gsl_ran_rayleigh_pdf (double @var{x}, double @var{sigma})
582 This function computes the probability density @math{p(x)} at @var{x}
583 for a Rayleigh distribution with scale parameter @var{sigma}, using the
589 \centerline{\input rand-rayleigh.tex}
592 @deftypefun double gsl_cdf_rayleigh_P (double @var{x}, double @var{sigma})
593 @deftypefunx double gsl_cdf_rayleigh_Q (double @var{x}, double @var{sigma})
594 @deftypefunx double gsl_cdf_rayleigh_Pinv (double @var{P}, double @var{sigma})
595 @deftypefunx double gsl_cdf_rayleigh_Qinv (double @var{Q}, double @var{sigma})
596 These functions compute the cumulative distribution functions
597 @math{P(x)}, @math{Q(x)} and their inverses for the Rayleigh
598 distribution with scale parameter @var{sigma}.
603 @node The Rayleigh Tail Distribution
604 @section The Rayleigh Tail Distribution
605 @deftypefun double gsl_ran_rayleigh_tail (const gsl_rng * @var{r}, double @var{a}, double @var{sigma})
606 @cindex Rayleigh Tail distribution
607 This function returns a random variate from the tail of the Rayleigh
608 distribution with scale parameter @var{sigma} and a lower limit of
609 @var{a}. The distribution is,
613 p(x) dx = {x \over \sigma^2} \exp ((a^2 - x^2) /(2 \sigma^2)) dx
620 p(x) dx = @{x \over \sigma^2@} \exp ((a^2 - x^2) /(2 \sigma^2)) dx
628 @deftypefun double gsl_ran_rayleigh_tail_pdf (double @var{x}, double @var{a}, double @var{sigma})
629 This function computes the probability density @math{p(x)} at @var{x}
630 for a Rayleigh tail distribution with scale parameter @var{sigma} and
631 lower limit @var{a}, using the formula given above.
636 \centerline{\input rand-rayleigh-tail.tex}
640 @node The Landau Distribution
641 @section The Landau Distribution
642 @deftypefun double gsl_ran_landau (const gsl_rng * @var{r})
643 @cindex Landau distribution
644 This function returns a random variate from the Landau distribution. The
645 probability distribution for Landau random variates is defined
646 analytically by the complex integral,
651 {1 \over {2 \pi i}} \int_{c-i\infty}^{c+i\infty} ds\, \exp(s \log(s) + x s)
658 p(x) = (1/(2 \pi i)) \int_@{c-i\infty@}^@{c+i\infty@} ds exp(s log(s) + x s)
661 For numerical purposes it is more convenient to use the following
662 equivalent form of the integral,
666 p(x) = (1/\pi) \int_0^\infty dt \exp(-t \log(t) - x t) \sin(\pi t).
673 p(x) = (1/\pi) \int_0^\infty dt \exp(-t \log(t) - x t) \sin(\pi t).
678 @deftypefun double gsl_ran_landau_pdf (double @var{x})
679 This function computes the probability density @math{p(x)} at @var{x}
680 for the Landau distribution using an approximation to the formula given
686 \centerline{\input rand-landau.tex}
690 @node The Levy alpha-Stable Distributions
691 @section The Levy alpha-Stable Distributions
692 @deftypefun double gsl_ran_levy (const gsl_rng * @var{r}, double @var{c}, double @var{alpha})
693 @cindex Levy distribution
694 This function returns a random variate from the Levy symmetric stable
695 distribution with scale @var{c} and exponent @var{alpha}. The symmetric
696 stable probability distribution is defined by a fourier transform,
700 p(x) = {1 \over 2 \pi} \int_{-\infty}^{+\infty} dt \exp(-it x - |c t|^\alpha)
707 p(x) = @{1 \over 2 \pi@} \int_@{-\infty@}^@{+\infty@} dt \exp(-it x - |c t|^alpha)
712 There is no explicit solution for the form of @math{p(x)} and the
713 library does not define a corresponding @code{pdf} function. For
714 @math{\alpha = 1} the distribution reduces to the Cauchy distribution. For
715 @math{\alpha = 2} it is a Gaussian distribution with @c{$\sigma = \sqrt{2} c$}
716 @math{\sigma = \sqrt@{2@} c}. For @math{\alpha < 1} the tails of the
717 distribution become extremely wide.
719 The algorithm only works for @c{$0 < \alpha \le 2$}
720 @math{0 < alpha <= 2}.
725 \centerline{\input rand-levy.tex}
729 @node The Levy skew alpha-Stable Distribution
730 @section The Levy skew alpha-Stable Distribution
732 @deftypefun double gsl_ran_levy_skew (const gsl_rng * @var{r}, double @var{c}, double @var{alpha}, double @var{beta})
733 @cindex Levy distribution, skew
734 @cindex Skew Levy distribution
735 This function returns a random variate from the Levy skew stable
736 distribution with scale @var{c}, exponent @var{alpha} and skewness
737 parameter @var{beta}. The skewness parameter must lie in the range
738 @math{[-1,1]}. The Levy skew stable probability distribution is defined
739 by a fourier transform,
743 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)))
750 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)))
755 When @math{\alpha = 1} the term @math{\tan(\pi \alpha/2)} is replaced by
756 @math{-(2/\pi)\log|t|}. There is no explicit solution for the form of
757 @math{p(x)} and the library does not define a corresponding @code{pdf}
758 function. For @math{\alpha = 2} the distribution reduces to a Gaussian
759 distribution with @c{$\sigma = \sqrt{2} c$}
760 @math{\sigma = \sqrt@{2@} c} and the skewness parameter has no effect.
761 For @math{\alpha < 1} the tails of the distribution become extremely
762 wide. The symmetric distribution corresponds to @math{\beta =
765 The algorithm only works for @c{$0 < \alpha \le 2$}
766 @math{0 < alpha <= 2}.
769 The Levy alpha-stable distributions have the property that if @math{N}
770 alpha-stable variates are drawn from the distribution @math{p(c, \alpha,
771 \beta)} then the sum @math{Y = X_1 + X_2 + \dots + X_N} will also be
772 distributed as an alpha-stable variate,
773 @c{$p(N^{1/\alpha} c, \alpha, \beta)$}
774 @math{p(N^(1/\alpha) c, \alpha, \beta)}.
776 @comment PDF not available because there is no analytic expression for it
778 @comment @deftypefun double gsl_ran_levy_pdf (double @var{x}, double @var{mu})
779 @comment This function computes the probability density @math{p(x)} at @var{x}
780 @comment for a symmetric Levy distribution with scale parameter @var{mu} and
781 @comment exponent @var{a}, using the formula given above.
782 @comment @end deftypefun
786 \centerline{\input rand-levyskew.tex}
790 @node The Gamma Distribution
791 @section The Gamma Distribution
792 @deftypefun double gsl_ran_gamma (const gsl_rng * @var{r}, double @var{a}, double @var{b})
793 @cindex Gamma distribution
794 This function returns a random variate from the gamma
795 distribution. The distribution function is,
799 p(x) dx = {1 \over \Gamma(a) b^a} x^{a-1} e^{-x/b} dx
806 p(x) dx = @{1 \over \Gamma(a) b^a@} x^@{a-1@} e^@{-x/b@} dx
812 @comment If @xmath{X} and @xmath{Y} are independent gamma-distributed random
813 @comment variables of order @xmath{a} and @xmath{b}, then @xmath{X+Y} has a gamma
814 @comment distribution of order @xmath{a+b}.
816 @cindex Erlang distribution
817 The gamma distribution with an integer parameter @var{a} is known as the Erlang distribution.
819 The variates are computed using the Marsaglia-Tsang fast gamma method.
820 This function for this method was previously called
821 @code{gsl_ran_gamma_mt} and can still be accessed using this name.
824 @deftypefun double gsl_ran_gamma_knuth (const gsl_rng * @var{r}, double @var{a}, double @var{b})
825 This function returns a gamma variate using the algorithms from Knuth (vol 2).
828 @deftypefun double gsl_ran_gamma_pdf (double @var{x}, double @var{a}, double @var{b})
829 This function computes the probability density @math{p(x)} at @var{x}
830 for a gamma distribution with parameters @var{a} and @var{b}, using the
836 \centerline{\input rand-gamma.tex}
839 @deftypefun double gsl_cdf_gamma_P (double @var{x}, double @var{a}, double @var{b})
840 @deftypefunx double gsl_cdf_gamma_Q (double @var{x}, double @var{a}, double @var{b})
841 @deftypefunx double gsl_cdf_gamma_Pinv (double @var{P}, double @var{a}, double @var{b})
842 @deftypefunx double gsl_cdf_gamma_Qinv (double @var{Q}, double @var{a}, double @var{b})
843 These functions compute the cumulative distribution functions
844 @math{P(x)}, @math{Q(x)} and their inverses for the gamma
845 distribution with parameters @var{a} and @var{b}.
849 @node The Flat (Uniform) Distribution
850 @section The Flat (Uniform) Distribution
851 @deftypefun double gsl_ran_flat (const gsl_rng * @var{r}, double @var{a}, double @var{b})
852 @cindex flat distribution
853 @cindex uniform distribution
854 This function returns a random variate from the flat (uniform)
855 distribution from @var{a} to @var{b}. The distribution is,
859 p(x) dx = {1 \over (b-a)} dx
866 p(x) dx = @{1 \over (b-a)@} dx
872 @math{a <= x < b} and 0 otherwise.
875 @deftypefun double gsl_ran_flat_pdf (double @var{x}, double @var{a}, double @var{b})
876 This function computes the probability density @math{p(x)} at @var{x}
877 for a uniform distribution from @var{a} to @var{b}, using the formula
883 \centerline{\input rand-flat.tex}
886 @deftypefun double gsl_cdf_flat_P (double @var{x}, double @var{a}, double @var{b})
887 @deftypefunx double gsl_cdf_flat_Q (double @var{x}, double @var{a}, double @var{b})
888 @deftypefunx double gsl_cdf_flat_Pinv (double @var{P}, double @var{a}, double @var{b})
889 @deftypefunx double gsl_cdf_flat_Qinv (double @var{Q}, double @var{a}, double @var{b})
890 These functions compute the cumulative distribution functions
891 @math{P(x)}, @math{Q(x)} and their inverses for a uniform distribution
892 from @var{a} to @var{b}.
897 @node The Lognormal Distribution
898 @section The Lognormal Distribution
899 @deftypefun double gsl_ran_lognormal (const gsl_rng * @var{r}, double @var{zeta}, double @var{sigma})
900 @cindex Lognormal distribution
901 This function returns a random variate from the lognormal
902 distribution. The distribution function is,
906 p(x) dx = {1 \over x \sqrt{2 \pi \sigma^2}} \exp(-(\ln(x) - \zeta)^2/2 \sigma^2) dx
913 p(x) dx = @{1 \over x \sqrt@{2 \pi \sigma^2@} @} \exp(-(\ln(x) - \zeta)^2/2 \sigma^2) dx
921 @deftypefun double gsl_ran_lognormal_pdf (double @var{x}, double @var{zeta}, double @var{sigma})
922 This function computes the probability density @math{p(x)} at @var{x}
923 for a lognormal distribution with parameters @var{zeta} and @var{sigma},
924 using the formula given above.
929 \centerline{\input rand-lognormal.tex}
932 @deftypefun double gsl_cdf_lognormal_P (double @var{x}, double @var{zeta}, double @var{sigma})
933 @deftypefunx double gsl_cdf_lognormal_Q (double @var{x}, double @var{zeta}, double @var{sigma})
934 @deftypefunx double gsl_cdf_lognormal_Pinv (double @var{P}, double @var{zeta}, double @var{sigma})
935 @deftypefunx double gsl_cdf_lognormal_Qinv (double @var{Q}, double @var{zeta}, double @var{sigma})
936 These functions compute the cumulative distribution functions
937 @math{P(x)}, @math{Q(x)} and their inverses for the lognormal
938 distribution with parameters @var{zeta} and @var{sigma}.
943 @node The Chi-squared Distribution
944 @section The Chi-squared Distribution
945 The chi-squared distribution arises in statistics. If @math{Y_i} are
946 @math{n} independent gaussian random variates with unit variance then the
963 has a chi-squared distribution with @math{n} degrees of freedom.
965 @deftypefun double gsl_ran_chisq (const gsl_rng * @var{r}, double @var{nu})
966 @cindex Chi-squared distribution
967 This function returns a random variate from the chi-squared distribution
968 with @var{nu} degrees of freedom. The distribution function is,
972 p(x) dx = {1 \over 2 \Gamma(\nu/2) } (x/2)^{\nu/2 - 1} \exp(-x/2) dx
979 p(x) dx = @{1 \over 2 \Gamma(\nu/2) @} (x/2)^@{\nu/2 - 1@} \exp(-x/2) dx
988 @deftypefun double gsl_ran_chisq_pdf (double @var{x}, double @var{nu})
989 This function computes the probability density @math{p(x)} at @var{x}
990 for a chi-squared distribution with @var{nu} degrees of freedom, using
991 the formula given above.
996 \centerline{\input rand-chisq.tex}
999 @deftypefun double gsl_cdf_chisq_P (double @var{x}, double @var{nu})
1000 @deftypefunx double gsl_cdf_chisq_Q (double @var{x}, double @var{nu})
1001 @deftypefunx double gsl_cdf_chisq_Pinv (double @var{P}, double @var{nu})
1002 @deftypefunx double gsl_cdf_chisq_Qinv (double @var{Q}, double @var{nu})
1003 These functions compute the cumulative distribution functions
1004 @math{P(x)}, @math{Q(x)} and their inverses for the chi-squared
1005 distribution with @var{nu} degrees of freedom.
1011 @node The F-distribution
1012 @section The F-distribution
1013 The F-distribution arises in statistics. If @math{Y_1} and @math{Y_2}
1014 are chi-squared deviates with @math{\nu_1} and @math{\nu_2} degrees of
1015 freedom then the ratio,
1019 X = { (Y_1 / \nu_1) \over (Y_2 / \nu_2) }
1026 X = @{ (Y_1 / \nu_1) \over (Y_2 / \nu_2) @}
1031 has an F-distribution @math{F(x;\nu_1,\nu_2)}.
1033 @deftypefun double gsl_ran_fdist (const gsl_rng * @var{r}, double @var{nu1}, double @var{nu2})
1034 @cindex F-distribution
1035 This function returns a random variate from the F-distribution with degrees of freedom @var{nu1} and @var{nu2}. The distribution function is,
1040 { \Gamma((\nu_1 + \nu_2)/2)
1041 \over \Gamma(\nu_1/2) \Gamma(\nu_2/2) }
1042 \nu_1^{\nu_1/2} \nu_2^{\nu_2/2}
1043 x^{\nu_1/2 - 1} (\nu_2 + \nu_1 x)^{-\nu_1/2 -\nu_2/2}
1051 @{ \Gamma((\nu_1 + \nu_2)/2)
1052 \over \Gamma(\nu_1/2) \Gamma(\nu_2/2) @}
1053 \nu_1^@{\nu_1/2@} \nu_2^@{\nu_2/2@}
1054 x^@{\nu_1/2 - 1@} (\nu_2 + \nu_1 x)^@{-\nu_1/2 -\nu_2/2@}
1063 @deftypefun double gsl_ran_fdist_pdf (double @var{x}, double @var{nu1}, double @var{nu2})
1064 This function computes the probability density @math{p(x)} at @var{x}
1065 for an F-distribution with @var{nu1} and @var{nu2} degrees of freedom,
1066 using the formula given above.
1071 \centerline{\input rand-fdist.tex}
1074 @deftypefun double gsl_cdf_fdist_P (double @var{x}, double @var{nu1}, double @var{nu2})
1075 @deftypefunx double gsl_cdf_fdist_Q (double @var{x}, double @var{nu1}, double @var{nu2})
1076 @deftypefunx double gsl_cdf_fdist_Pinv (double @var{P}, double @var{nu1}, double @var{nu2})
1077 @deftypefunx double gsl_cdf_fdist_Qinv (double @var{Q}, double @var{nu1}, double @var{nu2})
1078 These functions compute the cumulative distribution functions
1079 @math{P(x)}, @math{Q(x)} and their inverses for the F-distribution
1080 with @var{nu1} and @var{nu2} degrees of freedom.
1084 @node The t-distribution
1085 @section The t-distribution
1086 The t-distribution arises in statistics. If @math{Y_1} has a normal
1087 distribution and @math{Y_2} has a chi-squared distribution with
1088 @math{\nu} degrees of freedom then the ratio,
1092 X = { Y_1 \over \sqrt{Y_2 / \nu} }
1099 X = @{ Y_1 \over \sqrt@{Y_2 / \nu@} @}
1104 has a t-distribution @math{t(x;\nu)} with @math{\nu} degrees of freedom.
1106 @deftypefun double gsl_ran_tdist (const gsl_rng * @var{r}, double @var{nu})
1107 @cindex t-distribution
1108 @cindex Student t-distribution
1109 This function returns a random variate from the t-distribution. The
1110 distribution function is,
1114 p(x) dx = {\Gamma((\nu + 1)/2) \over \sqrt{\pi \nu} \Gamma(\nu/2)}
1115 (1 + x^2/\nu)^{-(\nu + 1)/2} dx
1122 p(x) dx = @{\Gamma((\nu + 1)/2) \over \sqrt@{\pi \nu@} \Gamma(\nu/2)@}
1123 (1 + x^2/\nu)^@{-(\nu + 1)/2@} dx
1128 for @math{-\infty < x < +\infty}.
1131 @deftypefun double gsl_ran_tdist_pdf (double @var{x}, double @var{nu})
1132 This function computes the probability density @math{p(x)} at @var{x}
1133 for a t-distribution with @var{nu} degrees of freedom, using the formula
1139 \centerline{\input rand-tdist.tex}
1142 @deftypefun double gsl_cdf_tdist_P (double @var{x}, double @var{nu})
1143 @deftypefunx double gsl_cdf_tdist_Q (double @var{x}, double @var{nu})
1144 @deftypefunx double gsl_cdf_tdist_Pinv (double @var{P}, double @var{nu})
1145 @deftypefunx double gsl_cdf_tdist_Qinv (double @var{Q}, double @var{nu})
1146 These functions compute the cumulative distribution functions
1147 @math{P(x)}, @math{Q(x)} and their inverses for the t-distribution
1148 with @var{nu} degrees of freedom.
1152 @node The Beta Distribution
1153 @section The Beta Distribution
1154 @deftypefun double gsl_ran_beta (const gsl_rng * @var{r}, double @var{a}, double @var{b})
1155 @cindex Beta distribution
1156 This function returns a random variate from the beta
1157 distribution. The distribution function is,
1161 p(x) dx = {\Gamma(a+b) \over \Gamma(a) \Gamma(b)} x^{a-1} (1-x)^{b-1} dx
1168 p(x) dx = @{\Gamma(a+b) \over \Gamma(a) \Gamma(b)@} x^@{a-1@} (1-x)^@{b-1@} dx
1173 for @c{$0 \le x \le 1$}
1177 @deftypefun double gsl_ran_beta_pdf (double @var{x}, double @var{a}, double @var{b})
1178 This function computes the probability density @math{p(x)} at @var{x}
1179 for a beta distribution with parameters @var{a} and @var{b}, using the
1180 formula given above.
1185 \centerline{\input rand-beta.tex}
1188 @deftypefun double gsl_cdf_beta_P (double @var{x}, double @var{a}, double @var{b})
1189 @deftypefunx double gsl_cdf_beta_Q (double @var{x}, double @var{a}, double @var{b})
1190 @deftypefunx double gsl_cdf_beta_Pinv (double @var{P}, double @var{a}, double @var{b})
1191 @deftypefunx double gsl_cdf_beta_Qinv (double @var{Q}, double @var{a}, double @var{b})
1192 These functions compute the cumulative distribution functions
1193 @math{P(x)}, @math{Q(x)} and their inverses for the beta
1194 distribution with parameters @var{a} and @var{b}.
1198 @node The Logistic Distribution
1199 @section The Logistic Distribution
1201 @deftypefun double gsl_ran_logistic (const gsl_rng * @var{r}, double @var{a})
1202 @cindex Logistic distribution
1203 This function returns a random variate from the logistic
1204 distribution. The distribution function is,
1208 p(x) dx = { \exp(-x/a) \over a (1 + \exp(-x/a))^2 } dx
1215 p(x) dx = @{ \exp(-x/a) \over a (1 + \exp(-x/a))^2 @} dx
1220 for @math{-\infty < x < +\infty}.
1223 @deftypefun double gsl_ran_logistic_pdf (double @var{x}, double @var{a})
1224 This function computes the probability density @math{p(x)} at @var{x}
1225 for a logistic distribution with scale parameter @var{a}, using the
1226 formula given above.
1231 \centerline{\input rand-logistic.tex}
1234 @deftypefun double gsl_cdf_logistic_P (double @var{x}, double @var{a})
1235 @deftypefunx double gsl_cdf_logistic_Q (double @var{x}, double @var{a})
1236 @deftypefunx double gsl_cdf_logistic_Pinv (double @var{P}, double @var{a})
1237 @deftypefunx double gsl_cdf_logistic_Qinv (double @var{Q}, double @var{a})
1238 These functions compute the cumulative distribution functions
1239 @math{P(x)}, @math{Q(x)} and their inverses for the logistic
1240 distribution with scale parameter @var{a}.
1244 @node The Pareto Distribution
1245 @section The Pareto Distribution
1246 @deftypefun double gsl_ran_pareto (const gsl_rng * @var{r}, double @var{a}, double @var{b})
1247 @cindex Pareto distribution
1248 This function returns a random variate from the Pareto distribution of
1249 order @var{a}. The distribution function is,
1253 p(x) dx = (a/b) / (x/b)^{a+1} dx
1260 p(x) dx = (a/b) / (x/b)^@{a+1@} dx
1269 @deftypefun double gsl_ran_pareto_pdf (double @var{x}, double @var{a}, double @var{b})
1270 This function computes the probability density @math{p(x)} at @var{x}
1271 for a Pareto distribution with exponent @var{a} and scale @var{b}, using
1272 the formula given above.
1277 \centerline{\input rand-pareto.tex}
1280 @deftypefun double gsl_cdf_pareto_P (double @var{x}, double @var{a}, double @var{b})
1281 @deftypefunx double gsl_cdf_pareto_Q (double @var{x}, double @var{a}, double @var{b})
1282 @deftypefunx double gsl_cdf_pareto_Pinv (double @var{P}, double @var{a}, double @var{b})
1283 @deftypefunx double gsl_cdf_pareto_Qinv (double @var{Q}, double @var{a}, double @var{b})
1284 These functions compute the cumulative distribution functions
1285 @math{P(x)}, @math{Q(x)} and their inverses for the Pareto
1286 distribution with exponent @var{a} and scale @var{b}.
1290 @node Spherical Vector Distributions
1291 @section Spherical Vector Distributions
1293 The spherical distributions generate random vectors, located on a
1294 spherical surface. They can be used as random directions, for example in
1295 the steps of a random walk.
1297 @deftypefun void gsl_ran_dir_2d (const gsl_rng * @var{r}, double * @var{x}, double * @var{y})
1298 @deftypefunx void gsl_ran_dir_2d_trig_method (const gsl_rng * @var{r}, double * @var{x}, double * @var{y})
1299 @cindex 2D random direction vector
1300 @cindex direction vector, random 2D
1301 @cindex spherical random variates, 2D
1302 This function returns a random direction vector @math{v} =
1303 (@var{x},@var{y}) in two dimensions. The vector is normalized such that
1304 @math{|v|^2 = x^2 + y^2 = 1}. The obvious way to do this is to take a
1305 uniform random number between 0 and @math{2\pi} and let @var{x} and
1306 @var{y} be the sine and cosine respectively. Two trig functions would
1307 have been expensive in the old days, but with modern hardware
1308 implementations, this is sometimes the fastest way to go. This is the
1309 case for the Pentium (but not the case for the Sun Sparcstation).
1310 One can avoid the trig evaluations by choosing @var{x} and
1311 @var{y} in the interior of a unit circle (choose them at random from the
1312 interior of the enclosing square, and then reject those that are outside
1313 the unit circle), and then dividing by @c{$\sqrt{x^2 + y^2}$}
1314 @math{\sqrt@{x^2 + y^2@}}.
1315 A much cleverer approach, attributed to von Neumann (See Knuth, v2, 3rd
1316 ed, p140, exercise 23), requires neither trig nor a square root. In
1317 this approach, @var{u} and @var{v} are chosen at random from the
1318 interior of a unit circle, and then @math{x=(u^2-v^2)/(u^2+v^2)} and
1319 @math{y=2uv/(u^2+v^2)}.
1322 @deftypefun void gsl_ran_dir_3d (const gsl_rng * @var{r}, double * @var{x}, double * @var{y}, double * @var{z})
1323 @cindex 3D random direction vector
1324 @cindex direction vector, random 3D
1325 @cindex spherical random variates, 3D
1326 This function returns a random direction vector @math{v} =
1327 (@var{x},@var{y},@var{z}) in three dimensions. The vector is normalized
1328 such that @math{|v|^2 = x^2 + y^2 + z^2 = 1}. The method employed is
1329 due to Robert E. Knop (CACM 13, 326 (1970)), and explained in Knuth, v2,
1330 3rd ed, p136. It uses the surprising fact that the distribution
1331 projected along any axis is actually uniform (this is only true for 3
1335 @deftypefun void gsl_ran_dir_nd (const gsl_rng * @var{r}, size_t @var{n}, double * @var{x})
1336 @cindex N-dimensional random direction vector
1337 @cindex direction vector, random N-dimensional
1338 @cindex spherical random variates, N-dimensional
1340 This function returns a random direction vector
1341 @c{$v = (x_1,x_2,\ldots,x_n)$}
1342 @math{v = (x_1,x_2,...,x_n)} in @var{n} dimensions. The vector is normalized
1344 @c{$|v|^2 = x_1^2 + x_2^2 + \cdots + x_n^2 = 1$}
1345 @math{|v|^2 = x_1^2 + x_2^2 + ... + x_n^2 = 1}. The method
1346 uses the fact that a multivariate gaussian distribution is spherically
1347 symmetric. Each component is generated to have a gaussian distribution,
1348 and then the components are normalized. The method is described by
1349 Knuth, v2, 3rd ed, p135--136, and attributed to G. W. Brown, Modern
1350 Mathematics for the Engineer (1956).
1354 @node The Weibull Distribution
1355 @section The Weibull Distribution
1356 @deftypefun double gsl_ran_weibull (const gsl_rng * @var{r}, double @var{a}, double @var{b})
1357 @cindex Weibull distribution
1358 This function returns a random variate from the Weibull distribution. The
1359 distribution function is,
1363 p(x) dx = {b \over a^b} x^{b-1} \exp(-(x/a)^b) dx
1370 p(x) dx = @{b \over a^b@} x^@{b-1@} \exp(-(x/a)^b) dx
1379 @deftypefun double gsl_ran_weibull_pdf (double @var{x}, double @var{a}, double @var{b})
1380 This function computes the probability density @math{p(x)} at @var{x}
1381 for a Weibull distribution with scale @var{a} and exponent @var{b},
1382 using the formula given above.
1387 \centerline{\input rand-weibull.tex}
1390 @deftypefun double gsl_cdf_weibull_P (double @var{x}, double @var{a}, double @var{b})
1391 @deftypefunx double gsl_cdf_weibull_Q (double @var{x}, double @var{a}, double @var{b})
1392 @deftypefunx double gsl_cdf_weibull_Pinv (double @var{P}, double @var{a}, double @var{b})
1393 @deftypefunx double gsl_cdf_weibull_Qinv (double @var{Q}, double @var{a}, double @var{b})
1394 These functions compute the cumulative distribution functions
1395 @math{P(x)}, @math{Q(x)} and their inverses for the Weibull
1396 distribution with scale @var{a} and exponent @var{b}.
1401 @node The Type-1 Gumbel Distribution
1402 @section The Type-1 Gumbel Distribution
1403 @deftypefun double gsl_ran_gumbel1 (const gsl_rng * @var{r}, double @var{a}, double @var{b})
1404 @cindex Gumbel distribution (Type 1)
1405 @cindex Type 1 Gumbel distribution, random variates
1406 This function returns a random variate from the Type-1 Gumbel
1407 distribution. The Type-1 Gumbel distribution function is,
1411 p(x) dx = a b \exp(-(b \exp(-ax) + ax)) dx
1418 p(x) dx = a b \exp(-(b \exp(-ax) + ax)) dx
1423 for @math{-\infty < x < \infty}.
1426 @deftypefun double gsl_ran_gumbel1_pdf (double @var{x}, double @var{a}, double @var{b})
1427 This function computes the probability density @math{p(x)} at @var{x}
1428 for a Type-1 Gumbel distribution with parameters @var{a} and @var{b},
1429 using the formula given above.
1434 \centerline{\input rand-gumbel1.tex}
1437 @deftypefun double gsl_cdf_gumbel1_P (double @var{x}, double @var{a}, double @var{b})
1438 @deftypefunx double gsl_cdf_gumbel1_Q (double @var{x}, double @var{a}, double @var{b})
1439 @deftypefunx double gsl_cdf_gumbel1_Pinv (double @var{P}, double @var{a}, double @var{b})
1440 @deftypefunx double gsl_cdf_gumbel1_Qinv (double @var{Q}, double @var{a}, double @var{b})
1441 These functions compute the cumulative distribution functions
1442 @math{P(x)}, @math{Q(x)} and their inverses for the Type-1 Gumbel
1443 distribution with parameters @var{a} and @var{b}.
1448 @node The Type-2 Gumbel Distribution
1449 @section The Type-2 Gumbel Distribution
1450 @deftypefun double gsl_ran_gumbel2 (const gsl_rng * @var{r}, double @var{a}, double @var{b})
1451 @cindex Gumbel distribution (Type 2)
1452 @cindex Type 2 Gumbel distribution
1453 This function returns a random variate from the Type-2 Gumbel
1454 distribution. The Type-2 Gumbel distribution function is,
1458 p(x) dx = a b x^{-a-1} \exp(-b x^{-a}) dx
1465 p(x) dx = a b x^@{-a-1@} \exp(-b x^@{-a@}) dx
1470 for @math{0 < x < \infty}.
1473 @deftypefun double gsl_ran_gumbel2_pdf (double @var{x}, double @var{a}, double @var{b})
1474 This function computes the probability density @math{p(x)} at @var{x}
1475 for a Type-2 Gumbel distribution with parameters @var{a} and @var{b},
1476 using the formula given above.
1481 \centerline{\input rand-gumbel2.tex}
1484 @deftypefun double gsl_cdf_gumbel2_P (double @var{x}, double @var{a}, double @var{b})
1485 @deftypefunx double gsl_cdf_gumbel2_Q (double @var{x}, double @var{a}, double @var{b})
1486 @deftypefunx double gsl_cdf_gumbel2_Pinv (double @var{P}, double @var{a}, double @var{b})
1487 @deftypefunx double gsl_cdf_gumbel2_Qinv (double @var{Q}, double @var{a}, double @var{b})
1488 These functions compute the cumulative distribution functions
1489 @math{P(x)}, @math{Q(x)} and their inverses for the Type-2 Gumbel
1490 distribution with parameters @var{a} and @var{b}.
1495 @node The Dirichlet Distribution
1496 @section The Dirichlet Distribution
1497 @deftypefun void gsl_ran_dirichlet (const gsl_rng * @var{r}, size_t @var{K}, const double @var{alpha}[], double @var{theta}[])
1498 @cindex Dirichlet distribution
1499 This function returns an array of @var{K} random variates from a Dirichlet
1500 distribution of order @var{K}-1. The distribution function is
1504 p(\theta_1,\ldots,\theta_K) \, d\theta_1 \cdots d\theta_K =
1505 {1 \over Z} \prod_{i=1}^{K} \theta_i^{\alpha_i - 1}
1506 \; \delta(1 -\sum_{i=1}^K \theta_i) d\theta_1 \cdots d\theta_K
1513 p(\theta_1, ..., \theta_K) d\theta_1 ... d\theta_K =
1514 (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
1519 for @c{$\theta_i \ge 0$}
1521 and @c{$\alpha_i \ge 0$}
1522 @math{alpha_i >= 0}. The delta function ensures that @math{\sum \theta_i = 1}.
1523 The normalization factor @math{Z} is
1527 Z = {\prod_{i=1}^K \Gamma(\alpha_i) \over \Gamma( \sum_{i=1}^K \alpha_i)}
1534 Z = @{\prod_@{i=1@}^K \Gamma(\alpha_i)@} / @{\Gamma( \sum_@{i=1@}^K \alpha_i)@}
1538 The random variates are generated by sampling @var{K} values
1539 from gamma distributions with parameters
1540 @c{$a=\alpha_i$, $b=1$}
1541 @math{a=alpha_i, b=1},
1543 See A.M. Law, W.D. Kelton, @cite{Simulation Modeling and Analysis} (1991).
1546 @deftypefun double gsl_ran_dirichlet_pdf (size_t @var{K}, const double @var{alpha}[], const double @var{theta}[])
1547 This function computes the probability density
1548 @c{$p(\theta_1, \ldots , \theta_K)$}
1549 @math{p(\theta_1, ... , \theta_K)}
1550 at @var{theta}[@var{K}] for a Dirichlet distribution with parameters
1551 @var{alpha}[@var{K}], using the formula given above.
1554 @deftypefun double gsl_ran_dirichlet_lnpdf (size_t @var{K}, const double @var{alpha}[], const double @var{theta}[])
1555 This function computes the logarithm of the probability density
1556 @c{$p(\theta_1, \ldots , \theta_K)$}
1557 @math{p(\theta_1, ... , \theta_K)}
1558 for a Dirichlet distribution with parameters
1559 @var{alpha}[@var{K}].
1563 @node General Discrete Distributions
1564 @section General Discrete Distributions
1566 Given @math{K} discrete events with different probabilities @math{P[k]},
1567 produce a random value @math{k} consistent with its probability.
1569 The obvious way to do this is to preprocess the probability list by
1570 generating a cumulative probability array with @math{K+1} elements:
1576 C[k+1] &= C[k]+P[k].
1590 Note that this construction produces @math{C[K]=1}. Now choose a
1591 uniform deviate @math{u} between 0 and 1, and find the value of @math{k}
1592 such that @c{$C[k] \le u < C[k+1]$}
1593 @math{C[k] <= u < C[k+1]}.
1594 Although this in principle requires of order @math{\log K} steps per
1595 random number generation, they are fast steps, and if you use something
1596 like @math{\lfloor uK \rfloor} as a starting point, you can often do
1599 But faster methods have been devised. Again, the idea is to preprocess
1600 the probability list, and save the result in some form of lookup table;
1601 then the individual calls for a random discrete event can go rapidly.
1602 An approach invented by G. Marsaglia (Generating discrete random numbers
1603 in a computer, Comm ACM 6, 37--38 (1963)) is very clever, and readers
1604 interested in examples of good algorithm design are directed to this
1605 short and well-written paper. Unfortunately, for large @math{K},
1606 Marsaglia's lookup table can be quite large.
1608 A much better approach is due to Alastair J. Walker (An efficient method
1609 for generating discrete random variables with general distributions, ACM
1610 Trans on Mathematical Software 3, 253--256 (1977); see also Knuth, v2,
1611 3rd ed, p120--121,139). This requires two lookup tables, one floating
1612 point and one integer, but both only of size @math{K}. After
1613 preprocessing, the random numbers are generated in O(1) time, even for
1614 large @math{K}. The preprocessing suggested by Walker requires
1615 @math{O(K^2)} effort, but that is not actually necessary, and the
1616 implementation provided here only takes @math{O(K)} effort. In general,
1617 more preprocessing leads to faster generation of the individual random
1618 numbers, but a diminishing return is reached pretty early. Knuth points
1619 out that the optimal preprocessing is combinatorially difficult for
1622 This method can be used to speed up some of the discrete random number
1623 generators below, such as the binomial distribution. To use it for
1624 something like the Poisson Distribution, a modification would have to
1625 be made, since it only takes a finite set of @math{K} outcomes.
1627 @deftypefun {gsl_ran_discrete_t *} gsl_ran_discrete_preproc (size_t @var{K}, const double * @var{P})
1628 @cindex Discrete random numbers
1629 @cindex Discrete random numbers, preprocessing
1630 This function returns a pointer to a structure that contains the lookup
1631 table for the discrete random number generator. The array @var{P}[] contains
1632 the probabilities of the discrete events; these array elements must all be
1633 positive, but they needn't add up to one (so you can think of them more
1634 generally as ``weights'')---the preprocessor will normalize appropriately.
1635 This return value is used
1636 as an argument for the @code{gsl_ran_discrete} function below.
1639 @deftypefun {size_t} gsl_ran_discrete (const gsl_rng * @var{r}, const gsl_ran_discrete_t * @var{g})
1640 @cindex Discrete random numbers
1641 After the preprocessor, above, has been called, you use this function to
1642 get the discrete random numbers.
1645 @deftypefun {double} gsl_ran_discrete_pdf (size_t @var{k}, const gsl_ran_discrete_t * @var{g})
1646 @cindex Discrete random numbers
1647 Returns the probability @math{P[k]} of observing the variable @var{k}.
1648 Since @math{P[k]} is not stored as part of the lookup table, it must be
1649 recomputed; this computation takes @math{O(K)}, so if @var{K} is large
1650 and you care about the original array @math{P[k]} used to create the
1651 lookup table, then you should just keep this original array @math{P[k]}
1655 @deftypefun {void} gsl_ran_discrete_free (gsl_ran_discrete_t * @var{g})
1656 @cindex Discrete random numbers
1657 De-allocates the lookup table pointed to by @var{g}.
1661 @node The Poisson Distribution
1662 @section The Poisson Distribution
1663 @deftypefun {unsigned int} gsl_ran_poisson (const gsl_rng * @var{r}, double @var{mu})
1664 @cindex Poisson random numbers
1665 This function returns a random integer from the Poisson distribution
1666 with mean @var{mu}. The probability distribution for Poisson variates is,
1670 p(k) = {\mu^k \over k!} \exp(-\mu)
1677 p(k) = @{\mu^k \over k!@} \exp(-\mu)
1686 @deftypefun double gsl_ran_poisson_pdf (unsigned int @var{k}, double @var{mu})
1687 This function computes the probability @math{p(k)} of obtaining @var{k}
1688 from a Poisson distribution with mean @var{mu}, using the formula
1694 \centerline{\input rand-poisson.tex}
1697 @deftypefun double gsl_cdf_poisson_P (unsigned int @var{k}, double @var{mu})
1698 @deftypefunx double gsl_cdf_poisson_Q (unsigned int @var{k}, double @var{mu})
1699 These functions compute the cumulative distribution functions
1700 @math{P(k)}, @math{Q(k)} for the Poisson distribution with parameter
1706 @node The Bernoulli Distribution
1707 @section The Bernoulli Distribution
1708 @deftypefun {unsigned int} gsl_ran_bernoulli (const gsl_rng * @var{r}, double @var{p})
1709 @cindex Bernoulli trial, random variates
1710 This function returns either 0 or 1, the result of a Bernoulli trial
1711 with probability @var{p}. The probability distribution for a Bernoulli
1733 @deftypefun double gsl_ran_bernoulli_pdf (unsigned int @var{k}, double @var{p})
1734 This function computes the probability @math{p(k)} of obtaining
1735 @var{k} from a Bernoulli distribution with probability parameter
1736 @var{p}, using the formula given above.
1741 \centerline{\input rand-bernoulli.tex}
1745 @node The Binomial Distribution
1746 @section The Binomial Distribution
1747 @deftypefun {unsigned int} gsl_ran_binomial (const gsl_rng * @var{r}, double @var{p}, unsigned int @var{n})
1748 @cindex Binomial random variates
1749 This function returns a random integer from the binomial distribution,
1750 the number of successes in @var{n} independent trials with probability
1751 @var{p}. The probability distribution for binomial variates is,
1755 p(k) = {n! \over k! (n-k)!} p^k (1-p)^{n-k}
1762 p(k) = @{n! \over k! (n-k)! @} p^k (1-p)^@{n-k@}
1767 for @c{$0 \le k \le n$}
1771 @deftypefun double gsl_ran_binomial_pdf (unsigned int @var{k}, double @var{p}, unsigned int @var{n})
1772 This function computes the probability @math{p(k)} of obtaining @var{k}
1773 from a binomial distribution with parameters @var{p} and @var{n}, using
1774 the formula given above.
1779 \centerline{\input rand-binomial.tex}
1782 @deftypefun double gsl_cdf_binomial_P (unsigned int @var{k}, double @var{p}, unsigned int @var{n})
1783 @deftypefunx double gsl_cdf_binomial_Q (unsigned int @var{k}, double @var{p}, unsigned int @var{n})
1784 These functions compute the cumulative distribution functions
1785 @math{P(k)}, @math{Q(k)} for the binomial
1786 distribution with parameters @var{p} and @var{n}.
1791 @node The Multinomial Distribution
1792 @section The Multinomial Distribution
1793 @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}[])
1794 @cindex Multinomial distribution
1796 This function computes a random sample @var{n}[] from the multinomial
1797 distribution formed by @var{N} trials from an underlying distribution
1798 @var{p}[@var{K}]. The distribution function for @var{n}[] is,
1802 P(n_1, n_2,\cdots, n_K) = {{ N!}\over{n_1 ! n_2 ! \cdots n_K !}} \,
1803 p_1^{n_1} p_2^{n_2} \cdots p_K^{n_K}
1810 P(n_1, n_2, ..., n_K) =
1811 (N!/(n_1! n_2! ... n_K!)) p_1^n_1 p_2^n_2 ... p_K^n_K
1816 where @c{($n_1$, $n_2$, $\ldots$, $n_K$)}
1817 @math{(n_1, n_2, ..., n_K)}
1818 are nonnegative integers with
1819 @c{$\sum_{k=1}^{K} n_k =N$}
1820 @math{sum_@{k=1@}^K n_k = N},
1822 @c{$(p_1, p_2, \ldots, p_K)$}
1823 @math{(p_1, p_2, ..., p_K)}
1824 is a probability distribution with @math{\sum p_i = 1}.
1825 If the array @var{p}[@var{K}] is not normalized then its entries will be
1826 treated as weights and normalized appropriately. The arrays @var{n}[]
1827 and @var{p}[] must both be of length @var{K}.
1829 Random variates are generated using the conditional binomial method (see
1830 C.S. David, @cite{The computer generation of multinomial random
1831 variates}, Comp. Stat. Data Anal. 16 (1993) 205--217 for details).
1834 @deftypefun double gsl_ran_multinomial_pdf (size_t @var{K}, const double @var{p}[], const unsigned int @var{n}[])
1835 This function computes the probability
1836 @c{$P(n_1, n_2, \ldots, n_K)$}
1837 @math{P(n_1, n_2, ..., n_K)}
1838 of sampling @var{n}[@var{K}] from a multinomial distribution
1839 with parameters @var{p}[@var{K}], using the formula given above.
1842 @deftypefun double gsl_ran_multinomial_lnpdf (size_t @var{K}, const double @var{p}[], const unsigned int @var{n}[])
1843 This function returns the logarithm of the probability for the
1844 multinomial distribution @c{$P(n_1, n_2, \ldots, n_K)$}
1845 @math{P(n_1, n_2, ..., n_K)} with parameters @var{p}[@var{K}].
1849 @node The Negative Binomial Distribution
1850 @section The Negative Binomial Distribution
1851 @deftypefun {unsigned int} gsl_ran_negative_binomial (const gsl_rng * @var{r}, double @var{p}, double @var{n})
1852 @cindex Negative Binomial distribution, random variates
1853 This function returns a random integer from the negative binomial
1854 distribution, the number of failures occurring before @var{n} successes
1855 in independent trials with probability @var{p} of success. The
1856 probability distribution for negative binomial variates is,
1860 p(k) = {\Gamma(n + k) \over \Gamma(k+1) \Gamma(n) } p^n (1-p)^k
1867 p(k) = @{\Gamma(n + k) \over \Gamma(k+1) \Gamma(n) @} p^n (1-p)^k
1872 Note that @math{n} is not required to be an integer.
1875 @deftypefun double gsl_ran_negative_binomial_pdf (unsigned int @var{k}, double @var{p}, double @var{n})
1876 This function computes the probability @math{p(k)} of obtaining @var{k}
1877 from a negative binomial distribution with parameters @var{p} and
1878 @var{n}, using the formula given above.
1883 \centerline{\input rand-nbinomial.tex}
1886 @deftypefun double gsl_cdf_negative_binomial_P (unsigned int @var{k}, double @var{p}, double @var{n})
1887 @deftypefunx double gsl_cdf_negative_binomial_Q (unsigned int @var{k}, double @var{p}, double @var{n})
1888 These functions compute the cumulative distribution functions
1889 @math{P(k)}, @math{Q(k)} for the negative binomial distribution with
1890 parameters @var{p} and @var{n}.
1894 @node The Pascal Distribution
1895 @section The Pascal Distribution
1897 @deftypefun {unsigned int} gsl_ran_pascal (const gsl_rng * @var{r}, double @var{p}, unsigned int @var{n})
1898 This function returns a random integer from the Pascal distribution. The
1899 Pascal distribution is simply a negative binomial distribution with an
1900 integer value of @math{n}.
1904 p(k) = {(n + k - 1)! \over k! (n - 1)! } p^n (1-p)^k
1911 p(k) = @{(n + k - 1)! \over k! (n - 1)! @} p^n (1-p)^k
1920 @deftypefun double gsl_ran_pascal_pdf (unsigned int @var{k}, double @var{p}, unsigned int @var{n})
1921 This function computes the probability @math{p(k)} of obtaining @var{k}
1922 from a Pascal distribution with parameters @var{p} and
1923 @var{n}, using the formula given above.
1928 \centerline{\input rand-pascal.tex}
1931 @deftypefun double gsl_cdf_pascal_P (unsigned int @var{k}, double @var{p}, unsigned int @var{n})
1932 @deftypefunx double gsl_cdf_pascal_Q (unsigned int @var{k}, double @var{p}, unsigned int @var{n})
1933 These functions compute the cumulative distribution functions
1934 @math{P(k)}, @math{Q(k)} for the Pascal distribution with
1935 parameters @var{p} and @var{n}.
1939 @node The Geometric Distribution
1940 @section The Geometric Distribution
1941 @deftypefun {unsigned int} gsl_ran_geometric (const gsl_rng * @var{r}, double @var{p})
1942 @cindex Geometric random variates
1943 This function returns a random integer from the geometric distribution,
1944 the number of independent trials with probability @var{p} until the
1945 first success. The probability distribution for geometric variates
1950 p(k) = p (1-p)^{k-1}
1957 p(k) = p (1-p)^(k-1)
1963 @math{k >= 1}. Note that the distribution begins with @math{k=1} with this
1964 definition. There is another convention in which the exponent @math{k-1}
1965 is replaced by @math{k}.
1968 @deftypefun double gsl_ran_geometric_pdf (unsigned int @var{k}, double @var{p})
1969 This function computes the probability @math{p(k)} of obtaining @var{k}
1970 from a geometric distribution with probability parameter @var{p}, using
1971 the formula given above.
1976 \centerline{\input rand-geometric.tex}
1979 @deftypefun double gsl_cdf_geometric_P (unsigned int @var{k}, double @var{p})
1980 @deftypefunx double gsl_cdf_geometric_Q (unsigned int @var{k}, double @var{p})
1981 These functions compute the cumulative distribution functions
1982 @math{P(k)}, @math{Q(k)} for the geometric distribution with parameter
1988 @node The Hypergeometric Distribution
1989 @section The Hypergeometric Distribution
1990 @cindex hypergeometric random variates
1991 @deftypefun {unsigned int} gsl_ran_hypergeometric (const gsl_rng * @var{r}, unsigned int @var{n1}, unsigned int @var{n2}, unsigned int @var{t})
1992 @cindex Geometric random variates
1993 This function returns a random integer from the hypergeometric
1994 distribution. The probability distribution for hypergeometric
1999 p(k) = C(n_1, k) C(n_2, t - k) / C(n_1 + n_2, t)
2006 p(k) = C(n_1, k) C(n_2, t - k) / C(n_1 + n_2, t)
2011 where @math{C(a,b) = a!/(b!(a-b)!)} and
2012 @c{$t \leq n_1 + n_2$}
2013 @math{t <= n_1 + n_2}. The domain of @math{k} is
2014 @c{$\hbox{max}(0,t-n_2), \ldots, \hbox{min}(t,n_1)$}
2015 @math{max(0,t-n_2), ..., min(t,n_1)}.
2017 If a population contains @math{n_1} elements of ``type 1'' and
2018 @math{n_2} elements of ``type 2'' then the hypergeometric
2019 distribution gives the probability of obtaining @math{k} elements of
2020 ``type 1'' in @math{t} samples from the population without
2024 @deftypefun double gsl_ran_hypergeometric_pdf (unsigned int @var{k}, unsigned int @var{n1}, unsigned int @var{n2}, unsigned int @var{t})
2025 This function computes the probability @math{p(k)} of obtaining @var{k}
2026 from a hypergeometric distribution with parameters @var{n1}, @var{n2},
2027 @var{t}, using the formula given above.
2032 \centerline{\input rand-hypergeometric.tex}
2035 @deftypefun double gsl_cdf_hypergeometric_P (unsigned int @var{k}, unsigned int @var{n1}, unsigned int @var{n2}, unsigned int @var{t})
2036 @deftypefunx double gsl_cdf_hypergeometric_Q (unsigned int @var{k}, unsigned int @var{n1}, unsigned int @var{n2}, unsigned int @var{t})
2037 These functions compute the cumulative distribution functions
2038 @math{P(k)}, @math{Q(k)} for the hypergeometric distribution with
2039 parameters @var{n1}, @var{n2} and @var{t}.
2044 @node The Logarithmic Distribution
2045 @section The Logarithmic Distribution
2046 @deftypefun {unsigned int} gsl_ran_logarithmic (const gsl_rng * @var{r}, double @var{p})
2047 @cindex Logarithmic random variates
2048 This function returns a random integer from the logarithmic
2049 distribution. The probability distribution for logarithmic random variates
2054 p(k) = {-1 \over \log(1-p)} {\left( p^k \over k \right)}
2061 p(k) = @{-1 \over \log(1-p)@} @{(p^k \over k)@}
2070 @deftypefun double gsl_ran_logarithmic_pdf (unsigned int @var{k}, double @var{p})
2071 This function computes the probability @math{p(k)} of obtaining @var{k}
2072 from a logarithmic distribution with probability parameter @var{p},
2073 using the formula given above.
2078 \centerline{\input rand-logarithmic.tex}
2082 @node Shuffling and Sampling
2083 @section Shuffling and Sampling
2085 The following functions allow the shuffling and sampling of a set of
2086 objects. The algorithms rely on a random number generator as a source of
2087 randomness and a poor quality generator can lead to correlations in the
2088 output. In particular it is important to avoid generators with a short
2089 period. For more information see Knuth, v2, 3rd ed, Section 3.4.2,
2090 ``Random Sampling and Shuffling''.
2092 @deftypefun void gsl_ran_shuffle (const gsl_rng * @var{r}, void * @var{base}, size_t @var{n}, size_t @var{size})
2094 This function randomly shuffles the order of @var{n} objects, each of
2095 size @var{size}, stored in the array @var{base}[0..@var{n}-1]. The
2096 output of the random number generator @var{r} is used to produce the
2097 permutation. The algorithm generates all possible @math{n!}
2098 permutations with equal probability, assuming a perfect source of random
2101 The following code shows how to shuffle the numbers from 0 to 51,
2106 for (i = 0; i < 52; i++)
2111 gsl_ran_shuffle (r, a, 52, sizeof (int));
2116 @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})
2117 This function fills the array @var{dest}[k] with @var{k} objects taken
2118 randomly from the @var{n} elements of the array
2119 @var{src}[0..@var{n}-1]. The objects are each of size @var{size}. The
2120 output of the random number generator @var{r} is used to make the
2121 selection. The algorithm ensures all possible samples are equally
2122 likely, assuming a perfect source of randomness.
2124 The objects are sampled @emph{without} replacement, thus each object can
2125 only appear once in @var{dest}[k]. It is required that @var{k} be less
2126 than or equal to @code{n}. The objects in @var{dest} will be in the
2127 same relative order as those in @var{src}. You will need to call
2128 @code{gsl_ran_shuffle(r, dest, n, size)} if you want to randomize the
2131 The following code shows how to select a random sample of three unique
2132 numbers from the set 0 to 99,
2135 double a[3], b[100];
2137 for (i = 0; i < 100; i++)
2142 gsl_ran_choose (r, a, 3, b, 100, sizeof (double));
2147 @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})
2148 This function is like @code{gsl_ran_choose} but samples @var{k} items
2149 from the original array of @var{n} items @var{src} with replacement, so
2150 the same object can appear more than once in the output sequence
2151 @var{dest}. There is no requirement that @var{k} be less than @var{n}
2156 @node Random Number Distribution Examples
2159 The following program demonstrates the use of a random number generator
2160 to produce variates from a distribution. It prints 10 samples from the
2161 Poisson distribution with a mean of 3.
2164 @verbatiminclude examples/randpoisson.c
2168 If the library and header files are installed under @file{/usr/local}
2169 (the default location) then the program can be compiled with these
2173 $ gcc -Wall demo.c -lgsl -lgslcblas -lm
2177 Here is the output of the program,
2181 @verbatiminclude examples/randpoisson.out
2185 The variates depend on the seed used by the generator. The seed for the
2186 default generator type @code{gsl_rng_default} can be changed with the
2187 @code{GSL_RNG_SEED} environment variable to produce a different stream
2191 $ GSL_RNG_SEED=123 ./a.out
2192 @verbatiminclude examples/randpoisson.2.out
2196 The following program generates a random walk in two dimensions.
2199 @verbatiminclude examples/randwalk.c
2203 Here is the output from the program, three 10-step random walks from the origin,
2206 \centerline{\input random-walk.tex}
2209 The following program computes the upper and lower cumulative
2210 distribution functions for the standard normal distribution at
2214 @verbatiminclude examples/cdf.c
2218 Here is the output of the program,
2221 @verbatiminclude examples/cdf.out
2224 @node Random Number Distribution References and Further Reading
2225 @section References and Further Reading
2227 For an encyclopaedic coverage of the subject readers are advised to
2228 consult the book @cite{Non-Uniform Random Variate Generation} by Luc
2229 Devroye. It covers every imaginable distribution and provides hundreds
2234 Luc Devroye, @cite{Non-Uniform Random Variate Generation},
2235 Springer-Verlag, ISBN 0-387-96305-7.
2239 The subject of random variate generation is also reviewed by Knuth, who
2240 describes algorithms for all the major distributions.
2244 Donald E. Knuth, @cite{The Art of Computer Programming: Seminumerical
2245 Algorithms} (Vol 2, 3rd Ed, 1997), Addison-Wesley, ISBN 0201896842.
2249 The Particle Data Group provides a short review of techniques for
2250 generating distributions of random numbers in the ``Monte Carlo''
2251 section of its Annual Review of Particle Physics.
2255 @cite{Review of Particle Properties}
2256 R.M. Barnett et al., Physical Review D54, 1 (1996)
2257 @uref{http://pdg.lbl.gov/}.
2261 The Review of Particle Physics is available online in postscript and pdf
2265 An overview of methods used to compute cumulative distribution functions
2266 can be found in @cite{Statistical Computing} by W.J. Kennedy and
2267 J.E. Gentle. Another general reference is @cite{Elements of Statistical
2268 Computing} by R.A. Thisted.
2272 William E. Kennedy and James E. Gentle, @cite{Statistical Computing} (1980),
2273 Marcel Dekker, ISBN 0-8247-6898-1.
2278 Ronald A. Thisted, @cite{Elements of Statistical Computing} (1988),
2279 Chapman & Hall, ISBN 0-412-01371-1.
2283 The cumulative distribution functions for the Gaussian distribution
2284 are based on the following papers,
2288 @cite{Rational Chebyshev Approximations Using Linear Equations},
2289 W.J. Cody, W. Fraser, J.F. Hart. Numerische Mathematik 12, 242--251 (1968).
2294 @cite{Rational Chebyshev Approximations for the Error Function},
2295 W.J. Cody. Mathematics of Computation 23, n107, 631--637 (July 1969).