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

1、v1.0可編輯可修改統(tǒng)計計算課程設(shè)計題 目基于蒙特卡洛方法求數(shù)值積分中文摘要蒙特卡洛方法，又稱隨機抽樣或統(tǒng)計試驗方法，屬于計算數(shù)學的一個分支，它是在上世紀四十年代中期為了適應(yīng)當時原子能事業(yè)的發(fā)展而發(fā)展起來的。傳統(tǒng)的經(jīng)驗方法由于不能逼近真實,物理過程，很難得到滿意的結(jié)果，而蒙特卡羅方法由于能夠真實地模擬實際物理過程，故解決_可題與實際非常符合，可以得到很圓滿的結(jié)果。v1.0可編輯可修改1用R軟件求e xdx ,04e xdx 和201 x2dx數(shù)值積分。計算以上各種估計的方差,利用隨機投點法，平均值法，重要性采樣法，分層抽樣法，控制變量法，對偶變量法，運給出精度與樣本量的關(guān)系，比較各種方法的效率

2、,關(guān)鍵字蒙特卡洛隨機投點法平均值法R軟件v1.0可編輯可修改1緒論蒙特卡洛的基本思想是，當所求解問題是某種隨機事件出現(xiàn)的概率，或者是某個隨機變量的 期望值時，通過某種“實驗”的方法，以這種事件出現(xiàn)的頻率估計這一隨機事件的概率，或者得 到這個隨機變量的某些數(shù)字特征，并將其作為問題的解。蒙特卡洛方法解題過程的三個主要步驟：(1)構(gòu)造或描述概率過程對于本身就具有隨機性質(zhì)的問題，如粒子輸運問題，主要是正確描述和模擬這個概率過程,對于本來不是隨機性質(zhì)的確定性問題，比如計算定積分，就必須事先構(gòu)造一個人為的概率過程， 它的某些參量正好是所要求問題的解。即要將不具有隨機性質(zhì)的問題轉(zhuǎn)化為隨機性質(zhì)的問題。(2)實

3、現(xiàn)從已知概率分布抽樣構(gòu)造了概率模型以后，由于各種概率模型都可以看作是由各種各樣的概率分布構(gòu)成的，因此 產(chǎn)生已知概率分布的隨機變量(或隨機向量)，就成為實現(xiàn)蒙特卡洛方法模擬實驗的基本手段， 這也是蒙特卡洛方法被稱為隨機抽樣的原因。最簡單、最基本、最重要的一個概率分布是(0, 1)上的均勻分布(或稱矩形分布)。隨機數(shù)就是具有這種均勻分布的隨機變量。隨機數(shù)序列就是具 有這種分布的總體的一個簡單子樣，也就是一個具有這種分布的相互獨立的隨機變數(shù)序列。產(chǎn)生 隨機數(shù)的問題，就是從這個分布的抽樣問題。在計算機上，可以用物理方法產(chǎn)生隨機數(shù)，但價格 昂貴，不能重復，使用不便。另一種方法是用數(shù)學遞推公式產(chǎn)生。這樣產(chǎn)

4、生的序列，與真正的隨 機數(shù)序列不同，所以稱為偽隨機數(shù)，或偽隨機數(shù)序列。不過，經(jīng)過多種統(tǒng)計檢驗表明，它與真正 的隨機數(shù)，或隨機數(shù)序列具有相近的性質(zhì)，因此可把它作為真正的隨機數(shù)來使用。由已知分布隨 機抽樣有各種方法，與從(0,1)上均勻分布抽樣不同，這些方法都是借助于隨機序列來實現(xiàn)的， 也就是說，都是以產(chǎn)生隨機數(shù)為前提的。由此可見，隨機數(shù)是我們實現(xiàn)蒙特卡洛模擬的基本工具。(3)建立各種估計量一般說來，構(gòu)造了概率模型并能從中抽樣后，即實現(xiàn)模擬實驗后，我們就要確定一個隨機變 量，作為所要求的問題的解，我們稱它為無偏估計。建立各種估計量，相當于對模擬實驗的結(jié)果 進行考察和登記，從中得到問題的解。v1.0

