randist.texi

开放gsl矩阵运算

\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 deftypefn

@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

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

@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

@page
@node The Pascal Distribution
@section The Pascal Distribution

@deftypefn Random {unsigned int} gsl_ran_pascal (const gsl_rng * @var{r}, double @var{p}, unsigned int @var{k})
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 deftypefn

@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

@page
@node The Geometric Distribution
@section The Geometric Distribution
@deftypefn Random {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}.
@end deftypefn

@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

@page
@node The Hypergeometric Distribution
@section The Hypergeometric Distribution
@deftypefn Random {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,k)
$$
\afterdisplay
@end tex
@ifinfo
@example
p(k) =  C(n_1,k) C(n_2, t-k) / C(n_1 + n_2,k)
@end example
@end ifinfo
@noindent
where @math{C(a,b) = a!/(b!(a-b)!)}.  The domain of @math{k} is
@c{$\hbox{max}(0\,t-n_2), \ldots, \hbox{max}(t,n_1)$}
@math{max(0,t-n_2), ..., max(t,n_1)}.
@end deftypefn

@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{n3}, using the formula given above.
@end deftypefun

@sp 1
@tex
\centerline{\input rand-hypergeometric.tex}
@end tex

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

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

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

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

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


@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
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

int
main (void)
@{
  const gsl_rng_type * T;
  gsl_rng * r;

  int i, n = 10;
  double mu = 3.0;

  /* create a generator chosen by the 
     environment variable GSL_RNG_TYPE */

  gsl_rng_env_setup();

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  /* print n random variates chosen from 
     the poisson distribution with mean 
     parameter mu */

  for (i = 0; i < n; i++) 
    @{
      unsigned int k = gsl_ran_poisson (r, mu);
      printf(" %u", k);
    @}

  printf("\n");
  return 0;
@}
@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 demo.c -lgsl -lgslcblas -lm
@end example

@noindent
Here is the output of the program,

@example
$ ./a.out 
 4 2 3 3 1 3 4 1 3 5
@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 
GSL_RNG_SEED=123
 1 1 2 1 2 6 2 1 8 7
@end example

@noindent
The following program generates a random walk in two dimensions.

@example
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

int
main (void)
@{
  const gsl_rng_type * T;
  gsl_rng * r;

  gsl_rng_env_setup();
  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  int i;
  double x = 0, y = 0, dx, dy;

  printf("%g %g\n", x, y);

  for (i = 0; i < 10; i++)
    @{
      gsl_ran_dir_2d (r, &dx, &dy);
      x += dx; y += dy; 
      printf("%g %g\n", x, y);
    @}
  return 0;
@}
@end example

@noindent
Example output from the program, three 10-step random walks from the origin.

@sp 1
@tex
\centerline{\input random-walk.tex}
@end tex

@node Random Number Distribution References and Further Reading
@section References and Further Reading
@noindent
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)
@url{http://pdg.lbl.gov/}.
@end itemize

@noindent
The Review of Particle Physics is available online in postscript and pdf
format.