montecarlo.texi

开放gsl矩阵运算

@end example
@end ifinfo
@noindent
with a corresponding variance,

@tex
\beforedisplay
$$
Var_g(f; N) = Var(f/g; N)
$$
\afterdisplay
@end tex
@ifinfo
@example
Var_g(f; N) = Var(f/g; N)
@end example
@end ifinfo
@noindent
If the probability distribution is chosen as @math{g = |f|/I(|f|)} then
it can be shown that the variance @math{V_g(f; N)} vanishes, and the
error in the estimate will be zero.  In practice it is not possible to
sample from the exact distribution @math{g} for an arbitrary function, so
importance sampling algorithms aim to produce efficient approximations
to the desired distribution.

The @sc{vegas} algorithm approximates the exact distribution by making a
number of passes over the integration region while histogramming the
function @math{f}. Each histogram is used to define a sampling
distribution for the next pass.  Asymptotically this procedure converges
to the desired distribution. In order
to avoid the number of histogram bins growing like @math{K^d} the
probability distribution is approximated by a separable function:
@c{$g(x_1, x_2,\ldots) = g_1(x_1) g_2(x_2)\ldots$} 
@math{g(x_1, x_2, ...) = g_1(x_1) g_2(x_2) ...}  
so that the number of bins required is only @math{Kd}.     
This is equivalent to locating the peaks of the function from the
projections of the integrand onto the coordinate axes.  The efficiency
of @sc{vegas} depends on the validity of this assumption.  It is most
efficient when the peaks of the integrand are well-localized.  If an
integrand can be rewritten in a form which is approximately separable
this will increase the efficiency of integration with @sc{vegas}.

@sc{vegas} incorporates a number of additional features, and combines both
stratified sampling and importance sampling.  The integration region is
divided into a number of ``boxes'', with each box getting in fixed
number of points (the goal is 2).  Each box can then have a fractional
number of bins, but if bins/box is less than two, Vegas switches to a
kind variance reduction (rather than importance sampling).


@deftypefun {gsl_monte_vegas_state *} gsl_monte_vegas_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_vegas_init (gsl_monte_vegas_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_vegas_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_vegas_state * @var{s}, double * @var{result}, double * @var{abserr})
This routines uses the @sc{vegas} 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}.  The result
and its error estimate are based on a weighted average of independent
samples. The chi-squared per degree of freedom for the weighted average
is returned via the state struct component, @var{s->chisq}, and must be
consistent with 1 for the weighted average to be reliable.
@end deftypefun

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

The @sc{vegas} algorithm computes a number of independent estimates of the
integral internally, according to the @code{iterations} parameter
described below, and returns their weighted average.  Random sampling of
the integrand can occasionally produce an estimate where the error is
zero, particularly if the function is constant in some regions. An
estimate with zero error causes the weighted average to break down and
must be handled separately. In the original Fortran implementations of
@sc{vegas} the error estimate is made non-zero by substituting a small
value (typically @code{1e-30}).  The implementation in GSL differs from
this and avoids the use of an arbitrary constant -- it either assigns
the value a weight which is the average weight of the preceding
estimates or discards it according to the following procedure,

@table @asis
@item current estimate has zero error, weighted average has finite error

The current estimate is assigned a weight which is the average weight of
the preceding estimates.

@item current estimate has finite error, previous estimates had zero error

The previous estimates are discarded and the weighted averaging
procedure begins with the current estimate.

@item current estimate has zero error, previous estimates had zero error

The estimates are averaged using the arithmetic mean, but no error is computed.
@end table

The @sc{vegas} algorithm is highly configurable. The following variables
can be accessed through the @code{gsl_monte_vegas_state} struct,

@deftypevar double result
@deftypevarx double sigma
These parameters contain the raw value of the integral @var{result} and
its error @var{sigma} from the last iteration of the algorithm.
@end deftypevar

@deftypevar double chisq
This parameter gives the chi-squared per degree of freedom for the
weighted estimate of the integral.  The value of @var{chisq} should be
close to 1.  A value of @var{chisq} which differs significantly from 1
indicates that the values from different iterations are inconsistent.
In this case the weighted error will be under-estimated, and further
iterations of the algorithm are needed to obtain reliable results.
@end deftypevar

@deftypevar double alpha
The parameter @code{alpha} controls the stiffness of the rebinning
algorithm.  It is typically set between one and two. A value of zero
prevents rebinning of the grid.  The default value is 1.5.
@end deftypevar

