Sampling from PDFs II CS295, Spring 2017 Shuang Zhao Computer - - PowerPoint PPT Presentation

Sampling from PDFs II

CS295, Spring 2017 Shuang Zhao

Computer Science Department University of California, Irvine

CS295, Spring 2017 Shuang Zhao 1

Announcements

• Additional readings on the course website
• Homework 1 due this Thursday (Apr 20)
• Programming Assignment 1 will be out on the

same day

CS295, Spring 2017 Shuang Zhao 2

Last Lecture

• Monte Carlo integration II
• Convergence properties
• Integrals over higher-dimensional domains
• Sampling from PDFs
• Inversion method

CS295, Spring 2017 Shuang Zhao 3

Today’s Lecture

• Sampling from PDFs II
• General methods
• Rejection sampling
• Metropolis-Hasting algorithm
• Alias method
• Sampling specific distributions
• Exponential distribution
• Normal distribution
• Distributions of directions (i.e., unit vectors)

CS295, Spring 2017 Shuang Zhao 4

Recap: Inversion Method

• Given a 1D distribution with CDF F, then

follows this given distribution

• Problem: some distrb.

(e.g., normal) have no closed-form CDFs

CS295, Spring 2017 Shuang Zhao 5

Rejection Sampling

• Consider a distribution over with PDF f
• Assume f is bounded so that
• Basic rejection sampling:
• 1. Draw x from U(Ω)
• 2. Draw y from U(0, M]
• 3. If

, return x

• 4. Otherwise, go to 1.
• Why valid?
• Accepted samples (x, y)

distribute uniformly over the subgraph of f(x)

CS295, Spring 2017 Shuang Zhao 6

Rejection Sampling

• To improve acceptance rate, use an envelop

distribution (that can be easily sampled)

• Given

, let M be a constant with . The sampling algorithm then becomes:

• 1. Draw x from
• 2. Draw y from U(0, M]
• 3. If

, return x

• 4. Otherwise, go to 1.
• can be made piecewise

for better performance

CS295, Spring 2017 Shuang Zhao 7

Rejection Sampling

• Generalize naturally to higher dimensions
• Assume pdf

and envelop over sample space :

• 1. Draw x from
• 2. Draw y from U(0, M]
• 3. If

, return x

• 4. Otherwise, go to 1.

CS295, Spring 2017 Shuang Zhao 8

Rejection Sampling

• Pros: only need the PDF f
• No CDF or its inverse
• Cons: not all generated samples are accepted
• Lower performance, especially for high-dimensional

problems (aka “curse of dimensionality”)

• In practice, only use rejection sampling when

there is no better option

CS295, Spring 2017 Shuang Zhao 9

Rejection Sampling

• The idea of rejecting generated samples can

be used in many other ways

• Example:
• To sample a conditional density

, one can generate samples according to and reject those not satisfying the condition y

• Metropolis-Hasting Algorithm

CS295, Spring 2017 Shuang Zhao 10

Metropolis-Hasting Algorithm

• A Markov-Chain Monte Carlo (MCMC) method
• Given a non-negative function f, generate a

chain of (correlated) samples X1, X2, X3, … that follow a probability density proportional to f

• Main advantage: f does not have to be a PDF

(i.e., unnormalized)

CS295, Spring 2017 Shuang Zhao 11

Metropolis-Hasting Algorithm

• Input
• Non-negative function f
• Probability density

suggesting a candidate for the next sample value x, given the previous sample value y

• The algorithm: given current sample Xi
• 1. Sample X’ from
• 2. Let

and draw

• 3. If

, set Xi+1 to X’; otherwise, set Xi+1 to Xi

• Start with arbitrary initial state X0

CS295, Spring 2017 Shuang Zhao 12

Metropolis-Hasting: Example

CS295, Spring 2017 Shuang Zhao 13

d = np.random.normal(scale=sigma) ang = np.pi*np.random.rand() X1 = X + np.array([d*np.cos(ang), d*np.sin(ang)]) a = f(X1)/f(X) # a only depends on f as g is symmetric if np.random.rand() < a: X = X1

Metropolis-Hasting: Example

CS295, Spring 2017 Shuang Zhao 14

Metropolis-Hasting: Issues

• The samples are correlated
• To obtain independent samples, have to take every

n-th sample: Xn, X2n, X3n, … for some n (determined

by examining autocorrelation of adjacent samples)

• Initial samples may follow a different

distribution

• Use a “burn-in” period by discarding the first few

samples (e.g., first 1000)

CS295, Spring 2017 Shuang Zhao 15

Metropolis Light Transport

• Applying the Metropolis-Hasting algorithm to

