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

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

《基于蒙特卡洛方法求数值积分与.docx》由会员分享，可在线阅读，更多相关《基于蒙特卡洛方法求数值积分与.docx（55页珍藏版）》请在冰点文库上搜索。

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

统计计算课程设计

题目

基于蒙特卡洛方法求数值积分

中文摘要

蒙特卡洛方法，又称随机抽样或统计试验方法，属于计算数学的一个分支，它是在上世纪四十年代中期为了适应当时原子能事业的发展而发展起来的。

传统的经验方法由于不能逼近真实的物理过程，很难得到满意的结果，而蒙特卡罗方法由于能够真实地模拟实际物理过程，故解决问题与实际非常符合，可以得到很圆满的结果。

利用随机投点法，平均值法，重要性采样法，分层抽样法，控制变量法，对偶变量法，运用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）随机投点的办法，给出如下MonteCarlo近似求积算法：

（1）产生服从给定分布的随机变量值,i=1,2,…,N;

（2）检查是否落入积分区间。

如果条件满足，则记录落入积分区间一次。

假设在N次实验以后，落入积分区间的总次数为n，那么用

作为概率积分的近似值，即

2.2平均值法

任取一组相互独立、同分布的随机变量，在[a,b]内服从分布率p（x）,令,则也是一组相互独立、同分布的随机变量，而且

由强大数定理

若记,则依概率1收敛到I，平均值法就是用作为I的近似值。

假如所需计算积分为，其中被积函数在[a,b]内可积，任意选择一个有简单办法可以进行抽样的概率密度函数p（x），使其满足条件：

●

●

记

则所求积分为

因而MonteCarlo近似求积算法为：

（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

j=m/n

var<-1/n*var（y

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

j

g=s*j+c*（b-a）

var<-1/n*var（y

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

j=m/n

var<-1/n*var（y

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

模