@deftypevar size_t iterations
The number of iterations to perform for each call to the routine. The
default value is 5 iterations.
@end deftypevar

@deftypevar int stage
Setting this determines the @dfn{stage} of the calculation.  Normally,
@code{stage = 0} which begins with a new uniform grid and empty weighted
average.  Calling vegas with @code{stage = 1} retains the grid from the
previous run but discards the weighted average, so that one can ``tune''
the grid using a relatively small number of points and then do a large
run with @code{stage = 1} on the optimized grid.  Setting @code{stage =
2} keeps the grid and the weighted average from the previous run, but
may increase (or decrease) the number of histogram bins in the grid
depending on the number of calls available.  Choosing @code{stage = 3}
enters at the main loop, so that nothing is changed, and is equivalent
to performing additional iterations in a previous call.
@end deftypevar

@deftypevar int mode
The possible choices are @code{GSL_VEGAS_MODE_IMPORTANCE},
@code{GSL_VEGAS_MODE_STRATIFIED}, @code{GSL_VEGAS_MODE_IMPORTANCE_ONLY}.
This determines whether @sc{vegas} will use importance sampling or
stratified sampling, or whether it can pick on its own.  In low
dimensions @sc{vegas} uses strict stratified sampling (more precisely,
stratified sampling is chosen if there are fewer than 2 bins per box).
@end deftypevar

@deftypevar int verbose
@deftypevarx {FILE *} ostream
These parameters set the level of information printed by @sc{vegas}. All
information is written to the stream @var{ostream}.  The default setting
of @var{verbose} is @code{-1}, which turns off all output.  A
@var{verbose} value of @code{0} prints summary information about the
weighted average and final result, while a value of @code{1} also
displays the grid coordinates.  A value of @code{2} prints information
from the rebinning procedure for each iteration.
@end deftypevar

@node Monte Carlo Examples
@section Examples

The example program below uses the Monte Carlo routines to estimate the
value of the following 3-dimensional integral from the theory of random
walks,

@tex
\beforedisplay
$$
I = \int_{-\pi}^{+\pi} {dk_x \over 2\pi} 
    \int_{-\pi}^{+\pi} {dk_y \over 2\pi} 
    \int_{-\pi}^{+\pi} {dk_z \over 2\pi} 
     { 1 \over (1 - \cos(k_x)\cos(k_y)\cos(k_z))}
$$
\afterdisplay
@end tex
@ifinfo
@example
I = \int_@{-pi@}^@{+pi@} @{dk_x/(2 pi)@} 
    \int_@{-pi@}^@{+pi@} @{dk_y/(2 pi)@} 
    \int_@{-pi@}^@{+pi@} @{dk_z/(2 pi)@} 
     1 / (1 - cos(k_x)cos(k_y)cos(k_z))
@end example
@end ifinfo
@noindent
The analytic value of this integral can be shown to be @math{I =
\Gamma(1/4)^4/(4 \pi^3) = 1.393203929685676859...}.  The integral gives
the mean time spent at the origin by a random walk on a body-centered
cubic lattice in three dimensions.

For simplicity we will compute the integral over the region
@math{(0,0,0)} to @math{(\pi,\pi,\pi)} and multiply by 8 to obtain the
full result.  The integral is slowly varying in the middle of the region
but has integrable singularities at the corners @math{(0,0,0)},
@math{(0,\pi,\pi)}, @math{(\pi,0,\pi)} and @math{(\pi,\pi,0)}.  The
Monte Carlo routines only select points which are strictly within the
integration region and so no special measures are needed to avoid these
singularities.

@smallexample
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_vegas.h>

/* Computation of the integral,

      I = int (dx dy dz)/(2pi)^3  1/(1-cos(x)cos(y)cos(z))

   over (-pi,-pi,-pi) to (+pi, +pi, +pi).  The exact answer
   is Gamma(1/4)^4/(4 pi^3).  This example is taken from
   C.Itzykson, J.M.Drouffe, "Statistical Field Theory -
   Volume 1", Section 1.1, p21, which cites the original
   paper M.L.Glasser, I.J.Zucker, Proc.Natl.Acad.Sci.USA 74
   1800 (1977) */

/* For simplicity we compute the integral over the region 
   (0,0,0) -> (pi,pi,pi) and multiply by 8 */

double exact = 1.3932039296856768591842462603255;

double
g (double *k, size_t dim, void *params)
@{
  double A = 1.0 / (M_PI * M_PI * M_PI);
  return A / (1.0 - cos (k[0]) * cos (k[1]) * cos (k[2]));
@}