5、可編輯可修改2方法介紹隨機投點法隨機投點法是進行n次試驗，當n充分大的時候，以隨機變量 k/n作為期望值 E(X)的近似估計值，即E(X) p k/n其中k是n次實驗中成功的次數(shù)。若一次投點試驗的成功概率為 p,并以1,表明試驗成功Xi0,表明試驗失敗則一次試驗成功的均值與方差為E(Xi) 1 p 0 (1 p) pVar(Xi) 12 p 02 (1 p) p2 p(1 p)若進行n次試驗，其中k次試驗成功，則k為具有參數(shù)為(n,p)的二項分布， 此時，隨機變量k的估計為p k/n顯然，隨機變量的均值和方差滿足_ k 1E p E Ek p n nVar p -p-1 ddv1.0可編輯可修

6、改設(shè)計算的定積分為I bf xdx,其中a,b為有限數(shù)，被積函數(shù)f(x)是連續(xù)隨機 a變量的概率密度函數(shù)，因此f(x)滿足如下條件：f x 非負，且 f (x)dx 1顯然I是一個概率積分，其積分值等于概率 P(a b)下面按給定分布f(x)隨機投點的辦法，給出如下Monte Carlo近似求積算法：(1)產(chǎn)生服從給定分布的隨機變量值，i=1 ,2,N;(2)檢查k是否落入積分區(qū)問。如果條件a xi b滿足，則記錄x落入積分區(qū)問 一次。假設(shè)在N實驗以后，落入積分區(qū)間的總次數(shù)為 n,那么用I -N作為概率積分的近似值，即I -N平均值法任取一組相互獨立、同分布的隨機變量i , i在a,b內(nèi)服從分

7、布率p(x),令,則也是一組相互獨立、同分布的隨機變量,而且b*E p i pax p x dxbf xdxa由強大數(shù)定理p lim N nNp*i 1v1.0可編輯可修改1 N_ 若記 L p* i I,則依概率1收斂到I,平均值法就是用I作為I的近似值。N i ib假如所需計算積分為I f xdx,其中被積函數(shù)在a,b內(nèi)可積，任意選擇一個a有簡單辦法可以進行抽樣的概率密度函數(shù)p(x),使其滿足條件:p x 0,當 f x 0時 a x bb p(x)dx 1af x、一 *, p x 0.一.記p x p x則所求積分為0, p x 0bI p* x p x dxa因而Monte Carl

8、o近似求積算法為：(1)產(chǎn)生服從分布率p(x)的隨機數(shù)xj 1,2,., N計算均值I -Np* i ,即有I重要性采樣法從數(shù)學角度上看定積分可以看成10g其中g(shù)(x)是某個隨機變量 X的密度函數(shù)，因此積分值I可看成隨機變量v1.0可編輯可修改Z=f(x)/g(x)的數(shù)學期望值NN11fxizi-N i i N i i g Xi為了減少模擬實驗的方差應(yīng)適當選取 g(x),使Var(I)盡可能小，如果被積函數(shù) f(x)0,可取g(x)=cf(x), 當c=1/I時就有Var(I)=0. 一般應(yīng)選取和f(x)相似的密度 函數(shù)g(x),使f(x)/g(x) 接近于常數(shù)，故而Var(I)接近于0,以達

9、到降低模擬實驗的 方差，這種減少方差的模擬試驗法為重要抽樣法。分層抽樣法分層抽樣法是利用貢獻率大小來降低估計方差的方法。它首先把樣本空間D分成m一些不交的小區(qū)間D U,然后在各個小區(qū)間內(nèi)的抽樣數(shù)由其貢獻大小決定。即，定 i 1義Pif xdx ,則Di內(nèi)的抽樣數(shù)ni應(yīng)與r成正比。Di1考慮積分 f xdx將0,1分成m個小區(qū)間： 0 1m aim0 a0 a1 . am 1f xdx f x dx Ii0i 1 ai 1i 1,在每個小區(qū)間上的積分值記li aia為第i個小區(qū)間的長度，i=1,2,.,m可用均值法估計出來，然后將其相加即可給出的一個估計。具體步驟為:1)獨立產(chǎn)生U(0, 1)隨

