Editing Monte Carlo integration (section)

== Recursive stratified sampling ==
{{see also|Stratified sampling}}
[[Image:Strata.png|thumb|right|An illustration of Recursive Stratified Sampling. In this example, the function:
<math>f(x,y) = \begin{cases}1 & x^2+y^2<1 \\0 & x^2+y^2 \ge 1 \end{cases}</math> <br>
from the above illustration was integrated within a unit square using the suggested algorithm. The sampled points were recorded and plotted. Clearly stratified sampling algorithm concentrates the points in the regions where the variation of the function is largest.]]

'''Recursive stratified sampling''' is a generalization of one-dimensional [[adaptive quadrature]]s to multi-dimensional integrals. On each recursion step the integral and the error are estimated using a plain Monte Carlo algorithm. If the error estimate is larger than the required accuracy the integration volume is divided into sub-volumes and the procedure is recursively applied to sub-volumes.

The ordinary 'dividing by two' strategy does not work for multi-dimensions as the number of sub-volumes grows far too quickly to keep track. Instead one estimates along which dimension a subdivision should bring the most dividends and only subdivides the volume along this dimension.

The stratified sampling algorithm concentrates the sampling points in the regions where the variance of the function is largest thus reducing the grand variance and making the sampling more effective, as shown on the illustration.

The popular MISER routine implements a similar algorithm.

=== MISER Monte Carlo ===
The MISER algorithm 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.<ref>{{harvnb|Press|Farrar|1990|pp=190-195}}</ref>

The idea of stratified sampling begins with the observation that for two [[Disjoint sets|disjoint]] regions ''a'' and ''b'' with Monte Carlo estimates of the integral <math>E_a(f)</math> and <math>E_b(f)</math> and variances <math>\sigma_a^2(f)</math> and <math>\sigma_b^2(f)</math>, the variance Var(''f'') of the combined estimate
<math display="block">E(f) = \tfrac{1}{2} \left (E_a(f) + E_b(f) \right )</math>
is given by,
<math display="block">\mathrm{Var}(f) = \frac{\sigma_a^2(f)}{4 N_a} + \frac{\sigma_b^2(f)}{4 N_b}</math>

It can be shown that this variance is minimized by distributing the points such that,
<math display="block">\frac{N_a}{N_a + N_b} = \frac{\sigma_a}{\sigma_a + \sigma_b}</math>

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<sub>a</sub>'' and ''N<sub>b</sub>''. 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.