
Go to the first, previous, next, last section, table of contents.
Monte Carlo IntegrationThis chapter describes routines for multidimensional Monte Carlo integration. These include the traditional Monte Carlo method and adaptive algorithms such as VEGAS and MISER which use importance sampling and stratified sampling techniques. Each algorithm computes an estimate of a multidimensional definite integral of the form, I = \int_xl^xu dx \int_yl^yu dy ... f(x, y, ...) over a hypercubic region ((x_l,x_u), (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,
InterfaceAll 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 `gsl_monte.h'.
Here is an example for a quadratic function in two dimensions, f(x,y) = a x^2 + b x y + c y^2
with a = 3, b = 2, c = 1. The following code
defines a 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 = ¶ms; The function f(x) can be evaluated using the following macro, #define GSL_MONTE_FN_EVAL(F,x) (*((F)>function))(x,(F)>dim,(F)>params) PLAIN Monte CarloThe 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 E(f; N) for N randomly distributed points x_i is given by, E(f; N) = = V <f> = (V / N) \sum_i^N f(x_i). where V is the volume of the integration region. The error on this estimate \sigma(E;N) is calculated from the estimated variance of the mean, \sigma^2 (E; N) = (V / N) \sum_i^N (f(x_i)  <f>)^2 For large N this variance decreases asymptotically as var(f)/N, where var(f) is the true variance of the function over the integration region. The error estimate itself should decrease as \sigma(f)/\sqrt{N}. The familiar law of errors decreasing as 1/\sqrt{N} applies  to reduce the error by a factor of 10 requires a 100fold increase in the number of sample points. The functions described in this section are declared in the header file `gsl_monte_plain.h'.
MISERThe 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 a and b with Monte Carlo estimates of the integral E_a(f) and E_b(f) and variances \sigma_a^2(f) and \sigma_b^2(f), the variance Var(f) of the combined estimate E(f) = (1/2) (E_a(f) + E_b(f)) is given by, Var(f) = (\sigma_a^2(f) / 4 N_a) + (\sigma_b^2(f) / 4 N_b) It can be shown that this variance is minimized by distributing the points such that, N_a / (N_a + N_b) = \sigma_a / (\sigma_a + \sigma_b) Hence the smallest error estimate is obtained by allocating sample points in proportion to the standard deviation of the function in each subregion. The MISER algorithm proceeds by bisecting the integration region along one coordinate axis to give two subregions at each step. The direction is chosen by examining all d possible bisections and selecting the one which will minimize the combined variance of the two subregions. The variance in the subregions 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 halfspaces from the best bisection. The remaining sample points are allocated to the subregions using the formula for N_a and N_b. This recursive allocation of integration points continues down to a userspecified depth where each subregion 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 `gsl_monte_miser.h'.
The MISER algorithm has several configurable parameters. The
following variables can be accessed through the
VEGASThe VEGAS algorithm of Lepage is based on importance sampling. It samples points from the probability distribution described by the function 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 f is sampled with points distributed according to a probability distribution described by the function g, we obtain an estimate E_g(f; N), E_g(f; N) = E(f/g; N) with a corresponding variance, Var_g(f; N) = Var(f/g; N) If the probability distribution is chosen as g = f/I(f) then it can be shown that the variance 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 g for an arbitrary function, so importance sampling algorithms aim to produce efficient approximations to the desired distribution. The VEGAS algorithm approximates the exact distribution by making a number of passes over the integration region while histogramming the function 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 K^d the probability distribution is approximated by a separable function: g(x_1, x_2, ...) = g_1(x_1) g_2(x_2) ... so that the number of bins required is only Kd. This is equivalent to locating the peaks of the function from the projections of the integrand onto the coordinate axes. The efficiency of VEGAS depends on the validity of this assumption. It is most efficient when the peaks of the integrand are welllocalized. If an integrand can be rewritten in a form which is approximately separable this will increase the efficiency of integration with VEGAS. 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).
The VEGAS algorithm computes a number of independent estimates of the
integral internally, according to the
The VEGAS algorithm is highly configurable. The following variables
can be accessed through the
ExamplesThe example program below uses the Monte Carlo routines to estimate the value of the following 3dimensional integral from the theory of random walks, 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)) The analytic value of this integral can be shown to be 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 bodycentered cubic lattice in three dimensions. For simplicity we will compute the integral over the region (0,0,0) to (\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 (0,0,0), (0,\pi,\pi), (\pi,0,\pi) and (\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. #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/(1cos(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 warmup", 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; }
With 500,000 function calls the plain Monte Carlo algorithm achieves a
fractional error of 0.6%. The estimated error plain ================== result = 1.385867 sigma = 0.007938 exact = 1.393204 error = 0.007337 = 0.9 sigma The MISER algorithm reduces the error by a factor of two, and also correctly estimates the error, miser ================== result = 1.390656 sigma = 0.003743 exact = 1.393204 error = 0.002548 = 0.7 sigma In the case of the VEGAS algorithm the program uses an initial warmup 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 chisquared 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. vegas warmup ================== 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
If the value of References and Further ReadingThe MISER algorithm is described in the following article,
The VEGAS algorithm is described in the following papers,
Go to the first, previous, next, last section, table of contents. 