This chapter describes the Monte Carlo integration routines in the library.  At the moment the algorithms implemented are: plain non-adaptive, Monte Carlo, VEGAS (following Peter Lepage) and MISER (following Numerical recipes).  We will dispense with an introduction and simply describe the algorithms and interfaces.  Note that the documentation (as well as the code) is still evolving - the code at a faster rate.

Algorithms
==========

PLAIN (or Simple) Monte Carlo
-----------------------------

We begin by establishing some notation.  Let I(f) denote the integral of the function f (for convenience, we take I(1)=1) We will write the Monte Carlo estimate of the integral as E(f; N) for the

     E(f; N) = {1 \over N } \sum_i^N f(x_i).

for N randomly distributed points $x_i$ x_i Similarly, we will denote the variance of f by

     \sigma^2(f) = I(f^2) - I(f)^2

and the variance of the estimate by

     Var(f; N) = E(f^2; N) - (E(f; N))^2.

The fundamental point of Monte Carlo integration is that for large N,

     E(f; N) - I(f) \approx {\sigma^2(f)\over N }

and furthermore, by the same argument,

     \sigma^2(f) \approx Var(f; N)

That's all there is to simple Monte Carlo.

MISER
-----

The starting point of all stratified sampling techniques is the observation that for two disjoint regions a and b with Monte Carlo estimates of the integral E_a(f) and E_b(f) of the integral of f in those regions and variances \sigma_a^2(f) and \sigma_b^2(f), the variance Var(f) of the combined estimate $E(f) = 1\over 2 E_a(f) + 1\over 2 E_b(f)$ E(f) = {1\over 2} E_a(f) + {1\over 2} E_b(f) is given by

     Var(f) = {\sigma_a^2(f) \over 4 N_a} + {\sigma_b^2(f) \over 4 N_b}

This variance is minimized (subject to the constrain that N_a + N_b is fixed) by choosing (assuming one can) N_a s.t.

     {N_a \over N_a+N_b } = { \sigma_a \over \sigma_a + \sigma_b }

For such a choice, the variance is given by

     Var(f) =  { (\sigma_a + \sigma_b)^2 \over 4 N }

In words, the variance of the estimate (and hence the error) is minimized by concentrating points in regions where the variance of the function f is large.

The most straightforward stratified sampling routine divides a hyper-cubic integration region into sub-cubes along the coordinate axes. This is simple to do, but gets out of hand for large numbers of dimensions d because the number of sub-cubes grows like K^d (where K is the number of sub-divisions along the axis).

What MISER does (and other recursive stratified samplers) is to bisect the hypercube along one coordinate axis.  The direction is chosen by examining all possible bisections (d of them) and picking the one that gives the best combined variance.  The same procedure is then done for each of the two half-spaces (choosing N_a and N_b as described above), and so on.  Since by assumption one does not know the variance Var_L(f) in the various sub-regions, it is estimated using some fraction of the total number of points alloted to the estimate.

VEGAS
-----

Vegas takes another approach to obtaining good results, namely it tries to sample points from the probability distribution described by the function f.  More precisely, suppose we estimate the integral of f with points distributed according to a probability distribution described by the function g.  If we call the new estimate E_g(f; N) we have

     E_g(f; N) = E(f/g; N)

and similarly

     Var_g(f; N) = Var(f/g; N)

It is clear from this that if g = |f|/I(|f|) then the variance V_g(f; N) vanishes. This it was what Vegas attempts to do.  It makes several passes, histograming f, and each time using the histogram to define the sampling distribution for the next pass.  This basic idea again has the problem that the number of histogram bins grows like $K^d$ K^d.  The compromise that Vegas (ie, Peter Lepage) adopts is to assume that g factors: $g(x_1, x_2, \ldots) = g_1(x_1) g_2(x_2) \ldots$ g(x_1, x_2, ...) = g_1(x_1) g_2(x_2) ... so that the number of bins required is only $d \cdot K$ d.K.  In this case, it is not hard to show that the optimal distribution is

     g_1(x_1) = [\int dx_1 ... dx_n {f^2(x_1,...,x_n)\over g_2(x_2)... g_n(x_n) } ]^{(1/2)}

Lepage's Vegas is actually a bit fancier than this, combining 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).

Interface
=========

