13.3 Sampling Random Variables

In order to evaluate the Monte Carlo estimator in Equation (13.3), it is necessary to be able to draw random samples from the chosen probability distribution. This section will introduce the basics of this process and demonstrate it with some straightforward examples. The next two sections will introduce more complex approaches to sampling before Section 13.6 develops the approach for the general multidimensional case. In Chapters 14, 15, and 16, we’ll then see how to use these techniques to generate samples from the distributions defined by BSDFs, light sources, cameras, and scattering media.

13.3.1 The Inversion Method

The inversion method uses one or more uniform random variables and maps them to random variables from the desired distribution. To explain how this process works in general, we will start with a simple discrete example. Suppose we have a process with four possible outcomes. The probabilities of each of the four outcomes are given by p 1 , p 2 , p 3 , and p 4 , respectively, with the requirement that sigma-summation Underscript i equals 1 Overscript 4 Endscripts p Subscript i Baseline equals 1 . The corresponding PDF is shown in Figure 13.1.

Figure 13.1: A Discrete PDF for Four Events, Each with a Probability p Subscript i . The sum of their probabilities sigma-summation Underscript i Endscripts p Subscript i is necessarily 1.

In order to draw a sample from this distribution, we first find the CDF upper P left-parenthesis x right-parenthesis . In the continuous case, upper P is the indefinite integral of p . In the discrete case, we can directly construct the CDF by stacking the bars on top of each other, starting at the left. This idea is shown in Figure 13.2. Notice that the height of the rightmost bar must be 1 because of the requirement that all probabilities sum to 1.

Figure 13.2: A Discrete CDF, Corresponding to the PDF in Figure 13.1. Each column’s height is given by the PDF for the event that it represents plus the sum of the PDFs for the previous events, upper P Subscript i Baseline equals sigma-summation Underscript j equals 1 Overscript i Endscripts p Subscript j .

To draw a sample from the distribution, we then take a uniform random number xi Subscript and use it to select one of the possible outcomes using the CDF, doing so in a way that chooses a particular outcome with probability equal to the outcome’s own probability. This idea is illustrated in Figure 13.3, where the events’ probabilities are projected onto the vertical axis and a random variable xi Subscript selects among them. It should be clear that this draws from the correct distribution—the probability of the uniform sample hitting any particular bar is exactly equal to the height of that bar.

Figure 13.3: To use the inversion method to draw a sample from the distribution described by the PDF in Figure 13.1, a canonical uniform random variable is plotted on the vertical axis. By construction, the horizontal extension of xi Subscript will intersect the box representing the i th outcome with probability p Subscript i . If the corresponding event is chosen for a set of random variables xi Subscript , then the resulting distribution of events will be distributed according to the PDF.

In order to generalize this technique to continuous distributions, consider what happens as the number of discrete possibilities approaches infinity. The PDF from Figure 13.1 becomes a smooth curve, and the CDF from Figure 13.2 becomes its integral. The projection process is still the same, although if the function is continuous, the projection has a convenient mathematical interpretation—it represents inverting the CDF and evaluating the inverse at xi Subscript . This technique is thus called the inversion method.

More precisely, we can draw a sample upper X Subscript i from an arbitrary PDF p left-parenthesis x right-parenthesis with the following steps:

  1. Compute the CDF upper P left-parenthesis x right-parenthesis equals integral Subscript 0 Superscript x Baseline p left-parenthesis x prime right-parenthesis normal d x Superscript prime Baseline .
  2. Compute the inverse upper P Superscript negative 1 Baseline left-parenthesis x right-parenthesis .
  3. Obtain a uniformly distributed random number xi Subscript .
  4. Compute upper X Subscript i Baseline equals upper P Superscript negative 1 Baseline left-parenthesis xi Subscript Baseline right-parenthesis .

Example: Power Distribution

As an example of how this procedure works, consider the task of drawing samples from a power distribution, p left-parenthesis x right-parenthesis proportional-to x Superscript n . The PDF of the power distribution is

p left-parenthesis x right-parenthesis equals c x Superscript n Baseline comma

for the constant c that normalizes the PDF. The first task to tackle is to find the PDF. In most cases, this simply involves computing the value of the proportionality constant  c , which can be found using the constraint that integral p left-parenthesis x right-parenthesis normal d x equals 1 :

