randist.texi

该文件为c++的数学函数库!是一个非常有用的编程工具.它含有各种数学函数,为科学计算、工程应用等程序编写提供方便！

@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 3d).
@end deftypefn

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

@page
@node The Weibull Distribution
@section The Weibull Distribution
@deftypefn Random 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 deftypefn

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

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

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

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

@deftypefn Random {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 deftypefn

@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
@deftypefn Random {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 deftypefn

@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

@page
@node The Bernoulli Distribution
@section The Bernoulli Distribution
@deftypefn Random {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 deftypefn

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