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

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

2、，当所求解问题是*种随机事件出现的概率，或者是*个随机变量的期望值时，通过*种实验的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的*些数字特征，并将其作为问题的解。蒙特卡洛方法解题过程的三个主要步骤：1构造或描述概率过程对于本身就具有随机性质的问题，如粒子输运问题，主要是正确描述和模拟这个概率过程，对于本来不是随机性质确实定性问题，比方计算定积分，就必须事先构造一个人为的概率过程，它的*些参量正好是所要求问题的解。即要将不具有随机性质的问题转化为随机性质的问题。2实现从概率分布抽样构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生概率

3、分布的随机变量或随机向量，就成为实现蒙特卡洛方法模拟实验的根本手段，这也是蒙特卡洛方法被称为随机抽样的原因。最简单、最根本、最重要的一个概率分布是0，1上的均匀分布或称矩形分布。随机数就是具有这种均匀分布的随机变量。随机数序列就是具有这种分布的总体的一个简单子样，也就是一个具有这种分布的相互独立的随机变数序列。产生随机数的问题，就是从这个分布的抽样问题。在计算机上，可以用物理方法产生随机数，但价格昂贵，不能重复，使用不便。另一种方法是用数学递推公式产生。这样产生的序列，与真正的随机数序列不同，所以称为伪随机数，或伪随机数序列。不过，经过多种统计检验说明，它与真正的随机数，或随机数序列具有相近的

4、性质，因此可把它作为真正的随机数来使用。由分布随机抽样有各种方法，与从(0,1)上均匀分布抽样不同，这些方法都是借助于随机序列来实现的，也就是说，都是以产生随机数为前提的。由此可见，随机数是我们实现蒙特卡洛模拟的根本工具。3建立各种估计量一般说来，构造了概率模型并能从中抽样后，即实现模拟实验后，我们就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计。建立各种估计量，相当于对模拟实验的结果进展考察和登记，从中得到问题的解。2方法介绍2.1随机投点法随机投点法是进展n次试验，当n充分大的时候，以随机变量k/n作为期望值E(*)的近似估计值，即其中k是n次实验中成功的次数。假设一次投点试

5、验的成功概率为p，并以则一次试验成功的均值与方差为假设进展n次试验，其中k次试验成功，则k为具有参数为n,p的二项分布，此时，随机变量k的估计为显然，随机变量的均值和方差满足dd设计算的定积分为，其中a,b为有限数，被积函数f(*)是连续随机变量的概率密度函数，因此f(*)满足如下条件:显然I是一个概率积分，其积分值等于概率。下面按给定分布f(*)随机投点的方法，给出如下Monte Carlo近似求积算法：(1)产生服从给定分布的随机变量值,i=1,2,N;(2)检查是否落入积分区间。如果条件满足，则记录落入积分区间一次。假设在N次实验以后，落入积分区间的总次数为n，则用作为概率积分的近似值，

6、即2.2平均值法任取一组相互独立、同分布的随机变量，在a,b服从分布率p(*),令,则也是一组相互独立、同分布的随机变量，而且由强大数定理假设记,则依概率1收敛到I，平均值法就是用作为I的近似值。假设所需计算积分为，其中被积函数在a,b可积，任意选择一个有简单方法可以进展抽样的概率密度函数p(*)，使其满足条件：记则所求积分为因而Monte Carlo近似求积算法为：(1) 产生服从分布率p(*)的随机数(2) 计算均值,即有2.3重要性采样法从数学角度上看定积分可以看成其中g(*)是*个随机变量*的密度函数，因此积分值I可看成随机变量Z=f(*)/g(*)的数学期望值为了减少模拟实验的方差应

7、适中选取g(*),使VarI尽可能小，如果被积函数f(*)0,可取g(*)=cf(*),当c=1/I时就有Var(I)=0.一般应选取和f(*)相似的密度函数g(*)，使f(*)/g(*)接近于常数，故而Var(I)接近于0，以到达降低模拟实验的方差，这种减少方差的模拟试验法为重要抽样法。2.4分层抽样法分层抽样法是利用奉献率大小来降低估计方差的方法。它首先把样本空间D分成一些不交的小区间，然后在各个小区间的抽样数由其奉献大小决定。即，定义，则的抽样数应与成正比。考虑积分将0,1分成m个小区间：则记为第i个小区间的长度，i=1,2,.,m，在每个小区间上的积分值可用均值法估计出来，然后将其相加