StartLayout 1st Row 1st Column integral Subscript 0 Superscript 1 Baseline c x Superscript n normal d x 2nd Column equals 1 2nd Row 1st Column c StartFraction x Superscript n plus 1 Baseline Over n plus 1 EndFraction vertical-bar Subscript 0 Superscript 1 2nd Column equals 1 3rd Row 1st Column StartFraction c Over n plus 1 EndFraction 2nd Column equals 1 4th Row 1st Column c 2nd Column equals n plus 1 period EndLayout

Therefore, p left-parenthesis x right-parenthesis equals left-parenthesis n plus 1 right-parenthesis x Superscript n . We can integrate this PDF to get the CDF:

upper P left-parenthesis x right-parenthesis equals integral Subscript 0 Superscript x Baseline p left-parenthesis x Superscript prime Baseline right-parenthesis normal d x Superscript prime Baseline equals x Superscript n plus 1 Baseline comma

and inversion is simple: upper P Superscript negative 1 Baseline left-parenthesis x right-parenthesis equals RootIndex n plus 1 StartRoot x EndRoot . Therefore, given a uniform random variable xi Subscript , samples can be drawn from the power distribution as

upper X equals RootIndex n plus 1 StartRoot xi Subscript Baseline EndRoot period

Another approach is to use a sampling trick that works only for the power distribution, selecting upper X equals max left-parenthesis xi 1 comma xi 2 comma ellipsis comma xi Subscript n plus 1 Baseline right-parenthesis . This random variable is distributed according to the power distribution as well. To see why, note that upper P r left-brace upper X less-than x right-brace is the probability that all the xi Subscript i Baseline less-than x . But the xi Subscript i are independent, so

upper P r left-brace upper X less-than x right-brace equals product Underscript i equals 1 Overscript n plus 1 Endscripts upper P r left-brace xi Subscript i Baseline less-than x right-brace equals x Superscript n plus 1 Baseline comma

which is exactly the desired CDF. Depending on the speed of your random number generator, this technique can be faster than the inversion method for small values of n .

Example: Exponential Distribution

When rendering images with participating media, it is frequently useful to draw samples from a distribution p left-parenthesis x right-parenthesis proportional-to normal e Superscript minus a x . As before, the first step is to normalize this distribution so that it integrates to one. In this case, we’ll assume for now that the range of values x we’d like the generated samples to cover is left-bracket 0 comma normal infinity right-parenthesis rather than left-bracket 0 comma 1 right-bracket , so

integral Subscript 0 Superscript normal infinity Baseline c normal e Superscript minus a x Baseline normal d x equals minus StartFraction c Over a EndFraction normal e Superscript minus a x Baseline vertical-bar Subscript 0 Superscript normal infinity Baseline equals StartFraction c Over a EndFraction equals 1 period

Thus we know that c equals a , and our PDF is p left-parenthesis x right-parenthesis equals a normal e Superscript minus a x . Now, we integrate to find upper P left-parenthesis x right-parenthesis :

upper P left-parenthesis x right-parenthesis equals integral Subscript 0 Superscript x Baseline a normal e Superscript minus a x Super Superscript prime Superscript Baseline normal d x Superscript prime Baseline equals 1 minus normal e Superscript minus a x Baseline period

This function is easy to invert:

upper P Superscript negative 1 Baseline left-parenthesis x right-parenthesis equals minus StartFraction ln left-parenthesis 1 minus x right-parenthesis Over a EndFraction comma

and we can draw samples thus:

upper X equals minus StartFraction ln left-parenthesis 1 minus xi Subscript Baseline right-parenthesis Over a EndFraction period

It may be tempting to simplify the log term from ln left-parenthesis 1 minus xi Subscript Baseline right-parenthesis to ln xi Subscript , under the theory that because xi Subscript Baseline element-of left-bracket 0 comma 1 right-parenthesis , these are effectively the same and a subtraction can thus be saved. The problem with this idea is that xi Subscript may have the value 0 but never has the value 1. With the simplification, it’s possible that we’d try to take the logarithm of 0, which is undefined; this danger is avoided with the first formulation. While a xi Subscript value of 0 may seem very unlikely, it does happen (especially in the world of floating-point arithmetic, rather than the real numbers). Some of the low-discrepancy sampling patterns introduced in Chapter 7 are particularly prone to generating the value 0.

