Monte Carlo integration

Relative error as a function of the number of samples, showing the scaling ${\displaystyle {\tfrac {1}{\sqrt {N}}}}$

A paradigmatic example of a Monte Carlo integration is the estimation of π. Consider the function

${\displaystyle H\left(x,y\right)={\begin{cases}1&{\text{if }}x^{2}+y^{2}\leq 1\\0&{\text{else}}\end{cases}}}$

and the set Ω = [−1,1] × [−1,1] with V = 4. Notice that

${\displaystyle I_{\pi }=\int _{\Omega }H(x,y)dxdy=\pi .}$

Thus, a crude way of calculating the value of π with Monte Carlo integration is to pick N random numbers on Ω and compute

${\displaystyle Q_{N}=4{\frac {1}{N}}\sum _{i=1}^{N}H(x_{i},y_{i})}$

In the figure on the right, the relative error ${\displaystyle {\tfrac {Q_{N}-\pi }{\pi }}}$ is measured as a function of N, confirming the ${\displaystyle {\tfrac {1}{\sqrt {N}}}}$.

Keep in mind that a true random number generator should be used.

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;

Made in 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)

The code below describes a process of integrating the function

${\displaystyle f(x)={\frac {1}{1+\sinh(2x)\log(x)^{2}}}}$

from ${\displaystyle 0.8<x<3}$ using the Monte-Carlo 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*)

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 ${\displaystyle E_{a}(f)}$ and ${\displaystyle E_{b}(f)}$ and variances ${\displaystyle \sigma _{a}^{2}(f)}$ and ${\displaystyle \sigma _{b}^{2}(f)}$, the variance Var(f) of the combined estimate

${\displaystyle E(f)={\tfrac {1}{2}}\left(E_{a}(f)+E_{b}(f)\right)}$

is given by,

${\displaystyle \mathrm {Var} (f)={\frac {\sigma _{a}^{2}(f)}{4N_{a}}}+{\frac {\sigma _{b}^{2}(f)}{4N_{b}}}}$

It can be shown that this variance is minimized by distributing the points such that,

${\displaystyle {\frac {N_{a}}{N_{a}+N_{b}}}={\frac {\sigma _{a}}{\sigma _{a}+\sigma _{b}}}}$

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_a and N_b. 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 provides a very important tool to perform Monte-Carlo integration.[3][8] The main result of importance sampling to this method is that the uniform sampling of ${\displaystyle {\overline {\mathbf {x} }}}$ is a particular case of a more generic choice, on which the samples are drawn from any distribution ${\displaystyle p({\overline {\mathbf {x} }})}$. The idea is that ${\displaystyle p({\overline {\mathbf {x} }})}$ 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

${\displaystyle p({\overline {\mathbf {x} }}):\qquad {\overline {\mathbf {x} }}_{1},\cdots ,{\overline {\mathbf {x} }}_{N}\in V,}$

the estimator for I is given by[3]

${\displaystyle Q_{N}\equiv {\frac {1}{N}}\sum _{i=1}^{N}{\frac {f({\overline {\mathbf {x} }}_{i})}{p({\overline {\mathbf {x} }}_{i})}}}$

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

The Metropolis–Hastings algorithm is one of the most used algorithms to generate ${\displaystyle {\overline {\mathbf {x} }}}$ from ${\displaystyle p({\overline {\mathbf {x} }})}$,[3] thus providing an efficient way of computing integrals.

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.[9] In order to avoid the number of histogram bins growing like K^d, the probability distribution is approximated by a separable function:

${\displaystyle g(x_{1},x_{2},\ldots )=g_{1}(x_{1})g_{2}(x_{2})\ldots }$

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.[9]