All of the integration routines use the same interface.  There is an allocator to allocate memory to hold control variables and workspace, a routine to initialize those control variables, the integrator itself, and of course a function to free the space when done.  For an integrator algorithm, call it COOL (substitute any of the currently existing algorithms) we then have
 -- Function: gsl_monte_COOL_state* gsl_monte_COOL_alloc (size_t DIM)
 -- Function: int gsl_monte_COOL_init (gsl_monte_COOL_state* S)
 -- Function: int gsl_monte_COOL_free (gsl_monte_COOL_state* S),
 -- Function: int gsl_monte_COOL_validate (gsl_monte_COOL_state* S,
     double *X_LOWER, double *X_UPPER, unsigned long DIMENSION, unsigned long FUNCTION_CALLS)
 -- Function: int gsl_monte_COOL_integrate (gsl_monte_COOL_state* S,
     double* X_LOWER, double* X_UPPER, unsigned long DIMENSION, unsigned long FUNCTION_CALLS, double* RESULT, double* ERROR, ...)

Notice the ellipses in the last argument to the actual integration routine. This is because the vegas algorithm (and perhaps others in the future) returns extra information -- in the case of vegas, it is the \chi^2 of the result.

In addition to the common function interface, the routines also share the state variables
 -- Control: gsl_rng* ranf
     which determines which random number generator will be used.

and
 -- Control: int verbose
     which says whether to print information about the calculation (though the actual use depends on the algorithm).

In the near future, we expect that there will be a common interface for selecting a ``log'' stream for reporting various information about the performance of the particular algorithm.

We describe the algorithm-specific variables.

PLAIN: None

MISER:
 -- Control: double alpha
     alpha controls how the variances for the two sub-regions are combined. The Numerical Recipes gang argue that for recursive sampling there is no reason to expect that the variance should scale like 1/N and so they allow the scaling to depend on \alpha

          Var(f) = {\sigma_a \over N_a^\alpha} + {\sigma_b \over N_b^\alpha}

 -- Control: double dither
     Rather than exactly bisection the integration region, dither allows the user to introduce a bit of fuzz.  This helps in the case when the function to be integrated has some symmetry, say if it is peaked in the center of the hypercube.  If needed, dither should be around 0.1.

VEGAS:
 -- Control: double alpha
     For Vegas, alpha controls the stiffness of the rebinning algorithm: alpha = 0 means never rebin.  It is typically set between one and two.

 -- Control: double acc
     Setting acc allows vegas to terminate when the desired accuracy has been achieved, rather than after a certain number of function calls. Setting it to a negative value disables this feature.

 -- Control: long int max_it_num
     The maximum number of iterations to perform.

 -- Control: int stage
     Setting this determines the "stage" of the calculation.  Normally, stage = 0.  Calling vegas with stage = 1 retains the grid (but not the answers) from the previous run, so that one can ``tune'' the grid using a relatively small number of points and then do a large run with stage = 1 on the optimized grid.  Setting stage = 2 keeps the grid and the answers from the previous run and stage = 3 enters at the main loop, so that nothing is changed -- this is like deciding to change max_it_num during the run.

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

Example
=======

Here we provide s simple example of using the integration routines.  Vegas was chosen at random.

     #include <math.h>
     #include <stdio.h>

     #include <gsl_math.h>
     #include <gsl_monte_vegas.h>

     double f0(double x[])
     {
       int i;
       double product = 1.0;

       for (i = 0; i < 10; i++)
         product *= x[i];

       return product;
     }

     main ()
     {
       double xl[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
       double xu[10] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1};

       double res = 0;
       double err = 0;
       double chisq = 0;
       int status = 0;
       unsigned long calls = 10000;
       unsigned long dimension = 10;

       gsl_monte_vegas_state* s = gsl_monte_vegas_alloc(10);
       gsl_monte_vegas_init(s);

       s->alpha = 1.5; /* default */
       s->verbose = 0; /* default, we prefer to remain ignorant! */
       s->stage = 0; /* default.  Start at the beginning */
       s->acc = -1.0; /* default. don't terminate when accuracy is reached */
       s->max_it_num = 5; /* default.  Do five iterations before quitting */

       status = gsl_monte_vegas_integrate(s, f0, xl, xu, dimension, calls,
                                         &res, &err, &chisq);
       if (status)
         printf("oops!\n");
       else
         printf("vegas(f0) = %f +- %f with \"chisq\" %f\n", res, err, chisq);

     }

The Future
==========

In the future, the author of the Monte Carlo routines intends to add more algorithms, greater control over the current ones,  more error handling, and a more consistent interface.  Probably something like the rng interface will evolve, so that something like

       integrator = gsl_monte_alloc(gsl_monte_vegas, dimension)

will be possible.  The old interface will probably stay around for a while (though, since this is version < 1, the author does not want to be held to such as statement).