montecarlo.texi

开放gsl矩阵运算

@cindex Monte Carlo integration

This chapter describes routines for multidimensional Monte Carlo
integration.  These include the traditional Monte Carlo method and
adaptive algorithms such as @sc{vegas} and @sc{miser} which use
importance sampling and stratified sampling techniques. Each algorithm
computes an estimate of a multidimensional definite integral of the
form,

@tex
\beforedisplay
$$
I = \int_{x_l}^{x_u} dx\,\int_{y_l}^{y_u}dy\,... f(x,y,...)
$$
\afterdisplay
@end tex
@ifinfo
@example
I = \int_xl^xu dx \int_yl^yu  dy ...  f(x, y, ...)
@end example
@end ifinfo
@noindent
over a hypercubic region @math{((x_l,x_u)}, @math{(y_l,y_u), ...)} using
a fixed number of function calls.  The routines also provide a
statistical estimate of the error on the result.  This error estimate
should be taken as a guide rather than as a strict error bound ---
random sampling of the region may not uncover all the important features
of the function, resulting in an underestimate of the error.

The functions are defined in separate header files for each routine,
@code{gsl_monte_plain.h}, @file{gsl_monte_miser.h} and
@file{gsl_monte_vegas.h}.

@menu
* Monte Carlo Interface::       
* PLAIN Monte Carlo::  
* MISER::                       
* VEGAS::                       
* Monte Carlo Examples::        
* Monte Carlo Integration References and Further Reading::  
@end menu

@node Monte Carlo Interface
@section Interface
All of the Monte Carlo integration routines use the same interface.
There is an allocator to allocate memory for control variables and
workspace, a routine to initialize those control variables, the
integrator itself, and a function to free the space when done.

Each integration function requires a random number generator to be
supplied, and returns an estimate of the integral and its standard
deviation.  The accuracy of the result is determined by the number of
function calls specified by the user.  If a known level of accuracy is
required this can be achieved by calling the integrator several times
and averaging the individual results until the desired accuracy is
obtained.  

Random sample points used within the Monte Carlo routines are always
chosen strictly within the integration region, so that endpoint
singularities are automatically avoided.

The function to be integrated has its own datatype, defined in the
header file @file{gsl_monte.h}.

@deftp {Data Type} gsl_monte_function 

This data type defines a general function with parameters for Monte
Carlo integration.

@table @code
@item double (* @var{function}) (double * @var{x}, size_t @var{dim}, void * @var{params})
this function should return the value
@c{$f(x,\hbox{\it params})$}
@math{f(x,params)} for argument @var{x} and parameters @var{params},
where @var{x} is an array of size @var{dim} giving the coordinates of
the point where the function is to be evaluated.

@item size_t @var{dim}
the number of dimensions for @var{x}

@item void * @var{params}
a pointer to the parameters of the function
@end table
@end deftp
@noindent
Here is an example for a quadratic function in two dimensions,

@tex
\beforedisplay
$$
f(x,y) = a x^2 + b x y + c y^2
$$
\afterdisplay
@end tex
@ifinfo
@example
f(x,y) = a x^2 + b x y + c y^2
@end example
@end ifinfo
@noindent
with @math{a = 3}, @math{b = 2}, @math{c = 1}.  The following code
defines a @code{gsl_monte_function} @code{F} which you could pass to an
integrator:

@example
struct my_f_params @{ double a; double b; double c; @};

double
my_f (double x, size_t dim, void * p) @{
   struct my_f_params * fp = (struct my_f_params *)p;

   if (dim != 2)
      @{
        fprintf(stderr, "error: dim != 2");
        abort();
      @}

   return  fp->a * x[0] * x[0] 
             + fp->b * x[0] * x[1] 
               + fp->c * x[1] * x[1];
@}

gsl_monte_function F;
struct my_f_params params = @{ 3.0, 2.0, 1.0 @};

F.function = &my_f;
F.dim = 2;
F.params = &params;
@end example
@noindent
The function @math{f(x)} can be evaluated using the following macro,

@example
#define GSL_MONTE_FN_EVAL(F,x) 
    (*((F)->function))(x,(F)->dim,(F)->params)
@end example

@node PLAIN Monte Carlo
@section PLAIN Monte Carlo
@cindex plain monte carlo
The plain Monte Carlo algorithm samples points randomly from the
integration region to estimate the integral and its error.  Using this
algorithm the estimate of the integral @math{E(f; N)} for @math{N}
randomly distributed points @math{x_i} is given by,

@tex
\beforedisplay
$$
E(f; N) =  V \langle f \rangle = {V \over N} \sum_i^N f(x_i).
$$
\afterdisplay
@end tex
@ifinfo
@example
E(f; N) = =  V <f> = (V / N) \sum_i^N f(x_i).
@end example
@end ifinfo
@noindent
where @math{V} is the volume of the integration region.  The error on
this estimate @math{\sigma(E;N)} is calculated from the estimated
variance of the mean,

