第5講 蒙特卡洛方法的應(yīng)用.ppt

1、實(shí)驗(yàn)?zāi)康?實(shí)驗(yàn)內(nèi)容,學(xué)習(xí)如何應(yīng)用蒙特卡洛方法解決實(shí)際問題,1、起源和發(fā)展 2、原理 3、計(jì)算機(jī)模擬應(yīng)用實(shí)例 4、實(shí)驗(yàn)作業(yè),蒙特卡洛方法的應(yīng)用,一、MC 的起源和發(fā)展,隨機(jī)模擬方法，也稱為Monte Carlo方法，是一種基于“隨機(jī)數(shù)”的計(jì)算方法。這一方法源于美國(guó)在第一次世界大戰(zhàn)進(jìn)行的研制原子彈的“曼哈頓計(jì)劃”。該計(jì)劃的主持人之一、數(shù)學(xué)家馮諾伊曼用馳名世界的賭城摩納哥的Monte Carlo來命名這種方法，為它蒙上了一層神秘色彩。馮諾伊曼是公理化方法和計(jì)算機(jī)體系的領(lǐng)袖人物，Monte Carlo方法也是他的功勞。,事實(shí)上，Monte Carlo方法的基本思想很早以前就被人們所發(fā)現(xiàn)和利用。早在17

2、世紀(jì)，人們就知道用事件發(fā)生的“頻率”來決定事件的“概率”。18世紀(jì)下半葉的法國(guó)學(xué)者Buffon提出用投針試驗(yàn)的方法來確定圓周率的值。這個(gè)著名的Buffon試驗(yàn)是Monte Carlo方法的最早的嘗試！,歷史上曾有幾位學(xué)者相繼做過這樣的試驗(yàn)。不過呢，他們的試驗(yàn)是費(fèi)時(shí)費(fèi)力的，同時(shí)精度不夠高，實(shí)施起來也很困難。然而，隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展，人們不需要具體實(shí)施這些試驗(yàn)，而只要在計(jì)算機(jī)上進(jìn)行大量的、快速的模擬試驗(yàn)就可以了。,Monte Carlo方法是現(xiàn)代計(jì)算技術(shù)的最為杰出的成果之一，它在工程領(lǐng)域的作用是不可比擬的。,Buffon試驗(yàn),假設(shè)平面上有無數(shù)條距離為1的等距平行線，現(xiàn)向該平面隨機(jī)投擲一根長(zhǎng)度

3、為 的針( )， 則我們可計(jì)算該針與任一平行線相交的概率。這里，隨機(jī)投針指的是：針的中心點(diǎn)與最近的平行線間的距離 均勻的分布在區(qū)間 上，針與平行線的夾角 （不管相交與否）均勻的分布在區(qū)間 上。 因此，針與線相交的充要條件是,Buffon試驗(yàn),從而針線相交的概率為 根據(jù)上式，若我們做大量的投針試驗(yàn)并記錄針與線相交的次數(shù)，則由大數(shù)定理可以估計(jì)出針線相交的概率 ，從而得到 的估計(jì)值。 針與線的位置關(guān)系：,function piguji=buffon(llength,mm) %llength 是針的長(zhǎng)度 %mm 是隨機(jī)實(shí)驗(yàn)次數(shù) frq=0; xrandnum = unifrnd(0,0.5,1,mm)

