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

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

2、，當(dāng)所求解問題是*種隨機(jī)事件出現(xiàn)的概率，或者是*個(gè)隨機(jī)變量的期望值時(shí)，通過*種實(shí)驗(yàn)的方法，以這種事件出現(xiàn)的頻率估計(jì)這一隨機(jī)事件的概率，或者得到這個(gè)隨機(jī)變量的*些數(shù)字特征，并將其作為問題的解。蒙特卡洛方法解題過程的三個(gè)主要步驟：1構(gòu)造或描述概率過程對(duì)于本身就具有隨機(jī)性質(zhì)的問題，如粒子輸運(yùn)問題，主要是正確描述和模擬這個(gè)概率過程，對(duì)于本來不是隨機(jī)性質(zhì)確實(shí)定性問題，比方計(jì)算定積分，就必須事先構(gòu)造一個(gè)人為的概率過程，它的*些參量正好是所要求問題的解。即要將不具有隨機(jī)性質(zhì)的問題轉(zhuǎn)化為隨機(jī)性質(zhì)的問題。2實(shí)現(xiàn)從概率分布抽樣構(gòu)造了概率模型以后，由于各種概率模型都可以看作是由各種各樣的概率分布構(gòu)成的，因此產(chǎn)生概率

3、分布的隨機(jī)變量或隨機(jī)向量，就成為實(shí)現(xiàn)蒙特卡洛方法模擬實(shí)驗(yàn)的根本手段，這也是蒙特卡洛方法被稱為隨機(jī)抽樣的原因。最簡單、最根本、最重要的一個(gè)概率分布是0，1上的均勻分布或稱矩形分布。隨機(jī)數(shù)就是具有這種均勻分布的隨機(jī)變量。隨機(jī)數(shù)序列就是具有這種分布的總體的一個(gè)簡單子樣，也就是一個(gè)具有這種分布的相互獨(dú)立的隨機(jī)變數(shù)序列。產(chǎn)生隨機(jī)數(shù)的問題，就是從這個(gè)分布的抽樣問題。在計(jì)算機(jī)上，可以用物理方法產(chǎn)生隨機(jī)數(shù)，但價(jià)格昂貴，不能重復(fù)，使用不便。另一種方法是用數(shù)學(xué)遞推公式產(chǎn)生。這樣產(chǎn)生的序列，與真正的隨機(jī)數(shù)序列不同，所以稱為偽隨機(jī)數(shù)，或偽隨機(jī)數(shù)序列。不過，經(jīng)過多種統(tǒng)計(jì)檢驗(yàn)說明，它與真正的隨機(jī)數(shù)，或隨機(jī)數(shù)序列具有相近的

4、性質(zhì)，因此可把它作為真正的隨機(jī)數(shù)來使用。由分布隨機(jī)抽樣有各種方法，與從(0,1)上均勻分布抽樣不同，這些方法都是借助于隨機(jī)序列來實(shí)現(xiàn)的，也就是說，都是以產(chǎn)生隨機(jī)數(shù)為前提的。由此可見，隨機(jī)數(shù)是我們實(shí)現(xiàn)蒙特卡洛模擬的根本工具。3建立各種估計(jì)量一般說來，構(gòu)造了概率模型并能從中抽樣后，即實(shí)現(xiàn)模擬實(shí)驗(yàn)后，我們就要確定一個(gè)隨機(jī)變量，作為所要求的問題的解，我們稱它為無偏估計(jì)。建立各種估計(jì)量，相當(dāng)于對(duì)模擬實(shí)驗(yàn)的結(jié)果進(jìn)展考察和登記，從中得到問題的解。2方法介紹2.1隨機(jī)投點(diǎn)法隨機(jī)投點(diǎn)法是進(jìn)展n次試驗(yàn)，當(dāng)n充分大的時(shí)候，以隨機(jī)變量k/n作為期望值E(*)的近似估計(jì)值，即其中k是n次實(shí)驗(yàn)中成功的次數(shù)。假設(shè)一次投點(diǎn)試

