2.4 Transforming between Distributions

In describing the inversion method, we introduced a technique that generates samples according to some distribution by transforming canonical uniform random variables in a particular manner. Here, we will investigate the more general question of which distribution results when we transform samples from an arbitrary distribution to some other distribution with a function  f . Understanding the effect of such transformations is useful for a few reasons, though here we will focus on how they allow us to derive multidimensional sampling algorithms.

Suppose we are given a random variable upper X drawn from some PDF p left-parenthesis x right-parenthesis with CDF upper P left-parenthesis x right-parenthesis . Given a function f left-parenthesis x right-parenthesis with y equals f left-parenthesis x right-parenthesis , if we compute upper Y equals f left-parenthesis upper X right-parenthesis , we would like to find the distribution of the new random variable  upper Y . In this case, the function f left-parenthesis x right-parenthesis must be a one-to-one transformation; if multiple values of x mapped to the same y value, then it would be impossible to unambiguously describe the probability density of a particular y value. A direct consequence of f being one-to-one is that its derivative must either be strictly greater than 0 or strictly less than 0, which implies that for a given x ,

probability left-brace upper Y less-than-or-equal-to f left-parenthesis x right-parenthesis right-brace equals upper P r left-brace upper X less-than-or-equal-to x right-brace period

From the definition of the CDF, Equation (2.3), we can see that

upper P Subscript f Baseline left-parenthesis y right-parenthesis equals upper P Subscript f Baseline left-parenthesis f left-parenthesis x right-parenthesis right-parenthesis equals upper P left-parenthesis x right-parenthesis period

This relationship between CDFs leads directly to the relationship between their PDFs. If we assume that f ’s derivative is greater than 0, differentiating gives

p Subscript f Baseline left-parenthesis y right-parenthesis StartFraction normal d f Over normal d x EndFraction equals p left-parenthesis x right-parenthesis comma

and so

p Subscript f Baseline left-parenthesis y right-parenthesis equals left-parenthesis StartFraction normal d f Over normal d x EndFraction right-parenthesis Superscript negative 1 Baseline p left-parenthesis x right-parenthesis period

In general, f ’s derivative is either strictly positive or strictly negative, and the relationship between the densities is

p Subscript f Baseline left-parenthesis y right-parenthesis equals StartAbsoluteValue StartFraction normal d f Over normal d x EndFraction EndAbsoluteValue Superscript negative 1 Baseline p left-parenthesis x right-parenthesis period

How can we use this formula? Suppose that p left-parenthesis x right-parenthesis equals 2 x over the domain left-bracket 0 comma 1 right-bracket , and let f left-parenthesis x right-parenthesis equals sine x . What is the PDF of the random variable upper Y equals f left-parenthesis upper X right-parenthesis ? Because we know that normal d f slash normal d x equals cosine x ,

p Subscript f Baseline left-parenthesis y right-parenthesis equals StartFraction p left-parenthesis x right-parenthesis Over StartAbsoluteValue cosine x EndAbsoluteValue EndFraction equals StartFraction 2 x Over cosine x EndFraction equals StartFraction 2 arc sine y Over StartRoot 1 minus y squared EndRoot EndFraction period

This procedure may seem backward—usually we have some PDF that we want to sample from, not a given transformation. For example, we might have upper X drawn from some p left-parenthesis x right-parenthesis and would like to compute upper Y from some distribution p Subscript f Baseline left-parenthesis y right-parenthesis . What transformation should we use? All we need is for the CDFs to be equal, or upper P Subscript f Baseline left-parenthesis y right-parenthesis equals upper P left-parenthesis x right-parenthesis , which immediately gives the transformation

f left-parenthesis x right-parenthesis equals upper P Subscript f Superscript negative 1 Baseline left-parenthesis upper P left-parenthesis x right-parenthesis right-parenthesis period

This is a generalization of the inversion method, since if upper X were uniformly distributed over left-bracket 0 comma 1 right-parenthesis then upper P left-parenthesis x right-parenthesis equals x , and we have the same procedure as was introduced previously.

2.4.1 Transformation in Multiple Dimensions

In the general d -dimensional case, a similar derivation gives the analogous relationship between different densities. We will not show the derivation here; it follows the same form as the 1D case. Suppose we have a d -dimensional random variable upper X with density function p left-parenthesis x right-parenthesis . Now let upper Y equals upper T left-parenthesis upper X right-parenthesis , where upper T is a bijection. In this case, the densities are related by

p Subscript upper T Baseline left-parenthesis y right-parenthesis equals p Subscript upper T Baseline left-parenthesis upper T left-parenthesis x right-parenthesis right-parenthesis equals StartFraction p left-parenthesis x right-parenthesis Over StartAbsoluteValue upper J Subscript upper T Baseline left-parenthesis x right-parenthesis EndAbsoluteValue EndFraction comma

