randist.texi

开放gsl矩阵运算

for @c{$0 \le x \le 1$}
@math{0 <= x <= 1}.
@end deftypefn

@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

@page
@node The Logistic Distribution
@section The Logistic Distribution

@deftypefn Random double gsl_ran_logistic (const gsl_rng * @var{r}, double @var{a})
@cindex Logistic random variates
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 deftypefn

@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

@page
@node The Pareto Distribution
@section The Pareto Distribution
@deftypefn Random double gsl_ran_pareto (const gsl_rng * @var{r}, double @var{a}, double @var{b})
@cindex Pareto random variates
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 deftypefn

@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

@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 my home Pentium (but not the case for my Sun Sparcstation 20 at
work).  Once 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}, int @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 random variates
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

@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), random variates
@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

@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), random variates
@cindex Type 2 Gumbel distribution, random variate
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

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

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

@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

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

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

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

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