5、驗(yàn)的成功概率為p，并以則一次試驗(yàn)成功的均值與方差為假設(shè)進(jìn)展n次試驗(yàn)，其中k次試驗(yàn)成功，則k為具有參數(shù)為n,p的二項(xiàng)分布，此時(shí)，隨機(jī)變量k的估計(jì)為顯然，隨機(jī)變量的均值和方差滿足dd設(shè)計(jì)算的定積分為，其中a,b為有限數(shù)，被積函數(shù)f(*)是連續(xù)隨機(jī)變量的概率密度函數(shù)，因此f(*)滿足如下條件:顯然I是一個(gè)概率積分，其積分值等于概率。下面按給定分布f(*)隨機(jī)投點(diǎn)的方法，給出如下Monte Carlo近似求積算法：(1)產(chǎn)生服從給定分布的隨機(jī)變量值,i=1,2,N;(2)檢查是否落入積分區(qū)間。如果條件滿足，則記錄落入積分區(qū)間一次。假設(shè)在N次實(shí)驗(yàn)以后，落入積分區(qū)間的總次數(shù)為n，則用作為概率積分的近似值，

6、即2.2平均值法任取一組相互獨(dú)立、同分布的隨機(jī)變量，在a,b服從分布率p(*),令,則也是一組相互獨(dú)立、同分布的隨機(jī)變量，而且由強(qiáng)大數(shù)定理假設(shè)記,則依概率1收斂到I，平均值法就是用作為I的近似值。假設(shè)所需計(jì)算積分為，其中被積函數(shù)在a,b可積，任意選擇一個(gè)有簡單方法可以進(jìn)展抽樣的概率密度函數(shù)p(*)，使其滿足條件：記則所求積分為因而Monte Carlo近似求積算法為：(1) 產(chǎn)生服從分布率p(*)的隨機(jī)數(shù)(2) 計(jì)算均值,即有2.3重要性采樣法從數(shù)學(xué)角度上看定積分可以看成其中g(shù)(*)是*個(gè)隨機(jī)變量*的密度函數(shù)，因此積分值I可看成隨機(jī)變量Z=f(*)/g(*)的數(shù)學(xué)期望值為了減少模擬實(shí)驗(yàn)的方差應(yīng)

7、適中選取g(*),使VarI盡可能小，如果被積函數(shù)f(*)0,可取g(*)=cf(*),當(dāng)c=1/I時(shí)就有Var(I)=0.一般應(yīng)選取和f(*)相似的密度函數(shù)g(*)，使f(*)/g(*)接近于常數(shù)，故而Var(I)接近于0，以到達(dá)降低模擬實(shí)驗(yàn)的方差，這種減少方差的模擬試驗(yàn)法為重要抽樣法。2.4分層抽樣法分層抽樣法是利用奉獻(xiàn)率大小來降低估計(jì)方差的方法。它首先把樣本空間D分成一些不交的小區(qū)間，然后在各個(gè)小區(qū)間的抽樣數(shù)由其奉獻(xiàn)大小決定。即，定義，則的抽樣數(shù)應(yīng)與成正比。考慮積分將0,1分成m個(gè)小區(qū)間：則記為第i個(gè)小區(qū)間的長度，i=1,2,.,m，在每個(gè)小區(qū)間上的積分值可用均值法估計(jì)出來，然后將其相加

8、即可給出的一個(gè)估計(jì)。具體步驟為：獨(dú)立產(chǎn)生U(0，1)隨機(jī)數(shù)計(jì)算計(jì)算于是的估計(jì)為，其方差為其中，2.5對(duì)偶變量法控制變量法利用數(shù)學(xué)上積分運(yùn)算的線性特性:選擇函數(shù)g(*)時(shí)要考慮到：g(*)在整個(gè)積分區(qū)間都是容易準(zhǔn)確算出,并且在上式右邊第一項(xiàng)的運(yùn)算中對(duì)(f-g)積分的方差應(yīng)當(dāng)要比第二項(xiàng)對(duì)f積分的方差小。在應(yīng)用這種方法時(shí)，在重要抽樣法中所遇到的，當(dāng)g(*)趨于零時(shí)，被積函數(shù)(f-g)趨于無窮大的困難就不再存在，因而計(jì)算出的結(jié)果穩(wěn)定性比擬好。該方法也不需要從分布密度函數(shù)g(*)，解析求出分布函數(shù)G(*)。由此我們可以看出選擇g(*)所受到的限制比重要抽樣法要小些。模擬過程：獨(dú)立產(chǎn)生U(0，1)隨機(jī)數(shù)計(jì)

