蒙特卡洛方法的应用2

f ( xi ) 2 1 e ˆ n i 1 g ( xi ) n i 1 1 xi
n
n
xi
g(x)的随机数对应分布函数为
0, 1 1 2 Fg ( x) ( x x ) 1 / 4, 2 2 1, x0 1 x 1 x 1
估计步骤：
重要抽样法 关键因素在于g(x)的选取，使得估计的方差较 小。 重要抽样法的基本思想，就是通过选取与f(x) 形状接近的密度函数g(x)来降低估计的方差。
例 利用Monte Carlo方法计算一个简单的积分

1 x e dx 0
(e 1)
(1) 先考虑样本平均值法:

xi n f (x ) n 1 3 e i （ 3） ˆ n i 1 g ( xi ) 2n i 1 1 xi
function result=zycy2( mm) %积分函数 %mm 是随机实验次数 sum=0; u = unifrnd(0,1,1,mm); xrandnum = -1+sqrt(1+3.* u); for ii=1:mm sum=sum+exp (xrandnum(1,ii))/(1+ xrandnum(1,ii)); end result=1.5*sum/mm; result=zycy2(0,1, 1000)=1.7222 result=zycy1(0,1, 10000)=1.7174 result=zycy1(0,1, 100000)=1.7185 True value is 1.7183
例 利用Monte Carlo方法计算一个简单的积分

1 x e dx 0
(e 1)
(2) 重要抽样法:
由重要抽样法思想，要选择一个与ex相似的 密度函数. 我们知道，ex的Taylor展开为
2 k x x x e 1 x ... ... 2! k!
利用线性近似，取（0，1）上密度函数
练习:用重要抽样法计算
I e dx
x 1
1
x x e 1 x ... ... 2! k!
x
2
k
1 g ( x) (1 x) 2

