2.1 Monte Carlo: Basics

Because Monte Carlo integration is based on randomization, we will start this chapter with a brief review of ideas from probability and statistics that provide the foundations of the approach. Doing so will allow us to introduce the basic Monte Carlo algorithm as well as mathematical tools for evaluating its error.

2.1.1 Background and Probability Review

We will start by defining some terms and reviewing basic ideas from probability. We assume that the reader is already familiar with basic probability concepts; readers needing a more complete introduction to this topic should consult a textbook such as Sheldon Ross’s Introduction to Probability Models (2002).

A random variable upper X is a value chosen by some random process. We will generally use capital letters to denote random variables, with exceptions made for a few Greek symbols that represent special random variables. Random variables are always drawn from some domain, which can be either discrete (e.g., a fixed, finite set of possibilities) or continuous (e.g., the real numbers bold upper R Superscript ). Applying a function  f to a random variable  upper X results in a new random variable upper Y equals f left-parenthesis upper X right-parenthesis .

For example, the result of a roll of a die is a discrete random variable sampled from the set of events upper X Subscript i Baseline element-of StartSet 1 comma 2 comma 3 comma 4 comma 5 comma 6 EndSet . Each event has a probability p Subscript i Baseline equals one-sixth , and the sum of probabilities sigma-summation p Subscript i is necessarily one. A random variable like this one that has the same probability for all potential values of it is said to be uniform. A function p left-parenthesis upper X right-parenthesis that gives a discrete random variable’s probability is termed a probability mass function (PMF), and so we could equivalently write p left-parenthesis upper X right-parenthesis equals one-sixth in this case.

Two random variables are independent if the probability of one does not affect the probability of the other. In this case, the joint probability p left-parenthesis upper X comma upper Y right-parenthesis of two random variables is given by the product of their probabilities:

p left-parenthesis upper X comma upper Y right-parenthesis equals p left-parenthesis upper X right-parenthesis p left-parenthesis upper Y right-parenthesis period

For example, two random variables representing random samples of the six sides of a die are independent.

For dependent random variables, one’s probability affects the other’s. Consider a bag filled with some number of black balls and some number of white balls. If we randomly choose two balls from the bag, the probability of the second ball being white is affected by the color of the first ball since its choice changes the number of balls of one type left in the bag. We will say that the second ball’s probability is conditioned on the choice of the first one. In this case, the joint probability for choosing two balls upper X and upper Y is given by

p left-parenthesis upper X comma upper Y right-parenthesis equals p left-parenthesis upper X right-parenthesis p left-parenthesis upper Y vertical-bar upper X right-parenthesis comma
(2.1)

where p left-parenthesis upper Y vertical-bar upper X right-parenthesis is the conditional probability of upper Y given a value of upper X .

In the following, it will often be the case that a random variable’s probability is conditioned on many values; for example, when choosing a light source from which to sample illumination, the BVHLightSampler in Section 12.6.3 considers the 3D position of the receiving point and its surface normal, and so the choice of light is conditioned on them. However, we will often omit the variables that a random variable is conditioned on in cases where there are many of them and where enumerating them would obscure notation.

A particularly important random variable is the canonical uniform random variable, which we will write as  xi Subscript . This variable takes on all values in its domain left-bracket 0 comma 1 right-parenthesis independently and with uniform probability. This particular variable is important for two reasons. First, it is easy to generate a variable with this distribution in software—most runtime libraries have a pseudo-random number generator that does just that. Second, we can take the canonical uniform random variable xi Subscript and map it to a discrete random variable, choosing upper X Subscript i if

sigma-summation Underscript j equals 1 Overscript i minus 1 Endscripts p Subscript j Baseline less-than-or-equal-to xi Subscript Baseline less-than sigma-summation Underscript j equals 1 Overscript i Endscripts p Subscript j Baseline period
(2.2)

For lighting applications, we might want to define the probability of sampling illumination from each light in the scene based on its power normal upper Phi Subscript i relative to the total power from all sources:

p Subscript i Baseline equals StartFraction normal upper Phi Subscript i Baseline Over sigma-summation Underscript j Endscripts normal upper Phi Subscript j Baseline EndFraction period

Notice that these p Subscript i values also sum to 1. Given such per-light probabilities, xi Subscript could be used to select a light source from which to sample illumination.

The cumulative distribution function (CDF) upper P left-parenthesis x right-parenthesis of a random variable is the probability that a value from the variable’s distribution is less than or equal to some value x :