where StartAbsoluteValue upper J Subscript upper T Baseline EndAbsoluteValue is the absolute value of the determinant of upper T ’s Jacobian matrix, which is

Start 3 By 3 Matrix 1st Row 1st Column partial-differential upper T 1 slash partial-differential x 1 2nd Column midline-horizontal-ellipsis 3rd Column partial-differential upper T 1 slash partial-differential x Subscript d Baseline 2nd Row 1st Column vertical-ellipsis 2nd Column down-right-diagonal-ellipsis 3rd Column vertical-ellipsis 3rd Row 1st Column partial-differential upper T Subscript d Baseline slash partial-differential x 1 2nd Column midline-horizontal-ellipsis 3rd Column partial-differential upper T Subscript d Baseline slash partial-differential x Subscript d Baseline EndMatrix comma

where subscripts index dimensions of upper T left-parenthesis x right-parenthesis and x .

For a 2D example of the use of Equation (2.21), the polar transformation relates Cartesian left-parenthesis x comma y right-parenthesis coordinates to a polar radius and angle,

StartLayout 1st Row 1st Column x 2nd Column equals r cosine theta 2nd Row 1st Column y 2nd Column equals r sine theta period EndLayout

Suppose we draw samples from some density p left-parenthesis r comma theta right-parenthesis . What is the corresponding density p left-parenthesis x comma y right-parenthesis ? The Jacobian of this transformation is

upper J Subscript upper T Baseline equals Start 2 By 2 Matrix 1st Row 1st Column StartFraction partial-differential x Over partial-differential r EndFraction 2nd Column StartFraction partial-differential x Over partial-differential theta EndFraction 2nd Row 1st Column StartFraction partial-differential y Over partial-differential r EndFraction 2nd Column StartFraction partial-differential y Over partial-differential theta EndFraction EndMatrix equals Start 2 By 2 Matrix 1st Row 1st Column cosine theta 2nd Column minus r sine theta 2nd Row 1st Column sine theta 2nd Column r cosine theta EndMatrix comma

and the determinant is r left-parenthesis cosine squared theta plus sine squared theta right-parenthesis equals r . So, p left-parenthesis x comma y right-parenthesis equals p left-parenthesis r comma theta right-parenthesis slash r . Of course, this is backward from what we usually want—typically we start with a sampling strategy in Cartesian coordinates and want to transform it to one in polar coordinates. In that case, we would have

p left-parenthesis r comma theta right-parenthesis equals r p left-parenthesis x comma y right-parenthesis period

In 3D, given the spherical coordinate representation of directions, Equation (3.7), the Jacobian of this transformation has determinant StartAbsoluteValue upper J Subscript upper T Baseline EndAbsoluteValue equals r squared sine theta , so the corresponding density function is

p left-parenthesis r comma theta comma phi right-parenthesis equals r squared sine theta p left-parenthesis x comma y comma z right-parenthesis period

This transformation is important since it helps us represent directions as points left-parenthesis x comma y comma z right-parenthesis on the unit sphere.

2.4.2 Sampling with Multidimensional Transformations

Suppose we have a 2D joint density function p left-parenthesis x comma y right-parenthesis that we wish to draw samples left-parenthesis upper X comma upper Y right-parenthesis from. If the densities are independent, they can be expressed as the product of 1D densities

p left-parenthesis x comma y right-parenthesis equals p Subscript x Baseline left-parenthesis x right-parenthesis p Subscript y Baseline left-parenthesis y right-parenthesis comma

and random variables left-parenthesis upper X comma upper Y right-parenthesis can be found by independently sampling upper X from p Subscript x and upper Y from p Subscript y . Many useful densities are not separable, however, so we will introduce the theory of how to sample from multidimensional distributions in the general case.

Given a 2D density function, the marginal density function p left-parenthesis x right-parenthesis is obtained by “integrating out” one of the dimensions:

p left-parenthesis x right-parenthesis equals integral p left-parenthesis x comma y right-parenthesis normal d y period

This can be thought of as the density function for upper X alone. More precisely, it is the average density for a particular x over all possible y values.

If we can draw a sample upper X tilde p left-parenthesis x right-parenthesis , then—using Equation (2.1)—we can see that in order to sample upper Y , we need to sample from the conditional probability density, upper Y tilde p left-parenthesis y vertical-bar x right-parenthesis , which is given by:

p left-parenthesis y vertical-bar x right-parenthesis equals StartFraction p left-parenthesis x comma y right-parenthesis Over integral p left-parenthesis x comma y right-parenthesis normal d y EndFraction period

Sampling from higher-dimensional distributions can be performed in a similar fashion, integrating out all but one of the dimensions, sampling that one, and then applying the same technique to the remaining conditional distribution, which has one fewer dimension.

