1 @cindex Monte Carlo integration
2 @cindex stratified sampling in monte carlo integration
3 This chapter describes routines for multidimensional Monte Carlo
4 integration. These include the traditional Monte Carlo method and
5 adaptive algorithms such as @sc{vegas} and @sc{miser} which use
6 importance sampling and stratified sampling techniques. Each algorithm
7 computes an estimate of a multidimensional definite integral of the
12 I = \int_{x_l}^{x_u} dx\,\int_{y_l}^{y_u}dy\,... f(x,y,...)
19 I = \int_xl^xu dx \int_yl^yu dy ... f(x, y, ...)
24 over a hypercubic region @math{((x_l,x_u)}, @math{(y_l,y_u), ...)} using
25 a fixed number of function calls. The routines also provide a
26 statistical estimate of the error on the result. This error estimate
27 should be taken as a guide rather than as a strict error bound---random
28 sampling of the region may not uncover all the important features
29 of the function, resulting in an underestimate of the error.
31 The functions are defined in separate header files for each routine,
32 @code{gsl_monte_plain.h}, @file{gsl_monte_miser.h} and
33 @file{gsl_monte_vegas.h}.
36 * Monte Carlo Interface::
40 * Monte Carlo Examples::
41 * Monte Carlo Integration References and Further Reading::
44 @node Monte Carlo Interface
46 All of the Monte Carlo integration routines use the same general form of
47 interface. There is an allocator to allocate memory for control
48 variables and workspace, a routine to initialize those control
49 variables, the integrator itself, and a function to free the space when
52 Each integration function requires a random number generator to be
53 supplied, and returns an estimate of the integral and its standard
54 deviation. The accuracy of the result is determined by the number of
55 function calls specified by the user. If a known level of accuracy is
56 required this can be achieved by calling the integrator several times
57 and averaging the individual results until the desired accuracy is
60 Random sample points used within the Monte Carlo routines are always
61 chosen strictly within the integration region, so that endpoint
62 singularities are automatically avoided.
64 The function to be integrated has its own datatype, defined in the
65 header file @file{gsl_monte.h}.
67 @deftp {Data Type} gsl_monte_function
69 This data type defines a general function with parameters for Monte
73 @item double (* f) (double * @var{x}, size_t @var{dim}, void * @var{params})
74 this function should return the value
75 @c{$f(x,\hbox{\it params})$}
76 @math{f(x,params)} for the argument @var{x} and parameters @var{params},
77 where @var{x} is an array of size @var{dim} giving the coordinates of
78 the point where the function is to be evaluated.
81 the number of dimensions for @var{x}.
84 a pointer to the parameters of the function.
89 Here is an example for a quadratic function in two dimensions,
93 f(x,y) = a x^2 + b x y + c y^2
100 f(x,y) = a x^2 + b x y + c y^2
105 with @math{a = 3}, @math{b = 2}, @math{c = 1}. The following code
106 defines a @code{gsl_monte_function} @code{F} which you could pass to an
110 struct my_f_params @{ double a; double b; double c; @};
113 my_f (double x[], size_t dim, void * p) @{
114 struct my_f_params * fp = (struct my_f_params *)p;
118 fprintf (stderr, "error: dim != 2");
122 return fp->a * x[0] * x[0]
123 + fp->b * x[0] * x[1]
124 + fp->c * x[1] * x[1];
127 gsl_monte_function F;
128 struct my_f_params params = @{ 3.0, 2.0, 1.0 @};
136 The function @math{f(x)} can be evaluated using the following macro,
139 #define GSL_MONTE_FN_EVAL(F,x)
140 (*((F)->f))(x,(F)->dim,(F)->params)
143 @node PLAIN Monte Carlo
144 @section PLAIN Monte Carlo
145 @cindex plain monte carlo
146 The plain Monte Carlo algorithm samples points randomly from the
147 integration region to estimate the integral and its error. Using this
148 algorithm the estimate of the integral @math{E(f; N)} for @math{N}
149 randomly distributed points @math{x_i} is given by,
153 E(f; N) = V \langle f \rangle = {V \over N} \sum_i^N f(x_i)
160 E(f; N) = = V <f> = (V / N) \sum_i^N f(x_i)
165 where @math{V} is the volume of the integration region. The error on
166 this estimate @math{\sigma(E;N)} is calculated from the estimated
167 variance of the mean,
171 \sigma^2 (E; N) = {V \over N } \sum_i^N (f(x_i) - \langle f \rangle)^2.
178 \sigma^2 (E; N) = (V / N) \sum_i^N (f(x_i) - <f>)^2.
183 For large @math{N} this variance decreases asymptotically as
184 @math{\Var(f)/N}, where @math{\Var(f)} is the true variance of the
185 function over the integration region. The error estimate itself should
186 decrease as @c{$\sigma(f)/\sqrt{N}$}
187 @math{\sigma(f)/\sqrt@{N@}}. The familiar law of errors
188 decreasing as @c{$1/\sqrt{N}$}
189 @math{1/\sqrt@{N@}} applies---to reduce the error by a
190 factor of 10 requires a 100-fold increase in the number of sample
193 The functions described in this section are declared in the header file
194 @file{gsl_monte_plain.h}.
196 @deftypefun {gsl_monte_plain_state *} gsl_monte_plain_alloc (size_t @var{dim})
197 This function allocates and initializes a workspace for Monte Carlo
198 integration in @var{dim} dimensions.
201 @deftypefun int gsl_monte_plain_init (gsl_monte_plain_state* @var{s})
202 This function initializes a previously allocated integration state.
203 This allows an existing workspace to be reused for different
207 @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})
208 This routines uses the plain Monte Carlo algorithm to integrate the
209 function @var{f} over the @var{dim}-dimensional hypercubic region
210 defined by the lower and upper limits in the arrays @var{xl} and
211 @var{xu}, each of size @var{dim}. The integration uses a fixed number
212 of function calls @var{calls}, and obtains random sampling points using
213 the random number generator @var{r}. A previously allocated workspace
214 @var{s} must be supplied. The result of the integration is returned in
215 @var{result}, with an estimated absolute error @var{abserr}.
218 @deftypefun void gsl_monte_plain_free (gsl_monte_plain_state * @var{s})
219 This function frees the memory associated with the integrator state
225 @cindex MISER monte carlo integration
226 @cindex recursive stratified sampling, MISER
228 The @sc{miser} algorithm of Press and Farrar is based on recursive
229 stratified sampling. This technique aims to reduce the overall
230 integration error by concentrating integration points in the regions of
233 The idea of stratified sampling begins with the observation that for two
234 disjoint regions @math{a} and @math{b} with Monte Carlo estimates of the
235 integral @math{E_a(f)} and @math{E_b(f)} and variances
236 @math{\sigma_a^2(f)} and @math{\sigma_b^2(f)}, the variance
237 @math{\Var(f)} of the combined estimate
238 @c{$E(f) = {1\over 2} (E_a(f) + E_b(f))$}
239 @math{E(f) = (1/2) (E_a(f) + E_b(f))}
244 \Var(f) = {\sigma_a^2(f) \over 4 N_a} + {\sigma_b^2(f) \over 4 N_b}.
251 \Var(f) = (\sigma_a^2(f) / 4 N_a) + (\sigma_b^2(f) / 4 N_b).
256 It can be shown that this variance is minimized by distributing the
261 {N_a \over N_a+N_b} = {\sigma_a \over \sigma_a + \sigma_b}.
268 N_a / (N_a + N_b) = \sigma_a / (\sigma_a + \sigma_b).
273 Hence the smallest error estimate is obtained by allocating sample
274 points in proportion to the standard deviation of the function in each
277 The @sc{miser} algorithm proceeds by bisecting the integration region
278 along one coordinate axis to give two sub-regions at each step. The
279 direction is chosen by examining all @math{d} possible bisections and
280 selecting the one which will minimize the combined variance of the two
281 sub-regions. The variance in the sub-regions is estimated by sampling
282 with a fraction of the total number of points available to the current
283 step. The same procedure is then repeated recursively for each of the
284 two half-spaces from the best bisection. The remaining sample points are
285 allocated to the sub-regions using the formula for @math{N_a} and
286 @math{N_b}. This recursive allocation of integration points continues
287 down to a user-specified depth where each sub-region is integrated using
288 a plain Monte Carlo estimate. These individual values and their error
289 estimates are then combined upwards to give an overall result and an
290 estimate of its error.
292 The functions described in this section are declared in the header file
293 @file{gsl_monte_miser.h}.
295 @deftypefun {gsl_monte_miser_state *} gsl_monte_miser_alloc (size_t @var{dim})
296 This function allocates and initializes a workspace for Monte Carlo
297 integration in @var{dim} dimensions. The workspace is used to maintain
298 the state of the integration.
301 @deftypefun int gsl_monte_miser_init (gsl_monte_miser_state* @var{s})
302 This function initializes a previously allocated integration state.
303 This allows an existing workspace to be reused for different
307 @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})
308 This routines uses the @sc{miser} Monte Carlo algorithm to integrate the
309 function @var{f} over the @var{dim}-dimensional hypercubic region
310 defined by the lower and upper limits in the arrays @var{xl} and
311 @var{xu}, each of size @var{dim}. The integration uses a fixed number
312 of function calls @var{calls}, and obtains random sampling points using
313 the random number generator @var{r}. A previously allocated workspace
314 @var{s} must be supplied. The result of the integration is returned in
315 @var{result}, with an estimated absolute error @var{abserr}.
318 @deftypefun void gsl_monte_miser_free (gsl_monte_miser_state * @var{s})
319 This function frees the memory associated with the integrator state
323 The @sc{miser} algorithm has several configurable parameters. The
324 following variables can be accessed through the
325 @code{gsl_monte_miser_state} struct,
327 @deftypevar double estimate_frac
328 This parameter specifies the fraction of the currently available number of
329 function calls which are allocated to estimating the variance at each
330 recursive step. The default value is 0.1.
333 @deftypevar size_t min_calls
334 This parameter specifies the minimum number of function calls required
335 for each estimate of the variance. If the number of function calls
336 allocated to the estimate using @var{estimate_frac} falls below
337 @var{min_calls} then @var{min_calls} are used instead. This ensures
338 that each estimate maintains a reasonable level of accuracy. The
339 default value of @var{min_calls} is @code{16 * dim}.
342 @deftypevar size_t min_calls_per_bisection
343 This parameter specifies the minimum number of function calls required
344 to proceed with a bisection step. When a recursive step has fewer calls
345 available than @var{min_calls_per_bisection} it performs a plain Monte
346 Carlo estimate of the current sub-region and terminates its branch of
347 the recursion. The default value of this parameter is @code{32 *
351 @deftypevar double alpha
352 This parameter controls how the estimated variances for the two
353 sub-regions of a bisection are combined when allocating points. With
354 recursive sampling the overall variance should scale better than
355 @math{1/N}, since the values from the sub-regions will be obtained using
356 a procedure which explicitly minimizes their variance. To accommodate
357 this behavior the @sc{miser} algorithm allows the total variance to
358 depend on a scaling parameter @math{\alpha},
362 \Var(f) = {\sigma_a \over N_a^\alpha} + {\sigma_b \over N_b^\alpha}.
369 \Var(f) = @{\sigma_a \over N_a^\alpha@} + @{\sigma_b \over N_b^\alpha@}.
374 The authors of the original paper describing @sc{miser} recommend the
375 value @math{\alpha = 2} as a good choice, obtained from numerical
376 experiments, and this is used as the default value in this
380 @deftypevar double dither
381 This parameter introduces a random fractional variation of size
382 @var{dither} into each bisection, which can be used to break the
383 symmetry of integrands which are concentrated near the exact center of
384 the hypercubic integration region. The default value of dither is zero,
385 so no variation is introduced. If needed, a typical value of
391 @cindex VEGAS monte carlo integration
392 @cindex importance sampling, VEGAS
394 The @sc{vegas} algorithm of Lepage is based on importance sampling. It
395 samples points from the probability distribution described by the
396 function @math{|f|}, so that the points are concentrated in the regions
397 that make the largest contribution to the integral.
399 In general, if the Monte Carlo integral of @math{f} is sampled with
400 points distributed according to a probability distribution described by
401 the function @math{g}, we obtain an estimate @math{E_g(f; N)},
405 E_g(f; N) = E(f/g; N)
412 E_g(f; N) = E(f/g; N)
417 with a corresponding variance,
421 \Var_g(f; N) = \Var(f/g; N).
428 \Var_g(f; N) = \Var(f/g; N).
433 If the probability distribution is chosen as @math{g = |f|/I(|f|)} then
434 it can be shown that the variance @math{V_g(f; N)} vanishes, and the
435 error in the estimate will be zero. In practice it is not possible to
436 sample from the exact distribution @math{g} for an arbitrary function, so
437 importance sampling algorithms aim to produce efficient approximations
438 to the desired distribution.
440 The @sc{vegas} algorithm approximates the exact distribution by making a
441 number of passes over the integration region while histogramming the
442 function @math{f}. Each histogram is used to define a sampling
443 distribution for the next pass. Asymptotically this procedure converges
444 to the desired distribution. In order
445 to avoid the number of histogram bins growing like @math{K^d} the
446 probability distribution is approximated by a separable function:
447 @c{$g(x_1, x_2,\ldots) = g_1(x_1) g_2(x_2)\ldots$}
448 @math{g(x_1, x_2, ...) = g_1(x_1) g_2(x_2) ...}
449 so that the number of bins required is only @math{Kd}.
450 This is equivalent to locating the peaks of the function from the
451 projections of the integrand onto the coordinate axes. The efficiency
452 of @sc{vegas} depends on the validity of this assumption. It is most
453 efficient when the peaks of the integrand are well-localized. If an
454 integrand can be rewritten in a form which is approximately separable
455 this will increase the efficiency of integration with @sc{vegas}.
457 @sc{vegas} incorporates a number of additional features, and combines both
458 stratified sampling and importance sampling. The integration region is
459 divided into a number of ``boxes'', with each box getting a fixed
460 number of points (the goal is 2). Each box can then have a fractional
461 number of bins, but if the ratio of bins-per-box is less than two, Vegas switches to a
462 kind variance reduction (rather than importance sampling).
465 @deftypefun {gsl_monte_vegas_state *} gsl_monte_vegas_alloc (size_t @var{dim})
466 This function allocates and initializes a workspace for Monte Carlo
467 integration in @var{dim} dimensions. The workspace is used to maintain
468 the state of the integration.
471 @deftypefun int gsl_monte_vegas_init (gsl_monte_vegas_state* @var{s})
472 This function initializes a previously allocated integration state.
473 This allows an existing workspace to be reused for different
477 @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})
478 This routines uses the @sc{vegas} Monte Carlo algorithm to integrate the
479 function @var{f} over the @var{dim}-dimensional hypercubic region
480 defined by the lower and upper limits in the arrays @var{xl} and
481 @var{xu}, each of size @var{dim}. The integration uses a fixed number
482 of function calls @var{calls}, and obtains random sampling points using
483 the random number generator @var{r}. A previously allocated workspace
484 @var{s} must be supplied. The result of the integration is returned in
485 @var{result}, with an estimated absolute error @var{abserr}. The result
486 and its error estimate are based on a weighted average of independent
487 samples. The chi-squared per degree of freedom for the weighted average
488 is returned via the state struct component, @var{s->chisq}, and must be
489 consistent with 1 for the weighted average to be reliable.
492 @deftypefun void gsl_monte_vegas_free (gsl_monte_vegas_state * @var{s})
493 This function frees the memory associated with the integrator state
497 The @sc{vegas} algorithm computes a number of independent estimates of the
498 integral internally, according to the @code{iterations} parameter
499 described below, and returns their weighted average. Random sampling of
500 the integrand can occasionally produce an estimate where the error is
501 zero, particularly if the function is constant in some regions. An
502 estimate with zero error causes the weighted average to break down and
503 must be handled separately. In the original Fortran implementations of
504 @sc{vegas} the error estimate is made non-zero by substituting a small
505 value (typically @code{1e-30}). The implementation in GSL differs from
506 this and avoids the use of an arbitrary constant---it either assigns
507 the value a weight which is the average weight of the preceding
508 estimates or discards it according to the following procedure,
511 @item current estimate has zero error, weighted average has finite error
513 The current estimate is assigned a weight which is the average weight of
514 the preceding estimates.
516 @item current estimate has finite error, previous estimates had zero error
518 The previous estimates are discarded and the weighted averaging
519 procedure begins with the current estimate.
521 @item current estimate has zero error, previous estimates had zero error
523 The estimates are averaged using the arithmetic mean, but no error is computed.
526 The @sc{vegas} algorithm is highly configurable. The following variables
527 can be accessed through the @code{gsl_monte_vegas_state} struct,
529 @deftypevar double result
530 @deftypevarx double sigma
531 These parameters contain the raw value of the integral @var{result} and
532 its error @var{sigma} from the last iteration of the algorithm.
535 @deftypevar double chisq
536 This parameter gives the chi-squared per degree of freedom for the
537 weighted estimate of the integral. The value of @var{chisq} should be
538 close to 1. A value of @var{chisq} which differs significantly from 1
539 indicates that the values from different iterations are inconsistent.
540 In this case the weighted error will be under-estimated, and further
541 iterations of the algorithm are needed to obtain reliable results.
544 @deftypevar double alpha
545 The parameter @code{alpha} controls the stiffness of the rebinning
546 algorithm. It is typically set between one and two. A value of zero
547 prevents rebinning of the grid. The default value is 1.5.
550 @deftypevar size_t iterations
551 The number of iterations to perform for each call to the routine. The
552 default value is 5 iterations.
555 @deftypevar int stage
556 Setting this determines the @dfn{stage} of the calculation. Normally,
557 @code{stage = 0} which begins with a new uniform grid and empty weighted
558 average. Calling vegas with @code{stage = 1} retains the grid from the
559 previous run but discards the weighted average, so that one can ``tune''
560 the grid using a relatively small number of points and then do a large
561 run with @code{stage = 1} on the optimized grid. Setting @code{stage =
562 2} keeps the grid and the weighted average from the previous run, but
563 may increase (or decrease) the number of histogram bins in the grid
564 depending on the number of calls available. Choosing @code{stage = 3}
565 enters at the main loop, so that nothing is changed, and is equivalent
566 to performing additional iterations in a previous call.
570 The possible choices are @code{GSL_VEGAS_MODE_IMPORTANCE},
571 @code{GSL_VEGAS_MODE_STRATIFIED}, @code{GSL_VEGAS_MODE_IMPORTANCE_ONLY}.
572 This determines whether @sc{vegas} will use importance sampling or
573 stratified sampling, or whether it can pick on its own. In low
574 dimensions @sc{vegas} uses strict stratified sampling (more precisely,
575 stratified sampling is chosen if there are fewer than 2 bins per box).
578 @deftypevar int verbose
579 @deftypevarx {FILE *} ostream
580 These parameters set the level of information printed by @sc{vegas}. All
581 information is written to the stream @var{ostream}. The default setting
582 of @var{verbose} is @code{-1}, which turns off all output. A
583 @var{verbose} value of @code{0} prints summary information about the
584 weighted average and final result, while a value of @code{1} also
585 displays the grid coordinates. A value of @code{2} prints information
586 from the rebinning procedure for each iteration.
589 @node Monte Carlo Examples
592 The example program below uses the Monte Carlo routines to estimate the
593 value of the following 3-dimensional integral from the theory of random
598 I = \int_{-\pi}^{+\pi} {dk_x \over 2\pi}
599 \int_{-\pi}^{+\pi} {dk_y \over 2\pi}
600 \int_{-\pi}^{+\pi} {dk_z \over 2\pi}
601 { 1 \over (1 - \cos(k_x)\cos(k_y)\cos(k_z))}.
608 I = \int_@{-pi@}^@{+pi@} @{dk_x/(2 pi)@}
609 \int_@{-pi@}^@{+pi@} @{dk_y/(2 pi)@}
610 \int_@{-pi@}^@{+pi@} @{dk_z/(2 pi)@}
611 1 / (1 - cos(k_x)cos(k_y)cos(k_z)).
616 The analytic value of this integral can be shown to be @math{I =
617 \Gamma(1/4)^4/(4 \pi^3) = 1.393203929685676859...}. The integral gives
618 the mean time spent at the origin by a random walk on a body-centered
619 cubic lattice in three dimensions.
621 For simplicity we will compute the integral over the region
622 @math{(0,0,0)} to @math{(\pi,\pi,\pi)} and multiply by 8 to obtain the
623 full result. The integral is slowly varying in the middle of the region
624 but has integrable singularities at the corners @math{(0,0,0)},
625 @math{(0,\pi,\pi)}, @math{(\pi,0,\pi)} and @math{(\pi,\pi,0)}. The
626 Monte Carlo routines only select points which are strictly within the
627 integration region and so no special measures are needed to avoid these
631 @verbatiminclude examples/monte.c
635 With 500,000 function calls the plain Monte Carlo algorithm achieves a
636 fractional error of 0.6%. The estimated error @code{sigma} is
637 consistent with the actual error, and the computed result differs from
638 the true result by about one standard deviation,
641 plain ==================
645 error = -0.007337 = 0.9 sigma
649 The @sc{miser} algorithm reduces the error by a factor of two, and also
650 correctly estimates the error,
653 miser ==================
657 error = -0.002548 = 0.7 sigma
661 In the case of the @sc{vegas} algorithm the program uses an initial
662 warm-up run of 10,000 function calls to prepare, or ``warm up'', the grid.
663 This is followed by a main run with five iterations of 100,000 function
664 calls. The chi-squared per degree of freedom for the five iterations are
665 checked for consistency with 1, and the run is repeated if the results
666 have not converged. In this case the estimates are consistent on the
670 vegas warm-up ==================
674 error = -0.006278 = 2 sigma
676 result = 1.392957 sigma = 0.000452 chisq/dof = 1.1
677 vegas final ==================
681 error = -0.000247 = 0.5 sigma
685 If the value of @code{chisq} had differed significantly from 1 it would
686 indicate inconsistent results, with a correspondingly underestimated
687 error. The final estimate from @sc{vegas} (using a similar number of
688 function calls) is significantly more accurate than the other two
691 @node Monte Carlo Integration References and Further Reading
692 @section References and Further Reading
694 The @sc{miser} algorithm is described in the following article by Press
699 W.H. Press, G.R. Farrar, @cite{Recursive Stratified Sampling for
700 Multidimensional Monte Carlo Integration},
701 Computers in Physics, v4 (1990), pp190--195.
705 The @sc{vegas} algorithm is described in the following papers,
710 @cite{A New Algorithm for Adaptive Multidimensional Integration},
711 Journal of Computational Physics 27, 192--203, (1978)
715 @cite{VEGAS: An Adaptive Multi-dimensional Integration Program},
716 Cornell preprint CLNS 80-447, March 1980