upper P left-parenthesis x right-parenthesis equals probability left-brace upper X less-than-or-equal-to x right-brace period
(2.3)

For the die example, upper P left-parenthesis 2 right-parenthesis equals one-third , since two of the six possibilities are less than or equal to 2.

Continuous random variables take on values over ranges of continuous domains (e.g., the real numbers, directions on the unit sphere, or the surfaces of shapes in the scene). Beyond xi Subscript , another example of a continuous random variable is the random variable that ranges over the real numbers between 0 and 2, where the probability of its taking on any particular value x is proportional to the value 2 minus x : it is twice as likely for this random variable to take on a value around 0 as it is to take one around 1, and so forth.

The probability density function (PDF) formalizes this idea: it describes the relative probability of a random variable taking on a particular value and is the continuous analog of the PMF. The PDF p left-parenthesis x right-parenthesis is the derivative of the random variable’s CDF,

p left-parenthesis x right-parenthesis equals StartFraction normal d upper P left-parenthesis x right-parenthesis Over normal d x EndFraction period

For uniform random variables, p left-parenthesis x right-parenthesis is a constant; this is a direct consequence of uniformity. For xi Subscript we have

p left-parenthesis x right-parenthesis equals StartLayout Enlarged left-brace 1st Row 1st Column 1 2nd Column x element-of left-bracket 0 comma 1 right-parenthesis 2nd Row 1st Column 0 2nd Column otherwise period EndLayout

PDFs are necessarily nonnegative and always integrate to 1 over their domains. Note that their value at a point x is not necessarily less than 1, however.

Given an interval left-bracket a comma b right-bracket in the domain, integrating the PDF gives the probability that a random variable lies inside the interval:

probability left-brace x element-of left-bracket a comma b right-bracket right-brace equals integral Subscript a Superscript b Baseline p left-parenthesis x right-parenthesis normal d x equals upper P left-parenthesis b right-parenthesis minus upper P left-parenthesis a right-parenthesis period

This follows directly from the first fundamental theorem of calculus and the definition of the PDF.

2.1.2 Expected Values

The expected value upper E Subscript p Baseline left-bracket f left-parenthesis x right-parenthesis right-bracket of a function f is defined as the average value of the function over some distribution of values p left-parenthesis x right-parenthesis over its domain upper D . It is defined as

upper E Subscript p Baseline left-bracket f left-parenthesis x right-parenthesis right-bracket equals integral Underscript upper D Endscripts f left-parenthesis x right-parenthesis p left-parenthesis x right-parenthesis normal d x period
(2.4)

As an example, consider finding the expected value of the cosine function between 0 and pi , where p is uniform. Because the PDF p left-parenthesis x right-parenthesis must integrate to 1 over the domain, p left-parenthesis x right-parenthesis equals 1 slash pi , so

upper E left-bracket cosine x right-bracket equals integral Subscript 0 Superscript pi Baseline StartFraction cosine x Over pi EndFraction normal d x equals StartFraction 1 Over pi EndFraction left-parenthesis sine pi minus sine 0 right-parenthesis equals 0 comma

which is precisely the expected result. (Consider the graph of cosine x over left-bracket 0 comma pi right-bracket to see why this is so.)

The expected value has a few useful properties that follow from its definition:

StartLayout 1st Row 1st Column upper E left-bracket a f left-parenthesis x right-parenthesis right-bracket 2nd Column equals a upper E left-bracket f left-parenthesis x right-parenthesis right-bracket 2nd Row 1st Column upper E left-bracket sigma-summation Underscript i equals 1 Overscript n Endscripts f left-parenthesis upper X Subscript i Baseline right-parenthesis right-bracket 2nd Column equals sigma-summation Underscript i equals 1 Overscript n Endscripts upper E left-bracket f left-parenthesis upper X Subscript i Baseline right-parenthesis right-bracket period EndLayout
(2.5)

We will repeatedly use these properties in derivations in the following sections.

2.1.3 The Monte Carlo Estimator

We can now define the Monte Carlo estimator, which approximates the value of an arbitrary integral. Suppose that we want to evaluate a 1D integral integral Subscript a Superscript b Baseline f left-parenthesis x right-parenthesis normal d x . Given a supply of independent uniform random variables upper X Subscript i Baseline element-of left-bracket a comma b right-bracket , the Monte Carlo estimator says that the expected value of the estimator

