In mathematics, Monte Carlo integration is a technique for numerical integration using random numbers. It is a particular Monte Carlo method that numerically computes a definite integral. While other algorithms usually evaluate the integrand at a regular grid,^{[1]} Monte Carlo randomly choose points at which the integrand is evaluated.^{[2]} This method is particularly useful for higherdimensional integrals.^{[3]}
There are different methods to perform a Monte Carlo integration, such as uniform sampling, stratified sampling, importance sampling, sequential Monte Carlo (also known as a particle filter), and mean field particle methods.
YouTube Encyclopedic

1/5Views:44 37313 535111 5782 597655

✪ Monte Carlo integration

✪ 1 Monte Carlo Simulation and Integration

✪ What is Monte Carlo?

✪ R Tutorial 6: Monte Carlo Integration

✪ Basic Monte Carlo integration with Matlab
Transcription
Contents
Overview
In numerical integration, methods such as the trapezoidal rule use a deterministic approach. Monte Carlo integration, on the other hand, employs a nondeterministic approach: each realization provides a different outcome. In Monte Carlo, the final outcome is an approximation of the correct value with respective error bars, and the correct value is likely to be within those error bars.
The problem Monte Carlo integration addresses is the computation of a multidimensional definite integral
where Ω, a subset of R^{m}, has volume
The naive Monte Carlo approach is to sample points uniformly on Ω:^{[4]} given N uniform samples,
I can be approximated by
 .
This is because the law of large numbers ensures that
 .
Given the estimation of I from Q_{N}, the error bars of Q_{N} can be estimated by the sample variance using the unbiased estimate of the variance.
which leads to
 .
