Monte-Carlo Integration

Often we are faced with integrals which cannot be done analytically. Especially in the case of multidimensional integrals the simplest methods of discretisation can become prohibitively expensive. For example, the error in a trapezium rule calculation of a $d$-dimensional integral falls as $N^{-2/d}$, where $N$ is the number of different values of the integrand used. In a Monte-Carlo calculation the error falls as $N^{-1/2}$ independently of the dimension. Hence for $d > 4$ Monte-Carlo integration will usually converge faster.

We consider the expression for the average of a statistic, $f(x)$ when $x$ is a random number distributed according to a distribution $p(x)$, then

$$\left<f(x)\right> = \int p(x) f(x) \d x,$$ (4.6)

which is just a generalisation of the well known results for (e.g. ) $\left<x\right>$ or $\left<x^2\right>$, where we are using the notation $<>$ to denote averaging. Now consider an integral of the sort which might arise while using Laplace transforms.

$$\int_0^\infty \exp(-x) f(x) \d x .$$ (4.7)

This integral can be evaluated by generating a set of $N$ random numbers, $\{x_1,\ldots,x_N\}$, from a Poisson distribution, $p(x) = \exp(-x)$, and calculating the mean of $f(x)$ as

$${1\over N}\sum_{i=1}^N f(x_i).$$ (4.8)

The error in this mean is evaluated as usual by considering the corresponding standard error of the mean

$$\sigma^2 = {1\over N-1}\left(\left<f^2(x)\right> - \left<f(x)\right>^2\right).$$ (4.9)

As mentioned earlier, Monte-Carlo integration can be particularly efficient in the case of multi-dimensional integrals. However this case is particularly susceptible to the flaws in random number generators. It is a common feature that when a random set of coordinates in a $d$-dimensional space is generated (i.e. a set of $d$ random numbers), the resulting distribution contains (hyper)planes on which the probability is either significantly higher or lower than expected.