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

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

2、投点法平均值法R软件运用R软件求=J】e-xdx,0给出精度与样本量的关系，0二J4e-xdx和91dx数值积分。计算以上各种估计的方差,201+x2比较各种方法的效率，如果您需要使用本文档，请点击下载按钮下载!如果您需要使用本文档，请点击下载按钮下载!1绪论蒙特卡洛的基本思想是，当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。蒙特卡洛方法解题过程的三个主要步骤：（1）构造或描述概率过程对于本身就具有随机性质的问题，如粒子输运问题，主要是正确描述和模拟

3、这个概率过程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。即要将不具有随机性质的问题转化为随机性质的问题。（2）实现从已知概率分布抽样构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量（或随机向量），就成为实现蒙特卡洛方法模拟实验的基本手段，这也是蒙特卡洛方法被称为随机抽样的原因。最简单、最基本、最重要的一个概率分布是（0，1）上的均匀分布（或称矩形分布）。随机数就是具有这种均匀分布的随机变量。随机数序列就是具有这种分布的总体的一个简单子样，也就是一个具有这种分布的

4、相互独立的随机变数序列。产生随机数的问题，就是从这个分布的抽样问题。在计算机上，可以用物理方法产生随机数，但价格昂贵，不能重复，使用不便。另一种方法是用数学递推公式产生。这样产生的序列，与真正的随机数序列不同，所以称为伪随机数，或伪随机数序列。不过，经过多种统计检验表明，它与真正的随机数，或随机数序列具有相近的性质，因此可把它作为真正的随机数来使用。由已知分布随机抽样有各种方法，与从（0,1）上均匀分布抽样不同，这些方法都是借助于随机序列来实现的，也就是说，都是以产生随机数为前提的。由此可见，随机数是我们实现蒙特卡洛模拟的基本工具。（3）建立各种估计量一般说来，构造了概率模型并能从中抽样后，即

5、实现模拟实验后，我们就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计。建立各种估计量，相当于对模拟实验的结果进行考察和登记，从中得到问题的解。2方法介绍2.1随机投点法随机投点法是进行n次试验，当n充分大的时候，以随机变量k/n作为期望值E(X)的近似估计值，即E(X)沁p=k/n其中k是n次实验中成功的次数。若一次投点试验的成功概率为P，并以v(1,表明试验成功Aio,表明试验失败则一次试验成功的均值与方差为E(X.)-1-p+0-(1-p)-piVar(X)-12-p+02(1-p)-p2-p(1-p)i若进行n次试验，其中k次试验成功，则k为具有参数为(n,p)的二项分布，

6、此时，随机变量k的估计为pk/n显然，随机变量的均值和方差满足E(p)-E-1E(k)-pVn丿nVar(p)-ddn设计算的定积分为I-f血fdx，其中a,b为有限数，被积函数f(x)是连续随机a变量g的概率密度函数，因此f(x)满足如下条件：f(x)非负，且Jf(x)dx=1-g显然I是一个概率积分，其积分值等于概率P(ab)。F面按给定分布f(x)随机投点的办法，给出如下MonteCarlo近似求积算法：产生服从给定分布的随机变量值，i=l,2,N;检查x是否落入积分区间。如果条件axb满足，则记录x落入积分区间iii一次。假设在N次实验以后，落入积分区间的总次数为n,那么用作为概率积分

7、的近似值，即I2.2平均值法任取一组相互独立、同分布的随机变量，g在a,b内服从分布率p(x)，令,ii则也是一组相互独立、同分布的随机变量，而且E百*()=fp*(x)p(x)dx=Jiaa由强大数定理plim兰p*(g)=I=1NT8N.i若记丄兰p吨)=I，则依概率1收敛到I,平均值法就是用门乍为I的近似值。Nii=1假如所需计算积分为I=Jf(xx，其中被积函数在a,b内可积，任意选择一个a有简单办法可以进行抽样的概率密度函数p(x)，使其满足条件：p(x)h0,当f(x)h0时(ax0,可取g(x)=cf(x),当c=1/I时就有Var(I)=0.般应选取和f(x)相似的密度函数g(

8、x)，使f(x)/g(x)接近于常数，故而Var(I)接近于0，以达到降低模拟实验的方差，这种减少方差的模拟试验法为重要抽样法。2.4分层抽样法分层抽样法是利用贡献率大小来降低估计方差的方法。它首先把样本空间D分成一些不交的小区间D=U，然后在各个小区间内的抽样数由其贡献大小决定。即，定义pf(xx，则D内的抽样数n应与p成正比。iiiiDi考虑积分0=/f(xx将0,1分成m个小区间：0则0=Jf(x=ff(x=瓦Ii0i=1ai=1i-1记l=a-a为第i个小区间的长度，i=1,2,.,m，在每个小区间上的积分值iii-1可用均值法估计出来，然后将其相加即可给出0的一个估计。具体步骤为：1

9、)独立产生U(0，1)随机数u,j=1,.ij2)计算x=a+1u,j=1,.iji-1iij3)计算i=丄iLf(x)nijij=1于是0的估计为0八=i,ii=1其方差为Var(0)=1空n其中，i-1I2ili2.5对偶变量法控制变量法利用数学上积分运算的线性特性：Jf(xI/x=Jf(x)-g(x)Lx+fg(x选择函数g(x)时要考虑到：g(x)在整个积分区间都是容易精确算出，并且在上式右边第一项的运算中对(f-g)积分的方差应当要比第二项对f积分的方差小。如果您需要使用本文档，请点击下载按钮下载!如果您需要使用本文档，请点击下载按钮下载!如果您需要使用本文档，请点击下载按钮下载!s

10、1(107)s1(107)i=1在应用这种方法时，在重要抽样法中所遇到的，当g(x)趋于零时，被积函数(f-g)趋于无穷大的困难就不再存在，因而计算出的结果稳定性比较好。该方法也不需要从分布密度函数g(x)，解析求出分布函数G(x)。由此我们可以看出选择g(x)所受到的限制比重要抽样法要小些。模拟过程：1)独立产生U(0,1)随机数u,j=1,.ij0=-工f(1-X)5nii=12)计算0二丄工f(x)5nii=13)计算0=.050=-ii2n2i=12.6控制变量法通常在蒙特卡洛计算中采用互相独立的随机点来进行计算。对偶变量法中却使用相关联的点来进行计算。它利用相关点间的关系可以是正关联

11、的，也可以是负关联的这个特点。两个函数值fi和f2之和的方差为Vf+f=Vf+Vf+2E(f-EfXf-Ef)12121122如果我们选择一些点，它们使f和f是负关联的。这样就可以使上12式所示的方差减小。当然这需要对具体的函数f和f有充分的了解121)独立产生U(0,1)随机数u,j=1,.ij2)计算1=n为f(xi)，找皿他)是相关的，且Eg(x)=卩i=13)计算0=0+九1(g(x)一卩)62ni3程序及实现结果9=i1e-xdx的求解0随机投点法先利用R软件产生服从0,1上均匀分布的随机数X,Y,，计算yf(x)的个数,即事件发生的频数，求出频率，即为积分的近似值。R程序s1-fu