10、機數(shù)Uj,j 1,.2)計算 xjai 1 Lu , j 1,.v1.0可編輯可修改3）計算I? Lnilif Xijj 1于是的估計為？mI?,其方差為Var ?L2 2h iniai乜Xli對偶變量法控制變量法利用數(shù)學上積分運算的線性特性：f x dx f x g x dx g xdx選擇函數(shù) g(x)時要考慮到：g(x)在整個積分區(qū)間都是容易精確算出，并且在上式右邊第一項的運算中對(f-g)積分的方差應(yīng)當要比第二項對f積分的方差小。在應(yīng)用這種方法時，在重要抽樣法中所遇到的，當 g(x)趨于零時，被積 函數(shù)(f-g)趨于無窮大的困難就不再存在，因而計算出的結(jié)果穩(wěn)定性比較 好。該方法也不需要

11、從分布密度函數(shù)g(x),解析求出分布函數(shù)G(x)。由此我們可以看出選擇g(x)所受到的限制比重要抽樣法要小些。模擬過程：1)獨立產(chǎn)生U(0, 1)隨機數(shù)Uj,j 1,.n n2)計算？ 一f (xi)5 f 1 xin i 1n i 1v1.0可編輯可修改1n f X fixi2控制變量法通常在蒙特卡洛計算中采用互相獨立的隨機點來進行計算。對偶變量法中卻使用相關(guān)聯(lián)的點來進行計算。它利用相關(guān)點間的關(guān)系可以是正關(guān)聯(lián)的，也可以是負關(guān)聯(lián)的這個特點。兩個函數(shù)值 fi和f2之和的方差為Vf1f2Vf1Vf22Ef1Ef1f2 E f2如果我們選擇一些點，它們使 片和f2是負關(guān)聯(lián)的。這樣就可以使上式所示的方