Example: Piecewise-Constant 1D Functions

An interesting exercise is to work out how to sample from 1D piecewise-constant functions (step functions). Without loss of generality, we will just consider piecewise-constant functions defined over left-bracket 0 comma 1 right-bracket .

Assume that the 1D function’s domain is split into upper N equal-sized pieces of size normal upper Delta equals 1 slash upper N . These regions start and end at points x Subscript i Baseline equals i normal upper Delta , where i ranges from 0 to upper N , inclusive. Within each region, the value of the function f left-parenthesis x right-parenthesis is a constant (Figure 13.4(a)).

Figure 13.4: (a) Probability density function for a piecewise-constant 1D function and (b) cumulative distribution function defined by this PDF.

The value of f left-parenthesis x right-parenthesis  is

f left-parenthesis x right-parenthesis equals StartLayout Enlarged left-brace 1st Row 1st Column v 0 2nd Column x 0 less-than-or-equal-to x less-than x 1 2nd Row 1st Column v 1 2nd Column x 1 less-than-or-equal-to x less-than x 2 3rd Row 1st Column vertical-ellipsis EndLayout period

The integral integral f left-parenthesis x right-parenthesis normal d x is

c equals integral Subscript 0 Superscript 1 Baseline f left-parenthesis x right-parenthesis normal d x equals sigma-summation Underscript i equals 0 Overscript upper N minus 1 Endscripts normal upper Delta v Subscript i Baseline equals sigma-summation Underscript i equals 0 Overscript upper N minus 1 Endscripts StartFraction v Subscript i Baseline Over upper N EndFraction comma

and so it is easy to construct the PDF p left-parenthesis x right-parenthesis for f left-parenthesis x right-parenthesis as f left-parenthesis x right-parenthesis slash c . By direct application of the relevant formulae, the CDF upper P left-parenthesis x right-parenthesis is a piecewise linear function defined at points x Subscript i  by

StartLayout 1st Row 1st Column upper P left-parenthesis x 0 right-parenthesis 2nd Column equals 0 2nd Row 1st Column upper P left-parenthesis x 1 right-parenthesis 2nd Column equals integral Subscript x 0 Superscript x 1 Baseline p left-parenthesis x right-parenthesis normal d x equals StartFraction v 0 Over upper N c EndFraction equals upper P left-parenthesis x 0 right-parenthesis plus StartFraction v 0 Over upper N c EndFraction 3rd Row 1st Column upper P left-parenthesis x 2 right-parenthesis 2nd Column equals integral Subscript x 0 Superscript x 2 Baseline p left-parenthesis x right-parenthesis normal d x equals integral Subscript x 0 Superscript x 1 Baseline p left-parenthesis x right-parenthesis normal d x plus integral Subscript x 1 Superscript x 2 Baseline p left-parenthesis x right-parenthesis normal d x equals upper P left-parenthesis x 1 right-parenthesis plus StartFraction v 1 Over upper N c EndFraction 4th Row 1st Column upper P left-parenthesis x Subscript i Baseline right-parenthesis 2nd Column equals upper P left-parenthesis x Subscript i minus 1 Baseline right-parenthesis plus StartFraction v Subscript i minus 1 Baseline Over upper N c EndFraction period EndLayout

Between two points x Subscript i and x Subscript i plus 1 , the CDF is linearly increasing with slope v Subscript i Baseline slash c .

Recall that in order to sample f left-parenthesis x right-parenthesis we need to invert the CDF to find the value x such that

xi Subscript Baseline equals integral Subscript 0 Superscript x Baseline p left-parenthesis x Superscript prime Baseline right-parenthesis normal d x Superscript prime Baseline equals upper P left-parenthesis x right-parenthesis period

Because the CDF is monotonically increasing, the value of x must be between the x Subscript i and x Subscript i plus 1 such that upper P left-parenthesis x Subscript i Baseline right-parenthesis less-than-or-equal-to xi Subscript and xi Subscript Baseline less-than-or-equal-to upper P left-parenthesis x Subscript i plus 1 Baseline right-parenthesis . Given an array of CDF values, this pair of upper P left-parenthesis x Subscript i Baseline right-parenthesis values can be efficiently found with a binary search.

