基于蒙特卡洛方法求数值积分与R.doc

统计计算课程设计
题 目
基于蒙特卡洛方法求数值积分
中文摘要
蒙特卡洛方法，又称随机抽样或统计试验方法，属于计算数学的一个分支，它是在上世纪四十年代中期为了适应当时原子能事业的发展而发展起来的.传统的经验方法由于不能逼近真实的物理过程，很难得到满意的结果,而蒙特卡罗方法由于能够真实地模拟实际物理过程，故解决问题与实际非常符合,可以得到很圆满的结果。
利用随机投点法,平均值法，重要性采样法,分层抽样法,控制变量法，对偶变量法，运用R软件求,和数值积分.计算以上各种估计的方差，给出精度与样本量的关系，比较各种方法的效率,
关键字 蒙特卡洛 随机投点法 平均值法 R软件

1 绪论
蒙特卡洛的基本思想是,当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。
蒙特卡洛方法解题过程的三个主要步骤：
(1）构造或描述概率过程
对于本身就具有随机性质的问题,如粒子输运问题，主要是正确描述和模拟这个概率过 程，对于本来不是随机性质的确定性问题，比如计算定积分,就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解.即要将不具有随机性质的问题转化为随机性质的问题。
（2）实现从已知概率分布抽样
构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量(或随机向量）,就成为实现蒙特卡洛方法模拟实验的基本手段,这也是蒙特卡洛方法被称为随机抽样的原因。最简单、最基本、最重要的一个概率分布是(0，1)上的均匀分布(或称矩形分布）。随机数就是具有这种均匀分布的随机变量。随机数序列就是具有这种分布的总体的一个简单子样，也就是一个具有这种分布的相互独立的随机变数序列。产生随机数的问题,就是从这个分布的抽样问题。在计算机上,可以用物理方法产生随机数,但价格昂贵，不能重复,使用不便.另一种方法是用数学递推公式产生.这样产生的序列，与真正的随机数序列不同，所以称为伪随机数，或伪随机数序列。不过，经过多种统计检验表明，它与真正的随机数，或随机数序列具有相近的性质,因此可把它作为真正的随机数来使用。由已知分布随机抽样有各种方法，与从(0，1）上均匀分布抽样不同，这些方法都是借助于随机序列来实现的,也就是说,都是以产生随机数为前提的.由此可见，随机数是我们实现蒙特卡洛模拟的基本工具。
（3）建立各种估计量
一般说来，构造了概率模型并能从中抽样后，即实现模拟实验后，我们就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计.建立各种估计量，相当于对模拟实验的结果进行考察和登记,从中得到问题的解。
2方法介绍
2。1随机投点法
随机投点法是进行n次试验，当n充分大的时候，以随机变量k/n作为期望值E（X)的近似估计值，即
其中k是n次实验中成功的次数.
若一次投点试验的成功概率为p,并以
则一次试验成功的均值与方差为
若进行n次试验，其中k次试验成功，则k为具有参数为（n,p）的二项分布，此时，随机变量k的估计为
显然，随机变量的均值和方差满足
dd
设计算的定积分为，其中a，b为有限数，被积函数f（x)是连续随机变量的概率密度函数，因此f（x)满足如下条件：
显然I是一个概率积分,其积分值等于概率.
下面按给定分布f（x)随机投点的办法,给出如下Monte Carlo近似求积算法：
（1)产生服从给定分布的随机变量值，i=1，2,…，N；
（2)检查是否落入积分区间。如果条件满足,则记录落入积分区间一次。
假设在N次实验以后，落入积分区间的总次数为n，那么用
作为概率积分的近似值，即
2。2平均值法
任取一组相互独立、同分布的随机变量，在[a，b］内服从分布率p（x），令，则也是一组相互独立、同分布的随机变量，而且
由强大数定理
若记 ，则依概率1收敛到I，平均值法就是用作为I的近似值。
假如所需计算积分为，其中被积函数在［a，b］内可积,任意选择一个有简单办法可以进行抽样的概率密度函数p(x）,使其满足条件：
记 则所求积分为
因而Monte Carlo近似求积算法为:
（1) 产生服从分布率p（x)的随机数
（2） 计算均值，即有
2。3重要性采样法
从数学角度上看定积分可以看成
其中g（x）是某个随机变量X的密度函数,因此积分值I可看成随机变量Z=f（x)/g（x)的数学期望值
为了减少模拟实验的方差应适当选取g（x)，使Var(I)尽可能小,如果被积函数f（x）>0,可取g（x）=cf（x)，当c=1/I时就有Var(I)=0。一般应选取和f（x）相似的密度函数g（x)，使f(x）/g(x）接近于常数，故而Var（I）接近于0,以达到降低模拟实验的方差,这种减少方差的模拟试验法为重要抽样法。
2。4分层抽样法
分层抽样法是利用贡献率大小来降低估计方差的方法。它首先把样本空间D分成一些不交的小区间，然后在各个小区间内的抽样数由其贡献大小决定。即，定义，则内的抽样数应与成正比。
考虑积分将[0,1]分成m个小区间:
则
记为第i个小区间的长度，i=1，2,。。.，m,在每个小区间上的积分值可用均值法估计出来，然后将其相加即可给出的一个估计。具体步骤为：
1） 独立产生U（0，1)随机数
2） 计算
3） 计算
于是的估计为,其方差为
其中,
2。5对偶变量法
控制变量法利用数学上积分运算的线性特性:
选择函数g（x）时要考虑到：g（x)在整个积分区间都是容易精确算出， 并且在上式右边第一项的运算中对（f—g)积分的方差应当要比第二项对f积分的方差小。
在应用这种方法时，在重要抽样法中所遇到的,当g（x)趋于零时，被积函数(f-g）趋于无穷大的困难就不再存在，因而计算出的结果稳定性比较好。 该方法也不需要从分布密度函数g（x)， 解析求出分布函数G(x)。 由此我们可以看出选择g(x)所受到的限制比重要抽样法要小些.
模拟过程：
1) 独立产生U(0，1)随机数
2) 计算
3) 计算
2.6 控制变量法
通常在蒙特卡洛计算中采用互相独立的随机点来进行计算。 对偶变量法中却使用相关联的点来进行计算.它利用相关点间的关系可以是正关联的，也可以是负关联的这个特点.
两个函数值和之和的方差为
如果我们选择一些点,它们使和是负关联的.这样就可以使上
式所示的方差减小。当然这需要对具体的函数和有充分的了解
1） 独立产生U（0，1）随机数
2） 计算，找g（x）,f(x）是相关的，且E［g（x)］=
3） 计算
3 程序及实现结果
3。1 的求解
3。1。1 随机投点法
先利用R 软件产生服从［0，1］上均匀分布的随机数 X，Y， ，计算的个数,即事件发生的频数，求出频率,即为积分的近似值。
R程序
s1<-function（n）
{
f〈-function(x) exp(-x)
a〈—0
b〈-1
x〈-runif（n）
y〈-runif(n)
m〈—sum（y〈f(x））
j=m/n
var〈—1/n*var(y〈f（x）)
lis〈-list（j,var）
return(lis）
｝
s1（10^4）
s1（10^5)
s1（10^6)
s1（10^7)
对模拟次数n调试了4次，分别为，得到精确值和模拟值。
表3。1.1 随机投点法的模拟次数和模拟值
n