Sampling the Bilinear Function

The bilinear function

f left-parenthesis x comma y right-parenthesis equals left-parenthesis 1 minus x right-parenthesis left-parenthesis 1 minus y right-parenthesis w 0 plus x left-parenthesis 1 minus y right-parenthesis w 1 plus y left-parenthesis 1 minus x right-parenthesis w 2 plus x y w 3

interpolates between four values w Subscript i at the four corners of left-bracket 0 comma 1 right-bracket squared . ( w 0 is at left-parenthesis 0 comma 0 right-parenthesis , w 1 is at left-parenthesis 1 comma 0 right-parenthesis , w 2 at left-parenthesis 0 comma 1 right-parenthesis , and w 3 at left-parenthesis 1 comma 1 right-parenthesis .) After integration and normalization, we can find that its PDF is

p left-parenthesis x comma y right-parenthesis equals StartFraction 4 f left-parenthesis x comma y right-parenthesis Over w 0 plus w 1 plus w 2 plus w 3 EndFraction period

<<Sampling Inline Functions>>+=  
Float BilinearPDF(Point2f p, pstd::span<const Float> w) { if (p.x < 0 || p.x > 1 || p.y < 0 || p.y > 1) return 0; if (w[0] + w[1] + w[2] + w[3] == 0) return 1; return 4 * ((1 - p[0]) * (1 - p[1]) * w[0] + p[0] * (1 - p[1]) * w[1] + (1 - p[0]) * p[1] * w[2] + p[0] * p[1] * w[3]) / (w[0] + w[1] + w[2] + w[3]); }

The two dimensions of this function are not independent, so the sampling method samples a marginal distribution before sampling the resulting conditional distribution.

<<Sampling Inline Functions>>+=  
Point2f SampleBilinear(Point2f u, pstd::span<const Float> w) { Point2f p; <<Sample y for bilinear marginal distribution>> 
p.y = SampleLinear(u[1], w[0] + w[1], w[2] + w[3]);
<<Sample x for bilinear conditional distribution>> 
p.x = SampleLinear(u[0], Lerp(p.y, w[0], w[2]), Lerp(p.y, w[1], w[3]));
return p; }

We can choose either x or y to be the marginal distribution. If we choose y and integrate out x , we find that

StartLayout 1st Row 1st Column p left-parenthesis y right-parenthesis equals integral Subscript 0 Superscript 1 Baseline p left-parenthesis x comma y right-parenthesis normal d x 2nd Column equals 2 StartFraction left-parenthesis 1 minus y right-parenthesis left-parenthesis w 0 plus w 1 right-parenthesis plus y left-parenthesis w 2 plus w 3 right-parenthesis Over w 0 plus w 1 plus w 2 plus w 3 EndFraction 2nd Row 1st Column Blank 2nd Column proportional-to left-parenthesis 1 minus y right-parenthesis left-parenthesis w 0 plus w 1 right-parenthesis plus y left-parenthesis w 2 plus w 3 right-parenthesis period EndLayout

p left-parenthesis y right-parenthesis performs linear interpolation between two constant values, and so we can use SampleLinear() to sample from the simplified proportional function since it normalizes the associated PDF.

<<Sample y for bilinear marginal distribution>>= 
p.y = SampleLinear(u[1], w[0] + w[1], w[2] + w[3]);

Applying Equation (2.1) and again canceling out common factors, we have

p left-parenthesis x vertical-bar y right-parenthesis equals StartFraction p left-parenthesis x comma y right-parenthesis Over p left-parenthesis y right-parenthesis EndFraction proportional-to left-parenthesis 1 minus x right-parenthesis left-bracket left-parenthesis 1 minus y right-parenthesis w 0 plus y w 2 right-bracket plus x left-bracket left-parenthesis 1 minus y right-parenthesis w 1 plus y w 3 right-bracket comma

which can also be sampled in x using SampleLinear().

<<Sample x for bilinear conditional distribution>>= 
p.x = SampleLinear(u[0], Lerp(p.y, w[0], w[2]), Lerp(p.y, w[1], w[3]));

Because the bilinear sampling routine is based on the composition of two 1D linear sampling operations, it can be inverted by applying the inverses of those two operations in reverse order.

<<Sampling Inline Functions>>+=  
Point2f InvertBilinearSample(Point2f p, pstd::span<const Float> w) { return {InvertLinearSample(p.x, Lerp(p.y, w[0], w[2]), Lerp(p.y, w[1], w[3])), InvertLinearSample(p.y, w[0] + w[1], w[2] + w[3])}; }

See Section A.5 for further examples of multidimensional sampling algorithms, including techniques for sampling directions on the unit sphere and hemisphere, sampling unit disks, and other useful distributions for rendering.