View source for Monte Carlo integration

{{Short description|Numerical technique}}
[[Image:MonteCarloIntegrationCircle.svg|thumb|An illustration of Monte Carlo integration. In this example, the domain ''D'' is the inner circle and the domain E is the square. Because the square's area (4) can be easily calculated, the area of the circle (π*1.0<sup>2</sup>) can be estimated by the ratio (0.8) of the points inside the circle (40) to the total number of points (50), yielding an approximation for the circle's area of 4*0.8 = 3.2 ≈ π.]]

In [[Mathematics|mathematics]], '''Monte Carlo integration''' is a technique for numerical integration using [[Pseudorandomness|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,<ref>{{harvnb|Press|Teukolsky|Vetterling|Flannery|2007|loc=Chap. 4}}</ref> Monte Carlo randomly chooses points at which the integrand is evaluated.<ref>{{harvnb|Press|Teukolsky|Vetterling|Flannery|2007|loc=Chap. 7}}</ref> This method is particularly useful for higher-dimensional integrals.<ref name=newman1999ch2/>

There are different methods to perform a Monte Carlo integration, such as [[Uniform distribution (continuous)|uniform sampling]], [[Stratified sampling|stratified sampling]], [[Importance sampling|importance sampling]], [[Particle filter|sequential Monte Carlo]] (also known as a particle filter), and [[Mean-field particle methods|mean-field particle methods]].

== Overview ==
In numerical integration, methods such as the [[Trapezoidal rule|trapezoidal rule]] use a [[Deterministic algorithm|deterministic approach]]. Monte Carlo integration, on the other hand, employs a [[Stochastic|non-deterministic]] 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 [[Multiple integral|multidimensional definite integral]]
<math display="block">I = \int_{\Omega}f(\overline{\mathbf{x}}) \, d\overline{\mathbf{x}}</math>
where Ω, a subset of '''R'''<sup>''m''</sup>, has volume
<math display="block">V = \int_{\Omega}d\overline{\mathbf{x}}</math>

The naive Monte Carlo approach is to sample points uniformly on Ω:<ref name=newman1999ch1>{{harvnb|Newman|Barkema|1999|at=Chap. 1}}</ref> given ''N'' uniform samples,
<math display="block">\overline{\mathbf{x}}_1, \cdots, \overline{\mathbf{x}}_N\in \Omega,</math>

''I'' can be approximated by
<math display="block"> I \approx Q_N \equiv V \frac{1}{N} \sum_{i=1}^N f(\overline{\mathbf{x}}_i) = V \langle f\rangle.</math>

This is because the [[Law of large numbers|law of large numbers]] ensures that
<math display="block"> \lim_{N \to \infty} Q_N = I.</math>

Given the estimation of ''I'' from ''Q<sub>N</sub>'', the error bars of ''Q<sub>N</sub>'' can be estimated by the  [[Sample variance#Population variance and sample variance|sample variance]] using the [[Bias of an estimator#Sample variance|unbiased estimate of the variance]].

<math display="block"> \mathrm{Var}(f)\equiv\sigma_N^2 = \frac{1}{N-1} \sum_{i=1}^N \left (f(\overline{\mathbf{x}}_i) - \langle f \rangle \right )^2. </math>
which leads to
<math display="block"> \mathrm{Var}(Q_N) = \frac{V^2}{N^2} \sum_{i=1}^N \mathrm{Var}(f) = V^2\frac{\mathrm{Var}(f)}{N} = V^2\frac{\sigma_N^2}{N}.</math>

As long as the sequence
<math display="block"> \left \{ \sigma_1^2, \sigma_2^2, \sigma_3^2, \ldots \right \} </math>
is bounded, this variance decreases asymptotically to zero as 1/''N''. The estimation of the error of ''Q<sub>N</sub>'' is thus
<math display="block">\delta Q_N\approx\sqrt{\mathrm{Var}(Q_N)}=V\frac{\sigma_N}{\sqrt{N}},</math>
which decreases as <math>\tfrac{1}{\sqrt{N}}</math>. This is standard error of the mean multiplied with <math>V</math>. 
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.<ref>{{harvnb|Press|Teukolsky|Vetterling|Flannery|2007}}</ref> 
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 can only 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 higher-dimensional integrands are very localized and only small subspace notably contributes to the integral.<ref>{{harvnb|MacKay|2003|pages=284&ndash;292}}</ref>
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 sub-domains—and importance sampling—sampling from non-uniform distributions—are two examples of such techniques.

=== Example ===
[[File:Relative error of a Monte Carlo integration to calculate pi.svg|thumb|right|350px|Relative error as a function of the number of samples, showing the scaling <math>\tfrac{1}{\sqrt{N}}</math>]]

A paradigmatic example of a Monte Carlo integration is the estimation of π. Consider the function
<math display="block">H\left(x,y\right)=\begin{cases}
1 & \text{if }x^{2}+y^{2}\leq1\\
0 & \text{else}
\end{cases}</math>
and the set Ω = [−1,1] × [−1,1] with ''V'' = 4. Notice that
<math display="block">I_\pi = \int_\Omega H(x,y) dx dy = \pi.</math>

Thus, a crude way of calculating the value of π with Monte Carlo integration is to pick ''N'' random numbers on Ω and compute
<math display="block">Q_N = 4 \frac{1}{N}\sum_{i=1}^N H(x_{i},y_{i})</math>

In the figure on the right, the relative error <math>\tfrac{Q_N-\pi}{\pi}</math>  is measured as a function of ''N'', confirming the <math>\tfrac{1}{\sqrt{N}}</math>.

=== C example ===
Keep in mind that a true random number generator should be used.
<source lang="c">
int i, throws = 99999, insideCircle = 0;
double randX, randY, pi;

srand(time(NULL));

for (i = 0; i < throws; ++i) {
  randX = rand() / (double) RAND_MAX;
  randY = rand() / (double) RAND_MAX;
  if (randX * randX + randY * randY < 1) ++insideCircle;
}

pi = 4.0 * insideCircle / throws;
</source>

=== Python example ===
Made in [[Python (programming language)|Python]]. 
<source lang="python">
from numpy import random
import numpy as np

throws = 2000
inside_circle = 0
i = 0
radius = 1
while i < throws:
    # Choose random X and Y centered around 0,0
    x = random.uniform(-radius, radius)
    y = random.uniform(-radius, radius)
    # If the point is inside circle, increase variable
    if x**2 + y**2 <= radius**2:
        inside_circle += 1
    i += 1

# Calculate area and print; should be closer to Pi with increasing number of throws
area = (((2 * radius) ** 2) * inside_circle) / throws
print(area)
</source>

=== Wolfram Mathematica example ===
The code below describes a process of integrating the function
<math display="block">f(x) = \frac{1}{1+\sinh(2x)\log(x)^2}</math>
from <math>0.8<x<3</math> using the Monte-Carlo method in [[Software:Mathematica|Mathematica]]:

<source lang="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*)
</source>

== Recursive 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|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|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.

== Importance sampling ==
{{Main|Importance sampling}}

There are a variety of importance sampling algorithms, such as 
=== Importance sampling algorithm ===

Importance sampling provides a very important tool to perform Monte-Carlo integration.<ref name="newman1999ch2">{{harvnb|Newman|Barkema|1999|loc=Chap. 2}}</ref><ref name="kr11">{{harvnb|Kroese|Taimre|Botev|2011}}</ref> The main result of importance sampling to this method is that the uniform sampling of <math>\overline{\mathbf{x}}</math> is a particular case of a more generic choice, on which the samples are drawn from any distribution <math>p(\overline{\mathbf{x}})</math>. The idea is that <math>p(\overline{\mathbf{x}})</math> can be chosen to decrease the variance of the measurement ''Q<sub>N</sub>''.

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
<math display="block">p(\overline{\mathbf{x}}) : \qquad \overline{\mathbf{x}}_1, \cdots, \overline{\mathbf{x}}_N \in V, </math>
the estimator for ''I'' is given by<ref name="newman1999ch2" />
<math display="block"> Q_N \equiv \frac{1}{N} \sum_{i=1}^N \frac{f(\overline{\mathbf{x}}_i)}{p(\overline{\mathbf{x}}_i)}</math>

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 <math>p(\overline{\mathbf{x}})</math> is constant.

The [[Metropolis–Hastings algorithm]] is one of the most used algorithms to generate <math>\overline{\mathbf{x}}</math> from <math>p(\overline{\mathbf{x}})</math>,<ref name="newman1999ch2" /> thus providing an efficient way of computing integrals.

=== VEGAS Monte Carlo ===
{{Main|VEGAS algorithm}}

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.<ref name="Lepage, 1978">{{harvnb|Lepage|1978}}</ref> In order to avoid the number of histogram bins growing like ''K<sup>d</sup>'', the probability distribution is approximated by a separable function:
<math display="block">g(x_1, x_2, \ldots) = g_1(x_1) g_2(x_2) \ldots </math>
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 well-localized. 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.<ref name="Lepage, 1978"/>

== See also ==
* [[Quasi-Monte Carlo method]]
* Auxiliary field Monte Carlo
* [[Physics:Monte Carlo method in statistical physics|Monte Carlo method in statistical physics]]
* [[Monte Carlo method]]
* [[Variance reduction]]

== Notes ==
{{reflist|30em}}

== References ==
{{refbegin}}
* {{cite journal |first=R. E. |last=Caflisch |title=Monte Carlo and quasi-Monte Carlo methods |journal=[[Acta Numerica]] |volume=7 |year=1998 |pages=1–49 |doi=10.1017/S0962492900002804 |bibcode=1998AcNum...7....1C |s2cid=5708790 }}
* {{cite arXiv |first=S. |last=Weinzierl |title=Introduction to Monte Carlo methods |eprint=hep-ph/0006269 |year=2000}}
* {{cite journal |first1=W. H. |last1=Press |first2=G. R. |last2=Farrar |title=Recursive Stratified Sampling for Multidimensional Monte Carlo Integration |journal=Computers in Physics |volume=4 |issue=2 |pages=190 |year=1990 |doi=10.1063/1.4822899 |bibcode=1990ComPh...4..190P |doi-access=free }}
* {{cite journal |first=G. P. |last=Lepage |title=A New Algorithm for Adaptive Multidimensional Integration |journal=[[Physics:Journal of Computational Physics|Journal of Computational Physics]] |volume=27 |issue=2 |pages=192–203 |year=1978 |doi=10.1016/0021-9991(78)90004-9 |bibcode=1978JCoPh..27..192L }}
* {{cite journal |first=G. P. |last=Lepage |title=VEGAS: An Adaptive Multi-dimensional Integration Program |journal=Cornell Preprint CLNS 80-447 |year=1980 }}
* {{cite book |first1=J. M. |last1=Hammersley |first2=D. C. |last2=Handscomb |year=1964 |title=Monte Carlo Methods |publisher=Methuen |isbn=978-0-416-52340-9 }}
* {{cite book|last1 = Kroese|first1 = D. P.|last2 = Taimre|first2 = T.|last3 = Botev|first3 = Z. I. |title = Handbook of Monte Carlo Methods|year = 2011|publisher = John Wiley & Sons}}
* {{Cite book | last1=Press | first1=WH | last2=Teukolsky | first2=SA | last3=Vetterling | first3=WT | last4=Flannery | first4=BP | year=2007 | title=Numerical Recipes: The Art of Scientific Computing | edition=3rd | publisher=Cambridge University Press |  location=New York | isbn=978-0-521-88068-8}}
* {{Cite book|url=http://www.inference.phy.cam.ac.uk/mackay/itila/book.html|title=Information Theory, Inference and Learning Algorithms | last=MacKay | first=David | publisher=Cambridge University Press|year=2003|isbn=978-0-521-64298-9|chapter=chapter 4.4 Typicality & chapter 29.1| mr=2012999 | chapter-url=http://www.inference.org.uk/itprnn/book.pdf}}
* {{Cite book | last1=Newman | first1=MEJ | last2=Barkema | first2=GT | year=1999 | title=Monte Carlo Methods in Statistical Physics | publisher=Clarendon Press}}
* {{Cite book | last1=Robert | first1=CP | last2=Casella | first2=G | year=2004 | title=Monte Carlo Statistical Methods | publisher=Springer | edition=2nd | isbn=978-1-4419-1939-7}}
{{refend}}

== External links ==
* [https://web.archive.org/web/20140202110318/http://www.cafemath.fr/mathblog/article.php?page=MonteCarlo.php Café math : Monte Carlo Integration] : A blog article describing Monte Carlo integration (principle, hypothesis, confidence interval)
* [http://www.boost.org/doc/libs/release/libs/math/doc/html/math_toolkit/naive_monte_carlo.html Boost.Math : Naive Monte Carlo integration: Documentation for the C++ naive Monte-Carlo routines]
* [https://sites.google.com/view/chremos-group/applets/monte-carlo: Monte Carlo applet applied in statistical physics problems] 

[[Category:Monte Carlo methods]]

{{Sourceattribution|Monte Carlo integration}}