As long as the sequence
is bounded, this variance decreases asymptotically to zero as 1/N. The estimation of the error of Q_{N} is thus
which decreases as . This is standard error of the mean multiplied with . This result does not depend on the number of dimensions of the integral, which is the promised advantage of Monte Carlo integration against most deterministic methods that depend exponentially on the dimension.^{[5]} It is important to notice that, unlike in deterministic methods, the estimate of the error is not a strict error bound; random sampling may not uncover all the important features of the integrand that can result in an underestimate of the error.
While the naive Monte Carlo works for simple examples, an improvement over deterministic algorithms only can be accomplished with algorithms that use problem specific sampling distributions. With an appropriate sample distribution it is possible to exploit the fact that almost all higherdimensional integrands are very localized and only small subspace notably contributes to the integral^{[6]}. A large part of the Monte Carlo literature is dedicated in developing strategies to improve the error estimates. In particular, stratified sampling—dividing the region in subdomains—, and importance sampling—sampling from nonuniform distributions—are two of such techniques.
Example
A paradigmatic example of a Monte Carlo integration is the estimation of π. Consider the function
and the set Ω = [−1,1] × [−1,1] with V = 4. Notice that
Thus, a crude way of calculating the value of π with Monte Carlo integration is to pick N random numbers on Ω and compute
In the figure on the right, the relative error is measured as a function of N, confirming the .
C example
Keep in mind that a true random number generator should be used.
int i, throws = 99999, circleDarts = 0;
long double randX, randY, pi;
srand(time(NULL));
for (i = 0; i < throws; ++i) {
randX = rand() / (double)RAND_MAX;
randY = rand() / (double)RAND_MAX;
if (1 > ((randX*randX) + (randY*randY))) ++circleDarts;
}
pi = 4 * (circleDarts/throws);
Wolfram Mathematica example
The code below describes a process of integrating the function
from using the MonteCarlo method in Mathematica:
func[x_] := 1/(1 + Sinh[2*x]*(Log[x])^2);
(*Sample from truncated normal distribution to speed up convergence*)
Distrib[x_, average_, var_] := PDF[NormalDistribution[average, var], 1.1*x  0.1];
n = 10;
RV = RandomVariate[TruncatedDistribution[{0.8, 3}, NormalDistribution[1, 0.399]], n];
Int = 1/n Total[func[RV]/Distrib[RV, 1, 0.399]]*Integrate[Distrib[x, 1, 0.399], {x, 0.8, 3}]
NIntegrate[func[x], {x, 0.8, 3}] (*Compare with real answer*)
C++ example
Boost.Math (version 1.67 and later) provides a multithreaded MonteCarlo integration routine.
# include <boost/math/quadrature/naive_monte_carlo.hpp>
// Define a function to integrate:
auto g = [](std::vector<double> const & x)
{
constexpr const double A = 1.0 / (M_PI * M_PI * M_PI);
return A / (1.0  cos(x[0])*cos(x[1])*cos(x[2]));
};
std::vector<std::pair<double, double>> bounds{{0, M_PI}, {0, M_PI}, {0, M_PI}};
double error_goal = 0.001;
naive_monte_carlo<double, decltype(g)> mc(g, bounds, error_goal);
std::future<double> task = mc.integrate();
while (task.wait_for(std::chrono::seconds(1)) != std::future_status::ready) {
display_progress(mc);
}
double I = task.get();
Recursive stratified sampling
Recursive stratified sampling is a generalization of onedimensional adaptive quadratures to multidimensional 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 subvolumes and the procedure is recursively applied to subvolumes.
The ordinary 'dividing by two' strategy does not work for multidimensions as the number of subvolumes 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.^{[7]}
The idea of stratified sampling begins with the observation that for two disjoint regions a and b with Monte Carlo estimates of the integral and and variances and , the variance Var(f) of the combined estimate
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.
Importance sampling
VEGAS Monte Carlo
The VEGAS algorithm takes advantage of the information stored during the sampling, and uses it and importance sampling to efficiently estimate the integral I. 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:
with a corresponding variance,
If the probability distribution is chosen as
then it can be shown that the variance vanishes, and the error in the estimate will be zero.^{[citation needed]} 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 which creates the histogram of the function f. Each histogram is used to define a sampling distribution for the next pass. Asymptotically this procedure converges to the desired distribution.^{[8]} In order to avoid the number of histogram bins growing like K^{d}, the probability distribution is approximated by a separable function:
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.^{[8]} The integration region is divided into a number of "boxes", with each box getting a 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).
This routines uses the VEGAS 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. The result and its error estimate are based on a weighted average of independent samples.
The VEGAS algorithm computes a number of independent estimates of the integral internally, according to the iterations parameter described below, and returns their weighted average. Random sampling of the integrand can occasionally produce an estimate where the error is zero, particularly if the function is constant in some regions. An estimate with zero error causes the weighted average to break down and must be handled separately.
Importance sampling algorithm
Importance sampling provides a very important tool to perform MonteCarlo integration.^{[3]}^{[9]} The main result of importance sampling to this method is that the uniform sampling of is a particular case of a more generic choice, on which the samples are drawn from any distribution . The idea is that can be chosen to decrease the variance of the measurement Q_{N}.
Consider the following example where one would like to numerically integrate a gaussian function, centered at 0, with σ = 1, from −1000 to 1000. Naturally, if the samples are drawn uniformly on the interval [−1000, 1000], only a very small part of them would be significant to the integral. This can be improved by choosing a different distribution from where the samples are chosen, for instance by sampling according to a gaussian distribution centered at 0, with σ = 1. Of course the "right" choice strongly depends on the integrand.
Formally, given a set of samples chosen from a distribution
the estimator for I is given by^{[3]}
Intuitively, this says that if we pick a particular sample twice as much as other samples, we weight it half as much as the other samples. This estimator is naturally valid for uniform sampling, the case where is constant.
The MetropolisHastings algorithm is one of the most used algorithms to generate from ,^{[3]} thus providing an efficient way of computing integrals.
Multiple and adaptive importance sampling
When different proposal distributions, , are jointly used for drawing the samples different proper weighting functions can be employed (e.g., see ^{[10]}^{[11]}^{[12]}). In an adaptive setting, the proposal distributions, , and are updated each iteration of the adaptive importance sampling algorithm. Hence, since a population of proposal densities is used, several suitable combinations of sampling and weighting schemes can be employed.^{[13]}^{[14]}^{[15]}^{[16]}^{[17]}
See also
 Auxiliary field Monte Carlo
 Monte Carlo method in statistical physics
 Monte Carlo method
 Variance reduction