模拟值
0。6269
0.63235
0。632503
0。6319298
方差
2。32599e—05
2.32414e—06
2.32713e-07
2.32477e—08
精确值为0.6321206
3。1。2 平均值法
先用R 软件产生n个服从［0,1］上均匀分布的随机数,计算，再计算的平均值，即为定积分的近似值
R程序
p1〈—function（n）
｛
f〈—function（x） exp（—x)
a〈—0
b〈—1
x<—runif(n）
y〈—mean(f（x)）
var〈-1/n*var（f（x））
lis<—list（y,var）
return(lis）
}
p1（10^4）
p1(10^5)
p1（10^6）
p1（10^7)
对模拟次数n调试了4次,分别为,得到精确值和模拟值。
表3。1.2 平均值法的模拟次数和模拟值
n

模拟值
0。6310117
0.6322643
0。6318199
0.632079
方差
3.25430e—06
3.27162e—07
3。2805e—08
3。27536e-09
精确值为0。6321206
3。1。3 重要性抽样法
R程序
z1〈—function(n）
｛
x<- 1-(sqrt(1—r)）
f〈-function(x） exp(-x)
g〈—function（x） （2*（1—x))
r〈-runif（n）
s=mean(f（x）/g（x））
var<-1/n＊var（f（x)/g(x））
lis<-list（s,var)
return(lis）
}
z1(10^4)
z1(10^5）
z1(10^6）
z1（10^7)
对模拟次数n调试了4次，分别为，得到精确值和模拟值.
表3。1。3 重要性抽样法法的模拟次数和模拟值
n

模拟值
0。6134819
0.6134819
0。6134819
0。6134819
方差
7。06957e-06
7。06957e—07
7.06957e-08
7。06957e—09
精确值为0。6321206
3。1.4 分层抽样法
R程序
f1〈—function(n,m)
{
r1〈—runif(n,min<-0，max<—0。5）
r2〈-runif(m，min〈-0。5，max<—1）
c<—1/2*mean（exp（-r1））+1/2*mean（exp（—r2))
var<—var(exp(—r1)）/（4＊n)+var（exp（-r2））/(4＊m）
j〈—list(c，var）
return（j)
}
f1(10，20）
f1（100,200）
f1(1000，2000）
f1（10^4,2＊10^4）
得到精确值和模拟值
表3.1.4 分层抽样法的模拟次数和模拟值
取值
n=10,m=20
n=100，
m=200
n=1000
M=2000
n=10000
M=20000
模拟值
0。6582702
0.6291755
0。6341718
0.6327054
方差
0.000269207
2。75023e-05
3。23184e—06
3。18334e-07
精确值为0。6321206
3。1。5 对偶变量法
先用 R 软件产生服从［0，1］上均匀分布的随机数 X ，函数f(x),
计算，
计算
R程序
d1〈—function（n）
｛
f〈—function(x） exp（-x)
y<—function（x） exp（—(1-x)）
x<-runif(n)
m〈—sum(f（x）)
p<-sum(y（x))
j〈—(m/n+p/n）/2
var〈—1/4*（var（f（x)）+var(y（x））+2＊cov(f(x），y（x）))
lis<-list（j，var）
return（lis)
}
d1(10^4）
d1(10^5）
d1(10^6)
d1（10^7)
对模拟次数n调试了4次，分别为,得到精确值和模拟值.
表3。1。5 对偶变量法的模拟次数和模拟值
n

模拟值
0。6321979
0.6320659
0。6320957
0.6321183
方差
0。00052472
0。000532748
0。0005305741
0。0005291691
精确值为0。6321206
3。1。6 控制变量法
R程序
k1<—function（n）
｛
f〈—function(x） exp（—x)
r〈—runif(n）
g〈-function（x） 2＊（1—x)
u〈-mean(g（r）)
l<—(—cov(f（r），g（r）））/var（g(r))
j〈-mean（f(r））+l＊mean（g(r）—u)
var〈—1/n＊var（f（r)+g（r)）
lst<—list(j,var）
return（lst）
｝
k1(10^4）
k1（10^5）
k1（10^6)
k1(10^7）
对模拟次数n调试了4次，分别为，得到精确值和模拟值.
表3。1.6 控制变量法的模拟次数和模拟值
n

模拟值
0。634575
0.6320884
0。631849
0.6321651
方差
5。713907e-05
5.72338e-06
5。736774e-07
5。731522e-08
精确值为0。6321206
3。2 对积分求解
3.2。1 随机投点法
R程序
s2<—function（n）
｛
f<-function（x) exp（—x）
t=function（y） （f（a+（b—a)＊y)—c）/（d—c）
a<-2
b〈—4
c<—f（4）
d〈—f(2)
s〈—(b—a)＊（d—c）
x<—runif(n)
y〈—runif（n）
m<—sum（y<f(x）)
j<m/n
g=s*j+c*(b-a）
var〈-1/n＊var（y〈f（x））
lis〈—list（g,var）
return(lis）
｝
s2(10^4)
s2(10^5）
s2（10^6）
s2（10^7）
对模拟次数n调试了4次,分别为，得到精确值和模拟值.
表3。2。1 随机投点法的模拟次数和模拟值
n

模拟值
0.1868845
0.1868845
0.1868845
0。1868845
方差
2。33071e—05
2.32230e—06
2。32479e—07
2.32611e—08
精确值为 0.1170196
3.2。2 平均值法
R程序
p2<—function(n）
｛
f〈—function（r） exp（—x）
a〈—2
b<—4
r〈—runif（n）
h〈-(b—a）＊f（a+(b-a）*r）
y<—mean（h）
var〈—1/n＊var（h)
lis〈—list（y，var）
return（lis）
｝
p2（10^4）
p2（10^5）
p2（10^6)
p2(10^7）
对模拟次数n调试了4次,分别为,得到精确值和模拟值。
表3.2。2 平均值法的模拟次数和模拟值
n

模拟值
0.1173333
0。1169582
0.1170038
0.1170311
方差
1.30285e—05
1。30285e—06
1。30285e—07
1.30285e—08
精确值为 0。1170196
3.2。3 分层抽样法
R程序
f2〈—function(n，m)｛
r1〈—runif（n,min<—0,max<—0。5）
r2〈—runif（m，min<—0。5，max〈—1）
c〈—1/2*mean（2＊exp（-2—2*r1））+1/2＊mean(2*exp（-2—2＊r2)）
var〈—var（2＊exp（-2—2＊r1））/(4*n)+var(2*exp（-2—2*r2）/（4*m））
j〈—list（c,var）
return(j)
｝
f2（10，20)
f2（100，200）
f2(1000，2000）
f2（10^4，2＊10^4）
对模拟次数n调试了4次,分别为，得到精确值和模拟值。
表3.2。3 分层抽样法的模拟次数和模拟值
取值
n=10,m=20
n=100,
m=200
n=1000
M=2000
n=10000
M=20000
模拟值
0。1108132
0.1219668
0。1173769
0.1169783
方差
6.12767e—05
6。31761e—06
6.20002e—07
6。02697e—08
精确值为 0.1170196
3.2。4 对偶变量法
R程序
d2〈-function(n）
{
a〈-2
b〈—4
f〈—function(x） exp(—x）
r〈—runif(n)
c<-f(b)
d<-f（c）
s<—（b—a）＊（d-c）
p〈-function(u) 1/(d—c）＊（f(a+(b—a）＊u）-c)
j1〈-mean(p（r）)
j2〈—mean(p（r)）
j3<—（j1+j2）/2
j〈—s*j3+c＊（b—a)
return（j)
}
d2(10^4）
d2(10^5)
d2（10^6）
d2（10^7）
对模拟次数n调试了4次，分别为，得到精确值和模拟值.
表3.2。4 对偶变量法的模拟次数和模拟值
n

模拟值
0.1169043
0。116741
0.116943
0.1170233
精确值为 0。1170196
3。2。5 控制变量法
R程序
k2<—function（n)
｛
f〈-function（x） exp(—x）
r〈—runif(n）
a〈—2
b〈—4
c〈—f（4)
d<—f(2）
s〈-（b-a）＊（d-c)
q〈-function（x) 1/（d—c）＊(f(a+(b-a）＊x)—c）
g〈—function（x） 2＊（1-x)
u〈—mean（g（r）)
l〈-（-cov(f(r)，g(r）））/var（g(r））
p〈—mean(q（r）)+l＊mean（g（r）—u）
j〈-s＊p+c＊（b—a）
var<—1/n＊var（f(r)+g（r））
lst〈—list(j,var)
return（lst）
｝
k2（10^4）
k2（10^5)
k2（10^6)
k2（10^7）
对模拟次数n调试了4次，分别为,得到精确值和模拟值。
表3。2。5 控制变量法的模拟次数和模拟值
n

模拟值
0。1174713
0.1171151
0。11703161
0.1170211
方差
5。68197e—05
5。73480e—06
5.73725e—07
5.73387e—08
精确值为 0。1170196
3。2。6 重要性采样法
R程序
z2<—function（n）
｛
a=2
b=4
f=function(x) exp(-x)
c=f（4）
d=f(2)
s=（b—a）*（d—c);
r=runif(n）
x〈- 1—（sqrt（1—r))
p<-function（x）
1/(d-c)＊（f(a+（b—a)＊x）—c)
q〈-function（x） (2*(1—x）)
j1〈—mean(p(x）/q（x）)
j=s＊j1+c＊(b—a）
var<—1/n＊var（p（x）/q(x）)
lis〈—list(j，var)
return（lis）
｝
z2(10^4)
z2(10^5)
z2(10^6)
z2（10^7)
对模拟次数n调试了4次，分别为，得到精确值和模拟值.
表3。2。6 重要性采样法的模拟次数和模拟值
n

模拟值
0.1173522
0.1170267
0。1170114
0.1170123
方差
8。17541e-07
8.16598e—08
8.17952e—09
8.17912e—10
精确值为 0。1170196
3.3 积分求解
3。3。1 随机投点法
R程序
s3<—function(n）
{
f<—function（x) exp(-x)/（1+x^2）
a<—0
b<—1
x〈—runif(n）
y<—runif（n)
m<—sum（y<f(x））
j=m/n
var〈—1/n＊var（y<f（x））
lis〈—list(j,var)
return(lis)
｝
s3（10^4）
s3(10^5)
s3（10^6）
s3(10^7)
对模拟次数n调试了4次,分别为，得到精确值和模拟值.
表3。3.1 随机投点法的模拟次数和模拟值
n

模拟值
0.5308
0.52507
0。525351
0.5247777
方差
2。49625e-05
2.49312e—06
2。49423e—07
2。49371e-08
精确值为 0。5247971
3。3.2平均值法
R程序
p3〈—function（n）
｛
f〈-function（x) exp(-x）/（1+x^2）
a〈-0
b〈—1
x〈—runif（n）
y<—mean（f（x）)
var〈—1/n＊var（f(x))
lis〈—list(y，var）
return（lis）
｝
p3(10^4)
p3（10^5）
p3(10^6)
p3（10^7）
对模拟次数n调试了4次，分别为，得到精确值和模拟值.
表3.3。2 平均值法的模拟次数和模拟值
n

模拟值
0。5239707
0.5252655
0。5245715
0.5247576
方差
5。97263e-06
6。02206e-07
5。99868e—08
5。99926e—09
精确值为 0。5247971
3.3.3分层抽样法
R程序
f3〈-function（n，m）
｛
r1〈-runif（n,min〈-0，max〈-0。5)
r2〈-runif（m,min<—0。5，max<—1）
z〈—function（u） exp(-u)/（1+u^2)
j〈—1/2*mean(z（r1））+1/2＊mean（z（r2))
var<-var(z(-r1））/（4*n）+var(z（-r2））/（4＊m)
lis〈—list（j，var)
return(lis)
｝
f3（10，20）
f3（100,200)
f3（1000，2000)
f3（10^4，2＊10^4)
得到精确值和模拟值
表3。3。3 分层抽样法的模拟次数和模拟值
取值
n=10,m=20
n=100,
m=200
n=1000
M=2000
n=10000
M=20000
模拟值
0。5338504
0。5305947
0.5228803
0。5254361
方差
0。000321813
2。11722e—05
2。16594e-06
2。17694e-07
精确值为0。5247971
3.3。4 对偶变量法
R程序
d3〈-function（n）
｛
f<—function(x） exp（—x）/(1+x^2）
r〈—runif（n)
m<—mean(f(r))
p〈-mean（f（1—r））
j<—(m+p）/2
var〈—1/4*(var（f(r))+var（f(1-r））+2＊cov(f（r），f（1-r)））
lis〈-list(j，var）
return(lis)
}
d3（10^4）
d3（10^5）
d3(10^6)
d3（10^7）
对模拟次数n调试了4次，分别为,得到精确值和模拟值。
表3。3.4 对偶法的模拟次数和模拟值
n

模拟值
0.6319788
0。6321444
0。6321236
0。6321169
方差
0。000519126
0.0005301334
0。0005297658
0.0005295388
精确值为 0。5247971
3.3.5 控制变量法
R程序
k3<—function（n）
｛
f〈—function（x) exp（—x）/（1+x^2)
r〈-runif(n）
g<—function（x） 2*（1—x）
u<-mean(g（r))
l〈-（—cov（f（r）,g（r））)/var（g（r）） ＃定义lambda
j〈—mean(f(r))+l＊mean（g(r）—u）
var<—1/n*var（f(r)+g(r））
lst〈—list(j，var)
return（lst）
｝
k1（10^4)
k1（10^5）
k1（10^6)
k1(10^7)
对模拟次数n调试了4次，分别为,得到精确值和模拟值
表3.3。5 控制变量法的模拟次数和模拟值
n

模拟值
0。6297403
0.632438
0.6320956
0.6321038
方差
5。80168e-05
5。74223e-06
5。73553e-07
5.73648e—08
精确值为 0.5247971
3。3。6 重要性抽样法
R程序
z3〈—function（n）
{
x<— 1—(sqrt（1—r））
f<-function(x） exp（-x）/(1+x^2）
g<—function（x） (2＊（1-x）)
r<-runif（n）
s=mean（f（x）/g(x)）
var<-1/n＊var（f（x）/g(x）)
lis〈-list(s，var)
return（lis)
}
z3(10^4）
z3(10^5)
z3（10^6）
z3(10^7)
对模拟次数n调试了4次，分别为,得到精确值和模拟值。
表3.3.6 重要性采样法的模拟次数和模拟值
n

模拟值
0.515594
0。515594
0。515594
0。515594
方差
1。18804e-06
1。18804e—07
1。18804e-08
1。18804e-09
精确值为 0。5247971
4方法结果比较
由输出结果可以看出，随机投点法，平均值法，对偶变量法，重要性采集法,分层抽样法和控制变量法中由不同模拟次数得出的积分估计值和真实值比较可已看出,他们相差不大，接近真实值，但对偶变量法和控制变量法的模拟值有些许差别.如下表：
4。1 积分的精确值和估计值

真实值
随机
平均
重要
分层
对偶
控制

0.6321206
0.6269
0。631011
0。613481
0。632705
0。632197
0。634575

0。1170196
0。1868845
0.117333
0.11735
0。116978
0。117019
0。11747

0。5247971
0。5308
0。523970
0.515594
0.525436
0。631978
0.629740
由方差可得,随机投点法的波动最小，对偶抽样法的波动最大，而且对偶抽样法的模拟值和精确值差距最大，因此我们可以看出对偶变量法的模拟效果最差。如下表
4。2 各种随机模拟的方差
N=10^4
随机
平均
重要
分层
对偶
控制

2.33E-05
3.25E-06
7.07E-06
3。18E—07
5。24E—04
5.71E—05

2.33E—05
1。31E—05
8.17E—07
2。17E—07
5,13E-04
5。68E-05

2。50E—05
5.97E-06
1。19E-06
2。18E—07
5，19E-04
5。80E-05
5总结
本文主要运用蒙特卡洛的随机投点法、平均值法、重要性抽样法、分层抽样法、对偶变量法5和控制变量法6种方法求数值积分，利用R 软件编程，进行模拟试验。最后通过进行 6 种方法的比较,分析积分的模拟值和真实值，同时比较方差大小，得出结论。
由结果可以看出随机投点法,分层抽样法，重要性采样法和平均值法的模拟值和真实值最为接近，且方差较小,模拟效果好。对偶变量法和控制变量法的模拟值和真实值较为接近，但有差距，并且这两种方法的方差较大，因此模拟效果不好。
参考文献
[1］ 高惠璇 统计计算［M] 北京大学出版社,1995
［2] 薛毅 陈丽平 统计建模与R软件［M］ 清华大学出版社 2009
［3］ 肖枝洪 朱强 统计模拟及其R实现［M] 武汉大学出版社 2010
［4] 汤银才 R语言与统计分析［M］ 高等教育出版社 2008
25