Distribution1D is a small utility class that represents a piecewise-constant 1D function’s PDF and CDF and provides methods to perform this sampling efficiently.

<<Sampling Declarations>>= 
struct Distribution1D { <<Distribution1D Public Methods>> 
Distribution1D(const Float *f, int n) : func(f, f + n), cdf(n + 1) { <<Compute integral of step function at x Subscript i >> 
cdf[0] = 0; for (int i = 1; i < n + 1; ++i) cdf[i] = cdf[i - 1] + func[i - 1] / n;
<<Transform step function integral into CDF>> 
funcInt = cdf[n]; if (funcInt == 0) { for (int i = 1; i < n + 1; ++i) cdf[i] = Float(i) / Float(n); } else { for (int i = 1; i < n + 1; ++i) cdf[i] /= funcInt; }
} int Count() const { return func.size(); } Float SampleContinuous(Float u, Float *pdf, int *off = nullptr) const { <<Find surrounding CDF segments and offset>> 
int offset = FindInterval(cdf.size(), [&](int index) { return cdf[index] <= u; });
if (off) *off = offset; <<Compute offset along CDF segment>> 
Float du = u - cdf[offset]; if ((cdf[offset + 1] - cdf[offset]) > 0) du /= (cdf[offset + 1] - cdf[offset]);
<<Compute PDF for sampled offset>> 
if (pdf) *pdf = func[offset] / funcInt;
<<Return x element-of left-bracket 0 comma 1 right-parenthesis corresponding to sample>> 
return (offset + du) / Count();
} int SampleDiscrete(Float u, Float *pdf = nullptr, Float *uRemapped = nullptr) const { <<Find surrounding CDF segments and offset>> 
int offset = FindInterval(cdf.size(), [&](int index) { return cdf[index] <= u; });
if (pdf) *pdf = func[offset] / (funcInt * Count()); if (uRemapped) *uRemapped = (u - cdf[offset]) / (cdf[offset + 1] - cdf[offset]); return offset; } Float DiscretePDF(int index) const { return func[index] / (funcInt * Count()); }
<<Distribution1D Public Data>> 
std::vector<Float> func, cdf; Float funcInt;
};

The Distribution1D constructor takes n values of a piecewise-constant function f. It makes its own copy of the function values, computes the function’s CDF, and also stores the integral of the function, funcInt. Note that the constructor allocates n+1 Floats for the cdf array because if f left-parenthesis x right-parenthesis has upper N step values, then we need to store the value of the CDF at each of the upper N plus 1 values of x Subscript i . Storing the CDF value of 1 at the end of the array is redundant but simplifies the sampling code later.

<<Distribution1D Public Methods>>= 
Distribution1D(const Float *f, int n) : func(f, f + n), cdf(n + 1) { <<Compute integral of step function at x Subscript i >> 
cdf[0] = 0; for (int i = 1; i < n + 1; ++i) cdf[i] = cdf[i - 1] + func[i - 1] / n;
<<Transform step function integral into CDF>> 
funcInt = cdf[n]; if (funcInt == 0) { for (int i = 1; i < n + 1; ++i) cdf[i] = Float(i) / Float(n); } else { for (int i = 1; i < n + 1; ++i) cdf[i] /= funcInt; }
}

<<Distribution1D Public Data>>= 
std::vector<Float> func, cdf; Float funcInt;

<<Distribution1D Public Methods>>+=  
int Count() const { return func.size(); }

This constructor computes the integral of f left-parenthesis x right-parenthesis using Equation (13.5). It stores the result in the cdf array for now so that it doesn’t need to allocate additional temporary space for it.

<<Compute integral of step function at x Subscript i >>= 
cdf[0] = 0; for (int i = 1; i < n + 1; ++i) cdf[i] = cdf[i - 1] + func[i - 1] / n;

Now that the value of the integral over all of left-bracket 0 comma 1 right-bracket is stored in cdf[n], this value can be copied into funcInt and the CDF can be normalized by dividing through all entries by this value.

