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,
It can be shown that this variance is minimized by distributing the
points such that,
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
.
gsl_monte_miser_state * gsl_monte_miser_alloc (size_t dim)

Function 
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.

int gsl_monte_miser_init (gsl_monte_miser_state* s)

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

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

Function 
This routines uses the MISER Monte Carlo algorithm to integrate the
function f over the dimdimensional 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.

void gsl_monte_miser_free (gsl_monte_miser_state * s)

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

The MISER algorithm has several configurable parameters. The
following variables can be accessed through the
gsl_monte_miser_state
struct,
double estimate_frac

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

size_t min_calls

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

size_t min_calls_per_bisection

Variable 
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 subregion and terminates its branch of
the recursion. The default value of this parameter is 32 *
min_calls .

This parameter controls how the estimated variances for the two
subregions 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 subregions 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,
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.

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.