4、; phi= unifrnd(0,pi,1,mm); for ii=1:mm if (xrandnum(1,ii)=(llength*sin(phi(1,ii)/2) frq=frq+1; end end piguji=2*llength/(frq/mm), buffon(.6,1000) piguji = 3.1662 buffon(.6,10000) piguji = 3.1072 buffon(.6,100000) piguji = 3.1522 buffon(.6,1000000) piguji = 3.1386 buffon(.6,1000000) piguji = 3.1451 b

5、uffon(.6,1000000) piguji = 3.1418 buffon(.6,1000000) piguji = 3.1448 buffon(.6,1000000) piguji = 3.1405 buffon(.6,1000000) piguji = 3.1394,二、MC 的原理,應(yīng)用Monte Carlo方法求解工程技術(shù)問題可以分為兩類： 確定性問題 隨機(jī)性問題,思路,1、 針對(duì)實(shí)際問題建立一個(gè)簡(jiǎn)單且便于實(shí)現(xiàn)的概率統(tǒng)計(jì)模型，使問題的解對(duì)應(yīng)于該模型中隨機(jī)變量的概率分布或其某些數(shù)字特征，比如，均值和方差等。所構(gòu)造的模型在主要特征參量方面要與實(shí)際問題或系統(tǒng)相一致的。 2、 根據(jù)模型中

6、各個(gè)隨機(jī)變量的分布，在計(jì)算機(jī)上產(chǎn)生隨機(jī)數(shù)，實(shí)現(xiàn)一次模擬過程所需的足夠數(shù)量的隨機(jī)數(shù)。通常先產(chǎn)生均勻分布的隨機(jī)數(shù)，然后生成服從某一分布的隨機(jī)數(shù)，再進(jìn)行隨機(jī)模擬試驗(yàn)。,3、 根據(jù)概率模型的特點(diǎn)和隨機(jī)變量的分布特性，設(shè)計(jì)和選取合適的抽樣方法，并對(duì)每個(gè)隨機(jī)變量進(jìn)行抽樣（包括直接抽樣、分層抽樣、相關(guān)抽樣、重要抽樣等）。 4、 按照所建立的模型進(jìn)行仿真試驗(yàn)、計(jì)算，求出問題的隨機(jī)解。 5、 統(tǒng)計(jì)分析模擬試驗(yàn)結(jié)果，給出問題的估計(jì)以及其精度估計(jì)。 6、 必要時(shí)，還應(yīng)改進(jìn)模型以降低估計(jì)方差和減少試驗(yàn)費(fèi)用，提高模擬計(jì)算的效率。,收斂性: 由大數(shù)定律, Monte-Carlo模擬的收斂是以概率而言的 誤差: 用頻率估計(jì)

7、概率時(shí)誤差的估計(jì)，可由中心極限定理，給定置信水平 的條件下，有： 模擬次數(shù):由誤差公式得,三、MC的應(yīng)用舉例,1、定積分的MC計(jì)算 隨機(jī)投點(diǎn)法 樣本平均值法 幾種降低估計(jì)方差的MC方法 2、 系統(tǒng)的可靠性數(shù)值模擬計(jì)算問題,1、定積分的MC計(jì)算,事實(shí)上，不少的統(tǒng)計(jì)問題，如計(jì)算概率、各階距等，最后都?xì)w結(jié)為定積分的近似計(jì)算問題。 下面考慮一個(gè)簡(jiǎn)單的定積分 為了說明問題，我們首先介紹兩種求 的簡(jiǎn)單的MC方法，然后給出幾種較為復(fù)雜而更有效的MC方法。,在計(jì)算積分上，MC的實(shí)用場(chǎng)合是計(jì)算重積分 其中 是 維空間的點(diǎn)，當(dāng) 較大時(shí)，用MC方法比一般的數(shù)值方法有優(yōu)點(diǎn)，主要是它的誤差與維數(shù) 無關(guān)。,隨機(jī)投點(diǎn)法,方

8、法簡(jiǎn)述： 設(shè) ，有限， ， ，并設(shè) 是在 上均勻分布的二維隨機(jī)變量，其聯(lián)合密度函數(shù)為 。則易見 是 中 曲線下方的面積。假設(shè)我們向 中進(jìn)行隨機(jī)投點(diǎn)，若點(diǎn)落在 下方，(即 稱為中的，否則稱為不中，則點(diǎn)中的概率為 。若我們進(jìn)行了 次投點(diǎn)，其中 次中的，則用頻率來估計(jì)概率 。即 。,那么我們可以得到 的一個(gè)估計(jì) 具體試驗(yàn)步驟為,求解定積分的算例,例 計(jì)算定積分 事實(shí)上，其精確解為 用隨機(jī)投點(diǎn)法求解： 注 增加樣本數(shù)目，可提高計(jì)算精度，但計(jì)算時(shí)間也會(huì)提高。 sj(0,4,4,1000000) result = 7.2336,function result=sj(a,b,m,mm) fq=0; x= u

9、nifrnd(a,b,1,mm); y = unifrnd(0,m,1,mm); for ii=1:mm if (cos(x(1,ii)+2=y(1,ii) fq=fq+1; end end result=fq*m*(b-a)/mm,注1 隨機(jī)投點(diǎn)法的思想簡(jiǎn)單明了，且每次投點(diǎn)結(jié)果服從二項(xiàng)分布，故 ，其中 注2 可證 是 的無偏估計(jì)。若用估計(jì)的標(biāo)準(zhǔn)差來衡量其精度，則估計(jì) 的精度的階為 。 注3 這里，定積分的解，就對(duì)應(yīng)我們選定的隨機(jī)變量的概率值。,2020/7/22,例 的計(jì)算,1.單位圓的面積等于 2. 3.,用隨機(jī)投點(diǎn)法求的值, sjtdf_pi1(1000) piguji = 3.0520

10、 sjtdf_pi1(10000) piguji = 3.1204 sjtdf_pi1(100000) piguji = 3.1296,function piguji=sjtdf_pi1(mm) %mm 是隨機(jī)實(shí)驗(yàn)次數(shù) frq=0;xrandnum = unifrnd(0,1,1,mm); yrandnum = unifrnd(0,1,1,mm); for ii=1:mm if xrandnum(1,ii)2+yrandnum(1,ii)2=1 frq=frq+1; end end piguji=4*frq/mm,sjtdf_pi2(100) piguji = 3.2000 sjtdf_pi2

11、(1000) piguji = 3.2120 sjtdf_pi2(10000) piguji = 3.1260 sjtdf_pi2(100000) piguji = 3.1373,function piguji=sjtdf_pi2(mm) %mm 是隨機(jī)實(shí)驗(yàn)次數(shù) frq=0; xrandnum = unifrnd(0,1,1,mm); yrandnum = unifrnd(0,1,1,mm); for ii=1:mm if (sqrt(1-xrandnum(1,ii)2)=yrandnum(1,ii) frq=frq+1; end End,piguji=4*frq/mm,樣本平均值法,基本原理

12、：對(duì)積分 ，設(shè) 是 上的一個(gè)密度函數(shù)，改寫 可見，任一積分均可以表示為某個(gè)隨機(jī)變量（函數(shù)）的期望。由矩法，若有 個(gè)來自 的觀測(cè)值，則可給出 的一個(gè)矩估計(jì)。 最簡(jiǎn)單的，若 ，有限，可取 。,設(shè) 是來自 的隨機(jī)數(shù)，則 的一個(gè)估計(jì)為 具體步驟為 注 可證 是 的無偏估計(jì)。一般而言，樣本均值法要比隨機(jī)投點(diǎn)法更有效。,求解定積分的算例,例 計(jì)算定積分 事實(shí)上，其精確解為 樣本平均值法求解： 注 增加樣本數(shù)目，可提高計(jì)算精度，但計(jì)算時(shí)間也會(huì)提高。 ybjzf(0,4,4,1000) result = 7.3036 ybjzf(0,4,4,10000) result = 7.2970 ybjzf(0,4,4

13、,100000) result = 7.2578,function result=ybjzf(a,b,m,mm) %a是積分的下限 %b是積分的上限 %積分函數(shù)cos(x)+2 %mm 是隨機(jī)實(shí)驗(yàn)次數(shù) sum=0; xrandnum = unifrnd(a,b,1,mm); for ii=1:mm sum=sum+cos(xrandnum(1,ii)+2; end result=sum*(b-a)/mm,2020/7/22,例 的計(jì)算,1.單位圓的面積等于 2. 3.,用樣本平均值法求的值,function result=ybjzf1(a,b,m,mm) %a是積分的下限 %b是積分的上限 %

14、積分函數(shù) %mm 是隨機(jī)實(shí)驗(yàn)次數(shù) xrandnum = unifrnd(a,b,1,mm); sum=sum(sqrt(1-xrandnum.2); result=sum*(b-a)/mm; result=result*4,ybjzf1(0,1,1,100) result = 3.08745746887753 ybjzf1(0,1,1,1000) result = 3.15500646827616 ybjzf1(0,1,1,10000) result = 3.11911714237706 ybjzf1(0,1,1,100000) result = 3.14014654862983 ybjzf1

15、(0,1,1,1000000) result = 3.14093979612119,function result=ybjzf1(a,b,m,mm) %a是積分的下限 %b是積分的上限 %積分函數(shù) %mm 是隨機(jī)實(shí)驗(yàn)次數(shù) sum=0; xrandnum = unifrnd(a,b,1,mm); for ii=1:mm sum=sum+1/(1+xrandnum(1,ii)2); end result=sum*(b-a)/mm; result=result*4,ybjzf2(0,1,1,100) result = 3.04500162146030 ybjzf2(0,1,1,1000) resul

16、t = 3.14857120401090 ybjzf2(0,1,1,10000) result = 3.14530178564491 ybjzf2(0,1,1,100000) result = 3.14069143292954,幾種降低估計(jì)方差的MC方法,重要抽樣法 特點(diǎn)：相對(duì)樣本均值法而言，樣本均值法是由于假設(shè) 是均勻分布的概率密度，故采用的是均勻抽樣，各隨機(jī)數(shù) 是均勻分布的隨機(jī)數(shù)，各 對(duì) 的貢獻(xiàn)是不同， 大則貢獻(xiàn)大，但在抽樣時(shí)，這種差別未能體現(xiàn)出來。 而重要抽樣法，則希望貢獻(xiàn)率大的隨機(jī)數(shù)出現(xiàn)的概率大，貢獻(xiàn)小的隨機(jī)數(shù)出現(xiàn)概率小，從而提高抽樣的效。,幾種降低估計(jì)方差的MC方法,重要抽樣法 關(guān)鍵

17、因素在于 的選取，使得估計(jì)的方差較小。 重要抽樣法的基本思想，就是通過選取與 形狀接近的密度函數(shù) 來降低估計(jì)的方差。,幾種降低估計(jì)方差的MC方法,分層抽樣法 同樣是利用貢獻(xiàn)率大小來降低估計(jì)方差的方法。它首先是把樣本空間 分成一些小區(qū)間 ，且諸 不交， ，然后在各個(gè)小區(qū)間內(nèi)的抽樣數(shù)由其貢獻(xiàn)大小決定。對(duì) 貢獻(xiàn)大的 抽樣多，可提高抽樣效率。如果能夠提出較好抽樣區(qū)間的分配和各子區(qū)間內(nèi)抽樣次數(shù)的分配方案，分層抽樣法估計(jì)積分可以達(dá)到非常令人滿意的效果。,幾種降低估計(jì)方差的MC方法,關(guān)聯(lián)抽樣法 將需要估計(jì)的積分分解成兩個(gè)積分之差， 對(duì) 的估計(jì)轉(zhuǎn)化為對(duì) ， 的估計(jì)的差。則相應(yīng)的，其估計(jì)的方差的大小則與 ， 的

18、估計(jì)的正相關(guān)度有關(guān)，若兩者的相關(guān)程度越高，則 的估計(jì)方差越小。這便是關(guān)聯(lián)抽樣法的基本出發(fā)點(diǎn)。,一個(gè)元件(或系統(tǒng))能正常工作的概率稱為 元件(或系統(tǒng))的可靠性,系統(tǒng)由元件組成,常見的元件連接方式：,串聯(lián),并聯(lián),2、系統(tǒng)的可靠性計(jì)算問題,例 設(shè)兩系統(tǒng)都是由 4 個(gè)元件組成,每個(gè)元件 正常工作的概率為 p=0.5 , 每個(gè)元件是否正常工作相互獨(dú)立.兩系統(tǒng)的連接方式如下圖所示，比較兩系統(tǒng)的可靠性.,S1:,function Rguji=litiR01(0.5,0.5,0.5,0.5,mm) frq=0; randnuma1 = binornd(1,0.5,1,mm); randnuma2 = bino

19、rnd(1,0.5,1,mm); randnumb1 = binornd(1,0.5,1,mm) randnumb2 = binornd(1,0.5,1,mm); Rguji=frq/mm,function Rguji=litiR01(0.5,0.5,0.5,0.5,mm) frq=0;randnuma1 = binornd(1,0.5,1,mm); randnuma2 = binornd(1,0.5,1,mm); randnumb1 = binornd(1,0.5,1,mm) randnumb2 = binornd(1,0.5,1,mm); for ii=1:mm if (randnuma1

20、(1,ii)=1) end End,Rguji=frq/mm,例 設(shè)兩系統(tǒng)都是由 4 個(gè)元件組成,每個(gè)元件 正常工作的概率為 p , 每個(gè)元件是否正常工 作相互獨(dú)立.兩系統(tǒng)的連接方式如下圖所示， 比較兩系統(tǒng)的可靠性.,S1:,S2:,例 設(shè)兩系統(tǒng)都是由 4 個(gè)元件組成,每個(gè)元件 的壽命服從參數(shù)為的指數(shù)分布，每個(gè)元件 是否正常工作相互獨(dú)立.兩系統(tǒng)的連接方式如 下圖所示，求兩系統(tǒng)壽命大于T=100的概率.,S1:,例 設(shè)兩系統(tǒng)都是由 4 個(gè)元件組成,每個(gè)元件 的壽命服從參數(shù)為的指數(shù)分布，每個(gè)元件 是否正常工作相互獨(dú)立.兩系統(tǒng)的連接方式如 下圖所示，求兩系統(tǒng)壽命大于T=100的概率.,S1:,S2:

21、,function Rguji=litiR1(t,thetaa1,thetaa2,thetab1,thetab2,mm) %t 是要求系統(tǒng)生存的壽命%thetaa1 是元件A1的數(shù)學(xué)期望%thetaa2 是元件A2的數(shù)學(xué)期望%thetab1 是元件B1的數(shù)學(xué)期望 %thetab2 是元件B2的數(shù)學(xué)期望%mm 是隨機(jī)實(shí)驗(yàn)次數(shù) frq=0;randnuma1 = exprnd(thetaa1,1,mm); randnuma2 = exprnd(thetaa2,1,mm);randnumb1 = exprnd(thetab1,1,mm);randnumb2 = exprnd(thetab2,1,mm); for ii=1:mm if (randnuma1(1,ii)t) end End,Rguji=frq/mm,function Rguji=litiR2(t,thetaa1,thetaa2,thetab1,thetab2,mm) %t 是要求系統(tǒng)生存的壽命%thetaa1 是元件A1的數(shù)學(xué)期望%thetaa2 是元件A2的數(shù)學(xué)期望%thetab1 是元件B1的數(shù)學(xué)期望 %thetab2 是元件B2的數(shù)學(xué)期望%mm