8、即可给出的一个估计。具体步骤为：独立产生U(0，1)随机数计算计算于是的估计为，其方差为其中，2.5对偶变量法控制变量法利用数学上积分运算的线性特性:选择函数g(*)时要考虑到：g(*)在整个积分区间都是容易准确算出,并且在上式右边第一项的运算中对(f-g)积分的方差应当要比第二项对f积分的方差小。在应用这种方法时，在重要抽样法中所遇到的，当g(*)趋于零时，被积函数(f-g)趋于无穷大的困难就不再存在，因而计算出的结果稳定性比拟好。该方法也不需要从分布密度函数g(*)，解析求出分布函数G(*)。由此我们可以看出选择g(*)所受到的限制比重要抽样法要小些。模拟过程：独立产生U(0，1)随机数计

9、算计算2.6 控制变量法通常在蒙特卡洛计算中采用互相独立的随机点来进展计算。对偶变量法中却使用相关联的点来进展计算。它利用相关点间的关系可以是正关联的，也可以是负关联的这个特点。两个函数值和之和的方差为如果我们选择一些点，它们使和是负关联的。这样就可以使上式所示的方差减小。当然这需要对具体的函数和有充分的了解独立产生U(0，1)随机数计算，找g(*),f(*)是相关的，且Eg(*)=计算3程序及实现结果3.1 的求解3.1.1 随机投点法先利用R 软件产生服从0,1上均匀分布的随机数 *,Y, ，计算的个数，即事件发生的频数，求出频率，即为积分的近似值。R程序s1-function(n) f-

