25.3 MISER

The 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 sub-region.

The 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 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 N_a and 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 gsl_monte_miser.h.

Function: gsl_monte_miser_state * gsl_monte_miser_alloc (size_t dim)

This function allocates and initializes a workspace for Monte Carlo integration in dim dimensions. The workspace is used to maintain the state of the integration.

Function: int gsl_monte_miser_init (gsl_monte_miser_state* s)

This function initializes a previously allocated integration state. This allows an existing workspace to be reused for different integrations.

Function: int gsl_monte_miser_integrate (gsl_monte_function * f, const double xl[], const double xu[], size_t dim, size_t calls, gsl_rng * r, gsl_monte_miser_state * s, double * result, double * abserr)

This routines uses the MISER Monte Carlo algorithm to integrate the function f over the dim-dimensional hypercubic region defined by the lower and upper limits in the arrays xl and xu, each of size dim. The integration uses a fixed number of function calls calls, and obtains random sampling points using the random number generator r. A previously allocated workspace s must be supplied. The result of the integration is returned in result, with an estimated absolute error abserr.

Function: void gsl_monte_miser_free (gsl_monte_miser_state * s)

This function frees the memory associated with the integrator state s.

The MISER algorithm has several configurable parameters which can be changed using the following two functions.¹³

Function: void gsl_monte_miser_params_get (const gsl_monte_miser_state * s, gsl_monte_miser_params * params)

This function copies the parameters of the integrator state into the user-supplied params structure.

Function: void gsl_monte_miser_params_set (gsl_monte_miser_state * s, const gsl_monte_miser_params * params)

This function sets the integrator parameters based on values provided in the params structure.

Typically the values of the parameters are first read using gsl_monte_miser_params_get, the necessary changes are made to the fields of the params structure, and the values are copied back into the integrator state using gsl_monte_miser_params_set. The functions use the gsl_monte_miser_params structure which contains the following fields:

Variable: 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.

Variable: 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 estimate_frac falls below min_calls then min_calls are used instead. This ensures that each estimate maintains a reasonable level of accuracy. The default value of min_calls is 16 * dim.

Variable: 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 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 32 * min_calls.

Variable: 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 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 MISER algorithm allows the total variance to depend on a scaling parameter \alpha,

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

The authors of the original paper describing MISER recommend the value \alpha = 2 as a good choice, obtained from numerical experiments, and this is used as the default value in this implementation.

Variable: double dither

This parameter introduces a random fractional variation of size 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 dither is 0.1.