<<Transform step function integral into CDF>>= 
funcInt = cdf[n]; if (funcInt == 0) { for (int i = 1; i < n + 1; ++i) cdf[i] = Float(i) / Float(n); } else { for (int i = 1; i < n + 1; ++i) cdf[i] /= funcInt; }

The Distribution1D::SampleContinuous() method uses the given random sample u to sample from its distribution. It returns the corresponding value x element-of left-bracket 0 comma 1 right-parenthesis and the value of the PDF p left-parenthesis x right-parenthesis . If the optional off parameter is not nullptr, it returns the offset into the array of function values of the largest index where the CDF was less than or equal to u. (In other words, cdf[*off] <= u < cdf[*off+1].)

<<Distribution1D Public Methods>>+=  
Float SampleContinuous(Float u, Float *pdf, int *off = nullptr) const { <<Find surrounding CDF segments and offset>> 
int offset = FindInterval(cdf.size(), [&](int index) { return cdf[index] <= u; });
if (off) *off = offset; <<Compute offset along CDF segment>> 
Float du = u - cdf[offset]; if ((cdf[offset + 1] - cdf[offset]) > 0) du /= (cdf[offset + 1] - cdf[offset]);
<<Compute PDF for sampled offset>> 
if (pdf) *pdf = func[offset] / funcInt;
<<Return x element-of left-bracket 0 comma 1 right-parenthesis corresponding to sample>> 
return (offset + du) / Count();
}

Mapping u to an interval matching the above criterion is carried out using the efficient binary search implemented in FindInterval() (see Appendix A for details).

<<Find surrounding CDF segments and offset>>= 
int offset = FindInterval(cdf.size(), [&](int index) { return cdf[index] <= u; });

Given the pair of CDF values that straddle u, we can compute x . First, we determine how far u is between cdf[offset] and cdf[offset+1], du, where du is 0 if u == cdf[offset] and goes up to 1 if u == cdf[offset+1]. Because the CDF is piecewise linear, the sample value x is the same offset between x Subscript i and x Subscript i plus 1 (Figure 13.4(b)).

<<Compute offset along CDF segment>>= 
Float du = u - cdf[offset]; if ((cdf[offset + 1] - cdf[offset]) > 0) du /= (cdf[offset + 1] - cdf[offset]);

The PDF for this sample p left-parenthesis x right-parenthesis is easily computed since we have the function’s integral in funcInt. (Note that the offset offset into the CDF array has been computed in a way so that func[offset] gives the value of the function in the CDF range that the sample landed in.)

<<Compute PDF for sampled offset>>= 
if (pdf) *pdf = func[offset] / funcInt;

Finally, the appropriate value of x is computed and returned.

<<Return x element-of left-bracket 0 comma 1 right-parenthesis corresponding to sample>>= 
return (offset + du) / Count();

In a small overloading of semantics, Distribution1D can also be used for discrete 1D probability distributions where there are some number of buckets n , each with some weight, and we’d like to sample among the buckets with probability proportional to their relative weights. This functionality is used, for example, by some of the Integrators that compute a discrete distribution for the light sources in the scene with weights given by the lights’ powers. Sampling from the discrete distribution just requires figuring out which pair of CDF values the sample value lies between; the PDF is computed as the discrete probability of sampling the corresponding bucket.

<<Distribution1D Public Methods>>+=  
int SampleDiscrete(Float u, Float *pdf = nullptr, Float *uRemapped = nullptr) const { <<Find surrounding CDF segments and offset>> 
int offset = FindInterval(cdf.size(), [&](int index) { return cdf[index] <= u; });
if (pdf) *pdf = func[offset] / (funcInt * Count()); if (uRemapped) *uRemapped = (u - cdf[offset]) / (cdf[offset + 1] - cdf[offset]); return offset; }

It’s also useful to be able to compute the PDF for sampling a given value from the discrete PDF.

<<Distribution1D Public Methods>>+= 
Float DiscretePDF(int index) const { return func[index] / (funcInt * Count()); }

13.3.2 The Rejection Method

