Monte-Carlo integration

Consider a one-dimensional integral: $\int_{x_l}^{x_h} f(x)\,dx$. We can evaluate this integral numerically by dividing the interval $x_l$ to $x_h$ into $N$ identical subdivisions of width

$$h = \frac{x_h-x_l}{N}.$$ (326)

Let $x_i$ be the midpoint of the $i$th subdivision, and let $f_i = f(x_i)$. Our approximation to the integral takes the form

$$\int_{x_l}^{x_h} f(x)\,dx\simeq \sum_{i=1}^{N} f_i\,h$$ (327)

This integration method--which is known as the midpoint method--is not particularly accurate, but is very easy to generalize to multi-dimensional integrals.

What is the error associated with the midpoint method? Well, the error is the product of the error per subdivision, which is $O(h^2)$, and the number of subdivisions, which is $O(h^{-1})$. The error per subdivision follows from the linear variation of $f(x)$ within each subdivision. Thus, the overall error is $O(h^2)\times O(h^{-1}) = O(h)$. Since, $h\propto N^{-1}$, we can write

$$\int_{x_l}^{x_h} f(x)\,dx\simeq \sum_{i=1}^{N} f_i\,h + O(N^{-1}).$$ (328)

Let us now consider a two-dimensional integral. For instance, the area enclosed by a curve. We can evaluate such an integral by dividing space into identical squares of dimension $h$, and then counting the number of squares, $N$ (say), whose midpoints lie within the curve. Our approximation to the integral then takes the form

$$A \simeq N\,h^2.$$ (329)

This is the two-dimensional generalization of the midpoint method.

What is the error associated with the midpoint method in two-dimensions? Well, the error is generated by those squares which are intersected by the curve. These squares either contribute wholly or not at all to the integral, depending on whether their midpoints lie within the curve. In reality, only those parts of the intersected squares which lie within the curve should contribute to the integral. Thus, the error is the product of the area of a given square, which is $O(h^2)$, and the number of squares intersected by the curve, which is $O(h^{-1})$. Hence, the overall error is $O(h^2)\times O(h^{-1}) = O(h)= O(N^{-1/2})$. It follows that we can write

$$A = N\,h^2 + O(N^{-1/2}).$$ (330)

Let us now consider a three-dimensional integral. For instance, the volume enclosed by a surface. We can evaluate such an integral by dividing space into identical cubes of dimension $h$, and then counting the number of cubes, $N$ (say), whose midpoints lie within the surface. Our approximation to the integral then takes the form

$$V \simeq N\,h^3.$$ (331)

This is the three-dimensional generalization of the midpoint method.

What is the error associated with the midpoint method in three-dimensions? Well, the error is generated by those cubes which are intersected by the surface. These cubes either contribute wholly or not at all to the integral, depending on whether their midpoints lie within the surface. In reality, only those parts of the intersected cubes which lie within the surface should contribute to the integral. Thus, the error is the product of the volume of a given cube, which is $O(h^3)$, and the number of cubes intersected by the surface, which is $O(h^{-2})$. Hence, the overall error is $O(h^3)\times O(h^{-2}) = O(h)= O(N^{-1/3})$. It follows that we can write

$$V = N\,h^3 + O(N^{-1/3}).$$ (332)

Let us, finally, consider using the midpoint method to evaluate the volume, $V$, of a $d$-dimensional hypervolume enclosed by a $(d-1)$-dimensional hypersurface. It is clear, from the above examples, that

$$V = N\,h^d + O(N^{-1/d}),$$ (333)

where $N$ is the number of identical hypercubes into which the hypervolume is divided. Note the increasingly slow fall-off of the error with $N$ as the dimensionality, $d$, becomes greater. The explanation for this phenomenon is quite simple. Suppose that $N=10^6$. With $N=10^6$ we can divide a unit line into (identical) subdivisions whose linear extent is $10^{-6}$, but we can only divide a unit area into subdivisions whose linear extent is $10^{-3}$, and a unit volume into subdivisions whose linear extent is $10^{-2}$. Thus, for a fixed number of subdivisions the grid spacing (and, hence, the integration error) increases dramatically with increasing dimension.