void
display_results (char *title, double result, double error)
@{
  printf ("%s ==================\n", title);
  printf ("result = % .6f\n", result);
  printf ("sigma  = % .6f\n", error);
  printf ("exact  = % .6f\n", exact);
  printf ("error  = % .6f = %.1g sigma\n", result - exact,
          fabs (result - exact) / error);
@}

int
main (void)
@{
  double res, err;

  double xl[3] = @{ 0, 0, 0 @};
  double xu[3] = @{ M_PI, M_PI, M_PI @};

  const gsl_rng_type *T;
  gsl_rng *r;

  gsl_monte_function G = @{ &g, 3, 0 @};

  size_t calls = 500000;

  gsl_rng_env_setup ();

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  @{
    gsl_monte_plain_state *s = gsl_monte_plain_alloc (3);
    gsl_monte_plain_integrate (&G, xl, xu, 3, calls, r, s, 
                               &res, &err);
    gsl_monte_plain_free (s);

    display_results ("plain", res, err);
  @}

  @{
    gsl_monte_miser_state *s = gsl_monte_miser_alloc (3);
    gsl_monte_miser_integrate (&G, xl, xu, 3, calls, r, s,
                               &res, &err);
    gsl_monte_miser_free (s);

    display_results ("miser", res, err);
  @}

  @{
    gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (3);

    gsl_monte_vegas_integrate (&G, xl, xu, 3, 10000, r, s,
                               &res, &err);
    display_results ("vegas warm-up", res, err);

    printf ("converging...\n");

    do
      @{
        gsl_monte_vegas_integrate (&G, xl, xu, 3, calls/5, r, s,
                                   &res, &err);
        printf ("result = % .6f sigma = % .6f "
                "chisq/dof = %.1f\n", res, err, s->chisq);
      @}
    while (fabs (s->chisq - 1.0) > 0.5);

    display_results ("vegas final", res, err);

    gsl_monte_vegas_free (s);
  @}
  return 0;
@}
@end smallexample
@noindent
With 500,000 function calls the plain Monte Carlo algorithm achieves a
fractional error of 0.6%.  The estimated error @code{sigma} is
consistent with the actual error, and the computed result differs from
the true result by about one standard deviation,

@example
plain ==================
result =  1.385867
sigma  =  0.007938
exact  =  1.393204
error  = -0.007337 = 0.9 sigma
@end example
@noindent
The @sc{miser} algorithm reduces the error by a factor of two, and also
correctly estimates the error,

@example
miser ==================
result =  1.390656
sigma  =  0.003743
exact  =  1.393204
error  = -0.002548 = 0.7 sigma
@end example
@noindent
In the case of the @sc{vegas} algorithm the program uses an initial
warm-up run of 10,000 function calls to prepare, or "warm up", the grid.
This is followed by a main run with five iterations of 100,000 function
calls. The chi-squared per degree of freedom for the five iterations are
checked for consistency with 1, and the run is repeated if the results
have not converged. In this case the estimates are consistent on the
first pass.

@example
vegas warm-up ==================
result =  1.386925
sigma  =  0.002651
exact  =  1.393204
error  = -0.006278 = 2 sigma
converging...
result =  1.392957 sigma =  0.000452 chisq/dof = 1.1
vegas final ==================
result =  1.392957
sigma  =  0.000452
exact  =  1.393204
error  = -0.000247 = 0.5 sigma
@end example
@noindent
If the value of @code{chisq} had differed significantly from 1 it would
indicate inconsistent results, with a correspondingly underestimated
error.  The final estimate from @sc{vegas} (using a similar number of
function calls) is significantly more accurate than the other two
algorithms.

@node Monte Carlo Integration References and Further Reading
@section References and Further Reading

@noindent
The @sc{miser} algorithm is described in the following article,

@itemize @asis
@item
W.H. Press, G.R. Farrar, @cite{Recursive Stratified Sampling for
Multidimensional Monte Carlo Integration},
Computers in Physics, v4 (1990), pp190-195.
@end itemize
@noindent
The @sc{vegas} algorithm is described in the following papers,

@itemize @asis
@item
G.P. Lepage,
@cite{A New Algorithm for Adaptive Multidimensional Integration},
Journal of Computational Physics 27, 192-203, (1978)

@item
G.P. Lepage,
@cite{VEGAS: An Adaptive Multi-dimensional Integration Program},
Cornell preprint CLNS 80-447, March 1980
@end itemize