upper F Subscript n Baseline equals StartFraction b minus a Over n EndFraction sigma-summation Underscript i equals 1 Overscript n Endscripts f left-parenthesis upper X Subscript i Baseline right-parenthesis comma
(2.6)

upper E left-bracket upper F Subscript n Baseline right-bracket , is equal to the integral. This fact can be demonstrated with just a few steps. First, note that the PDF p left-parenthesis x right-parenthesis corresponding to the random variable upper X Subscript i must be equal to 1 slash left-parenthesis b minus a right-parenthesis , since p must not only be a constant but also integrate to 1 over the domain left-bracket a comma b right-bracket . Algebraic manipulation using the properties from Equations (2.4) and (2.5) then shows that

StartLayout 1st Row 1st Column upper E left-bracket upper F Subscript n Baseline right-bracket 2nd Column equals upper E left-bracket StartFraction b minus a Over n EndFraction sigma-summation Underscript i equals 1 Overscript n Endscripts f left-parenthesis upper X Subscript i Baseline right-parenthesis right-bracket 2nd Row 1st Column Blank 2nd Column equals StartFraction b minus a Over n EndFraction sigma-summation Underscript i equals 1 Overscript n Endscripts upper E left-bracket f left-parenthesis upper X Subscript i Baseline right-parenthesis right-bracket 3rd Row 1st Column Blank 2nd Column equals StartFraction b minus a Over n EndFraction sigma-summation Underscript i equals 1 Overscript n Endscripts integral Subscript a Superscript b Baseline f left-parenthesis x right-parenthesis p left-parenthesis x right-parenthesis normal d x 4th Row 1st Column Blank 2nd Column equals StartFraction 1 Over n EndFraction sigma-summation Underscript i equals 1 Overscript n Endscripts integral Subscript a Superscript b Baseline f left-parenthesis x right-parenthesis normal d x 5th Row 1st Column Blank 2nd Column equals integral Subscript a Superscript b Baseline f left-parenthesis x right-parenthesis normal d x period EndLayout

Extending this estimator to multiple dimensions or complex integration domains is straightforward: n independent samples upper X Subscript i are taken from a uniform multidimensional PDF, and the estimator is applied in the same way. For example, consider the 3D integral

integral Subscript z 0 Superscript z 1 Baseline integral Subscript y 0 Superscript y 1 Baseline integral Subscript x 0 Superscript x 1 Baseline f left-parenthesis x comma y comma z right-parenthesis normal d x normal d y normal d z period

If samples upper X Subscript i Baseline equals left-parenthesis x Subscript i Baseline comma y Subscript i Baseline comma z Subscript i Baseline right-parenthesis are chosen uniformly from the cube from left-bracket x 0 comma x 1 right-bracket times left-bracket y 0 comma y 1 right-bracket times left-bracket z 0 comma z 1 right-bracket , then the PDF p left-parenthesis upper X right-parenthesis is the constant value

StartFraction 1 Over left-parenthesis x 1 minus x 0 right-parenthesis EndFraction StartFraction 1 Over left-parenthesis y 1 minus y 0 right-parenthesis EndFraction StartFraction 1 Over left-parenthesis z 1 minus z 0 right-parenthesis EndFraction comma

and the estimator is

StartFraction left-parenthesis x 1 minus x 0 right-parenthesis left-parenthesis y 1 minus y 0 right-parenthesis left-parenthesis z 1 minus z 0 right-parenthesis Over n EndFraction sigma-summation Underscript i equals 1 Overscript n Endscripts f left-parenthesis upper X Subscript i Baseline right-parenthesis period

The restriction to uniform random variables can be relaxed with a small generalization. This is an important step, since carefully choosing the PDF from which samples are drawn leads to a key technique for reducing error in Monte Carlo that will be introduced in Section 2.2.2. If the random variables upper X Subscript i are drawn from a PDF p left-parenthesis x right-parenthesis , then the estimator

upper F Subscript n Baseline equals StartFraction 1 Over n EndFraction sigma-summation Underscript i equals 1 Overscript n Endscripts StartFraction f left-parenthesis upper X Subscript i Baseline right-parenthesis Over p left-parenthesis upper X Subscript i Baseline right-parenthesis EndFraction
(2.7)

can be used to estimate the integral instead. The only limitation on p left-parenthesis x right-parenthesis is that it must be nonzero for all  x where StartAbsoluteValue f left-parenthesis x right-parenthesis EndAbsoluteValue greater-than 0 .

