18-660: Numerical Methods for
Engineering Design and Optimization
Xin Li
Department of ECE
Carnegie Mellon University
Pittsburgh, PA 15213
Slide 1
Overview
Monte Carlo Analysis
Latin hypercube sampling
Importance sampling
Slide 2
Latin Hypercube Sampling (LHS)
A great number of samples are typically required in traditional
Monte Carlo to achieve good accuracy
Various techniques exist to improve Monte Carlo accuracy
Controlling sampling points is the key
Latin hypercube sampling is a widely-used method to generate
controlled random samples
The basic idea is to make sampling point distribution close to
probability density function (PDF)
M. Mckay, R. Beckman and W. Conover, A comparison of three methods
for selecting values of input variables in the analysis of output from a
computer code, Technometrics, vol. 21, no. 2, pp. 239-245, May. 1979
Slide 3
Latin Hypercube Sampling (LHS)
One dimensional Latin hypercube sampling
Evenly partition CDF into N regions
Randomly pick up one sampling point in each region
One random sample in
one region
CDF
Evenly partitioned
regions
x
Avoid the probability that all sampling points come from
the same local region
Slide 4
Latin Hypercube Sampling (LHS)
Two dimensional Latin hypercube sampling
and x2 must be independent
Generate one-dimensional LHS samples for x1
Generate one-dimensional LHS samples for x2
Randomly combine the LHS samples to two-dimensional pairs
1
Probability of x2
x1
0.8
One sample in each row and each column
Sampling is random in each grid
0.6
0.4
0.2
0
0.2 0.4 0.6 0.8
Probability of x1
Higher-dimensional LHS samples can be
similarly generated
1
Slide 5
Latin Hypercube Sampling (LHS)
Matlab code for LHS sampling of independent standard
Normal distributions
data = rand(NSample,NVar);
for i = 1:NVar
index = randperm(NSample);
prob = (index'-data(:,i))/NSample;
data(:,i) = sqrt(2)*erfinv(2*prob-1);
end;
NVar:
NSample:
data:
# of random variables
# of samples
LHS sampling points
Slide 6
Latin Hypercube Sampling (LHS)
Compare Monte Carlo accuracy for a simple example
x ~ N(0,1) (standard Normal distribution)
Repeatedly estimate the mean value by Monte Carlo analysis
Estimated Mean Value
Accuracy is improved by increasing
sampling #
Estimated Mean Value
0.5
-0.5
0
10
10
# of Samples
Random sampling
10
0.5
-0.5
0
10
10
# of Samples
10
Latin hypercube sampling
Monte Carlo is not deterministic
Slide 7
Importance Sampling
Even with Latin hypercube sampling, Monte Carlo analysis
requires a HUGE number of sampling points
Example: rare event estimation
x ~ N (0,1)
Estimate
Standard Normal distribution
P( x 5) = ???
The theoretical answer for P(x -5) is equal to 2.8710-7
~100M sampling points are required if we attempt to estimate
this probability by random sampling or LHS
Slide 8
Importance Sampling
Key idea:
Do not generate random samples from pdfx(t)
Instead, find a good distorted pdfy(t) to improve Monte Carlo
sampling accuracy
Example: if x ~ N(0,1), what is P(x -5)?
Intuitively, if we draw sampling points based on pdfy(t), more
samples will fall into the grey area
pdfy(t)
PDF
pdfx(t)
How do we calculate P(x -5)
when sampling pdfy(t)?
P(x-5)
0
t
Slide 9
Importance Sampling
Assume that we want to estimate the following expected value
E [ f ( x )] =
f (t ) pdf (t ) dt
x
Example: if we want to estimate P(x -5), then
(x 5)
(x > 5)
1
f (x ) =
0
5
E [ f ( x )] = 1 pdf x (t ) dt + 0 pdf x (t ) dt = P( x 5)
Slide 10
Importance Sampling
Estimate E[f(x)] where x ~ pdfx(t) by importance sampling
pdf y (t )
)
(
)
(
)
(
)
(
)
(
[
]
E f x = f t pdf x t dt = f t pdf x t
dt
pdf (t )
+
pdf x (t )
E [ f (x )] = f (t )
pdf y (t ) dt = g (t ) pdf y (t ) dt = E [g ( x )]
pdf y (t )
Slide 11
Importance Sampling
Estimate E[f(x)] where x ~ pdfx(t) by traditional sampling
Step 1: draw M random samples {t1,t2,...,tM} based on pdfx(t)
Step 2: calculate fm = f(tm) at each sampling point m = 1,2,...,M
Step 3: calculate E[f] (f1+f2+...+fM)/M
Estimate E[f(x)] where x ~ pdfx(t) by importance sampling
Step 1: draw M random samples {t1,t2,...,tM} based on pdfy(t)
Step 2: calculate gm = f(tm)pdfx(tm)/pdfy(tm) at each sampling
point m = 1,2,...,M
Step 3: calculate E[f] (g1+g2+...+gM)/M
How do we decide the optimal pdfy(t) to achieve minimal
Monte Carlo analysis error?
Slide 12
Importance Sampling
Determine optimal pdfy(t) for importance sampling
pdf x (t m )
1 M
(
)
f tm
E[ f ] f =
M m =1
pdf y (t m )
Estimator
The accuracy of an estimator can be quantitatively measured
by its variance
[ ]
Error ~ VAR f
To improve Monte Carlo analysis accuracy, we should
minimize VAR[f]
Slide 13
Importance Sampling
Determine optimal pdfy(t) for importance sampling
pdf x (t m )
1 M
(
)
f tm
E[ f ] f =
M m =1
pdf y (t m )
We achieve the minimal VAR[f] = 0 if
f (t )
pdf y (t ) =
pdf x (t )
pdf y (t )
=k
f (t ) pdf x (t )
k
(Constant)
f(tm)pdfx(tm)/pdfy(tm) is equal to
the same constant for all
sampling points
(Optimal PDF)
Slide 14
Importance Sampling
pdf y (t ) =
f (t ) pdf x (t )
How do we decide the value k?
K cannot be arbitrarily selected
pdfy(t) must be a valid PDF that satisfies the following condition
pdf y (t ) dt =
f (t ) pdf x (t )
k
k=
dt = 1
f (t ) pdf (t ) dt = E[ f ]
Finding the optimal pdfy(t)
requires to know E[f], i.e., the
answer of our Monte Carlo
analysis!!!
Slide 15
Importance Sampling
In practice, such an optimal pdfy(t) cannot be easily applied
Instead, we typically look for a sub-optimal solution that
satisfies the following constraints
Easy to construct we do not have to know k = E[f]
Easy to sample not all random distributions can be easily
sampled by a random number generator
Minimal estimator variance the sub-optimal pdfy(t) is close to
the optimal case as much as possible
Slide 16
Importance Sampling
Finding the right pdfy(t) is nontrivial for practical problems
No magic equation exists in general
Engineering approach is based on heuristics
Sometimes require a lot of human experience and numerical
optimization
The criterion to choose pdfy(t) is also application-dependent
Slide 17
Importance Sampling
Example: if x ~ N(0,1), what is P(x -5)?
Shifted Normal
PDF
pdfx(t)
Uniform
P(x-5)
0
Several possible choices for pdfy(t)
Slide 18
Summary
Monte Carlo analysis
Latin hypercube sampling
Importance sampling
Slide 19