12、nction(n)f-function(x)exp(-x)a-0b-1x-runif(n)y-runif(n)m-sum(yf(x)j=m/nvar-1/n*var(yf(x)lis-list(j,var)return(lis)s1(104)s1(105)s1(10飞)如果您需要使用本文档，请点击下载按钮下载!如果您需要使用本文档，请点击下载按钮下载!如果您需要使用本文档，请点击下载按钮下载!p1(104)对模拟次数n调试了4次，分别为n4,n5,n6,n7，得到精确值和模拟值。表3.1.1随机投点法的模拟次数和模拟值模拟值0.62690.632350.6325030.6319298方差2.3

13、2599e-052.32414e-062.32713e-072.32477e-08精确值为0.63212063.1.2平均值法先用R软件产生n个服从0,1上均匀分布的随机数x，计算f(x)，再ii计算f(x)的平均值，即为定积分的近似值iR程序p1-function(n)f-function(x)exp(-x)a-0b-1x-runif(n)y-mean(f(x)var-1/n*var(f(x)lis-list(y,var)return(lis)pl(105)pl(10飞)pl(lOS)对模拟次数n调试了4次，分别为n4,n5,n6,n7，得到精确值和模拟值。表3.1.2平均值法的模拟次数和模

14、拟值n104105106107模拟值0.63101170.63226430.63181990.632079方差3.25430e-063.27162e-073.2805e-083.27536e-09精确值为0.63212063.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)如果您需要使用本文档，请点击下载按钮下载!如果您需要使用本文档，请

15、点击下载按钮下载!如果您需要使用本文档，请点击下载按钮下载!psum(y(x)psum(y(x)f1(10,20)zl(104)zl(105)z1(10飞)zl(107)对模拟次数n调试了4次，分别为n4,n5,n6,n7，得到精确值和模拟值。表3.1.3重要性抽样法法的模拟次数和模拟值n104105106107模拟值0.61348190.61348190.61348190.6134819方差7.06957e-067.06957e-077.06957e-087.06957e-09精确值为0.63212063.1.4分层抽样法R程序f1-function(n,m)r1-runif(n,min-0

16、,max-0.5)r2-runif(m,min-0.5,maxT)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)fl(100,200)fl(1000,2000)fl(104,2*104)得到精确值和模拟值n=10,m=20n=100,n=1000n=10000为m=200M=2000M=20000值0.65827020.62917550.63417180.63270540.0002692072.75023e053.23184e063.1833

17、4e07表3.1.4分层抽样法的模拟次数和模拟值精确值0.63212063.1.5对偶变量法先用R软件产生服从0,1上均匀分布的随机数X，函数f(x).计算二1Yf(x)(二1Ef(1-x)TOC o 1-5 h z5ni5n1i=1，i=1sk+()ief(x)+f(i-x) HYPERLINK l bookmark52 V=54=乙ii计算2ni=12R程序d1-function(n)f-function(x)exp(-x)y-function(x)exp(1x)x-runif(n)msum(f(x)如果您需要使用本文档，请点击下载按钮下载!如果您需要使用本文档，请点击下载按钮下载!j-(

18、m/n+p/n)/2var-l/4*(var(f(x)+var(y(x)+2*cov(f(x),y(x)lis-list(j,var)return(lis)dl(104)dl(105)dl(10飞)dl(107)对模拟次数n调试了4次，分别为n4,n5,n6,n7，得到精确值和模拟值。表3.1.5对偶变量法的模拟次数和模拟值模拟值0.63219790.63206590.63209570.6321183方差0.000524720.0005327480.00053057410.0005291691精确值为0.63212063.1.6控制变量法R程序k1-function(n)f-function(

19、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)+1*mean(g(r)-u)var-l/n*var(f(r)+g(r)1st-list(j,var)return(1st)kl(104)kl(105)kl(10飞)kl(107)对模拟次数n调试了4次，分别为n4,n5,n6,n7，得到精确值和模拟值。表3.1.6控制变量法的模拟次数和模拟值104105106107模拟值0.6345750.63208840.6318490.6321651方差5.713907e-055.7

20、2338e-065.736774e-075.731522e-08精确值为0.63212063.2对0=J4e-xdx积分求解23.2.1随机投点法R程序s2function(n)ffunction(x)exp(x)t二function(y)(f(a+(b-a)*y)-c)/(d-c)a-2b-4c-f(4)d-f(2)s-(b-a)*(d-c)x-runif(n)y-runif(n)m-sum(yf(x)jm/ng=s*j+c*(b-a)var-1/n*var(yf(x)lis-list(g,var)return(lis)s2(104)s2(105)s2(10飞)s2(107)对模拟次数n调试

21、了4次，分别为n4,n5,n6,n7，得到精确值和模拟值。表3.2.1随机投点法的模拟次数和模拟值模拟值0.18688450.18688450.18688450.1868845方差2.33071e052.32230e062.32479e072.32611e08精确值为0.11701963.2.2平均值法R程序p2function(n)ffunction(r)exp(x)a2b4rrunif(n)h(ba)*f(a+(b-a)*r)ymean(h)var1/n*var(h)lislist(y,var)return(lis)p2(104)p2(105)p2(10飞)p2(107)对模拟次数n调试了

22、4次，分别为n4,n5,n6,n7，得到精确值和模拟值。表3.2.2平均值法的模拟次数和模拟值104105106107模拟值0.11733330.11695820.11700380.1170311方差1.30285e-051.30285e-061.30285e-071.30285e-08精确值为0.11701963.2.3分层抽样法R程序f2-function(n,m)r1-runif(n,min-0,max-0.5)r2-runif(m,min-0.5,maxT)c-1/2*mean(2*exp(-2-2*r1)+1/2*mean(2*exp(-2-2*r2)var-var(2*exp(-2

23、-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(104,2*104)对模拟次数n调试了4次，分别为n4,n5,n6,n7，得到精确值和模拟值。表3.2.3分层抽样法的模拟次数和模拟值n=10,m=20n=100,n=1000n=10000确值m=200M=2000M=20000值0.11081320.12196680.11737690.11697836.12767e-056.31761e-066.20002e-076.02697e-08精为0.11701

24、963.2.4对偶变量法R程序d2-function(n)a-2b-4f-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)/2j-s*j3+c*(b-a)return(j)d2(104)d2(105)如果您需要使用本文档，请点击下载按钮下载!如果您需要使用本文档，请点击下载按钮下载!c=f(4)d2(10飞)d2(107)对模拟次数n调试了4次，分别为n4,n5,n6,n7，得到精确值和模拟值。表

25、3.2.4对偶变量法的模拟次数和模拟值n104105106107模拟值0.11690430.1167410.1169430.1170233精确值为0.11701963.2.5控制变量法R程序k2-function(n)f-function(x)exp(-x)r-runif(n)a-2b-4c-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)+1*mean(g(r)-u)j-s*p+c*(b

26、-a)var-l/n*var(f(r)+g(r)1st-list(j,var)return(1st)k2(104)k2(105)k2(10飞)k2(107)对模拟次数n调试了4次，分别为n4,n5,n6,n7，得到精确值和模拟值。表3.2.5控制变量法的模拟次数和模拟值模拟值0.11747130.11711510.117031610.1170211方差5.68197e-055.73480e-065.73725e-075.73387e-08精确值为0.11701963.2.6重要性采样法R程序z2-function(n)a=2b=4f=function(x)exp(-x)d=f(2)s=(b-a

27、)*(d-c);r二runif(n)x-1-(sqrt(l-r)p-function(x)1/(d-c)*(f(a+(b-a)*x)-c)q-function(x)(2*(1-x)jl-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(104)z2(105)z2(10飞)z2(107)对模拟次数n调试了4次，分别为n4,n5,n6,n7，得到精确值和模拟值。表3.2.6重要性采样法的模拟次数和模拟值如果您需要使用本文档，请点击下载按钮下载!如果您需要使用本文档，请点击下载按钮下载!如果您需

28、要使用本文档，请点击下载按钮下载!p3(l04)p3(l04)s3(104)104105106107模拟值0.11735220.11702670.11701140.1170123方差8.17541e-078.16598e-088.17952e-098.17912e-10精确值为0.11701963.31二dx积分求解01+X23.3.1随机投点法R程序s3-function(n)f-function(x)exp(-x)/(1+x“2)a-0b-1x-runif(n)y-runif(n)m-sum(yf(x)j=m/nvar-1/n*var(yf(x)lis-list(j,var)return(

29、lis)s3(105)s3(10飞)s3(107)对模拟次数n调试了4次，分别为n4,n5,n6,n7，得到精确值和模拟值。表3.3.1随机投点法的模拟次数和模拟值104105106107模拟值0.53080.525070.5253510.5247777方差2.49625e-052.49312e-062.49423e-072.49371e-08精确值为0.52479713.3.2平均值法R程序p3-function(n)f-function(x)exp(-x)/(l+x“2)a-0b-lx-runif(n)y-mean(f(x)var-l/n*var(f(x)lis-list(y,var)re

30、turn(lis)如果您需要使用本文档，请点击下载按钮下载!如果您需要使用本文档，请点击下载按钮下载!d3(105)p3(105)p3(10飞)p3(107)对模拟次数n调试了4次，分别为n4,n5,n6,n7，得到精确值和模拟值。表3.3.2平均值法的模拟次数和模拟值n104105106107模拟值0.52397070.52526550.52457150.5247576方差5.97263e-066.02206e-075.99868e-085.99926e-09精确值为0.52479713.3.3分层抽样法R程序f3-function(n,m)r1-runif(n,min-0,max-0.5)

31、r2-runif(m,min-0.5,maxT)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(104,2*10八4)得到精确值和模拟值n=10,m=20n=100,n=1000n=10000为m=200M=2000M=20000值0.53385040.53059470.52288030.52543610.0003218132.1

32、1722e-052.16594e-062.17694e-07表3.3.3分层抽样法的模拟次数和模拟值精确值0.52479713.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)/2var-1/4*(var(f(r)+var(f(1r)+2*cov(f(r),f(1-r)lis-list(j,var)return(lis)d3(104)d3(10飞)d3(107)如果您需要使用本文档，请点击下载按钮下载!如果您需要使用本文档，请点击下载按钮下载!k1(10飞

33、)k1(10飞)对模拟次数n调试了4次，分别为n4,n5,n6,n7，得到精确值和模拟值。表3.3.4对偶法的模拟次数和模拟值104105106107模拟值0.63197880.63214440.63212360.6321169方差0.0005191260.00053013340.00052976580.0005295388精确值为0.52479713.3.5控制变量法R程序k3function(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)#定义

34、lambdaj-mean(f(r)+1*mean(g(r)-u)var-1/n*var(f(r)+g(r)1st-list(j,var)return(1st)k1(104)k1(105)k1(107)如果您需要使用本文档，请点击下载按钮下载!如果您需要使用本文档，请点击下载按钮下载!如果您需要使用本文档，请点击下载按钮下载!z3(105)对模拟次数n调试了4次，分别为n4,n5,n6,n7，得到精确值和模拟值表3.3.5控制变量法的模拟次数和模拟值n104105106107模拟值0.62974030.6324380.63209560.6321038方差5.80168e-055.74223e-065.73553e-075.73648e-08精确值为0.52479713.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