It is similarly not too hard to see that the expected value of this estimator is the desired integral of  f :

StartLayout 1st Row 1st Column upper E left-bracket upper F Subscript n Baseline right-bracket 2nd Column equals upper E left-bracket StartFraction 1 Over n EndFraction sigma-summation Underscript i equals 1 Overscript n Endscripts StartFraction f left-parenthesis upper X Subscript i Baseline right-parenthesis Over p left-parenthesis upper X Subscript i Baseline right-parenthesis EndFraction right-bracket 2nd Row 1st Column Blank 2nd Column equals StartFraction 1 Over n EndFraction sigma-summation Underscript i equals 1 Overscript n Endscripts integral Subscript a Superscript b Baseline StartFraction f left-parenthesis x right-parenthesis Over p left-parenthesis x right-parenthesis EndFraction p left-parenthesis x right-parenthesis normal d x 3rd Row 1st Column Blank 2nd Column equals StartFraction 1 Over n EndFraction sigma-summation Underscript i equals 1 Overscript n Endscripts integral Subscript a Superscript b Baseline f left-parenthesis x right-parenthesis normal d x 4th Row 1st Column Blank 2nd Column equals integral Subscript a Superscript b Baseline f left-parenthesis x right-parenthesis normal d x period EndLayout

We can now understand the factor of 1 slash left-parenthesis 4 pi right-parenthesis in the implementation of the RandomWalkIntegrator: directions are uniformly sampled over the unit sphere, which has area 4 pi . Because the PDF is normalized over the sampling domain, it must have the constant value 1 slash left-parenthesis 4 pi right-parenthesis . When the estimator of Equation (2.7) is applied, that value appears in the divisor.

With Monte Carlo, the number of samples n can be chosen arbitrarily, regardless of the dimensionality of the integrand. This is another important advantage of Monte Carlo over traditional deterministic quadrature techniques, which typically require a number of samples that is exponential in the dimension.

2.1.4 Error in Monte Carlo Estimators

Showing that the Monte Carlo estimator converges to the right answer is not enough to justify its use; its rate of convergence is important too. Variance, the expected squared deviation of a function from its expected value, is a useful way to characterize Monte Carlo estimators’ convergence. The variance of an estimator  upper F is defined as

upper V left-bracket upper F right-bracket equals upper E left-bracket left-parenthesis upper F minus upper E left-bracket upper F right-bracket right-parenthesis squared right-bracket comma
(2.8)

from which it follows that

upper V left-bracket a upper F right-bracket equals a squared upper V left-bracket upper F right-bracket period

This property and Equation (2.5) yield an alternative expression for the variance:

upper V left-bracket upper F right-bracket equals upper E left-bracket upper F squared right-bracket minus upper E left-bracket upper F right-bracket squared period
(2.9)

Thus, the variance is the expected value of the square minus the square of the expected value.

If the estimator is a sum of independent random variables (like the Monte Carlo estimator upper F Subscript n ), then the variance of the sum is the sum of the individual random variables’ variances:

upper V left-bracket sigma-summation Underscript i equals 1 Overscript n Endscripts upper X Subscript i Baseline right-bracket equals sigma-summation Underscript i equals 1 Overscript n Endscripts upper V left-bracket upper X Subscript i Baseline right-bracket period
(2.10)

From Equation (2.10) it is easy to show that variance decreases linearly with the number of samples n . Because variance is squared error, the error in a Monte Carlo estimate therefore only goes down at a rate of upper O left-parenthesis n Superscript negative 1 slash 2 Baseline right-parenthesis in the number of samples. Although standard quadrature techniques converge at a faster rate in one dimension, their performance becomes exponentially worse as the dimensionality of the integrand increases, while Monte Carlo’s convergence rate is independent of the dimension, making Monte Carlo the only practical numerical integration algorithm for high-dimensional integrals.

The upper O left-parenthesis n Superscript negative 1 slash 2 Baseline right-parenthesis characteristic of Monte Carlo’s rate of error reduction is apparent when watching a progressive rendering of a scene where additional samples are incrementally taken in all pixels. The image improves rapidly for the first few samples when doubling the number of samples is relatively little additional work. Later on, once tens or hundreds of samples have been taken, each additional sample doubling takes much longer and remaining error in the image takes a long time to disappear.

