第九章
Monte Carlo积分
第九章 Monte Carlo积分
Monte Carlo法的重要应用领域之一：计算积分和多重积分
适用于求解,
1,被积函数、积分边界复杂，难以用解析方法或一般的数
值方法求解；
2,被积函数的具体形式未知，只知道由模拟返回的函数值。
本章内容,
用 Monte Carlo法求定积分的几种方法,
均匀投点法、期望值估计法、重要抽样法、半解析
法,…
第九章 Monte Carlo积分
Goal,Evaluate an integral,
??
b
a
I dx g( x )
Why use random methods?
Computation by,deterministic quadrature”
can become expensive and inaccurate,
? grid points add up quickly in high dimensions
? bad choices of grid may misrepresent g(x)
第九章 Monte Carlo积分
?Monte Carlo method can be used to compute
integral of any dimension d (d-fold integrals)
?Error comparison of d-fold integrals
Simpson’s rule,…
dNE /1??
21?? NE
purely statistical,
not rely on the dimension !
Monte Carlo method WINS,when d >> 3
Monte Carlo method
approximating the integral of a
function f using quadratic polynomials
)]()(4)([31)()( 210000 xfxfxfhdxxfdxxf hxxxx ???? ?? ?
hxxxx ???? 1201
第九章 Monte Carlo积分
?Hit-or-Miss Method
?Sample Mean Method
?Variance Reduction Technique
?Variance Reduction using Rejection Technique
?Importance Sampling Method
Hit-or-Miss Method
?Evaluation of a definite integral
dxxI ba?? )(?
xxh a n y f o r )(??
a b
h
X
X
X
X
X X
O
O
O O
O
O
O ?Probability that a random point
reside inside the area
N
M
hab
Ip ?
?? )(
N
MhabI )( ??
N, Total number of points
M, points that reside inside
the region
Hit-or-Miss Method
Sample uniformly from the rectangular region
[a,b]x[0,h]
a)-h ( b
I, p ?
The probability that we are below the curve is
So,if we can estimate p,we can estimate I,
a)-h ( bpI ?? ?
where is our estimate of p p?
Hit-or-Miss Method
We can easily estimate p,
? throw N,uniform darts” at the rectangle
N
Mp ??? let
? let M be the number of times you end up under
the curve y=g(x)
Hit-or-Miss Method
a b
h
X
X
X
X
X X
O
O
O O
O
O
O
Start
Set N, large integer
M = 0
Choose a point x in [a,b]
Choose a point y in [0,h]
if [x,y] reside inside then M = M+1
I = (b-a) h (M/N)
End
Loop
N times auabx ??? 1)(
2huy ?
Hit-or-Miss Method
Error Analysis of the Hit-or-Miss Method
It is important to know how accurate the result of
simulations are
note that M is binomial(M,p)
)1()()( 2 PNpMNpME ??? ?
? ?
IpabhME
N
abh
abh
N
MEabhpEIE
?????
??
?
??
? ????
)()()(
)()(?)?( ? ?
N
ppabhM
N
abh
abh
N
MabhpI
)1()()()(
)()(?)?(
22
2
2
22
222
?????
??
?
??
? ????
?
???
NMpp ?? ? )1()()?( NMMN abhI ????
21)1()()?( ????? NN ppabhI?
第九章 Monte Carlo积分
?Hit-or-Miss Method
?Sample Mean Method
?Variance Reduction Technique
?Variance Reduction using Rejection Technique
?Importance Sampling Method
Sample Mean Method
Start
Set N, large integer
s1 = 0,s2 = 0
xn = (b-a) un + a
yn = ?(xn)
s1 = s1 + yn,s2 = s2 + yn2
Estimate mean m’=s1/N
Estimate variance V’ = s2/N – m’2
End
Loop
N times
N
VabI SM
e r r o r
'6 7 4 5.0)( ???
Sample Mean Method
dx g ( x ) I
b
a
??
Write this as,
dx
a-b
1 g ( x )a)-(b d x g ( x ) I b
a
b
a
?? ??
? ?g ( X )E a)-(b ? where X~unif(a,b)
Sample Mean Method
? ?g ( X )E a)-(b I ? where X~unif(a,b)
So,we will estimate I by estimating E[g(X)] with
?
?
?
n
1i
) ? ig (Xn1 [g (X )]E
where X1,X2,…,X n is a random sample from the
uniform(a,b) distribution,
Sample Mean Method
dx e I
3
0
x??
Example,
(we know that the answer is e3-1 19.08554)
? write this as
]3E [ e d x
3
1e3 I X3
0
x ?? ?
where X~unif(0,3)
Sample Mean Method
? write this as
]3E [ e d x
3
1e3 I X3
0
x ?? ?
where X~unif(0,3)
estimate this with
?
?
??
n
1i
xX ie
n
1 3][eE3 I ??
where X1,X2,…,X n are n independent unif(0,3)’s,
Sample Mean Method
Simulation Results,true = 19.08554,n=100,000
1 19.10724
I?
2 19.08260
3 18.97227
4 19.06814
5 19.13261
Simulation
Sample Mean Method
Don’t ever give an estimate without a
confidence interval!
This estimator is,unbiased”,
??
?
??
?? ?
?
n
1i
i )g ( Xn
1a)-(bE ]IE[ ?
?
?
?
n
1i
i )]E [ g (Xn
1a)-(b
??
b
a
dx
a-b
1 g ( x )
n
1a)-(b ?? b
a
d x g ( x ) I ?
Sample Mean Method
]IV a r [, σ 2I ?? ? ?
?
?
??
?? ?
?
n
1i
i )g ( X n
1 a)-(bV a r
??
?
??
?? ?
?
n
1i
i2
2
)g ( XV a rn a)-(b
( i n d e p ) )]V a r [ g ( Xn a)-(b
n
1i
i2
2
?
?
?
( i n d e n t ) )]V a r [ g ( X n a)-(b 1
2
?
? ?? ??
b
a
2
2
dx
a-b
1 E [ g ( X ) ]g ( x )
n
a)-(b
? ?????? ??
b
a
22
dx
a-b
1
a-b
Ig ( x )
n
a)-(b
Sample Mean Method
? an approximation
1n
ab
I
)g (x
n
a)-(b
s
n
1i
i2
2
I ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
Sample Mean Method
?
?
?
n
1i
i )g (Xn
1a)-(b I?
? X1,X2,…,X n iid -> g(X1),g(X2),…,g(X n) iid
? Let Yi=g(Xi) for i=1,2,…,n
Then
Y a)-(b I ??
and we can once again invoke the CLT,
Sample Mean Method
)N ( I,I 2I?? ??
For n,large enough” (n>30),
So,a confidence interval for I is roughly given by
I/2 z I ?? ???
but since we don’t know,we’ll have to be
content with the further approximation,I?
?
I/2 s z I ?? ??
Sample Mean Method
? ??
b
a
2x1 / 2 dx e x,I
By the way…
No one ever said that you have to use the uniform
distribution
Example,
?
?
??
0
b][ a,
2x1 / 2 dx ( x )I 2e x
2
1
? ?( X )I XE21 b][ a,1 / 2?
where X~exp(rate=2),
Sample Mean Method
Comparison of Hit-and-Miss and Sample Mean Monte Carlo
? Let be the hit-and-miss estimator of I
HMI?
)IV a r ( )IV a r ( SMHM ?? ?
Then
? Let be the sample mean estimator of I
SMI?
Sample Mean Method
Comparison of Hit-and-Miss and Sample Mean Monte Carlo
Sample mean Monte Carlo is generally preferred over Hit-
and-Miss Monte Carlo because,
? the estimator from SMMC has lower variance
? SMMC does not require a non-negative integrand
(or adjustments)
? H&M MC requires that you be able to put g(x) in a
,box”,so you need to figure out the max
value of g(x) over [a,b] and you need to be
integrating over a finite integral,
2.1 Variance Reduction Technique - Introduction
第九章 Monte Carlo积分
?Hit-or-Miss Method
?Sample Mean Method
?Variance Reduction Technique
?Variance Reduction using Rejection Technique
?Importance Sampling Method
Variance Reduction Technique
Introduction
Monte Carlo Method and Sampling Distribution
Monte Carlo Method, Take values from random sample
From central limit theorem,
3? rule
Most probable error
Important characteristics
N/ 22 ??mm ??
9973.0)33( ????? ?m?m XP
NE r r o r
?6745.0??
NE r r o r /1? ??Error
Variance Reduction Technique
Introduction
Reducing error
*100 samples reduces the error order of 10
Reducing variance ? Variance Reduction
Technique
The value of variance is closely related to how
samples are taken
Unbiased sampling
Biased sampling
More points are taken in important parts of the population
Variance Reduction Technique
Motivation
If we are using sample-mean Monte Carlo Method
Variance depends very much on the behavior of ?(x)
?(x) varies little ? variance is small
?(x) = const ? variance=0
Evaluation of a integral
Near minimum points ? contribute less to the summation
Near maximum points ? contribute more to the summation
More points are sampled near the peak ?”importance
sampling strategy”
?
?
???? N
n
nY xN
ababI
1
)()(' ?m
第九章 Monte Carlo积分
?1.2 Hit-or-Miss Method
?1.3 Sample Mean Method
?2.1 Variance Reduction Technique
?2.3 Variance Reduction using Rejection Technique
?2.4 Importance Sampling Method
Variance Reduction Technique
Variance Reduction for Hit-or-Miss method
?In the domain [a,b] choose a comparison function
dxxwA
dxxwxW
xxw
b
a
x
?
?
?
?
?
??
)(
)()(
)()( ?
)( xwAu ? )(1 AuWx ??
a b
w(x)
X
X
X
X X
O
O
O O
O
O
O
?(x)
Points are generated on the area under w(x)
function
Random variable that follows distribution w(x)
Variance Reduction Technique
Points lying above ?(x) is rejected
N
NAI '?
nnn uxwy )(?
???
?
?
??
)( if 0
)( if 1
nn
nn
n xy
xyq
?
?
q 1 0
P(q) r 1-r AIr /?
a b
w(x)
X
X
X
X X
O
O
O O
O
O
O
?(x)
Variance Reduction Technique
Error Analysis
)1()(,)( rrQVrQE ???
)( QAEArI ??
N
IAII RJ
e r r o r
)(67.0 ??
Hit or Miss method
Error reduction
habA )( ??
0E r r o r th e n ?? IA
Variance Reduction Technique
Start
Set N, large integer
N’ = 0
Generate u1,x= W-1(Au1)
Generate u2,y=u2 w(x)
If y<= f(x) accept value N’ = N’+1
Else, reject value
I = (b-a) h (N’/N)
End
Loop
N times
N
IAII RJ
e r r o r
)'('67.0 ??
a b
w(x)
X
X
X
X X
O
O
O O
O
O
O
?(x)
第九章 Monte Carlo积分
?1.2 Hit-or-Miss Method
?1.3 Sample Mean Method
?2.1 Variance Reduction Technique
?2.3 Variance Reduction using Rejection Technique
?2.4 Importance Sampling Method
Importance Sampling Method
Basic idea
Put more points near maximum
Put less points near minimum
F(x), transformation function (or
weight function_
? ??? x dxxfxF )()(
)(
)(
1 yFx
xFy
??
?
Importance Sampling Method
)(/)(/ xfdydxxfdxdy ???
?? ???????? baba dxxfxf xdyxf xI )()( )()( )( ??
???? baf dxxfx )()(??
)(
)()(
xf
xx ?? ?
if we choose f(x) = c?(x),then variance will be small
The magnitude of error depends on the choice of f(x)
f
b
a dxxfxI ???? ? ?? )()(
Importance Sampling Method
Estimate of error
?
?
???? N
n
nf xNI
1
)(1 ??
N
VI f
e r r o r
)(67.0 ??
22 )()( fffV ?????? ???
N
II fIS
e r r o r
22
67.0 ???? ?
Importance Sampling Method
Start
Set N, large integer
s1 = 0,s2 = 0
Generate xn according to f(x)
?n = ? (xn) / f ( xn)
Add ?n to s1
Add ?n to s2
I ‘ =s1/N,V’=s2/N-I’2
End
Loop
N times
N
VI IS
e rro r
'67.0?
Monte Carlo积分
第九章 Monte Carlo积分
Monte Carlo法的重要应用领域之一：计算积分和多重积分
适用于求解,
1,被积函数、积分边界复杂，难以用解析方法或一般的数
值方法求解；
2,被积函数的具体形式未知，只知道由模拟返回的函数值。
本章内容,
用 Monte Carlo法求定积分的几种方法,
均匀投点法、期望值估计法、重要抽样法、半解析
法,…
第九章 Monte Carlo积分
Goal,Evaluate an integral,
??
b
a
I dx g( x )
Why use random methods?
Computation by,deterministic quadrature”
can become expensive and inaccurate,
? grid points add up quickly in high dimensions
? bad choices of grid may misrepresent g(x)
第九章 Monte Carlo积分
?Monte Carlo method can be used to compute
integral of any dimension d (d-fold integrals)
?Error comparison of d-fold integrals
Simpson’s rule,…
dNE /1??
21?? NE
purely statistical,
not rely on the dimension !
Monte Carlo method WINS,when d >> 3
Monte Carlo method
approximating the integral of a
function f using quadratic polynomials
)]()(4)([31)()( 210000 xfxfxfhdxxfdxxf hxxxx ???? ?? ?
hxxxx ???? 1201
第九章 Monte Carlo积分
?Hit-or-Miss Method
?Sample Mean Method
?Variance Reduction Technique
?Variance Reduction using Rejection Technique
?Importance Sampling Method
Hit-or-Miss Method
?Evaluation of a definite integral
dxxI ba?? )(?
xxh a n y f o r )(??
a b
h
X
X
X
X
X X
O
O
O O
O
O
O ?Probability that a random point
reside inside the area
N
M
hab
Ip ?
?? )(
N
MhabI )( ??
N, Total number of points
M, points that reside inside
the region
Hit-or-Miss Method
Sample uniformly from the rectangular region
[a,b]x[0,h]
a)-h ( b
I, p ?
The probability that we are below the curve is
So,if we can estimate p,we can estimate I,
a)-h ( bpI ?? ?
where is our estimate of p p?
Hit-or-Miss Method
We can easily estimate p,
? throw N,uniform darts” at the rectangle
N
Mp ??? let
? let M be the number of times you end up under
the curve y=g(x)
Hit-or-Miss Method
a b
h
X
X
X
X
X X
O
O
O O
O
O
O
Start
Set N, large integer
M = 0
Choose a point x in [a,b]
Choose a point y in [0,h]
if [x,y] reside inside then M = M+1
I = (b-a) h (M/N)
End
Loop
N times auabx ??? 1)(
2huy ?
Hit-or-Miss Method
Error Analysis of the Hit-or-Miss Method
It is important to know how accurate the result of
simulations are
note that M is binomial(M,p)
)1()()( 2 PNpMNpME ??? ?
? ?
IpabhME
N
abh
abh
N
MEabhpEIE
?????
??
?
??
? ????
)()()(
)()(?)?( ? ?
N
ppabhM
N
abh
abh
N
MabhpI
)1()()()(
)()(?)?(
22
2
2
22
222
?????
??
?
??
? ????
?
???
NMpp ?? ? )1()()?( NMMN abhI ????
21)1()()?( ????? NN ppabhI?
第九章 Monte Carlo积分
?Hit-or-Miss Method
?Sample Mean Method
?Variance Reduction Technique
?Variance Reduction using Rejection Technique
?Importance Sampling Method
Sample Mean Method
Start
Set N, large integer
s1 = 0,s2 = 0
xn = (b-a) un + a
yn = ?(xn)
s1 = s1 + yn,s2 = s2 + yn2
Estimate mean m’=s1/N
Estimate variance V’ = s2/N – m’2
End
Loop
N times
N
VabI SM
e r r o r
'6 7 4 5.0)( ???
Sample Mean Method
dx g ( x ) I
b
a
??
Write this as,
dx
a-b
1 g ( x )a)-(b d x g ( x ) I b
a
b
a
?? ??
? ?g ( X )E a)-(b ? where X~unif(a,b)
Sample Mean Method
? ?g ( X )E a)-(b I ? where X~unif(a,b)
So,we will estimate I by estimating E[g(X)] with
?
?
?
n
1i
) ? ig (Xn1 [g (X )]E
where X1,X2,…,X n is a random sample from the
uniform(a,b) distribution,
Sample Mean Method
dx e I
3
0
x??
Example,
(we know that the answer is e3-1 19.08554)
? write this as
]3E [ e d x
3
1e3 I X3
0
x ?? ?
where X~unif(0,3)
Sample Mean Method
? write this as
]3E [ e d x
3
1e3 I X3
0
x ?? ?
where X~unif(0,3)
estimate this with
?
?
??
n
1i
xX ie
n
1 3][eE3 I ??
where X1,X2,…,X n are n independent unif(0,3)’s,
Sample Mean Method
Simulation Results,true = 19.08554,n=100,000
1 19.10724
I?
2 19.08260
3 18.97227
4 19.06814
5 19.13261
Simulation
Sample Mean Method
Don’t ever give an estimate without a
confidence interval!
This estimator is,unbiased”,
??
?
??
?? ?
?
n
1i
i )g ( Xn
1a)-(bE ]IE[ ?
?
?
?
n
1i
i )]E [ g (Xn
1a)-(b
??
b
a
dx
a-b
1 g ( x )
n
1a)-(b ?? b
a
d x g ( x ) I ?
Sample Mean Method
]IV a r [, σ 2I ?? ? ?
?
?
??
?? ?
?
n
1i
i )g ( X n
1 a)-(bV a r
??
?
??
?? ?
?
n
1i
i2
2
)g ( XV a rn a)-(b
( i n d e p ) )]V a r [ g ( Xn a)-(b
n
1i
i2
2
?
?
?
( i n d e n t ) )]V a r [ g ( X n a)-(b 1
2
?
? ?? ??
b
a
2
2
dx
a-b
1 E [ g ( X ) ]g ( x )
n
a)-(b
? ?????? ??
b
a
22
dx
a-b
1
a-b
Ig ( x )
n
a)-(b
Sample Mean Method
? an approximation
1n
ab
I
)g (x
n
a)-(b
s
n
1i
i2
2
I ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
Sample Mean Method
?
?
?
n
1i
i )g (Xn
1a)-(b I?
? X1,X2,…,X n iid -> g(X1),g(X2),…,g(X n) iid
? Let Yi=g(Xi) for i=1,2,…,n
Then
Y a)-(b I ??
and we can once again invoke the CLT,
Sample Mean Method
)N ( I,I 2I?? ??
For n,large enough” (n>30),
So,a confidence interval for I is roughly given by
I/2 z I ?? ???
but since we don’t know,we’ll have to be
content with the further approximation,I?
?
I/2 s z I ?? ??
Sample Mean Method
? ??
b
a
2x1 / 2 dx e x,I
By the way…
No one ever said that you have to use the uniform
distribution
Example,
?
?
??
0
b][ a,
2x1 / 2 dx ( x )I 2e x
2
1
? ?( X )I XE21 b][ a,1 / 2?
where X~exp(rate=2),
Sample Mean Method
Comparison of Hit-and-Miss and Sample Mean Monte Carlo
? Let be the hit-and-miss estimator of I
HMI?
)IV a r ( )IV a r ( SMHM ?? ?
Then
? Let be the sample mean estimator of I
SMI?
Sample Mean Method
Comparison of Hit-and-Miss and Sample Mean Monte Carlo
Sample mean Monte Carlo is generally preferred over Hit-
and-Miss Monte Carlo because,
? the estimator from SMMC has lower variance
? SMMC does not require a non-negative integrand
(or adjustments)
? H&M MC requires that you be able to put g(x) in a
,box”,so you need to figure out the max
value of g(x) over [a,b] and you need to be
integrating over a finite integral,
2.1 Variance Reduction Technique - Introduction
第九章 Monte Carlo积分
?Hit-or-Miss Method
?Sample Mean Method
?Variance Reduction Technique
?Variance Reduction using Rejection Technique
?Importance Sampling Method
Variance Reduction Technique
Introduction
Monte Carlo Method and Sampling Distribution
Monte Carlo Method, Take values from random sample
From central limit theorem,
3? rule
Most probable error
Important characteristics
N/ 22 ??mm ??
9973.0)33( ????? ?m?m XP
NE r r o r
?6745.0??
NE r r o r /1? ??Error
Variance Reduction Technique
Introduction
Reducing error
*100 samples reduces the error order of 10
Reducing variance ? Variance Reduction
Technique
The value of variance is closely related to how
samples are taken
Unbiased sampling
Biased sampling
More points are taken in important parts of the population
Variance Reduction Technique
Motivation
If we are using sample-mean Monte Carlo Method
Variance depends very much on the behavior of ?(x)
?(x) varies little ? variance is small
?(x) = const ? variance=0
Evaluation of a integral
Near minimum points ? contribute less to the summation
Near maximum points ? contribute more to the summation
More points are sampled near the peak ?”importance
sampling strategy”
?
?
???? N
n
nY xN
ababI
1
)()(' ?m
第九章 Monte Carlo积分
?1.2 Hit-or-Miss Method
?1.3 Sample Mean Method
?2.1 Variance Reduction Technique
?2.3 Variance Reduction using Rejection Technique
?2.4 Importance Sampling Method
Variance Reduction Technique
Variance Reduction for Hit-or-Miss method
?In the domain [a,b] choose a comparison function
dxxwA
dxxwxW
xxw
b
a
x
?
?
?
?
?
??
)(
)()(
)()( ?
)( xwAu ? )(1 AuWx ??
a b
w(x)
X
X
X
X X
O
O
O O
O
O
O
?(x)
Points are generated on the area under w(x)
function
Random variable that follows distribution w(x)
Variance Reduction Technique
Points lying above ?(x) is rejected
N
NAI '?
nnn uxwy )(?
???
?
?
??
)( if 0
)( if 1
nn
nn
n xy
xyq
?
?
q 1 0
P(q) r 1-r AIr /?
a b
w(x)
X
X
X
X X
O
O
O O
O
O
O
?(x)
Variance Reduction Technique
Error Analysis
)1()(,)( rrQVrQE ???
)( QAEArI ??
N
IAII RJ
e r r o r
)(67.0 ??
Hit or Miss method
Error reduction
habA )( ??
0E r r o r th e n ?? IA
Variance Reduction Technique
Start
Set N, large integer
N’ = 0
Generate u1,x= W-1(Au1)
Generate u2,y=u2 w(x)
If y<= f(x) accept value N’ = N’+1
Else, reject value
I = (b-a) h (N’/N)
End
Loop
N times
N
IAII RJ
e r r o r
)'('67.0 ??
a b
w(x)
X
X
X
X X
O
O
O O
O
O
O
?(x)
第九章 Monte Carlo积分
?1.2 Hit-or-Miss Method
?1.3 Sample Mean Method
?2.1 Variance Reduction Technique
?2.3 Variance Reduction using Rejection Technique
?2.4 Importance Sampling Method
Importance Sampling Method
Basic idea
Put more points near maximum
Put less points near minimum
F(x), transformation function (or
weight function_
? ??? x dxxfxF )()(
)(
)(
1 yFx
xFy
??
?
Importance Sampling Method
)(/)(/ xfdydxxfdxdy ???
?? ???????? baba dxxfxf xdyxf xI )()( )()( )( ??
???? baf dxxfx )()(??
)(
)()(
xf
xx ?? ?
if we choose f(x) = c?(x),then variance will be small
The magnitude of error depends on the choice of f(x)
f
b
a dxxfxI ???? ? ?? )()(
Importance Sampling Method
Estimate of error
?
?
???? N
n
nf xNI
1
)(1 ??
N
VI f
e r r o r
)(67.0 ??
22 )()( fffV ?????? ???
N
II fIS
e r r o r
22
67.0 ???? ?
Importance Sampling Method
Start
Set N, large integer
s1 = 0,s2 = 0
Generate xn according to f(x)
?n = ? (xn) / f ( xn)
Add ?n to s1
Add ?n to s2
I ‘ =s1/N,V’=s2/N-I’2
End
Loop
N times
N
VI IS
e rro r
'67.0?
