# Why Monte Carlo Integration? Introduction to Monte Carlo Method

```Introduction to Monte Carlo
Method
IRISA
Email: [email protected]
Continuous Probability
• A continuous random variable x is a
variable that randomly takes on a value
from its domain
• The behavior of x is completely described
by the distribution of values it take.
Why Monte Carlo Integration?
• To generate realistic looking images, we
need to solve integrals of 2 or higher
dimension
– Pixel filtering and lens simulation both involve
solving a 2 dimensional integral
– Combining pixel filtering and lens simulation
requires solving a 4 dimensional integral
• Normal quadrature algorithms don’t extend
well beyond 1 dimension
Continuous Probability
• The distribution of values that x takes on
is described by a probability distribution
function p
– We say that x is distributed according to p, or
x~p
– p(x) ≥ 0
∞
– ∫− ∞ p ( x ) dx = 1
• The probability that x takes on a value
in
b
the interval [a, b] is: P( x ∈ [a, b ]) = ∫a p ( x ) dx
Multi-Dimensional Random
Variables
Expected Value
• The expected value of x~p is defined as
E ( x) = ∫ xp( x)dx
– As a function of a random variable is itself a
random variable, the expected value of f(x) is
E ( f ( x)) = ∫ f ( x) p( x)dx
• For some space S, we can define a pdf p:S →ℜ
• If x is a random variable, and x~p, the probability
that x takes on a value in Si ⊂ S is:
P( x ∈ Si ) = ∫ p( x)dµ
Si
• The expected value of a sum of random
variables is the sum of the expected
values:
E(x + y) = E(x) + E(y)
Monte Carlo Integration
• Expected value of a real valued function f :S →ℜ
extends naturally to the multidimensional case:
E ( f ( x )) = ∫ f ( x ) p ( x ) dµ
S
Monte Carlo Integration
• Suppose we have a function f(x) defined over
the domain x є [a, b]
• We would like to evaluate the integral
• Finding the estimated value of the
estimator we get: E[ I ] = E ⎡⎢ 1 ∑ f ( x ) ⎤⎥
N
i
⎣ N i =1 p ( xi ) ⎦
1 N ⎡ f ( xi ) ⎤
= ∑ E⎢
⎥
N i =1 ⎣ p ( xi ) ⎦
f ( xi )
1
= N∫
p ( xi )dx
N
p ( xi )
b
I = ∫ f ( x)dx
a
• The Monte Carlo approach is to consider N
samples, selected randomly with pdf p(x), to
estimate the integral
1 N f ( xi )
• We get the following estimator: I = ∑
N
i =1
p ( xi )
= ∫ f ( xi )dx = I
• In other words:
1
lim
N →∞ N
N
f ( xi )
∑ p( x ) = I
i =1
i
Variance
Variance Reduction
• Increase number of samples
• Choose p such that f/p has low variance (f and p
should have similar shape)
• The variance of the estimator is
σ2 =
1
N
⎛ f ( xi )
∫ ⎜⎜⎝ p( x ) − I
i
2
⎞
⎟⎟ p ( xi )dx
⎠
– This is called importance sampling because if p is
large when f is large, and small when f is small, there
will be more sample in important regions
• As the error in the estimator is proportional
to σ, the error is proportional to 1
• Partition the domain of the integral into several
smaller regions and evaluate the integral as a
the sum of integrals over the smaller regions
N
– So, to halve the error we need to use four
times as many samples
– This is called stratified sampling
Variance Reduction: Stratified
sampling
• Region D divided into disjoint sub-regions
∑iPi=1
• Di=[xi,xi+1]: sub-region
•
I =
Variance Reduction: Stratified
sampling
f (x)
• Ni samples per sub-region i
•
∫
0
x
f ( x ) dx = ∫ ( f ( x ) / p ( x )) p ( x ) dx +
1
0
x2
∫
( f ( x ) / p ( x )) p ( x ) dx +
x1
... +
1
∫(
x
n −1
f ( x ) / p ( x )) p ( x ) dx
n sub-regions
• p(x)/Pi : pdf for a sub-region i
( x) / p ( x)) p ( x)dx = Pi
∫Di ( f ( x ) / Pip ( x )) p ( x ) dx = 1
∫Di ( f
1
f (x)
0
1
( x) / p ( x)) p ( x)dx = Pi
∫Di ( f ( x ) / p ( x ))( p ( x ) / Pi ) dx = 1
∫Di ( f
I=
n −1
∑ Pi
i =0
f ( x) f ( x)
∫Di p( x).
Pi
0
i = n −1
dx
I≈∑
i =0
1
Pi
Ni
f ( Xki )
∑ p( X
Ni
ki = 0
ki
)
Variance Reduction: Stratified
sampling
• 2D example
•
•
•
•
•
•
•
•
•
→ N2 samples
•
•
•
•
•
•
•
Example
• Sample the pdf p(x)=3x2/2, x є [-1, 1]:
– First we need to find P(x): P( x) = ∫ 3x 2 / 2dx = x 3 / 2 + C
Sampling Random Variables
• Given a pdf p(x), defined over the interval
[xmin, xmax], we can sample a random variable
x~p from a set of uniform random numbers ξi є
x
[0,1)
prob(α < x) = P ( x) = ∫
xmin
• To do this, we need the cumulative probability
distribution function:
– To get xi we transform ξi: xi = P-1(ξi)
– P-1 is guaranteed to exist for all valid pdfs
Sampling 2D Random Variables
• If we have a 2D random variable α=(αx, αy) with pdf p(αx,
αy), we need the two dimensional cumulative pdf:
prob(α x < x & α y < y ) = P( x, y ) = ∫
y
∫
x
y min xmin
– Choosing C such that P(-1)=0 and P(1)=1, we
get P(x)=(x3+1)/2 and
P-1(ξ)= 3 2ξ − 1
– Thus, given a random number ξi є [0,1), we
can ’warp’ it according to p(x) to get xi: xi = 3 2ξ i − 1
p ( µ ) dµ
p ( µ x , µ y )dµ x dµ y
– We can choose xi using the marginal distribution pG(x)
and then choose yi according to p(y|x), where
y max
p ( x, y )
pG ( x) = ∫
p ( x, y )dy and p ( y | x) =
y min
pG ( x)
– If p is separable, that is p(x, y)=q(x,)r(y), the one
dimensional technique can be used on each
2D example: Uniform Sampling of
Triangles
2D example: Uniform Sampling of
Triangles
• To uniformly sample a triangle, we use barycentric
coordinates in a parametric space
• Let A,B,C the 3 vertices of a triangle
• Then a point P on the triangle is expressed as: C
– P= αA + βB + γC
1
– α+ β+ γ= 1
– α = 1- γ - β
P
– P =(1- γ - β ) A + βB + γC
γ
0
A
2D example: Uniform Sampling of
Triangles
• To uniformly sample a triangle, we use
barycentric coordinates
• Integrating the constant 1 across the
1
1−γ
triangle gives ∫γ =0 ∫β =0 dβdγ =0.5
– Thus our pdf is p(β,γ)=2
• Since β depends on γ (or γ depends on β),
we use the marginal density for γ, pG(γ):
1−γ '
pG (γ ' ) =
∫
0
2dβ = 2 − 2γ '
• To uniformly sample a triangle, we use
barycentric coordinates in a parametric space
(0,1)
β
B
(0,0)
γ
(1,0)
β +γ =1
2D example: Uniform Sampling of
Triangles
• From pG(γ’), we find p(β’|γ’) = p(γ’, β’)/pG(γ’) =
2/(2-2γ’) = 1/(1-γ’)
• To find γ’ we look at the cummulative pdf for
γ:
γ'
γ'
0
0
ξ1 = PG (γ ' ) = ∫ pG (γ )dγ = ∫ 2 − 2γdγ = 2γ '−γ '2
• Solving for γ’ we get
• We then turn to β’:
γ ' = 1 − 1 − ξ1
β'
β'
0
0
ξ 2 = P( β ' | γ ' ) = ∫ p(β ' | γ ' )dβ = ∫
1
β'
dβ =
1− γ '
1− γ '
2D example: Uniform Sampling of
Triangles
• Solving for β’ we get:
β ' = ξ 2 (1 − γ ' ) = ξ 2 (1 − (1 − 1 − ξ1 )) = ξ 2 1 − ξ1
• Thus given a set of random numbers ξ1
and ξ2, we warp these to a set of
barycentric coordinates sampling a
triangle: ( β , γ ) = (ξ 2 1 − ξ1 ,1 − 1 − ξ1 )
2D example: Uniform Sampling of
Discs
2D example: Uniform Sampling of
Discs
• Disc
– Generate random point on unit disk with
probability density:
p(x) = 1/(π.R2)
ϕ ∈ [0,2π ] et r ∈ [0, R ]
dµ = dS = rdrdφ
F (r , ϕ ) = ∫0ϕ ∫0r (r ' /(π .R 2 )dr ' dϕ '
– 2D CDF:
2D example: Uniform Sampling of
Spheres
• Sphere
• Disc
– Generate random point on unit disk with
probability density:
p(x) = 1/(π.R2)
ϕ ∈ [0,2π ] et r ∈ [0, R ]
dµ = dS = rdrdφ
2
– 2D CDF: F (r , ϕ ) = ∫0ϕ ∫0r (2r ' / R )(1 / 2π )dr ' dϕ '
– Compute marginal pdf and associated CDF
– ζ 1 and ζ 2 ∈ [0,1] are uniform random numbers
ϕ = 2πζ 1 and r = R ζ 2
θin
π /2
dωΘ in
θin
0
0
ϕin
ϕ in
2π
p(Θ)= cosθ ⇒θ =acos( ξ1 ) etϕ =2πξ2
π
Θ = (θ , ϕ ) ⇒ dµ = dωΘ = sin θ ⋅ dθ ⋅ dϕ
p(Θ)= n+1 cosnθ ⇒θ =acos(ξ1n+1) etϕ =2πξ2
2π
1
Summary
• Given a function f(µ), defined over an n-dimensional
domain S, we can estimate the integral of f over S by a
sum:
1 N f (x i )
(
)
f
d
µ
µ
≈
∑
∫
S
N
i =1
p(x i )
where x~p is a random variable and xi are samples of x
selected according to p.
• To reduce the variance and get faster convergence we:
– Use importance sampling: p should have similar
shape as f
– Use Stratified sampling: Subdivide S into smaller
regions, evaluate the integral for each region and sum
these together
```