基于蒙特卡洛方法求數(shù)值積分與R

x0001.x0001.頁腳內(nèi)容頁腳內(nèi)容頁腳內(nèi)容頁腳內(nèi)容統(tǒng)計計算課程設(shè)計題目基于蒙特卡洛方法求數(shù)值積分中文摘要蒙特卡洛方法，又稱隨機抽樣或統(tǒng)計試驗方法，屬于計算數(shù)學(xué)的一個分支，它是在上世紀(jì)四十年代中期為了適應(yīng)當(dāng)時原子能事業(yè)的發(fā)展而發(fā)展起來的。傳統(tǒng)的經(jīng)驗方法由于不能逼近真實的物理過程，很難得到滿意的結(jié)果，而蒙特卡羅方法由于能夠真實地模擬實際物理過程，故解決問題與實際非常符合，可以得到很圓滿的結(jié)果。利用隨機投點法，平均值法，重要性采樣法，分層抽樣法，控制變量法，對偶變量法，運用氏軟件求0二J1e-xdx,00二J氏軟件求0二J1e-xdx,0201+x2出精度與樣本量的關(guān)系，比較各種方法的效率，關(guān)鍵字蒙特卡洛隨機投點法平均值法R軟件1緒論蒙特卡洛的基本思想是，當(dāng)所求解問題是某種隨機事件出現(xiàn)的概率，或者是某個隨機變量的期望值時，通過某種“實驗”的方法，以這種事件出現(xiàn)的頻率估計這一隨機事件的概率，或者得到這個隨機變量的某些數(shù)字特征，并將其作為問題的解。蒙特卡洛方法解題過程的三個主要步驟：（1）構(gòu)造或描述概率過程對于本身就具有隨機性質(zhì)的問題，如粒子輸運問題，主要是正確描述和模擬這個概率過程，對于本來不是隨機性質(zhì)的確定性問題，比如計算定積分，就必須事先構(gòu)造一個人為的概率過程，它的某些參量正好是所要求問題的解。即要將不具有隨機性質(zhì)的問題轉(zhuǎn)化為隨機性質(zhì)的問題。（2）實現(xiàn)從已知概率分布抽樣構(gòu)造了概率模型以后，由于各種概率模型都可以看作是由各種各樣的概率分布構(gòu)成的，因此產(chǎn)生已知概率分布的隨機變量（或隨機向量），就成為實現(xiàn)蒙特卡洛方法模擬實驗的基本手段，這也是蒙特卡洛方法被稱為隨機抽樣的原因。最簡單、最基本、最重要的一個概率分布是（，1）上的均勻分布（或稱矩形分布）。隨機數(shù)就是具有這種均勻分布的隨機變量。隨機數(shù)序列就是具有這種分布的總體的一個簡單子樣，也就是一個具有這種分布的相互獨立的隨機變數(shù)序列。產(chǎn)生隨機數(shù)的問題，就是從這個分布的抽樣問題。在計算機上，可以用物理方法產(chǎn)生隨機數(shù)，但價格昂貴，不能重復(fù)，使用不便。另一種方法是用數(shù)學(xué)遞推公式產(chǎn)生。這樣產(chǎn)生的序列，與真正的隨機數(shù)序列不同，所以稱為偽隨機數(shù)，或偽隨機數(shù)序列。不過，經(jīng)過多種統(tǒng)計檢驗表明，它與真正的隨機數(shù)，或隨機數(shù)序列具有相近的性質(zhì)，因此可把它作為真正的隨機數(shù)來使用。由已知分布隨機抽樣有各種方法，與從0,刀上均勻分布抽樣不同，這些方法都是借助于隨機序列來實現(xiàn)的，也就是說，都是以產(chǎn)生隨機數(shù)為前提的。由此可見，隨機數(shù)是我們實現(xiàn)蒙特卡洛模擬的基本工具。（3）建立各種估計量一般說來，構(gòu)造了概率模型并能從中抽樣后，即實現(xiàn)模擬實驗后，我們就要確定一個隨機變量，作為所要求的問題的解，我們稱它為無偏估計。建立各種估計量，相當(dāng)于對模擬實驗的結(jié)果進行考察和登記，從中得到問題的解。2方法介紹2.1隨機投點法隨機投點法是進行n次試驗，當(dāng)n充分大的時候，以隨機變量k/n作為期望值E(X)的近似估計值，即E(X)沁p=k/n其中k是n次實驗中成功的次數(shù)。若一次投點試驗的成功概率為p,并以v(1,表明試驗成功Ai[o,表明試驗失敗則一次試驗成功的均值與方差為E(X.)-1-p+0-(1-p)-piVar(X)-12-p+02-(1-p)-p2-p(1-p)i若進行n次試驗，其中化次試驗成功，則化為具有參數(shù)為(n,p)的二項分布，此時，隨機變量化的估計為p—k/n顯然，隨機變量的均值和方差滿足E(p)-E[-[-丄E(k)-pVn丿nVar(p)-p(-p)ddx0001.x0001.x0001.x0001.=1=1頁腳內(nèi)容=1=1頁腳內(nèi)容設(shè)計算的定積分為I=ibf(Xkx，其中a,b為有限數(shù)，被積函數(shù)yx是連續(xù)隨機變a量g的概率密度函數(shù)，因此yx滿足如下條件：f(x)非負，且Jf(x)dx=1顯然I是一個概率積分，其積分值等于概率P(a<g<b)。下面按給定分布yx隨機投點的辦法，給出如下％0沁Ca旳近似求積算法：1產(chǎn)生服從給定分布的隨機變量值,i=l,2,…N2檢查x是否落入積分區(qū)間。如果條件a<x.<b滿足，則記錄x落入積分區(qū)間一iii次。假設(shè)在N次實驗以后，落入積分區(qū)間的總次數(shù)為〃，那么用I=-N作為概率積分的近似值，即I沁-N2.2平均值法任取一組相互獨立、同分布的隨機變量4}，g在/a,b/內(nèi)服從分布率卩何令,則也是ii一組相互獨立、同分布的隨機變量，而且E百*()}=Jp*(x)p(x)dx=Jf(x》x=Iiaa由強大數(shù)定理lim丄￡p*(g)=INT8N.1x0001.x0001.x0001.x0001.頁腳內(nèi)容頁腳內(nèi)容頁腳內(nèi)容頁腳內(nèi)容若記Np弋)=1，則依概率1收斂到人平均值法就是用門乍為啲近似值。i=1假如所需計算積分為I=ff(x》x，其中被積函數(shù)在間6丿內(nèi)可積，任意選擇一個有a簡單辦法可以進行抽樣的概率密度函數(shù)pX，使其滿足條件:?p(x)h0,當(dāng)f(x)h0時(a<x<b)記p*(x)=lpC°則所求積分為10,p(x)=0I=fp*(x)p(x)dxa因而％o^e3也近似求積算法為:(1)產(chǎn)生服從分布率的隨機數(shù)x.,i=1,2,...,N(2)計算均值I=丄弋p*G),即有I沁INii=12.3重要性采樣法從數(shù)學(xué)角度上看定積分可以看成其中g(shù)(X是某個隨機變量X的密度函數(shù)，因此積分值I可看成隨機變量Z=f(x)/g(x)的數(shù)學(xué)期望值I=E(Z人丄迓z二丄迓蟲NiNgQ丿

i=1i=1i為了減少模擬實驗的方差應(yīng)適當(dāng)選取gX，使伽(刀盡可能小,如果被積函數(shù)f(X>0，可取當(dāng)c=i/i時就有伽池)=0.—般應(yīng)選取和fX相似的密度函數(shù)g(X，吏f(x)/g(x)接近于常數(shù)，故而伽厲接近于0,以達到降低模擬實驗的方差，這種減少方差的模擬試驗法為重要抽樣法。2.4分層抽樣法分層抽樣法是利用貢獻率大小來降低估計方差的方法。它首先把樣本空間D分成一些不交的小區(qū)間D=U，然后在各個小區(qū)間內(nèi)的抽樣數(shù)由其貢獻大小決定。即，定義p=If(x，則D內(nèi)的抽樣數(shù)n應(yīng)與p成正比。iiiiDi考慮積分0=If(x將0,1］分成m個小區(qū)間：0則0=If(x=XIf(x=瓦ITOC\o"1-5"\h\zi0i=1ai=1i-1記l=a-a為第i個小區(qū)間的長度,i=1,2,...,m，在每個小區(qū)間上的積分值可用均iii-1值法估計出來，然后將其相加即可給出0的一個估計。具體步驟為：1)獨立產(chǎn)生U0，1)隨機數(shù)｛uj,j=1,...2)計算x=a+1u,j=1,...iji-1iij3)計算I=丄為f(x)nijij=1

于是9于是9的估計為9二遲I,ii=1其方差為Var空n其中，i-1I2其中，i-1i2.5對偶變量法控制變量法利用數(shù)學(xué)上積分運算的線性特性??Jf(x=Jf(x)—g(x)Lx+fg(x選擇函數(shù)g(x時要考慮到：g(x在整個積分區(qū)間都是容易精確算出,并且在上式右邊第一項的運算中對(f-g)積分的方差應(yīng)當(dāng)要比第二項對f積分的方差小。在應(yīng)用這種方法時，在重要抽樣法中所遇到的，當(dāng)g(x｝趨于零時，被積函數(shù)f-g)趨于無窮大的困難就不再存在，因而計算出的結(jié)果穩(wěn)定性比較好。該方法也不需要從分布密度函數(shù)g(x，解析求出分布函數(shù)G(x)。由此我們可以看出選擇g^(j)所受到的限制比重要抽樣法要小些。模擬過程：1)獨立產(chǎn)生U(0,1)隨機數(shù)｛u.,j=1,...2計算9二丄工f(x)(~二丄工f(1—x)TOC\o"1-5"\h\z5ni5nii=1i=13)計算9=￡^?=1g—ii—i=12n2i=126控制變量法通常在蒙特卡洛計算中采用互相獨立的隨機點來進行計算。對偶變量法中卻使用相關(guān)聯(lián)的點來進行計算。它利用相關(guān)點間的關(guān)系可以是正關(guān)聯(lián)的，也可以是負關(guān)聯(lián)的這個特點。兩個函數(shù)值f和f之和的方差為12V{/+f}=V{f}+V{f}+2E《f-e{/Rf-e{/DI12121122如果我們選擇一些點，它們使f和f是負關(guān)聯(lián)的。這樣就可以使上12式所示的方差減小。當(dāng)然這需要對具體的函數(shù)f和f有充分的了解121)獨立產(chǎn)生U0，狎隨機數(shù)｛uij,j=1,...刀計算°2=AfU找翊^是相關(guān)的，且T卩i=13）計算e=e+九—g（g（x）—卩）62nii=13程序及實現(xiàn)結(jié)果e=11e-xdx的求解0隨機投點法先利用R軟件產(chǎn)生服從［0,1］上均勻分布的隨機數(shù)XY，計算y<f(x)的個數(shù)，即事件發(fā)生的頻數(shù)，求出頻率，即為積分的近似值。R程序s1<-function（n）f<-function(x)exp(-x)a<-0b<-1x<-runif(n)y<-runif(n)m<-sum(y<f(x))j=m/nvar<-1/n*var(y<f(x))lis<-list(j,var)return(lis)s1(10'4)s1(10A5)s1(10'6)s1(10'7)對模擬次數(shù)n調(diào)試了4次，分別為n4,n5,n6,n7，得到精確值和模擬值。表3.1.1隨機投點法的模擬次數(shù)和模擬值n104105106107模擬值^.6269^.632350.6325030.6319298方差2.32599e-052.32414e-062.32713e-072.32477e-08精確值為0.63212063.1.2平均值法先用R軟件產(chǎn)生n個服從0,1]上均勻分布的隨機數(shù)x.，計算f(x),再計算iif(x)的平均值，即為定積分的近似值iR程序p1<-function(n){f<-function(Xexp(-x)a<-0b<-1X<-runif(n)y<-meanf(x))var<-1/n*var(f(x))lis<-[ist(y,var)return(lis)pl(10"4)pl(10'5)p1(10，6)p1(10"7)對模擬次數(shù)n調(diào)試了4次，分別為n4,n5,n6,n7，得到精確值和模擬值。表3.1.2平均值法的模擬次數(shù)和模擬值n104105106107模擬值0.63101170.63226430.63181990.632079方差3.25430e-063.27162e-073.2805e-083.27536e-09精確值為0.63212063.1.3重要性抽樣法R程序z1<-function(n){-1_(sqrt(1_r))f<-function(Xexp(-x)g<-function(X(2*(1_勸r<_runf(n)s=mean(f(x)/g(x))var<-1/n*var(f(x)/g(xj)Hs<-list(s,var)return(lis)}zl(10"4)z1(10A5)zl(10'6)zl(10"7)對模擬次數(shù)n調(diào)試了4次，分別為n4,n5,n6,n7，得到精確值和模擬值。表3.1.3重要性抽樣法法的模擬次數(shù)和模擬值n104105106107模擬值0.61348190.61348190.61348190.6134819方差7.06957e-067.06957e-077.06957e-087.06957e-09精確值為0.6321206.X0OQ1.X0OQ1x0001.x0001..X0OO1..X0OO1.頁腳內(nèi)容頁腳內(nèi)容頁腳內(nèi)容頁腳內(nèi)容頁腳內(nèi)容頁腳內(nèi)容3.1.4分層抽樣法R程序f1<function(n,m)r1<-runf(n,min<-0,max<-0.5)r2<-runf(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分層抽樣法的模擬次數(shù)和模擬值

n=10,m=20n=100,n=1000n=10000m=200M=2000n=10,m=20n=100,n=1000n=10000m=200M=2000M=200000.65827020.62917550.63417180.63270540.0002692072.75023e-053.23184e-063.18334e-07為值精確值0.63212063.1.5對偶變量法先用R軟件產(chǎn)生服從［0,1］上均勻分布的隨機數(shù)X，函數(shù)f(x),計算『二丄￡f(x)0~=丄￡f(1-x)5n15n1i=1，i=1f(x)+fG-x)iif(x)+fG-x)ii—2U=—5A=—乙計算2ni=1R程序d1<function(n)f<-function(x)exp(-x)y<-function(x)exp(-(1-x))X<-runf(n)m<-sumf(x))p<^-sum^(y(x))j<-(m/n+p/n)/2var<-1/4^(yar(f(xj)+var(y(xj)+2*cov(f(x),y(x))j)Hs<-Hst(j,var)return(Us)}dl(10"4)dl(10'5)dl(10,6)dl(10"7)對模擬次數(shù)n調(diào)試了4次，分別為n4,n5,n6,n7，得到精確值和模擬值。表3.1.5對偶變量法的模擬次數(shù)和模擬值n104105106107模擬值0.63219790.6320659^.632^9570.6321183方差0.000524720.0005327480.00053057410.0005291691精確值為0.63212063.1.6控制變量法R程序kj<-function(n)f<function(x)exp(-x)r<-runf(n)g<function(x)2*(1-x)u<-mean(g(r))l<-(-cov(f(r),g(r)))/var(g(r))j<-meanf(r》+l*mean(g(r)-u)var<-1/n*varf(r)+g(r))lst<-list(j,var)return(lst)}k!(io,4)k!(10'5)k1(10'6)k1(10'7)對模擬次數(shù)n調(diào)試了4次，分別為n4,n5,n6,n7，得到精確值和模擬值。表3.1.6控制變量法的模擬次數(shù)和模擬值x0001.x0001.x0001.x0001.頁腳內(nèi)容頁腳內(nèi)容頁腳內(nèi)容頁腳內(nèi)容n104105106107模擬值0.6345750.63208840.6318490.6321651方差5.713907e-055.72338e-065.736774e-075.731522e-08精確值為0.63212063.2對0e-xdx積分求解23.2.1隨機投點法R程序s2<-function(n){f<-function(x)exp(-x)t=function(y)(f(a+(b_a)*y)-c)/(d-c)a<-2b<-4cfd<f2)s<-(b-a)*(d-c)X<-runf(n)y<_runf(n)m<-sum(y<f(x))j<m/ng=s*j+c*(b-a)var<-l/n*var(y<f(X)Hs<-Hst(g,var)return(Us)}s2(10'4)s2(10"5)s2(10"6)s2(10'7)對模擬次數(shù)n調(diào)試了4次，分別為n4,n5,n6,n7，得到精確值和模擬值。表3.21隨機投點法的模擬次數(shù)和模擬值n104105106107模擬值0.18688450.18688450.18688450.1868845方差2.33071e-052.32^3^e-^62.32479e-072.32611e-08精確值為0.11701963.2.2平均值法R程序p2<-function(n)f<function(r)exp(-x)b<-4r<-runf(n)h<-(b-ayf(a+(b-a)*r)y<-mean(h)var<-1/n*var(h)Hs<-Hst(y,var)return(lis)}p2(10'4)p2(10"5)p2(10，6)p2(10'7)對模擬次數(shù)n調(diào)試了4次，分別為n4,n5,n6,n7，得到精確值和模擬值。表3.22平均值法的模擬次數(shù)和模擬值n104105106107模擬值0.11733330.11695820.11700380.1170311方差1.30285e-051.3^^85e-^61.30285e-071.30285e-08精確值為0.11701963.2.3分層抽樣法R程序f2<function(n,m){r1<-runf(n,min<-0,max<-0.5)r2<-runf(m,min<-0.5,max<-l)c<-l/2*mean(2*exp(-2-2*rl))+l/2*mean(2*exp(-2-2*r2))var<-var(2*exp(-2_2*rl))/(4*n)+var(2*exp(-2_2*r2)/(4*m))j<-list(c,var)return(j)}f2(10,20)f2(100,200)f2(1O00,200O)f2(10'4,2*10'4)對模擬次數(shù)n調(diào)試了4次，分別為n4,n5,n6,n，得到精確值和模擬值。表3.23分層抽樣法的模擬次數(shù)和模擬值n=10,m=20n=100,n=1000n=10000確值m=200M=2000M=20000值0.11081320.12196680.11737690.11697836.12767e_056.31761e_066.2^^^2e_^76.^2697e_^80.11701963.2.4對偶變量法R程序d2<-function(n){a^<-2b<-4f<-function(x)exp(-x)r<-runf(n)dfs<_(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-o)return(j)}d2(10"4)d2(10'5)d2(10，6)d2(10'7)對模擬次數(shù)n調(diào)試了4次，分別為n4,n5,n6,n，得到精確值和模擬值。表3.2.A對偶變量法的模擬次數(shù)和模擬值n104105106107模擬值0.11690430.1167410.1169430.1170233精確值為0.1170196325控制變量法R程序k2<-function(n)xO001.xO001.xO001.xO001.x0001.x0001.c=f(c=f(4)頁腳內(nèi)容c=f(c=f(4)頁腳內(nèi)容k2(10Ak2(10A4)頁腳內(nèi)容f<function(x)exp(-x)r<-runf(n)a<-2b<-4c<-f(4)d今⑵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))[<-(-cov(f(r),g(r))^/var(g(r))p<-mean(q(r))+C*mean(g(r》u)j<-s*p+c*(b-a)var<-1/n*varf(r)+g(r))lst<-[ist(j,var)return(lSt)k?(io、5)k2(io*)k2(i0，7)對模擬次數(shù)n調(diào)試了4次，分別為n4,n5,n6,n7，得到精確值和模擬值。表3.25控制變量法的模擬次數(shù)和模擬值n104105106107模擬值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)^xp(-x)x0001.x0001.x0001.x0001.x0001.x0001.m<-sum(y<f(x))頁腳內(nèi)容m<-sum(y<f(x))頁腳內(nèi)容m<-sum(y<f(x))頁腳內(nèi)容m<-sum(y<f(x))頁腳內(nèi)容對模擬次數(shù)n對模擬次數(shù)n調(diào)試了4次，分別為n4,n5,n6,n，得到精確值和模擬值。頁腳內(nèi)容d=f(2)s=(b-a)*(d-c);r=runf(n)x<-l_(sqrt(l_r))p<function(x)1/(d-c)*^(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<-[ist(j,var)return(lis)}z2(10"4)z2(10'5)z2(10"6)z2(10'7)表3.26重要性采樣法的模擬次數(shù)和模擬值n104105106107模擬值0.11735220.11702670.11701140.1170123方差8.17541e-078.16598e-088.17952e-098.17912^-10精確值為0.11701963.30Jdx積分求解01+X23.3.1隨機投點法R程序s3<-function(n){f<function(Xexp(-x)/(1+x?2)a<-0b<-1x<-runif(n)y<-runif(n)x0001.x0001.x0001.x0001.x0001.x0001.頁腳內(nèi)容頁腳內(nèi)容頁腳內(nèi)容頁腳內(nèi)容p3<-function(n)頁腳內(nèi)容p3<-function(n)頁腳內(nèi)容j=m/nvar<-1/n*var(y<f(x))lis<-list(j,var)return(lis)}s3(10'4)s3(10'5)s3(10'6)s3(10'7)對模擬次數(shù)n調(diào)試了4次，分別為n4,n5,n6,n7，得到精確值和模擬值。表3.3.1隨機投點法的模擬次數(shù)和模擬值n104105106107模擬值C^.53C^8C.525C7C.525351C.5247777方差2.49625e-052.49312e-C62.49423e-C72.49371e-C8精確值為C.52479713.3.2平均值法R程序f<function(Xexp(-x)/(1+x?2)a<-0b<-1X<-runif(n)y<-meanf(x))var<-1/n*varf(X)Hs<-Hst(y,var)return(lis)}p3(10,4)p3(10"5)p3(10"6)p3(10'7)對模擬次數(shù)n調(diào)試了4次，分別為n4,n5,n6,n7，得到精確值和模擬值。表3.3.2平均值法的模擬次數(shù)和模擬值x0001.x0001.x0001.x0001.f3(100,200)頁腳內(nèi)容f3(100,200)頁腳內(nèi)容f3(100,200)頁腳內(nèi)容f3(100,200)頁腳內(nèi)容n104105106107模擬值^.52397^70.52526550.52457150.5247576方差5.97263e-066.02206e-075.99868e-^85.99926e-09精確值為0.52479713.3.3分層抽樣法R程序f3<function(n,m){r1<-runf(n,min<-0,max<-0.5)r2<-runf(m,min<-0.5,max<-l)z<function(u)exp(-u)/(l+u，2)j^-1/2*m^ea^n^(z(i^1)(^+1/2*m^e^n^(z(i^2))var<-var(z(-rl))/(4*n)+var(z(-r2))/(4*m)lis<-[ist(j,var)return([is)}f3(10,20)^3(1000,2000)x0001.x0001.x0001.x0001.x0001.x0001.頁腳內(nèi)容頁腳內(nèi)容頁腳內(nèi)容頁腳內(nèi)容var<_1/4'*(varf(r))+varf(1-r))+2'*covf(r),f(1-r)))頁腳內(nèi)容var<_1/4'*(varf(r))+varf(1-r))+2'*covf(r),f(1-r)))頁腳內(nèi)容f3(10'4,2*10'4)得到精確值和模擬值表3.3.3分層抽樣法的模擬次數(shù)和模擬值n=10,m=20為n=100,n=1000n=10,m=20為精確值m=200M=2000M=20000值0.53385040.53059470.52288030.52543610.0003218132.11722e-052.16594e-062.17694e-070.52479713.3.4對偶變量法R程序d3<function(n){f<function(Xexp(-x)^/(1+x^2)r<-runf(n)m<-meanf(r))p<-meanf(1-r))j<_(m+p)/2Us<_[ist(j,var)return(Rs)d3(10"4)d3(10A5)d3(10，6)d3(10'7)對模擬次數(shù)n調(diào)試了4次，分別為n4,n5,n6,n7，得到精確值和模擬值。表3.3.4對偶法的模擬次數(shù)和模擬值n104105106107模擬值0.63197880.63214440.63212360.6321169方差0.0005191260.00053013340.00052976580.0005295388精確值為0.52479713.3.5控制變量法R程序k}<-function(n){f<-function(x)exp(-x)/(1+x^2)r<-runf(n)g<function(x)2*(1-x)u<-mean(g(r》l<_(-covf(r),g(r)))/var(g(r))#定義lambdaj<-meanf(r》+C*mean(g(r)-u)var<_l/n*varf(r)+g(r))(st<_(ist(j,var)return(lSt)}k!(io,4)k!(10'5)k!(10'6)k!(10'7)對模擬次數(shù)n調(diào)試了4次，分別為n4,n5,n6,n7，得到精確值和模擬值表3.3.5控制變量法的模擬次數(shù)和模擬值n104105106107模擬值0.62974030.6324380.63209560.6321038方差5.80168e_055.74223e_065.73553e_075.73648e_08精確值為0.5247971x0001.x0001.x0001.x0001..X0OO1..X0OO1.頁腳內(nèi)容頁腳內(nèi)容頁腳內(nèi)容頁腳內(nèi)容頁腳內(nèi)容頁腳內(nèi)容3.3.6重要性抽樣法R程序z3<-function(n)x<-1-(sqrt(1-r))f<function(Xexp(-x)/(1+x^2)g<function(x)(2*(1-x))r<-runf(n)s=meanf