12、差減小。當然這需要對具體的函數(shù)f1和f2有充分的了解1)獨立產(chǎn)生U(0, 1)隨機數(shù)uj,j 1,.1 n2)計算？； f Xi ,找 g(x),f(x)是相關(guān)的，且 Eg(x)=n i 1 n3)計算 Z ?21g(xi)n i 13程序及實現(xiàn)結(jié)果1 e xdx的求解 0隨機投點法先利用R軟件產(chǎn)生服從0,1上均勻分布的隨機數(shù)X,Y,計算y f(x)的個數(shù),v1.0可編輯可修改即事件發(fā)生的頻數(shù)，求出頻率，即為積分的近似值R程序s1-function(n)(f-function(x) exp(-x)a-0b-1x-runif(n)y-runif(n)m-sum(yf(x)j=m/nvar-1/n

13、*var(yf(x)lis-list(j,var)return(lis)s1(10A4)s1(10A5)s1(10A6)v1.0可編輯可修改s1(10A7)對模擬次數(shù)n調(diào)試了 4次，分別為n4,n5,n6,n7,得到精確值和模擬值。表 隨機投點 法的模 擬次數(shù)和模 擬值n104105106107模擬值萬差平均值法先用R軟件產(chǎn)生n個服從0,1上均勻分布的隨機數(shù)x ,計算f(xi),再計算”為)的平均值，即為定積分的近似值R 程序p1-function(n)(f-function(x) exp(-x)a-0 |b-1x-runif(n)y-mean(f(x)v1.0可編輯可修改var-1/n*va

14、r(f(x)lis-list(y,var)return(lis)p1(10A4)p1(10A5)p1(10A6)p1(10A7)對模擬次數(shù)n調(diào)試了 4次，分別為n4,n5,n6,n7 ,得到精確值和模擬值。表 平均值法 的模擬 次數(shù)和模擬 值n104105106107模擬值、.、.廣. 力左精確值為重要性抽樣法R 程序z1-function(n)10v1.0可編輯可修改(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

15、(s,var)return(lis)z1(10A4)z1(10A5)z1(10A6)z1(10A7)對模擬次數(shù)n調(diào)試了 4次，分別為n4,n5,n6,n7 ,得到精確值和模擬值。表重要性抽樣法法的模擬次數(shù)和模擬值11v1.0可編輯可修改n104105106107模擬值萬差精確值為分層抽樣法R 程序f1-function(n,m)(r1-runif(n,min-0,maxr2-runif(m,min,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)r

16、eturn(j)f1(10,20)f1(100,200)12v1.0可編輯可修改f1(1000,2000)f1(10A4,2*10A4)得到精確值和模擬值表分層抽樣法的模擬次數(shù)和模擬值n=10,m=20n=100,n=1000n=10000精確值m=200M=2000M=20000值對偶變量法先用 R軟件產(chǎn)生服從0,1上均勻分布的隨機數(shù) X ,函數(shù)f(x), TOC o 1-5 h z 、一c 1n1n計算 2 f (Xi)5 - f 1 Xin i 1n i 1?51 n f X f 1 X計算.Hi12R 程序d1-function(n)(f-function(x) exp(-x)y-fu

17、nction(x) exp(-(1-x)13v1.0可編輯可修改x-runif(n)m-sum(f(x)p-sum(y(x)j-(m/n+p/n)/2var-1/4*(var(f(x)+var(y(x)+2*cov(f(x),y(x)lis-list(j,var)return(lis)d1(10A4)d1(10A5)d1(10A6)d1(10A7)對模擬次數(shù)n調(diào)試了 4次，分別為n4,n5,n6,n7 ,得到精確值和模擬值。表 對偶變量 法的模 擬次數(shù)和模 擬值n104105106107模擬值萬差14v1.0可編輯可修改精確值為控制變量法R程序k1-function(n)(f-function

18、(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(10A4)k1(10A5)k1(10A6)k1(10A7)15v1.0可編輯可修改對模擬次數(shù)n調(diào)試了 4次，分別為n4,n5,n6,n7 ,得到精確值和模擬值。表 控制變量 法的模 擬次數(shù)和模 擬值n104105106107模擬值.、.、.廣.萬差精確值為4對 e Xdx積分求解2隨

19、機投點法R程序s2-function(n)(f-function(x) exp(-x)t=function(y) (f(a+(b-a)*y)-c)/(d-c)a-2b-4c-f(4)16v1.0可編輯可修改d-fs-(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(10A4)s2(10A5)s2(10A6)s2(10A7)對模擬次數(shù)n調(diào)試了 4次，分別為n4,n5,n6,n7,得到精確值和模擬值。表 隨機投點 法的模 擬次數(shù)和模 擬值1

20、7v1.0可編輯可修改104105106107模擬值、.、.r. 萬差精確值為平均值法R程序p2-function(n)f-function(r) exp(-x)a-2b-4r-runif(n)h-(b-a)*f(a+(b-a)*r)y-mean(h)var-1/n*var(h)lis-list(y,var)return(lis)18v1.0可編輯可修改p2(10A4)p2(10A5)p2(10A6)p2(10A7)對模擬次數(shù) n調(diào)試了 4次，分別為n4,n5,n6,n7,得到精確值和模擬值。表 平均值法 的模擬 次數(shù)和模擬 值n104105106107模擬值萬差精確值為分層抽樣法R程序f2-

21、function(n,m)r1-runif(n,min-0,maxr2-runif(m,min,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)19v1.0可編輯可修改return(j)f2(10,20)f2(100,200)f2(1000,2000)f2(10A4,2*10A4)對模擬次數(shù)n調(diào)試了 4次，分別為n4,n5,n6,n7,得到精確值和模擬值。表分層抽樣法的模擬次數(shù)和模擬值n=10,m=20

22、n=100,n=1000n=10000精,m=200M=2000M=20000確值為值對偶變量法R程序d2-function(n)(a-2b-420v1.0可編輯可修改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)/2j-s*j3+c*(b-a)return(j)d2(10A4)d2(10A5)d2(10A6)d2(10A7)對模擬次數(shù)n調(diào)試了 4次，分別為n4,n5,n6,n7,得到精

23、確值和模擬值。表 對偶變量 法的模 擬次數(shù)和模 擬值21v1.0可編輯可修改n104105106107模擬值精確值為控制變量法R程序k2-function(n)(f-function(x) exp(-x)r-runif(n)a-2b-4c-f(4)d-fs-(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)22v1.0可編輯可修改p-mean(q(r)+l*mean(g(r)-u)j-s*p+c*(b-a)var-1/n*var(f

24、(r)+g(r)lst-list(j,var)return(lst)k2(10A4)k2(10A5)k2(10A6)k2(10A7)對模擬次數(shù) n調(diào)試了 4次，分別為n4,n5,n6,n7,得到精確值和模擬值。表 控制變量 法的模 擬次數(shù)和模 擬值n104105106107模擬值0.萬差精確值為23v1.0可編輯可修改重要性采樣法R程序z2-function(n)(a=2b=4f=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

25、-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)24v1.0可編輯可修改return(lis)z2(10A4)z2(10A5)z2(10A6)z2(10A7)對模擬次數(shù)n調(diào)試了 4次，分別為n4,n5,n6,n7,得到精確值和模擬值。表 重要性采 樣法的 模擬次數(shù)和 模擬值n104105106107模擬值萬差精確值為XJdx積分求解0 1 X隨機投點法R程序25v1.0可編輯可修改s3-function(n)(f-function(x) exp(-x)/(1+xA2

26、)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)s3(10A4)s3(10A5)s3(10A6)s3(10A7)對模擬次數(shù) n調(diào)試了 4次，分別為n4,n5,n6,n7,得到精確值和模擬值。26v1.0可編輯可修改表 隨機投點 法的模 擬次數(shù)和模 擬值n104105106107模擬值萬差精確值為平均值法R程序p3-function(n)(f-function(x) exp(-x)/(1+xA2)a-0 |b-1x-runif(n)y-mean(f(x)var-1/n*va

27、r(f(x)lis-list(y,var)return(lis)27v1.0可編輯可修改p3(10A4)p3(10A5)p3(10A6)p3(10A7)對模擬次數(shù) n調(diào)試了 4次，分別為n4,n5,n6,n7,得到精確值和模擬值。表 平均值法 的模擬 次數(shù)和模擬 值n104105106107模擬值萬差精確值為分層抽樣法R 程序f3-function(n,m)(r1-runif(n,min-0,maxr2-runif(m,min,max-1) z-function(u) exp(-u)/(1+uA2)28v1.0可編輯可修改j-1/2*mean(z(r1)+1/2*mean(z(r2)var-v

28、ar(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(10A4,2*10A4)得到精確值和模擬值表分層抽樣法的模擬次數(shù)和模擬值精確值n=10,m=20n=100,n=1000n=10000m=200M=2000M=20000偶變量程序d3-function(n)29v1.0可編輯可修改f-function(x) exp(-x)/(1+xA2)r-runif(n)m-mean(f(r)p-mean(f(1-r)j-(m+p)/2var-1/4*(var(f(r)+

29、var(f(1-r)+2*cov(f(r),f(1-r)lis-list(j,var)return(lis)d3(10A4)d3(10A5)d3(10A6)d3(10A7)對模擬次數(shù)n調(diào)試了 4次，分別為n4,n5,n6,n7,得到精確值和模擬值。表對偶法的模擬次數(shù)和模擬值n104105106107模擬值萬差30v1.0可編輯可修改精確值為控制變量法R程序k3-function(n)(f-function(x) exp(-x)/(1+xA2)r-runif(n)g-function(x) 2*(1-x)u-mean(g(r)l-(-cov(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(10A4)k1(10A5)k1(10A6)k1(10A7)31v1.0可編輯可修改對模擬次數(shù)n調(diào)試了 4次，分別為n4,n5,n6,n7 ,得到精確值和模擬值 表