Say It With Science

Monte Carlo Integration

In the course of theoretical and computational science, we sometimes encounter integrals that cannot be solved exactly, e.g. integrals without analytic solutions. One tool for numerically estimating these integrals uses the Monte Carlo method, and is called Monte Carlo (MC) integration.

Roughly, an MC integration finds the average value of a function over a region of integration and multiplies it by the hypervolume of that region. This is just an extension of the mean value theorem any first year calculus student learns about integrals over R. The first part is easy: to find the average value of the integrand f : Rⁿ→R, we simply select M points x_i uniformly distributed throughout the region of integration Ω. Then the average value ⟨f⟩ = (1/M) Σ_if (x_i). The more difficult problem: how do we find the volume of Ω?

Luckily, the MC method makes this easy as well, and thereby becomes a useful tool for computing integrals on domains of any finite dimension. For a given region of integration Ω, it is usually easy to define some hypercube H such that Ω ⊂ H. Now instead of only trying to generate M points inside of Ω, we generate N points inside H (a much easier task, since H is naturally uniform in each dimension). Some of them will fall within Ω anyway; we again define their count to be M. But some of them will not, such that M < N. Since we can easily compute the volume of H (call it V), the volume of Ω is roughly (M/N)⋅V.

In the image above, we’ve chosen the unit circle as our region Ω. Of the 25,000 points plotted, 19,610 fell inside the unit circle and 5,391 fell outside, so M/N = 19 610/(5391+19 610) ≈ 0.7844. Since our “hypercube” is a square of “hypervolume” (or, area) V = 4, (M/N)⋅V = 3.13747 ≈ π, as we’d expect.

Suppose we want to integrate a function f over the unit circle we’ve pictured above; at each point x_i inside the circle, we’d evaluate f(x_i), and then we’d let ⟨f⟩ be the average of all these various calculations. (Notice that we don’t even evaluate points outside the circle.) Now we have everything we need: the volume of the region and the average value of the function over that region. Their product approximates the integral.

The pseudocode for a simple MC integration might look something like this:

Program MonteCarlo(f, Ω, H, N):
    M ← 0
    F ← 0
    For i, 0 → N:
        p ← RandomPoint(H)
        If p ∈ Ω:
            F ← F + f(p)
            M ← M + 1
    avg ← F/M
    vol ← (M/N) ∗ Volume(H)
    Return avg ∗ vol

With a higher N, the precision of the integration will be higher. And the easiness of finding an n-dimensional hypercube volume makes Monte Carlo integration ideal for multidimensional integrals. But it has its drawbacks. Because the points are randomly distributed over H, integrals that have erratic behavior near specific points may not have those values proportionally represented, and highly oscillatory integrals may be difficult to properly average. Adaptive Monte Carlo methods can sometimes identify these points and recursively focus on the more extreme values of the function to get a better estimate. Finding a tight fit for H is also important; any large H will work in theory, but in practice the best results will come for a hypercube that fits tightly around the region of integration.