randist.texi

用于VC.net的gsl的lib库文件包

@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 The Spherical Distribution (2D & 3D)
@section The Spherical Distribution (2D & 3D)

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.

@deftypefn Random void gsl_ran_dir_2d (const gsl_rng * @var{r}, double *@var{x}, double *@var{y})
@deftypefnx Random 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=uv/(u^2+v^2)}.
@end deftypefn



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