b
a
f X f x g x dx E g x g X
设x1,…,xn是来自g(x)的随机数，则 的估计为
几种降低估计方差的MC方法
重要抽样法
特点：相对样本均值法而言，样本均值法是由
于假设g(x)是均匀分布的概率密度，故采用的 是均匀抽样，各随机数xi是均匀分布的随机数， 各xi 对 ˆ 的贡献是不同，f(xi) 大则贡献大，但 在抽样时，这种差别未能体现出来。 而重要抽样法，则希望贡献率大的随机数出现 的概率大，贡献小的随机数出现概率小，从而 提高抽样的效率。
function Rguji=litiR4(t,thetaa1,thetaa2,thetab1,thetab2,mm)
%t 是要求系统生存的寿命%thetaa1 是元件A1的数学期望%thetaa2 是元件A2的数学期望 %thetab1 是元件B1的数学期望 %thetab2 是元件B2的数学期望%mm 是随机实验次数
例 设系统 L 由相互独立的 n 个元件组成，连 接方式为 (1) 串联； (2) 并联； (3) 冷贮备(起初由一个元件工作，其它 n – 1 个元件做冷贮备，当工作元件失效时， 贮备的元件件 中有 k 个或 k 个以上的元件正常工作时， 系统 L 才正常工作)
(1)
X min{ X 1 , X 2 ,, X n }
FX ( x) 1 (1 FX ( x))
i 1
i
n
e , x 0, 1 FX ( x) x0 1,
i
x
ne , x 0 f X ( x) x0 0,
nx
(2) X max{ X 1 , X 2 ,, X n }
sum=0;
u = unifrnd(0,1,1,mm); xrandnum = 2*sqrt(u )-1; for ii=1:mm sum=sum+exp (xrandnum(1,ii))/(1+ xrandnum(1,ii)); end I=2*sum/mm;
function I=exp_3_2( mm) %重要抽样法 %mm 是随机实验次数 u = unifrnd(0,1,1,mm); xrandnum = 2*sqrt(u )-1; s=sum(exp (xrandnum)./(1+ xrandnum)); I=2*s/mm;
所以, 估计方差的大小与I1，I2 的估计的相关 度有关，若两者的正相关程度越高，则 的估计 方差越小。这便是关联抽样法的基本出发点。
2、系统的可靠性计算问题
一个元件(或系统)能正常工作的概率称为 元件(或系统)的可靠性 系统由元件组成,常见的元件连接方式： 串联
1 1 2
并联
2
例R3设两系统都是由 4 个元件组成,每个元件的 寿命服从参数为θ的指数分布，每个元件是否正 常工作相互独立.两系统的连接方式如下图所示， 求两系统寿命大于T=100的概率. A2 A1 S1: B1 B2
xi n f (x ) n 1 3 e i ˆ n i 1 g ( xi ) 2n i 1 1 xi
g(x)的随机数对应分布函数为
0, 1 Fg ( x) (2 x x 2 ), 3 1, x0 0 x 1 x 1
估计步骤：
（1）产生n个U(0,1)随机数u1,…,un, 则 （2）xi= 1 1 3ui
1 x e 0
f ( x)dx
f(x)=1, 0<x<1,为U(0,1)对应的概率密度. 由此 产生n个U(0,1)随机数x1,…,xn, 则
n 1 ˆ e xi n i 1
function result=zycy1(a,b, mm) %a是积分的下限 %b是积分的上限 %积分函数 %mm 是随机实验次数 sum=0; xrandnum = unifrnd(a,b,1,mm); for ii=1:mm sum=sum+exp (xrandnum(1,ii)); end result=sum/mm result=zycy1(0,1, 1000)=1.7267 result=zycy1(0,1, 10000)=1.7199 result=zycy1(0,1, 100000)=1.7171 True value is 1.7183
关联抽样法 将需要估计的积分分解成两个积分之差，
f x dx f1 x dx f 2 x dx I1 I 2
b b b a a a
对的估计转化为对I1，I2 的估计的差。即
ˆI ˆ I ˆ 1 2
由于
ˆ) D( I ˆ ) D( I ˆ ) 2 ˆ ˆ D( 1 2 I1I 2 DI1 DI 2
frq=0;randnuma1 = exprnd(thetaa1,1,mm); randnuma2 = exprnd(thetaa2,1,mm); randnumb1 = exprnd(thetab1,1,mm); randnumb2 = exprnd(thetab2,1,mm); for ii=1:mm if (randnuma1(1,ii)>t)|(randnumb1(1,ii)>t) pass1=1; else pass1=0; end if (randnuma2(1,ii)>t)|(randnumb2(1,ii)>t) pass2=1; else pass2=0; end if (pass1*pass2)==1 frq=frq+1; end end，Rguji=frq/mm
1 2 1 2
x
x t ( x t ) e e dt , x 0 0 0, x0
2 xe x , x 0 x0 0,
t
（1）产生n个U(0,1)随机数u1,…,un, 则 （2）xi= 2 u 1 i
xi n n f ( x ) 1 2 e i ˆ （ 3） n i 1 g ( xi ) n i 1 1 xi
function I=exp_3 ( mm) %重要抽样法 %mm 是随机实验次数
如果 n 个元件的寿命分别为 X 1 , X 2 ,, X n 且 X i ~ E ( ), i 1,2,, n 求在以上 4 种组成方式下，系统 L 的寿命 X 的密度函数. 解
e x , x 0 f Xi (x ) 其它 0,
1 e x , x 0 FX i ( x ) 其它 0,
2 g ( x) (1 x) 3
2 (1 x), g ( x) 3 0,

b
0 x 1 else
a
f X f x g x dx E g x gX
设x1,…,xn是来自g(x)的随机数，则 的估计为
FX ( x) FX ( x)
i 1
i
n
(1 e ) , x 0, 0, x0
ne f X ( x)
x
x n
(1 e 0,
x n 1
) , x0 x0
(3)
X X1 X 2 X n
n = 2 时，

f X X ( x) f X (t ) f X ( x t )dt
P( S1 ) P(( A1 A2 ) ( B1B2 )) P( A1 A2 ) P( B1B2 ) P( A1 A2 B1 B2 )
例R4 S2:
A1 B1
2