@tex
\beforedisplay
$$
\sigma^2 (E; N) = {V \over N } \sum_i^N (f(x_i) - \langle f \rangle)^2
$$
\afterdisplay
@end tex
@ifinfo
@example
\sigma^2 (E; N) = (V / N) \sum_i^N (f(x_i) -  <f>)^2
@end example
@end ifinfo
@noindent
For large @math{N} this variance decreases asymptotically as
@math{var(f)/N}, where @math{var(f)} is the true variance of the
function over the integration region.  The error estimate itself should
decrease as @c{$\sigma(f)/\sqrt{N}$}
@math{\sigma(f)/\sqrt@{N@}}.  The familiar law of errors
decreasing as @c{$1/\sqrt{N}$}
@math{1/\sqrt@{N@}} applies --- to reduce the error by a
factor of 10 requires a 100-fold increase in the number of sample
points.

The functions described in this section are declared in the header file
@file{gsl_monte_plain.h}.

@deftypefun {gsl_monte_plain_state *} gsl_monte_plain_alloc (size_t @var{dim})
This function allocates and initializes a workspace for Monte Carlo
integration in @var{dim} dimensions.  
@end deftypefun

@deftypefun int gsl_monte_plain_init (gsl_monte_plain_state* @var{s})
This function initializes a previously allocated integration state.
This allows an existing workspace to be reused for different
integrations.
@end deftypefun

@deftypefun int gsl_monte_plain_integrate (gsl_monte_function * @var{f}, double * @var{xl}, double * @var{xu}, size_t @var{dim}, size_t @var{calls}, gsl_rng * @var{r}, gsl_monte_plain_state * @var{s}, double * @var{result}, double * @var{abserr})
This routines uses the plain Monte Carlo algorithm to integrate the
function @var{f} over the @var{dim}-dimensional hypercubic region
defined by the lower and upper limits in the arrays @var{xl} and
@var{xu}, each of size @var{dim}.  The integration uses a fixed number
of function calls @var{calls}, and obtains random sampling points using
the random number generator @var{r}. A previously allocated workspace
@var{s} must be supplied.  The result of the integration is returned in
@var{result}, with an estimated absolute error @var{abserr}.
@end deftypefun

@deftypefun void gsl_monte_plain_free (gsl_monte_plain_state* @var{s}), 
This function frees the memory associated with the integrator state
@var{s}.
@end deftypefun

@node MISER
@section MISER
@cindex MISER monte carlo integration
@cindex recursive stratified sampling, MISER

The @sc{miser} algorithm of Press and Farrar is based on recursive
stratified sampling.  This technique aims to reduce the overall
integration error by concentrating integration points in the regions of
highest variance.

The idea of stratified sampling begins with the observation that for two
disjoint regions @math{a} and @math{b} with Monte Carlo estimates of the
integral @math{E_a(f)} and @math{E_b(f)} and variances
@math{\sigma_a^2(f)} and @math{\sigma_b^2(f)}, the variance
@math{Var(f)} of the combined estimate 
@c{$E(f) = {1\over 2} (E_a(f) + E_b(f))$} 
@math{E(f) = (1/2) (E_a(f) + E_b(f))} 
is given by,

@tex
\beforedisplay
$$
Var(f) = {\sigma_a^2(f) \over 4 N_a} + {\sigma_b^2(f) \over 4 N_b}
$$
\afterdisplay
@end tex
@ifinfo
@example
Var(f) = (\sigma_a^2(f) / 4 N_a) + (\sigma_b^2(f) / 4 N_b)
@end example
@end ifinfo
@noindent
It can be shown that this variance is minimized by distributing the
points such that,

@tex
\beforedisplay
$$
{N_a \over N_a+N_b} = {\sigma_a \over \sigma_a + \sigma_b}
$$
\afterdisplay
@end tex
@ifinfo
@example
N_a / (N_a + N_b) = \sigma_a / (\sigma_a + \sigma_b)
@end example
@end ifinfo
@noindent
Hence the smallest error estimate is obtained by allocating sample
points in proportion to the standard deviation of the function in each
sub-region.

The @sc{miser} algorithm proceeds by bisecting the integration region
along one coordinate axis to give two sub-regions at each step.  The
direction is chosen by examining all @math{d} possible bisections and
selecting the one which will minimize the combined variance of the two
sub-regions.  The variance in the sub-regions is estimated by sampling
with a fraction of the total number of points available to the current
step.  The same procedure is then repeated recursively for each of the
two half-spaces from the best bisection. The remaining sample points are
allocated to the sub-regions using the formula for @math{N_a} and
@math{N_b}.  This recursive allocation of integration points continues
down to a user-specified depth where each sub-region is integrated using
a plain Monte Carlo estimate.  These individual values and their error
estimates are then combined upwards to give an overall result and an
estimate of its error.