10、function(*) e*p(-*) a-0 b-1 *-runif(n)y-runif(n) m-sum(yf(*) j=m/nvar-1/n*var(yf(*) lis-list(j,var)return(lis)s1(104)s1(105)s1(106)s1(107)对模拟次数n调试了4次，分别为，得到准确值和模拟值。表3.1.1 随机投点法的模拟次数和模拟值n模拟值0.62690.632350.6325030.6319298方差2.32599e-052.32414e-062.32713e-072.32477e-08准确值为0.6321206 3.1.2 平均值法先用R 软件产生n个服

11、从0,1上均匀分布的随机数，计算,再计算的平均值，即为定积分的近似值 R程序p1-function(n) f-function(*) e*p(-*) a-0 b-1 *-runif(n) y-mean(f(*)var-1/n*var(f(*)lis-list(y,var) return(lis)p1(104)p1(105)p1(106)p1(107)对模拟次数n调试了4次，分别为，得到准确值和模拟值。表3.1.2 平均值法的模拟次数和模拟值n模拟值0.63101170.63226430.63181990.632079方差3.25430e-063.27162e-073.2805e-083.275

12、36e-09准确值为0.6321206 3.1.3 重要性抽样法 R程序z1-function(n) *- 1-(sqrt(1-r) f-function(*) e*p(-*) g-function(*) (2*(1-*) r-runif(n) s=mean(f(*)/g(*)var-1/n*var(f(*)/g(*)lis-list(s,var) return(lis)z1(104)z1(105)z1(106)z1(107)对模拟次数n调试了4次，分别为，得到准确值和模拟值。表3.1.3 重要性抽样法法的模拟次数和模拟值n模拟值0.61348190.61348190.61348190.613

13、4819方差7.06957e-067.06957e-077.06957e-087.06957e-09准确值为0.6321206 3.1.4 分层抽样法 R程序f1-function(n,m) r1-runif(n,min-0,ma*-0.5) r2-runif(m,min-0.5,ma*-1)c-1/2*mean(e*p(-r1)+1/2*mean(e*p(-r2)var-var(e*p(-r1)/(4*n)+var(e*p(-r2)/(4*m)j-list(c,var)return(j)f1(10,20)f1(100,200)f1(1000,2000)f1(104,2*104)得到准确值和模

14、拟值表3.1.4 分层抽样法的模拟次数和模拟值取值n=10,m=20n=100,m=200n=1000M=2000n=10000M=20000模拟值0.65827020.62917550.63417180.6327054方差0.0002692072.75023e-053.23184e-063.18334e-07准确值为0.6321206 3.1.5 对偶变量法先用 R 软件产生服从0,1上均匀分布的随机数 * ，函数f(*),计算，计算 R程序d1-function(n)f-function(*) e*p(-*)y-function(*) e*p(-(1-*)*-runif(n)m-sum(f

15、(*)p-sum(y(*)j-(m/n+p/n)/2var-1/4*(var(f(*)+var(y(*)+2*cov(f(*),y(*)lis-list(j,var)return(lis)d1(104)d1(105)d1(106)d1(107) 对模拟次数n调试了4次，分别为，得到准确值和模拟值。表3.1.5 对偶变量法的模拟次数和模拟值n模拟值0.63219790.63206590.63209570.6321183方差0.000524720.0005327480.00053057410.0005291691准确值为0.6321206 3.1.6 控制变量法R程序k1-function(n)f

16、-function(*) e*p(-*)r-runif(n)g-function(*) 2*(1-*)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(104)k1(105)k1(106)k1(107)对模拟次数n调试了4次，分别为，得到准确值和模拟值。表3.1.6 控制变量法的模拟次数和模拟值n模拟值0.6345750.63208840.6318490.6321651方差5.713907e-055.72338e-06

17、5.736774e-075.731522e-08准确值为0.6321206 3.2 对积分求解3.2.1 随机投点法R程序s2-function(n) f-function(*) e*p(-*)t=function(y) (f(a+(b-a)*y)-c)/(d-c) a-2 b-4c-f(4)d-f(2)s-(b-a)*(d-c)*-runif(n)y-runif(n)m-sum(yf(*)jm/ng=s*j+c*(b-a)var-1/n*var(yf(*) lis-list(g,var)return(lis)s2(104)s2(105)s2(106)s2(107)对模拟次数n调试了4次，分别

18、为，得到准确值和模拟值。表3.2.1 随机投点法的模拟次数和模拟值n模拟值0.18688450.18688450.18688450.1868845方差2.33071e-052.32230e-062.32479e-072.32611e-08准确值为 0.1170196 3.2.2 平均值法R程序p2-function(n) f-function(r) e*p(-*) 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(104)p2(105)p2(106)p2(1

19、07)对模拟次数n调试了4次，分别为，得到准确值和模拟值。表3.2.2 平均值法的模拟次数和模拟值n模拟值 0.11733330.11695820.1170038 0.1170311方差1.30285e-051.30285e-061.30285e-071.30285e-08准确值为 0.1170196 3.2.3 分层抽样法R程序f2-function(n,m) r1-runif(n,min-0,ma*-0.5) r2-runif(m,min-0.5,ma*-1) c-1/2*mean(2*e*p(-2-2*r1)+1/2*mean(2*e*p(-2-2*r2) var-var(2*e*p(-

20、2-2*r1)/(4*n)+var(2*e*p(-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次，分别为，得到准确值和模拟值。表3.2.3 分层抽样法的模拟次数和模拟值取值n=10,m=20n=100,m=200n=1000M=2000n=10000M=20000模拟值0.11081320.12196680.11737690.1169783方差6.12767e-056.31761e-066.20002e-076.02697e-08准确值为 0.1170196

21、 3.2.4 对偶变量法R程序d2-function(n)a-2b-4f-function(*) e*p(-*)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)d2(106)d2(107)对模拟次数n调试了4次，分别为，得到准确值和模拟值。表3.2.4 对偶变量法的模拟次数和模拟值n模拟值 0.11690430.1167410.116943 0.1170

22、233准确值为 0.11701963.2.5 控制变量法R程序k2-function(n)f-function(*) e*p(-*)r-runif(n)a-2b-4c-f(4) d-f(2)s-(b-a)*(d-c)q-function(*) 1/(d-c)*(f(a+(b-a)*)-c)g-function(*) 2*(1-*)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(104)k2(

23、105)k2(106)k2(107)对模拟次数n调试了4次，分别为，得到准确值和模拟值。表3.2.5 控制变量法的模拟次数和模拟值n模拟值 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(*) e*p(-*)c=f(4)d=f(2)s=(b-a)*(d-c);r=runif(n)*- 1-(sqrt(1-r)p-function(*) 1/(d-c)*(f(a

24、+(b-a)*)-c)q-function(*) (2*(1-*)j1-mean(p(*)/q(*) j=s*j1+c*(b-a) var-1/n*var(p(*)/q(*)lis-list(j,var) return(lis)z2(104)z2(105)z2(106)z2(107)对模拟次数n调试了4次，分别为，得到准确值和模拟值。表3.2.6 重要性采样法的模拟次数和模拟值n模拟值 0.11735220.11702670.11701140.1170123方差8.17541e-078.16598e-088.17952e-098.17912e-10准确值为 0.11701963.3 积分求解3

25、.3.1 随机投点法R程序s3-function(n) f-function(*) e*p(-*)/(1+*2) a-0 b-1 *-runif(n)y-runif(n) m-sum(yf(*) j=m/nvar-1/n*var(yf(*) lis-list(j,var)return(lis)s3(104)s3(105)s3(106)s3(107)对模拟次数n调试了4次，分别为，得到准确值和模拟值。表3.3.1 随机投点法的模拟次数和模拟值n模拟值 0.53080.52507 0.5253510.5247777方差2.49625e-052.49312e-062.49423e-072.49371

26、e-08准确值为 0.5247971 平均值法R程序p3-function(n) f-function(*) e*p(-*)/(1+*2) a-0 b-1 *-runif(n) y-mean(f(*)var-1/n*var(f(*)lis-list(y,var) return(lis)p3(104)p3(105)p3(106)p3(107)对模拟次数n调试了4次，分别为，得到准确值和模拟值。表3.3.2 平均值法的模拟次数和模拟值n模拟值0.52397070.5252655 0.5245715 0.5247576方差5.97263e-066.02206e-075.99868e-085.9992

27、6e-09准确值为 0.5247971分层抽样法 R程序f3-function(n,m) r1-runif(n,min-0,ma*-0.5) r2-runif(m,min-0.5,ma*-1) z-function(u) e*p(-u)/(1+u2)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*104)得到准确值和模拟值表3.3.3 分层抽样法的模拟次数和模拟值取值n=1

28、0,m=20n=100,m=200n=1000M=2000n=10000M=20000模拟值0.53385040.53059470.52288030.5254361方差0.0003218132.11722e-052.16594e-062.17694e-07准确值为0.5247971 3.3.4 对偶变量法R程序d3-function(n)f-function(*) e*p(-*)/(1+*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,

29、var)return(lis)d3(104)d3(105)d3(106)d3(107)对模拟次数n调试了4次，分别为，得到准确值和模拟值。表3.3.4 对偶法的模拟次数和模拟值n模拟值0.63197880.6321444 0.6321236 0.6321169方差0.0005191260.00053013340.00052976580.0005295388准确值为 0.5247971 3.3.5 控制变量法R程序k3-function(n)f-function(*) e*p(-*)/(1+*2)r-runif(n)g-function(*) 2*(1-*)u-mean(g(r)l-(-cov(

30、f(r),g(r)/var(g(r) #定义lambdaj-mean(f(r)+l*mean(g(r)-u)var-1/n*var(f(r)+g(r)lst-list(j,var)return(lst)k1(104)k1(105)k1(106)k1(107)对模拟次数n调试了4次，分别为，得到准确值和模拟值表3.3.5 控制变量法的模拟次数和模拟值n模拟值0.6297403 0.632438 0.6320956 0.6321038方差5.80168e-055.74223e-065.73553e-075.73648e-08准确值为 0.5247971 3.3.6 重要性抽样法R程序z3-function(n) *- 1-(sqrt(1-r) f-function(*) e*p(-*)/(1+*2) g-function(*) (2*(1-*) r-runif(n) s=mean(f(*)/g(*)var-1/n*var(f(*)/g(*)lis-list(s,var) return(lis)z3(104)z3(105)z3(106)z3(107)对模拟次数