9、算計(jì)算2.6 控制變量法通常在蒙特卡洛計(jì)算中采用互相獨(dú)立的隨機(jī)點(diǎn)來進(jìn)展計(jì)算。對(duì)偶變量法中卻使用相關(guān)聯(lián)的點(diǎn)來進(jìn)展計(jì)算。它利用相關(guān)點(diǎn)間的關(guān)系可以是正關(guān)聯(lián)的，也可以是負(fù)關(guān)聯(lián)的這個(gè)特點(diǎn)。兩個(gè)函數(shù)值和之和的方差為如果我們選擇一些點(diǎn)，它們使和是負(fù)關(guān)聯(lián)的。這樣就可以使上式所示的方差減小。當(dāng)然這需要對(duì)具體的函數(shù)和有充分的了解獨(dú)立產(chǎn)生U(0，1)隨機(jī)數(shù)計(jì)算，找g(*),f(*)是相關(guān)的，且Eg(*)=計(jì)算3程序及實(shí)現(xiàn)結(jié)果3.1 的求解3.1.1 隨機(jī)投點(diǎn)法先利用R 軟件產(chǎn)生服從0,1上均勻分布的隨機(jī)數(shù) *,Y, ，計(jì)算的個(gè)數(shù)，即事件發(fā)生的頻數(shù)，求出頻率，即為積分的近似值。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)對(duì)模擬次數(shù)n調(diào)試了4次，分別為，得到準(zhǔn)確值和模擬值。表3.1.1 隨機(jī)投點(diǎn)法的模擬次數(shù)和模擬值n模擬值0.62690.632350.6325030.6319298方差2.32599e-052.32414e-062.32713e-072.32477e-08準(zhǔn)確值為0.6321206 3.1.2 平均值法先用R 軟件產(chǎn)生n個(gè)服

11、從0,1上均勻分布的隨機(jī)數(shù)，計(jì)算,再計(jì)算的平均值，即為定積分的近似值 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)對(duì)模擬次數(shù)n調(diào)試了4次，分別為，得到準(zhǔn)確值和模擬值。表3.1.2 平均值法的模擬次數(shù)和模擬值n模擬值0.63101170.63226430.63181990.632079方差3.25430e-063.27162e-073.2805e-083.275

12、36e-09準(zhǔn)確值為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)對(duì)模擬次數(shù)n調(diào)試了4次，分別為，得到準(zhǔn)確值和模擬值。表3.1.3 重要性抽樣法法的模擬次數(shù)和模擬值n模擬值0.61348190.61348190.61348190.613

13、4819方差7.06957e-067.06957e-077.06957e-087.06957e-09準(zhǔn)確值為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)得到準(zhǔn)確值和模

14、擬值表3.1.4 分層抽樣法的模擬次數(shù)和模擬值取值n=10,m=20n=100,m=200n=1000M=2000n=10000M=20000模擬值0.65827020.62917550.63417180.6327054方差0.0002692072.75023e-053.23184e-063.18334e-07準(zhǔn)確值為0.6321206 3.1.5 對(duì)偶變量法先用 R 軟件產(chǎn)生服從0,1上均勻分布的隨機(jī)數(shù) * ，函數(shù)f(*),計(jì)算，計(jì)算 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) 對(duì)模擬次數(shù)n調(diào)試了4次，分別為，得到準(zhǔn)確值和模擬值。表3.1.5 對(duì)偶變量法的模擬次數(shù)和模擬值n模擬值0.63219790.63206590.63209570.6321183方差0.000524720.0005327480.00053057410.0005291691準(zhǔn)確值為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)對(duì)模擬次數(shù)n調(diào)試了4次，分別為，得到準(zhǔn)確值和模擬值。表3.1.6 控制變量法的模擬次數(shù)和模擬值n模擬值0.6345750.63208840.6318490.6321651方差5.713907e-055.72338e-06