Notes
 ^ Press et al, 2007, Chap. 4.
 ^ Press et al, 2007, Chap. 7.
 ^ ^{a} ^{b} ^{c} ^{d} Newman, 1999, Chap. 2.
 ^ Newman, 1999, Chap. 1.
 ^ Press et al, 2007
 ^ MacKay, David (2003). "chapter 4.4 Typicality & chapter 29.1" (PDF). Information Theory, Inference and Learning Algorithms. Cambridge University Press. pp. 284–292. ISBN 9780521642989. MR 2012999.
 ^ Press, 1990, pp 190195.
 ^ ^{a} ^{b} Lepage, 1978
 ^ Kroese, D. P.; Taimre, T.; Botev, Z. I. (2011). Handbook of Monte Carlo Methods. John Wiley & Sons.
 ^ Veach, Eric; Guibas, Leonidas J. (19950101). Optimally Combining Sampling Techniques for Monte Carlo Rendering10.1145/218380.218498 (PDF). Proceedings of the 22Nd Annual Conference on Computer Graphics and Interactive Techniques. SIGGRAPH '95. New York, NY, USA. pp. 419–428. CiteSeerX 10.1.1.127.8105. doi:10.1145/218380.218498. ISBN 9780897917018.
 ^ Owen, Art; Associate, Yi Zhou (20000301). "Safe and Effective Importance Sampling". Journal of the American Statistical Association. 95 (449): 135–143. CiteSeerX 10.1.1.36.4536. doi:10.1080/01621459.2000.10473909. ISSN 01621459.
 ^ Elvira, V.; Martino, L.; Luengo, D.; Bugallo, M.F. (20151001). "Efficient Multiple Importance Sampling Estimators". IEEE Signal Processing Letters. 22 (10): 1757–1761. arXiv:1505.05391. Bibcode:2015ISPL...22.1757E. doi:10.1109/LSP.2015.2432078. ISSN 10709908.
 ^ Cappé, O.; Guillin, A.; Marin, J. M.; Robert, C. P. (20041201). "Population Monte Carlo". Journal of Computational and Graphical Statistics. 13 (4): 907–929. doi:10.1198/106186004X12803. ISSN 10618600.
 ^ Cappé, Olivier; Douc, Randal; Guillin, Arnaud; Marin, JeanMichel; Robert, Christian P. (20080425). "Adaptive importance sampling in general mixture classes". Statistics and Computing. 18 (4): 447–459. arXiv:0710.4242. doi:10.1007/s112220089059x. ISSN 09603174.
 ^ Cornuet, JeanMarie; Marin, JeanMichel; Mira, Antonietta; Robert, Christian P. (20121201). "Adaptive Multiple Importance Sampling". Scandinavian Journal of Statistics. 39 (4): 798–812. arXiv:0907.1254. doi:10.1111/j.14679469.2011.00756.x. ISSN 14679469.
 ^ Martino, L.; Elvira, V.; Luengo, D.; Corander, J. (20150801). "An Adaptive Population Importance Sampler: Learning From Uncertainty". IEEE Transactions on Signal Processing. 63 (16): 4422–4437. Bibcode:2015ITSP...63.4422M. CiteSeerX 10.1.1.464.9395. doi:10.1109/TSP.2015.2440215. ISSN 1053587X.
 ^ Bugallo, Mónica F.; Martino, Luca; Corander, Jukka (20151201). "Adaptive importance sampling in signal processing". Digital Signal Processing. Special Issue in Honour of William J. (Bill) Fitzgerald. 47: 36–49. doi:10.1016/j.dsp.2015.05.014.
References
 Caflisch, R. E. (1998). "Monte Carlo and quasiMonte Carlo methods". Acta Numerica. 7: 1–49. Bibcode:1998AcNum...7....1C. doi:10.1017/S0962492900002804.
 Weinzierl, S. (2000). "Introduction to Monte Carlo methods". arXiv:hepph/0006269.
 Press, W. H.; Farrar, G. R. (1990). "Recursive Stratified Sampling for Multidimensional Monte Carlo Integration". Computers in Physics. 4 (2): 190. Bibcode:1990ComPh...4..190P. doi:10.1063/1.4822899.
 Lepage, G. P. (1978). "A New Algorithm for Adaptive Multidimensional Integration". Journal of Computational Physics. 27 (2): 192–203. Bibcode:1978JCoPh..27..192L. doi:10.1016/00219991(78)900049.
 Lepage, G. P. (1980). "VEGAS: An Adaptive Multidimensional Integration Program". Cornell Preprint CLNS 80447.
 Hammersley, J. M.; Handscomb, D. C. (1964). Monte Carlo Methods. Methuen. ISBN 9780416523409.
 Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). Numerical Recipes: The Art of Scientific Computing (3rd ed.). New York: Cambridge University Press. ISBN 9780521880688.
 Newman, MEJ; Barkema, GT (1999). Monte Carlo Methods in Statistical Physics. Clarendon Press.
 Robert, CP; Casella, G (2004). Monte Carlo Statistical Methods (2nd ed.). Springer. ISBN 9781441919397.
External links
 Café math : Monte Carlo Integration : A blog article describing Monte Carlo integration (principle, hypothesis, confidence interval)
 Boost.Math : Naive Monte Carlo integration: Documentation for the C++ naive MonteCarlo routines
 Monte Carlo applet applied in statistical physics problems