For some functions f left-parenthesis x right-parenthesis , it may not be possible to integrate them in order to find their PDFs, or it may not be possible to analytically invert their CDFs. The rejection method is a technique for generating samples according to a function’s distribution without needing to do either of these steps; it is essentially a dart-throwing approach. Assume that we want to draw samples from some such function f left-parenthesis x right-parenthesis but we do have a PDF p left-parenthesis x right-parenthesis that satisfies f left-parenthesis x right-parenthesis less-than c p left-parenthesis x right-parenthesis for some scalar constant c , and suppose that we do know how to sample from  p . The rejection method is then:

StartLayout 1st Row 1st Column Blank 2nd Column loop forever colon 2nd Row 1st Column Blank 2nd Column sample upper X from p prime s distribution 3rd Row 1st Column Blank 2nd Column if xi Subscript Baseline less-than f left-parenthesis upper X right-parenthesis slash left-parenthesis c p left-parenthesis upper X right-parenthesis right-parenthesis then 4th Row 1st Column Blank 2nd Column return upper X EndLayout

This procedure repeatedly chooses a pair of random variables left-parenthesis upper X comma xi Subscript Baseline right-parenthesis . If the point left-parenthesis upper X comma xi Subscript Baseline c p left-parenthesis upper X right-parenthesis right-parenthesis lies under f left-parenthesis upper X right-parenthesis , then the sample upper X is accepted. Otherwise, it is rejected and a new sample pair is chosen. This idea is illustrated in Figure 13.5. Without going into too much detail, it should be clear that the efficiency of this scheme depends on how tightly c p left-parenthesis x right-parenthesis bounds f left-parenthesis x right-parenthesis . This technique works in any number of dimensions.

Figure 13.5: Rejection sampling generates samples according to the distribution of an arbitrary function f left-parenthesis x right-parenthesis even if f ’s PDF is unknown or its CDF can’t be inverted. If some distribution p left-parenthesis x right-parenthesis and a scalar constant c are known such that f left-parenthesis x right-parenthesis less-than c p left-parenthesis x right-parenthesis , then samples can be drawn from p left-parenthesis x right-parenthesis and randomly accepted with the rejection method. The closer the fit of c p left-parenthesis x right-parenthesis to f left-parenthesis x right-parenthesis , the more efficient this process is.

Rejection sampling isn’t actually used in any of the Monte Carlo algorithms currently implemented in pbrt. We will normally prefer to find distributions that are similar to f left-parenthesis x right-parenthesis that can be sampled directly, so that well-distributed points on left-bracket 0 comma 1 right-parenthesis Superscript n can be mapped to sample points that are in turn well-distributed, as will be discussed in Section 13.8. Nevertheless, rejection sampling is an important technique to be aware of, particularly when debugging Monte Carlo implementations. For example, if one suspects the presence of a bug in code that draws samples from some distribution using the inversion method, then one can replace it with a straightforward implementation based on the rejection method and see if the Monte Carlo estimator computes the same result. Of course, it’s necessary to take many samples in situations like these, so that variance in the estimates doesn’t mask errors.

Example: Rejection Sampling a Unit Circle

Suppose we want to select a uniformly distributed point inside a unit circle. Using the rejection method, we simply select a random left-parenthesis x comma y right-parenthesis position inside the circumscribed square and return it if it falls inside the circle. This process is shown in Figure 13.6.

Figure 13.6: Rejection Sampling a Circle. One approach to finding uniform points in the unit circle is to sample uniform random points in the unit square and reject all that lie outside the circle (red points). The remaining points will be uniformly distributed within the circle.

The function RejectionSampleDisk() implements this algorithm. A similar approach will work to generate uniformly distributed samples on the inside of any complex shape as long as it has an inside–outside test.

<<Sampling Function Definitions>>+=  
Point2f RejectionSampleDisk(RNG &rng) { Point2f p; do { p.x = 1 - 2 * rng.UniformFloat(); p.y = 1 - 2 * rng.UniformFloat(); } while (p.x * p.x + p.y * p.y > 1); return p; }

In general, the efficiency of rejection sampling depends on the percentage of samples that are expected to be rejected. For the problem of finding uniform points in the 2D case, this is easy to compute. It is the area of the circle divided by the area of the square: StartFraction pi Over 4 EndFraction almost-equals 78.5 percent-sign . If the method is applied to generate samples in hyperspheres in the general n -dimensional case, however, the volume of an n -dimensional hypersphere actually goes to 0 as n increases, and this approach becomes increasingly inefficient.