The linear decrease in variance with increasing numbers of samples makes it easy to compare different Monte Carlo estimators. Consider two estimators, where the second has half the variance of the first but takes three times as long to compute an estimate; which of the two is better? In that case, the first is preferable: it could take three times as many samples in the time consumed by the second, in which case it would achieve a 3 times variance reduction. This concept can be encapsulated in the efficiency of an estimator upper F , which is defined as

epsilon left-bracket upper F right-bracket equals StartFraction 1 Over upper V left-bracket upper F right-bracket upper T left-bracket upper F right-bracket EndFraction comma

where upper V left-bracket upper F right-bracket is its variance and upper T left-bracket upper F right-bracket is the running time to compute its value.

Not all estimators of integrals have expected values that are equal to the value of the integral. Such estimators are said to be biased, where the difference

beta equals upper E left-bracket upper F right-bracket minus integral f left-parenthesis x right-parenthesis normal d x

is the amount of bias. Biased estimators may still be desirable if they are able to get close to the correct result more quickly than unbiased estimators. Kalos and Whitlock (1986, pp. 36–37) gave the following example: consider the problem of computing an estimate of the mean value of a uniform distribution upper X Subscript i Baseline tilde p over the interval from 0 to 1. One could use the estimator

StartFraction 1 Over n EndFraction sigma-summation Underscript i equals 1 Overscript n Endscripts upper X Subscript i Baseline comma

or one could use the biased estimator

one-half max left-parenthesis upper X 1 comma upper X 2 comma ellipsis comma upper X Subscript n Baseline right-parenthesis period

The first estimator is unbiased but has variance with order upper O left-parenthesis n Superscript negative 1 Baseline right-parenthesis . The second estimator’s expected value is

0.5 StartFraction n Over n plus 1 EndFraction not-equals 0.5 comma

so it is biased, although its variance is upper O left-parenthesis n Superscript negative 2 Baseline right-parenthesis , which is much better. This estimator has the useful property that its error goes to 0 in the limit as the number of samples n goes to infinity; such estimators are consistent. Most of the Monte Carlo estimators used in pbrt are unbiased, with the notable exception of the SPPMIntegrator, which implements a photon mapping algorithm.

Closely related to the variance is the mean squared error (MSE), which is defined as the expectation of the squared difference of an estimator and the true value,

upper M upper S upper E left-bracket upper F right-bracket equals upper E left-bracket left-parenthesis upper F minus integral f left-parenthesis x right-parenthesis normal d x right-parenthesis squared right-bracket period

For an unbiased estimator, MSE is equal to the variance; otherwise it is the sum of variance and the squared bias of the estimator.

It is possible to work out the variance and MSE of some simple estimators in closed form, but for most of the ones of interest in rendering, this is not possible. Yet it is still useful to be able to quantify these values. For this purpose, the sample variance can be computed using a set of independent random variables upper X Subscript i . Equation (2.8) points at one way to compute the sample variance for a set of n random variables upper X Subscript i . If the sample mean is computed as their average, upper X overbar equals left-parenthesis 1 slash n right-parenthesis sigma-summation upper X Subscript i , then the sample variance is

StartFraction 1 Over n minus 1 EndFraction sigma-summation Underscript i equals 1 Overscript n Endscripts left-parenthesis upper X Subscript i Baseline minus upper X overbar right-parenthesis squared period
(2.11)

The division by n minus 1 rather than n is Bessel’s correction, and ensures that the sample variance is an unbiased estimate of the variance. (See also Section B.2.11, where a numerically stable approach for computing the sample variance is introduced.)

The sample variance is itself an estimate of the variance, so it has variance itself. Consider, for example, a random variable that has a value of 1 99.99% of the time, and a value of one million 0.01% of the time. If we took ten random samples of it that all had the value 1, the sample variance would suggest that the random variable had zero variance even though its variance is actually much higher.

If an accurate estimate of the integral upper F overTilde almost-equals integral f left-parenthesis x right-parenthesis normal d x can be computed (for example, using a large number of samples), then the mean squared error can be estimated by

upper M upper S upper E left-bracket upper F right-bracket almost-equals StartFraction 1 Over n EndFraction sigma-summation Underscript i equals 1 Overscript n Endscripts left-parenthesis f left-parenthesis upper X Subscript i Baseline right-parenthesis minus upper F overTilde right-parenthesis squared period

The imgtool utility program that is provided in pbrt’s distribution can compute an image’s MSE with respect to a reference image via its diff option.