17、5.736774e-075.731522e-08準(zhǔn)確值為0.6321206 3.2 對(duì)積分求解3.2.1 隨機(jī)投點(diǎn)法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)對(duì)模擬次數(shù)n調(diào)試了4次，分別

18、為，得到準(zhǔn)確值和模擬值。表3.2.1 隨機(jī)投點(diǎn)法的模擬次數(shù)和模擬值n模擬值0.18688450.18688450.18688450.1868845方差2.33071e-052.32230e-062.32479e-072.32611e-08準(zhǔn)確值為 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)對(duì)模擬次數(shù)n調(diào)試了4次，分別為，得到準(zhǔn)確值和模擬值。表3.2.2 平均值法的模擬次數(shù)和模擬值n模擬值 0.11733330.11695820.1170038 0.1170311方差1.30285e-051.30285e-061.30285e-071.30285e-08準(zhǔn)確值為 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)對(duì)模擬次數(shù)n調(diào)試了4次，分別為，得到準(zhǔn)確值和模擬值。表3.2.3 分層抽樣法的模擬次數(shù)和模擬值取值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準(zhǔn)確值為 0.1170196

21、 3.2.4 對(duì)偶變量法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)對(duì)模擬次數(shù)n調(diào)試了4次，分別為，得到準(zhǔn)確值和模擬值。表3.2.4 對(duì)偶變量法的模擬次數(shù)和模擬值n模擬值 0.11690430.1167410.116943 0.1170

22、233準(zhǔn)確值為 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)對(duì)模擬次數(shù)n調(diào)試了4次，分別為，得到準(zhǔn)確值和模擬值。表3.2.5 控制變量法的模擬次數(shù)和模擬值n模擬值 0.11747130.11711510.117031610.1170211方差5.68197e-055.73480e-065.73725e-075.73387e-08準(zhǔn)確值為 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)對(duì)模擬次數(shù)n調(diào)試了4次，分別為，得到準(zhǔn)確值和模擬值。表3.2.6 重要性采樣法的模擬次數(shù)和模擬值n模擬值 0.11735220.11702670.11701140.1170123方差8.17541e-078.16598e-088.17952e-098.17912e-10準(zhǔn)確值為 0.11701963.3 積分求解3

25、.3.1 隨機(jī)投點(diǎn)法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)對(duì)模擬次數(shù)n調(diào)試了4次，分別為，得到準(zhǔn)確值和模擬值。表3.3.1 隨機(jī)投點(diǎn)法的模擬次數(shù)和模擬值n模擬值 0.53080.52507 0.5253510.5247777方差2.49625e-052.49312e-062.49423e-072.49371

26、e-08準(zhǔn)確值為 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)對(duì)模擬次數(shù)n調(diào)試了4次，分別為，得到準(zhǔn)確值和模擬值。表3.3.2 平均值法的模擬次數(shù)和模擬值n模擬值0.52397070.5252655 0.5245715 0.5247576方差5.97263e-066.02206e-075.99868e-085.9992

27、6e-09準(zhǔn)確值為 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)得到準(zhǔn)確值和模擬值表3.3.3 分層抽樣法的模擬次數(shù)和模擬值取值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準(zhǔn)確值為0.5247971 3.3.4 對(duì)偶變量法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)對(duì)模擬次數(shù)n調(diào)試了4次，分別為，得到準(zhǔn)確值和模擬值。表3.3.4 對(duì)偶法的模擬次數(shù)和模擬值n模擬值0.63197880.6321444 0.6321236 0.6321169方差0.0005191260.00053013340.00052976580.0005295388準(zhǔn)確值為 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)對(duì)模擬次數(shù)n調(diào)試了4次，分別為，得到準(zhǔn)確值和模擬值表3.3.5 控制變量法的模擬次數(shù)和模擬值n模擬值0.6297403 0.632438 0.6320956 0.6321038方差5.80168e-055.74223e-065.73553e-075.73648e-08準(zhǔn)確值為 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)對(duì)模擬次數(shù)