Monte Carlo Method for Integrals

We present here the simplest application of Monte Carlo techniques. We will use Monte Carlo to approximate the integral

$$\int_a^b f'(x') \, dx'.$$

Without loss of generality, we will be further assuming that the function $f'(x')$ can be transformed to $0 \leq f(x) \leq 1$, and that the limits can be mapped to the interval $0$ and $1$, respectively, so that the task is now converted to approximating the integral

$$I \equiv \int_0^1 f(x) \, dx.$$

Suppose we choose $N$ pairs of points $(x_i,y_i)$ with uniform distribution (see Figure 51) and define

Figure 51: Domain of $(x,y)$ pairs.

$$z_i = \begin{cases} 0 \qquad y_i > f(x_i) \\ 1 \qquad y_i \leq f(x_i). \end{cases}$$

See Figure 52. Then, if $n = \sum_i z_i$, we have $n/N \approx I$, or more precisely,

$$I = \frac{n}{N} + O(1/\sqrt N).$$

Figure 52: The value of $z_i$ is $0$ for $y_i > f(x_i)$, and $1$ if in the shaded region.

Can also consider $I$ as the mean value of the probability density function $f(\xi)$ with $\xi$ from a uniformly distributed density function. Then, the estimate of the mean value is

$$I \approx \frac{1}{N} \sum_{i=1}^N f(\xi_i).$$

The basic theorem for Monte Carlo estimates the integral in $\mathbf{R}^n$ space of $f$ as

$$\int_V f(y_1,y_2,\ldots,y_n) dy_1 dy_2 \ldots dy_n \approx V <f> \pm \sqrt{\frac{<f^2>-<f>^2}{N}},$$ (50)

where

$$\begin{eqnarray*} <f> &\equiv & \frac{1}{N} \sum_{i=1}^N f(y_i) \\ <f^2> &\equiv & \frac{1}{N} \sum_{i=1}^N f^2(y_i) \end{eqnarray*}$$

so we recognize the last term on the right hand side of (51) as the volume $V$ times the standard deviation.

This is not a rigorous bound and it assumes that the error is Gaussian distributed. In fact, it also needs to be pointed out that the the approximation also depends on the quality of the random number generator and the actual creation of good and fast random number generators is a matter of current research.

The important thing about (51) is that it says that the integral converges, but very slowly: $O(1/\sqrt N)$, where $N$ is the number of $(x_i,y_i)$ pairs. So it is typically impractical, i.e. computationally expensive for integrals in one or two dimensions. However, it is competitive for higher dimensional integrals. In particular, there are many applications in estimation theory in which the integral is a Feynman path integral, which is an integral over ALL possible paths.