physically-based rendering

• To be discussed later in this course!

CS295, Spring 2017 Shuang Zhao 16

[Veach & Guibas 1997] [Kelemen et al. 2002] [Jakob & Marschner 2012]

Alias Method

• Efficient way to sample finite discrete

distributions with finite sample spaces

• Assuming
• O(N) preprocessing
• O(1) sampling
• Inversion method
• O(N) preprocessing (for creating the CMF)
• O(logN) sampling

CS295, Spring 2017 Shuang Zhao 17

Alias Method

• Observation: sampling a uniform distribution

with N outcomes is easy:

• Can we “convert” a non-uniform one into a uniform
• ne?
• Yes!

CS295, Spring 2017 Shuang Zhao 18

Alias Method

• All bins have the probability 1/N to be chosen
• Each bin involves at most two outcomes
• Construction process: O(N)
• Repeatedly choose a bin with a probability < 1/N and “steal”

the missing portion from another bin with probability > 1/N to form a bin with exactly two outcomes and probability 1/N

CS295, Spring 2017 Shuang Zhao 19

Bin 1 Bin 3 Bin 2

Today’s Lecture

• Sampling from PDFs II
• General methods
• Rejection sampling
• Metropolis-Hasting algorithm
• Alias method
• Sampling specific distributions
• Exponential distribution
• Normal distribution
• Distributions of directions (i.e., unit vectors)

CS295, Spring 2017 Shuang Zhao 20

Exponential Distribution

• Sample space:
• PDF:

for some λ>0

• CDF:
• Sampling method (inversion)
• More on this distribution later!

CS295, Spring 2017 Shuang Zhao 21

Normal Distribution

• Sample space:
• PDF:
• CDF:
• Inversion sampling is possible but requires

numerically inverting the error function erf

CS295, Spring 2017 Shuang Zhao 22

Normal Distribution

• Box-Muller Transform
• Let

be drawn independently from U(0, 1]

• Then,

both have standard normal distribution and are independent

CS295, Spring 2017 Shuang Zhao 23

Distributions of Directions

• Uniform position on a 2D disc
• Uniform position on the surface of a unit

sphere in 3D (i.e., )

• Uniform position on the surface of a unit

sphere in (n+1)-dimensional space (i.e., )

CS295, Spring 2017 Shuang Zhao 24

Uniform Distribution on 2D Disc

Solution 1

• Sample space:
• Using the idea of rejecting samples:

uniformly sample points in the bounding square and reject those not in the disc

• 1. Draw

from U(0, 1) 2.

• 3. If

, return x

• 4. Otherwise, go to 1.

CS295, Spring 2017 Shuang Zhao 25

1

Uniform Distribution on 2D Disc

Solution 2

• Sample space:

(polar coordinates)

• Due to symmetry,

. Next we focus on sampling r

• Observation: the probability

for a sampled point to locate within a disc with radius a is for all

CS295, Spring 2017 Shuang Zhao 26

1 a

Uniform Distribution on 2D Disc

Solution 2

• Thus, r can be sampled with inversion method:
• Putting everything together:
• 1. Draw

from U[0, 1) 2.

(polar coordinates)

3.

(Cartesian coordinates)

CS295, Spring 2017 Shuang Zhao 27

Uniform Distribution on

Solution 1

• Sample space:
• Using the idea of rejecting samples:

uniformly sample points in the bounding cube and reject those not in the disc

• 1. Draw

from U(0,1) 2.

• 3. If

, return

• 4. Otherwise, go to 1.

CS295, Spring 2017 Shuang Zhao 28

Uniform Distribution on

Solution 2

• Sample space:

(spherical coordinates)

• Due to symmetry,

. Next we focus on sampling θ

• Observation: the probability for

a sampled point to locate within a spherical cap with angle θ’ is

CS295, Spring 2017 Shuang Zhao 29

Surface area: 2πr2(1 – cosθ)

Uniform Distribution on

Solution 2

• Thus, θ can be sampled with inversion method:
• Putting everything together:
• 1. Draw

from U[0, 1) 2.

(spherical coord.)

3.

(Cartesian coord.)

CS295, Spring 2017 Shuang Zhao 30

Uniform Distribution on

• Sample space:
• General approach

1. Draw from N(0, 1) (standard normal) 2. 3. Return

• Why is this valid?
• x has the multivariate normal distribution with identity

covariance matrix

• This distribution is rotationally symmetric around the
• rigin

CS295, Spring 2017 Shuang Zhao 31

Next Lecture

• The rendering equation
• Monte Carlo path tracing I

CS295, Spring 2017 Shuang Zhao 32