Monte Carlo Integration

A Simple Example

Consider the figure Figure 1. The square represents the domain of some function \(g\). That is, \(g\) takes some 2D point \(\vec{x}\) from the square and returns a value.

Within the subdomain \(A\), it turns out that \(g\) returns \(3\). Similarly, in subdomain \(B\), \(g(\vec{x})=4\) and in subdomain \(C\), \(g(\vec{x})=5\).

Say that the whole square (domain) is \(40\) units on a side. Thus, it has an area of \(40 \times 40=1600\). The dividing lines make it so that subdomain \(A\) has area \(20 \times 20=400\), subdomain \(B\) has area \(20 \times 40=800\), and subdomain \(C\) has area \(20 \times 20=400\).

def g(x):
	if   x[0]>=20.0: return 4.0 #Region B
	elif x[1]>=20.0: return 3.0 #Region A
	else           : return 5.0 #Region C

Let's suppose we're interested in the integral of \(g\) over this domain.

In this case, we can figure it out in closed-form, just integrating each subdomain separately as the value times the area:

We'll be able to use this to check our work, but what if we couldn't do this? There are many complicated functions that don't have elementary integrals, and the domains can be very complicated indeed. Both cases are common in graphics.

Solution Via Monte Carlo

Let's now add the Monte Carlo bit. The first key idea of Monte Carlo is estimating the value of the whole integral, but only testing a single randomly chosen point.

Say we pick our random point, and it happens to be within \(C\), so \(g\) is \(5\). That \(5\) somehow needs to tell us a guess for the entire integral. But how?

The second key idea of Monte Carlo integration is probability density. We know that the probability of all possible events sums to \(1\), right? We can now distribute (divide) that amount of probability over the domain to get a notion of probability per area—probability density. Here, the area is \(40 \times 40=1600\), so the probability density, or probability density function, is \(1/1600\). Note that we could do a more complicated distribution than this simple uniform division, but we'll keep it simple here.

The key statement of Monte Carlo integration is this:

That is, we can choose a random \(x_k\) from the domain, sampled according to \(\pdf(x)\), divide the integrand by the probability density of that choice \(\pdf(x_k)\), and on average ('in expectation'; the \(E[\cdots]\) is the expected value), the result will be the correct value of the integral.

For our example, our sample was \(5\) and the PDF was \(1/1600\). Thus, an estimate for the integral is \(5~/~(1/1600)=5\cdot 1600=8000\).

This estimate is too high, but if we repeat this procedure with other random samples, we'll get a bunch of different estimates—\(8000, 4800, 6400, 8000, 6400, 6400, 4800, 6400, 4800, 6400, 6400, 8000\), etc.—and on average those estimates will be the right answer! This is what we mean when we say 'expected value'.

The following simple Python code demonstrates this:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
import random

def g(x):
	if   x[0]>=20.0: return 4.0 #Region B
	elif x[1]>=20.0: return 3.0 #Region A
	else           : return 5.0 #Region C

def get_random_x():
	return random.random()*40.0, random.random()*40.0
def monte_carlo_estimate():
	x_k = get_random_x()
	pdf = 1.0 / 1600.0
	return g(x_k) / pdf

#Average 1,000,000 estimates together, showing that we get the right answer in expectation
N = 1_000_000
print(sum([ monte_carlo_estimate() for k in range(N)])/N)

Try it! This code simply averages a bunch of separate Monte Carlo estimates together, producing (for example) the output \(6399.4896\), which is just about right. In fact, as you increase the number of samples, the average gets closer and closer to correct.

Going Further

This toy example might seem slightly contrived, but Monte Carlo integration allows us to solve extremely complicated problems.

In pathtracing, we are interested in integrating the amount of light that flows from all lights to the camera: a set of paths called 'path-space'. Path-space is difficult to even imagine—it is the set of all possible paths, including those with an infinite number of reflections, with points varying infinitesimally over the surface of objects, scattering in any which way through volumes. There is no chance that we could integrate this with a traditional approach. However, the Monte Carlo method allows us to pick a single sample from this space—a single path—and use it to estimate the entire monster. Average enough estimates, and you can get the correct answer—and a photorealistic image.

I have a tutorial about applying Monte Carlo integration to the rendering equation.

It may also be interesting to check out what that looks like in a bare-bores teaching pathtracer.