Editing Monte Carlo integration (section)

== Overview ==
In numerical integration, methods such as the [[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 <math>\mathbb{R}^m</math>, 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]] 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)=\mathrm{E}(\sigma_N^2) \equiv \frac{1}{N-1} \sum_{i=1}^N \mathrm{E}\left [\left (f(\overline{\mathbf{x}}_i) - \langle f \rangle \right )^2\right ]. </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{\mathrm{E}(\sigma_N^2)}{N}.</math>

Since the sequence
<math display="block"> \left \{ \mathrm{E}(\sigma_1^2), \mathrm{E}(\sigma_2^2), \mathrm{E}(\sigma_3^2), \ldots \right \} </math>
is bounded due to being identically equal to ''Var(f)'', as long as this is assumed finite, 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{\sqrt{\mathrm{Var}(f)}}{\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/C++ example ===
<syntaxhighlight lang="c">
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
int main() {
    // Initialize the number of counts to 0, setting the total number to be 100000 in the loop.
    int throws = 99999, insideCircle = 0;
    double randX, randY, pi;

    srand(time(NULL));

    // Checks for each random pair of x and y if they are inside circle of radius 1.
    for (int i = 0; i < throws; i++) {
        randX = rand() / (double) RAND_MAX;
        randY = rand() / (double) RAND_MAX;
        if (randX * randX + randY * randY < 1) {
            insideCircle++;
        }
    }

    // Calculating pi and printing.
    pi = 4.0 * insideCircle / throws;
    printf("%lf\n", pi);
}
</syntaxhighlight>

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

rng = np.random.default_rng(0)

throws = 2000
radius = 1

# Choose random X and Y data centered around 0,0
x = rng.uniform(-radius, radius, throws)
y = rng.uniform(-radius, radius, throws)

# Count the times (x, y) is inside the circle,
# which happens when sqrt(x^2 + y^2) <= radius.
inside_circle = np.count_nonzero(np.hypot(x, y) <= radius)

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

=== 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 [[Mathematica]]:

<syntaxhighlight 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*)
</syntaxhighlight>