The functions described in this section are declared in the header file
@file{gsl_monte_miser.h}.

@deftypefun {gsl_monte_miser_state *} gsl_monte_miser_alloc (size_t @var{dim})
This function allocates and initializes a workspace for Monte Carlo
integration in @var{dim} dimensions.  The workspace is used to maintain
the state of the integration.
@end deftypefun

@deftypefun int gsl_monte_miser_init (gsl_monte_miser_state* @var{s})
This function initializes a previously allocated integration state.
This allows an existing workspace to be reused for different
integrations.
@end deftypefun

@deftypefun int gsl_monte_miser_integrate (gsl_monte_function * @var{f}, double * @var{xl}, double * @var{xu}, size_t @var{dim}, size_t @var{calls}, gsl_rng * @var{r}, gsl_monte_miser_state * @var{s}, double * @var{result}, double * @var{abserr})
This routines uses the @sc{miser} Monte Carlo algorithm to integrate the
function @var{f} over the @var{dim}-dimensional hypercubic region
defined by the lower and upper limits in the arrays @var{xl} and
@var{xu}, each of size @var{dim}.  The integration uses a fixed number
of function calls @var{calls}, and obtains random sampling points using
the random number generator @var{r}. A previously allocated workspace
@var{s} must be supplied.  The result of the integration is returned in
@var{result}, with an estimated absolute error @var{abserr}.
@end deftypefun

@deftypefun void gsl_monte_miser_free (gsl_monte_miser_state* @var{s}), 
This function frees the memory associated with the integrator state
@var{s}.
@end deftypefun

The @sc{miser} algorithm has several configurable parameters. The
following variables can be accessed through the
@code{gsl_monte_miser_state} struct,

@deftypevar double estimate_frac
This parameter specifies the fraction of the currently available number of
function calls which are allocated to estimating the variance at each
recursive step. The default value is 0.1.
@end deftypevar

@deftypevar size_t min_calls
This parameter specifies the minimum number of function calls required
for each estimate of the variance. If the number of function calls
allocated to the estimate using @var{estimate_frac} falls below
@var{min_calls} then @var{min_calls} are used instead.  This ensures
that each estimate maintains a reasonable level of accuracy.  The
default value of @var{min_calls} is @code{16 * dim}.
@end deftypevar

@deftypevar size_t min_calls_per_bisection
This parameter specifies the minimum number of function calls required
to proceed with a bisection step.  When a recursive step has fewer calls
available than @var{min_calls_per_bisection} it performs a plain Monte
Carlo estimate of the current sub-region and terminates its branch of
the recursion.  The default value of this parameter is @code{32 *
min_calls}.
@end deftypevar

@deftypevar double alpha
This parameter controls how the estimated variances for the two
sub-regions of a bisection are combined when allocating points.  With
recursive sampling the overall variance should scale better than
@math{1/N}, since the values from the sub-regions will be obtained using
a procedure which explicitly minimizes their variance.  To accommodate
this behavior the @sc{miser} algorithm allows the total variance to
depend on a scaling parameter @math{\alpha},

@tex
\beforedisplay
$$
Var(f) = {\sigma_a \over N_a^\alpha} + {\sigma_b \over N_b^\alpha}
$$
\afterdisplay
@end tex
@ifinfo
@example
Var(f) = @{\sigma_a \over N_a^\alpha@} + @{\sigma_b \over N_b^\alpha@}
@end example
@end ifinfo
@noindent
The authors of the original paper describing @sc{miser} recommend the
value @math{\alpha = 2} as a good choice, obtained from numerical
experiments, and this is used as the default value in this
implementation.
@end deftypevar

@deftypevar double dither
This parameter introduces a random fractional variation of size
@var{dither} into each bisection, which can be used to break the
symmetry of integrands which are concentrated near the exact center of
the hypercubic integration region.  The default value of dither is zero,
so no variation is introduced. If needed, a typical value of
@var{dither} is around 0.1.
@end deftypevar

@node VEGAS
@section VEGAS
@cindex VEGAS monte carlo integration
@cindex importance sampling, VEGAS

The @sc{vegas} algorithm of Lepage is based on importance sampling.  It
samples points from the probability distribution described by the
function @math{|f|}, so that the points are concentrated in the regions
that make the largest contribution to the integral.

In general, if the Monte Carlo integral of @math{f} is sampled with
points distributed according to a probability distribution described by
the function @math{g}, we obtain an estimate @math{E_g(f; N)},

@tex
\beforedisplay
$$
E_g(f; N) = E(f/g; N)
$$
\afterdisplay
@end tex
@ifinfo
@example
E_g(f; N) = E(f/g; N)