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