Let us now consider the so-called Monte-Carlo method for evaluating multi-dimensional integrals. Consider, for example, the evaluation of the area, $A$, enclosed by a curve, $C$. Suppose that the curve lies wholly within some simple domain of area $A'$, as illustrated in Fig. 97. Let us generate $N'$ points which are randomly distributed throughout $A'$. Suppose that $N$ of these points lie within curve $C$. Our estimate for the area enclosed by the curve is simply

$$A = \frac{N}{N'}\,A'.$$ (334)

Figure 97: The Monte-Carlo integration method.

What is the error associated with the Monte-Carlo integration method? Well, each point has a probability $p=A/A'$ of lying within the curve. Hence, the determination of whether a given point lies within the curve is like the measurement of a random variable $x$ which has two possible values: 1 (corresponding to the point being inside the curve) with probability $p$, and 0 (corresponding to the point being outside the curve) with probability $1-p$. If we make $N'$ measurements of $x$ (i.e., if we scatter $N'$ points throughout $A'$) then the number of points lying within the curve is

$$N = \sum_{i=1,N'} x_i,$$ (335)

where $x_i$ denotes the $i$th measurement of $x$. Now, the mean value of $N$ is

$$\bar{N} = \sum_{i=1,N'} \bar{x}= N'\,\bar{x},$$ (336)

where

$$\bar{x} = 1\times p + 0\times (1-p) = p.$$ (337)

Hence,

$$\bar{N} = N'\,p = N'\,\frac{A}{A'},$$ (338)

which is consistent with Eq. (334). We conclude that, on average, a measurement of $N$ leads to the correct answer. But, what is the scatter in such a measurement? Well, if $\sigma$ represents the standard deviation of $N$ then we have

$$\sigma^2 = \overline{(N-\bar{N})^2},$$ (339)

which can also be written

$$\sigma^2 = \sum_{i,j,=1,N'} \overline{(x_i-\bar{x})(x_j-\bar{x})}.$$ (340)

However, $\overline{(x_i-\bar{x})(x_j-\bar{x})}$ equals $\overline{(x-\bar{x})^2}$ if $i=j$, and equals zero, otherwise, since successive measurements of $x$ are uncorrelated. Hence,

$$\sigma^2 = N'\,\overline{(x-\bar{x})^2}.$$ (341)

Now,

$$\overline{(x-\bar{x})^2}=\overline{(x^2-2\,x\,\bar{x}+\bar{x}^2)} = \overline{x^2} - \bar{x}^2,$$ (342)

and

$$\overline{x^2} = 1^2\times p + 0^2\times (1-p) = p.$$ (343)

Thus,

$$\overline{(x-\bar{x})^2}= p -p^2 =p\,(1-p),$$ (344)

giving

$$\sigma = \sqrt{N'\,p\,(1-p)}.$$ (345)

Finally, since the likely values of $N$ lie in the range $N=\bar{N}\pm\sigma$, we can write

$$N = N'\,\frac{A}{A'} \pm \sqrt{N'\,\frac{A}{A'}\left(1-\frac{A}{A'}\right)}.$$ (346)

It follows from Eq. (334) that

$$A = A'\,\frac{N}{N'} \pm \frac{\sqrt{A\,(A'-A)}}{\sqrt{N'}}.$$ (347)

In other words, the error scales like $(N')^{-1/2}$.

The Monte-Carlo method generalizes immediately to $d$-dimensions. For instance, consider a $d$-dimensional hypervolume $V$ enclosed by a $(d-1)$-dimensional hypersurface $A$. Suppose that $A$ lies wholly within some simple hypervolume $V'$. We can generate $N'$ points randomly distributed throughout $V'$. Let $N$ be the number of these points which lie within $A$. It follows that our estimate for $V$ is simply

$$V = \frac{N}{N'}\,V'.$$ (348)

Now, there is nothing in our derivation of Eq. (347) which depends on the fact that the integral in question is two-dimensional. Hence, we can generalize this equation to give

$$V = V'\,\frac{N}{N'} \pm \frac{\sqrt{V\,(V'-V)}}{\sqrt{N'}}.$$ (349)

We conclude that the error associated with Monte-Carlo integration always scales like $(N')^{-1/2}$, irrespective of the dimensionality of the integral.

We are now in a position to compare and contrast the midpoint and Monte-Carlo methods for evaluating multi-dimensional integrals. In the midpoint method, we fill space with an evenly spaced mesh of $N$ (say) points (i.e., the midpoints of the subdivisions), and the overall error scales like $N^{-1/d}$, where $d$ is the dimensionality of the integral. In the Monte-Carlo method, we fill space with $N$ (say) randomly distributed points, and the overall error scales like $N^{-1/2}$, irrespective of the dimensionality of the integral. For a one-dimensional integral ($d=1$), the midpoint method is more efficient than the Monte-Carlo method, since in the former case the error scales like $N^{-1}$, whereas in the latter the error scales like $N^{-1/2}$. For a two-dimensional integral ($d=2$), the midpoint and Monte-Carlo methods are both equally efficient, since in both cases the error scales like $N^{-1/2}$. Finally, for a three-dimensional integral ($d=3$), the midpoint method is less efficient than the Monte-Carlo method, since in the former case the error scales like $N^{-1/3}$, whereas in the latter the error scales like $N^{-1/2}$. We conclude that for a sufficiently high dimension integral the Monte-Carlo method is always going to be more efficient than an integration method (such as the midpoint method) which relies on a uniform grid.

Figure 98: Example calculation: volume of unit-radius 2-dimensional sphere enclosed in a close-fitting 2-dimensional cube.

Up to now, we have only considered how the Monte-Carlo method can be employed to evaluate a rather special class of integrals in which the integrand function can only take the values 0 or 1. However, the Monte-Carlo method can easily be adapted to evaluate more general integrals. Suppose that we wish to evaluate $\int f\,dV$, where $f$ is a general function and the domain of integration is of arbitrary dimension. We proceed by randomly scattering $N$ points throughout the integration domain and calculating $f$ at each point. Let $x_i$ denote the $i$th point. The Monte-Carlo approximation to the integral is simply

$$\int f\,dV = \frac{1}{N}\sum_{i=1,N} f(x_i) + O\left(\frac{1}{\sqrt{N}}\right).$$ (350)

Figure 99: The integration error, $\epsilon$, versus the number of grid-points, $N$, for three integrals evaluated using the midpoint method. The integrals are the area of a unit-radius circle (solid curve), the volume of a unit-radius sphere (dotted curve), and the volume of a unit-radius 4-sphere (dashed curve).

We end this section with an example calculation. Let us evaluate the volume of a unit-radius $d$-dimensional sphere, where $d$ runs from 2 to 4, using both the midpoint and Monte-Carlo methods. For both methods, the domain of integration is a cube, centred on the sphere, which is such that the sphere just touches each face of the cube, as illustrated in Fig. 98.

Figure 100: The integration error, $\epsilon$, versus the number of points, $N$, for three integrals evaluated using the Monte-Carlo method. The integrals are the area of a unit-radius circle (solid curve), the volume of a unit-radius sphere (dotted curve), and the volume of a unit-radius 4-sphere (dashed curve).

Figure 99 shows the integration error associated with the midpoint method as a function of the number of grid-points, $N$. It can be seen that as the dimensionality of the integral increases the error falls off much less rapidly as $N$ increases.

Figure 100 shows the integration error associated with the Monte-Carlo method as a function of the number of points, $N$. It can be seen that there is very little change in the rate at which the error falls off with increasing $N$ as the dimensionality of the integral varies. Hence, as the dimensionality, $d$, increases the Monte-Carlo method eventually wins out over the midpoint method.