3.2 Monte Carlo methods

Monte Carlo methods were first developed as a method for estimating integrals that could not be evaluated analytically. Although many statistical techniques are now included in the category of ``Monte Carlo methods''[16,17], the method used in this thesis is principally Monte Carlo integration.

3.2.1 Monte Carlo integration

A straightforward application of Monte Carlo is the evaluation of definite integrals. Consider the one dimensional integral

$$E=\int_{a}^{b}f(x)dx \;\;\; .$$ (3.1)

By application of the mean value theorem of calculus, the integral may be approximated by

$$E_N=\frac{(b-a)}{N} \sum_{i=1}^{N} f(x_{i}) \;\;\;,$$ (3.2)

where the points $x_{i}$ fully cover the range of integration. In the limit of a large number of points $N$, $E_N$ tends to the exact value $E$.

A conventional choice for the points $x_{i}$ would be a uniform grid. More accurate methods of numeric quadrature, such as Simpson's rule and Gaussian quadrature[18] use a weighted average of the points:

$$E_N=\frac{(b-a)}{N} \frac{\sum_{i=1}^{N} w_i f(x_{i})}{\sum_{i=1}^{N} w_i} \;\;\;.$$ (3.3)

These methods are highly effective for low dimensional integrals. The computational cost of evaluating an integral, to a fixed accuracy, of dimensionality $d$ increases as $N^d$. An alternate and more efficient approach is to select the points $x_i$ randomly, from a given probability distribution by Monte Carlo methods.

If points are selected at random over the interval $[a,b]$, the Monte Carlo estimate of the integral, equation $\ref{eqn_simple_int}$, becomes

$$E_N=(b-a) \overline{f}=\frac{(b-a)}{N} \sum_{i=1}^{N} f(x_{i}) + O(1/\sqrt(N))\;\;\;,$$ (3.4)

where $\overline{f}$ denotes the mean value of $f$ over the set of sampled points $\{x_i\}$. By the central limit theorem, [19] the set of all possible sums over different $\{x_i\}$ will have a Gaussian distribution. The standard deviation $\sigma_N$ of the different values of $E_N$ is a measure of the uncertainty in the integral's value

$$\sigma_N=\sqrt{\frac{\frac{b-a}{N} \sum_{i=1}^{N} f^2(x_i)-E^{2}_{N}}{N-1}} \;\;\;.$$ (3.5)

The probability that $E$ is within $E_N\pm\sigma_N$ is $\approx 0.68$ and the probability of being with $2\sigma_N$ is $\approx 0.95$. This error decays as $1/\sqrt{N}$ independent of the dimensionality of the integral, unlike grid based methods which have a strong dimensional dependence.

3.2.2 Importance sampling

Simple schemes for Monte Carlo integration can suffer from low efficiency. Many functions of interest have significant weight in only a few regions. For example, most of the contributions to an integral of a simple Gaussian are located near the central peak. In a simple Monte Carlo integration scheme, points are sampled uniformly, wasting considerable effort sampling the tails of the Gaussian. Techniques for overcoming this problem act to increase the density of points in regions of interest and hence improve the overall efficiency. These techniques are called importance sampling.

Points are sampled from a non-uniform distribution $w(x)$, where $w(x)$ is always positive and chosen to approximate $f(x)$ over the region of interest. The integral may now be evaluated by selecting points from the probability distribution $p(x)$,

$$p(x)=\frac{w(x)}{\int_{a}^{b}w(x)dx} \;\;\;,$$ (3.6)

so that the integral may be evaluated

$$E=\int_{a}^{b}g(x)p(x)dx \approx \frac{1}{N}\sum_{i=1}^{N} g(x_i)\;\;\; ,$$ (3.7)

where $g(x)=f(x)/p(x)$ and the points $x_i$ are generated from the probability distribution $p(x)$.

Provided $f(x)$ is well described by $w(x)$, $g(x)$ remains close to unity, so that the standard deviation defined in equation 3.5 is greatly reduced. (The best possible choice for $p(x)$ is $\vert f(x)/g(x)\vert$.) The integration procedure is therefore more efficient, although the error still decays as $1/\sqrt{N}$.

3.2.3 The Metropolis algorithm

In this thesis, the Metropolis algorithm[20] is used in the importance sampling of multi-dimensional integrals. The Metropolis algorithm generates a random walk of points distributed according to a required probability distribution. From an initial ``position'' in phase or configuration space, a proposed ``move'' is generated and the move either accepted or rejected according to the Metropolis algorithm. By taking a sufficient number of trial steps all of phase space is explored and the Metropolis algorithm ensures that the points are distributed according to the required probability distribution.

The Metropolis algorithm may be derived by considering the probability density, $\rho$ at two points in configuration space $X$ and $X^\prime$. Configuration space is specified by the integration variables and limits. For one dimensional integration, $X=\{x\}$ for $x \in [a,b]$. Averaged over many trial steps, the average number of samples at $X$ and $X^\prime$ should be equal, $\rho(X)=\rho(X^\prime)$.

The probability for a trial move from $X$ to $X^\prime$ defined $T(X\rightarrow X^\prime)$ and we require

$$T(X\rightarrow X^\prime)=T(X^\prime\rightarrow X) \;\;\;.$$ (3.8)

If the probability of accepting a move from $X$ to $X^\prime$ is $P(X\rightarrow X^\prime)$ then the total probability of accepting a move from $X$ to $X^\prime$ is $P(X\rightarrow X^\prime)T(X\rightarrow X^\prime)$. Therefore, at equilibrium

$$\rho(X)P(X\rightarrow X^\prime)T(X\rightarrow X^\prime)= \rh... ...prime)P(X^\prime\rightarrow X)T(X^\prime\rightarrow X) \;\;\;.$$ (3.9)

The Metropolis algorithm corresponds to choosing

$$P(X\rightarrow X^\prime)=\mathrm{min}\left\{ 1, \frac{\rho(X^\prime)}{\rho(X)} \right\} \;\;\;.$$ (3.10)

For the Metropolis algorithm to be valid, it is essential that the random walk is ergodic, that is any point $X^\prime$ in configuration space may be reached from any other point $X$. In some applications of the Metropolis algorithm, parts of configuration space may be difficult to reach. Long simulations or a modification of the algorithm are then necessary. Many methods have been developed to cope with this ``slowing down'', [16] but for the applications presented here the Metropolis algorithm has been found both adequate and sufficient for the problems investigated.