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

2、nte Carlo方法也是他的功勞。 事實(shí)上，Monte Carlo方法的基本思想很早以前就被人們所發(fā)現(xiàn)和利用。早在17世紀(jì)，人們就知道用事件發(fā)生的“頻率”來(lái)決定事件的“概率”。18世紀(jì)下半葉的法國(guó)學(xué)者Buffon提出用投針試驗(yàn)的方法來(lái)確定圓周率的值。這個(gè)著名的Buffon試驗(yàn)是Monte Carlo方法的最早的嘗試！ 歷史上曾有幾位學(xué)者相繼做過(guò)這樣的試驗(yàn)。不過(guò)呢，他們的試驗(yàn)是費(fèi)時(shí)費(fèi)力的，同時(shí)精度不夠高，實(shí)施起來(lái)也很困難。然而，隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展，人們不需要具體實(shí)施這些試驗(yàn)，而只要在計(jì)算機(jī)上進(jìn)行大量的、快速的模擬試驗(yàn)就可以了。 Monte Carlo方法是現(xiàn)代計(jì)算技術(shù)的最為杰出的成果之一

3、，它在工程領(lǐng)域的作用是不可比擬的。 BuffonBuffon試驗(yàn)試驗(yàn) 假設(shè)平面上有無(wú)數(shù)條距離為1的等距平行線，現(xiàn)向該平面隨機(jī)投擲一根長(zhǎng)度為 的針( )， 則我們可計(jì)算該針與任一平行線相交的概率。這里，隨機(jī)投針指的是：針的中心點(diǎn)與最近的平行線間的距離 均勻的分布在區(qū)間 上，針與平行線的夾角 （不管相交與否）均勻的分布在區(qū)間 上。 因此，針與線相交的充要條件是 l1l21 , 0 x, 021sinxBuffonBuffon試驗(yàn)試驗(yàn) 從而針線相交的概率為 根據(jù)上式，若我們做大量的投針試驗(yàn)并記錄針與線相交的次數(shù)，則由大數(shù)定理可以估計(jì)出針線相交的概率 ，從而得到 的估計(jì)值。 針與線的位置關(guān)系：ldxd

4、lXPpl22sin20sin20 pfunction piguji=buffon(llength,mm)%llength 是針的長(zhǎng)度%mm 是隨機(jī)實(shí)驗(yàn)次數(shù)frq=0;xrandnum = unifrnd(0,0.5,1,mm);phi= unifrnd(0,pi,1,mm);for ii=1:mm if (xrandnum(1,ii) buffon(.6,1000) piguji = 3.1662 buffon(.6,10000) piguji = 3.1072 buffon(.6,100000) piguji = 3.1522 buffon(.6,1000000) piguji = 3.1

5、386 buffon(.6,1000000) piguji = 3.1451 buffon(.6,1000000) piguji = 3.1418 buffon(.6,1000000) piguji = 3.1448 buffon(.6,1000000) piguji = 3.1405 buffon(.6,1000000) piguji = 3.1394二、二、MC MC 的原理的原理 應(yīng)用Monte Carlo方法求解工程技術(shù)問(wèn)題可以分為兩類： 確定性問(wèn)題 隨機(jī)性問(wèn)題思路思路1、 針對(duì)實(shí)際問(wèn)題建立一個(gè)簡(jiǎn)單且便于實(shí)現(xiàn)的概率統(tǒng)計(jì)模型，使問(wèn)題的解對(duì)應(yīng)于該模型中隨機(jī)變量的概率分布或其某些數(shù)字特征，比

6、如，均值和方差等。所構(gòu)造的模型在主要特征參量方面要與實(shí)際問(wèn)題或系統(tǒng)相一致的。2、 根據(jù)模型中各個(gè)隨機(jī)變量的分布，在計(jì)算機(jī)上產(chǎn)生隨機(jī)數(shù)，實(shí)現(xiàn)一次模擬過(guò)程所需的足夠數(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ì)算，求出問(wèn)題的隨機(jī)解。5、 統(tǒng)計(jì)分析模擬試驗(yàn)結(jié)果，給出問(wèn)題的估計(jì)以及其精度估計(jì)。6、 必要時(shí)，還應(yīng)改進(jìn)模型以降低估計(jì)方差和減少試驗(yàn)費(fèi)用，提高模擬計(jì)算的效率。 收斂

7、性收斂性: : 由大數(shù)定律, Monte-Carlo模擬的收斂是以概率而言的 誤差誤差: : 用頻率估計(jì)概率時(shí)誤差的估計(jì)，可由中心極限定理，給定置信水平 的條件下，有： 模擬次數(shù)模擬次數(shù): :由誤差公式得NU2/1|)(XgVar22/1)(UN三、三、MCMC的應(yīng)用舉例的應(yīng)用舉例 1 1、定積分的、定積分的MCMC計(jì)算計(jì)算 隨機(jī)投點(diǎn)法隨機(jī)投點(diǎn)法 樣本平均值法樣本平均值法 幾種降低估計(jì)方差的幾種降低估計(jì)方差的MCMC方法方法 2 2、 系統(tǒng)的可靠性數(shù)值模擬計(jì)算問(wèn)題系統(tǒng)的可靠性數(shù)值模擬計(jì)算問(wèn)題1 1、定積分的、定積分的MCMC計(jì)算計(jì)算 事實(shí)上，不少的統(tǒng)計(jì)問(wèn)題，如計(jì)算概率、各階距等，最后都?xì)w結(jié)為定

8、積分的近似計(jì)算問(wèn)題。 下面考慮一個(gè)簡(jiǎn)單的定積分 為了說(shuō)明問(wèn)題，我們首先介紹兩種求 的簡(jiǎn)單的MC方法，然后給出幾種較為復(fù)雜而更有效的MC方法。 dxxfba 在計(jì)算積分上，MC的實(shí)用場(chǎng)合是計(jì)算重積分 其中 是 維空間的點(diǎn)，當(dāng) 較大時(shí)，用MC方法比一般的數(shù)值方法有優(yōu)點(diǎn)，主要是它的誤差與維數(shù) 無(wú)關(guān)。 dPPgIkPmmm隨機(jī)投點(diǎn)法隨機(jī)投點(diǎn)法 方法簡(jiǎn)述方法簡(jiǎn)述：設(shè) ，有限， ， ，并設(shè) 是在 上均勻分布的二維隨機(jī)變量，其聯(lián)合密度函數(shù)為 。則易見(jiàn) 是 中 曲線下方的面積。假設(shè)我們向 中進(jìn)行隨機(jī)投點(diǎn)，若點(diǎn)落在 下方，(即 稱為中的，否則稱為不中，則點(diǎn)中的概率為 。若我們進(jìn)行了 次投點(diǎn)，其中 次中的，則用頻

9、率來(lái)估計(jì)概率 。即 。ab Mxf0Mybxayx0 ,:,YX,MybxaIabM0 ,1 dxxfba xfy xfy xfy abMpn0npnnp0 那么我們可以得到 的一個(gè)估計(jì) 具體試驗(yàn)步驟為nnabM01 用上式來(lái)估計(jì)的個(gè)數(shù)統(tǒng)計(jì)，和，計(jì)算，隨機(jī)數(shù)，個(gè)獨(dú)立地產(chǎn)生0,.,11 , 02nyxfxfMvyabuaxnivuUniiiiiiiii求解定積分的算例求解定積分的算例例例 計(jì)算定積分事實(shí)上，其精確解為用隨機(jī)投點(diǎn)法求解：注注 增加樣本數(shù)目，可提高計(jì)算精度，但計(jì)算時(shí)間也會(huì)提高。sjtdf(0,4,4,1000000) result = 7.23367.2432 = 0 . 4sin0

10、 . 8dxx400 . 2cosfunction result=sjtdf(a,b,m,mm)%a是積分的下限%b是積分的上限%m是函數(shù)的上界%mm 是隨機(jī)實(shí)驗(yàn)次數(shù)frq=0;xrandnum = unifrnd(a,b,1,mm);yrandnum = unifrnd(0,m,1,mm);for ii=1:mm if (cos(xrandnum(1,ii)+2=yrandnum(1,ii) frq=frq+1; end end result=frq*m*(b-a)/mm 注注1 1 隨機(jī)投點(diǎn)法的思想簡(jiǎn)單明了，且每次投點(diǎn)結(jié)果服從二項(xiàng)分布，故 ，其中 注注2 2 可證 是 的無(wú)偏估計(jì)。若用估計(jì)

11、的標(biāo)準(zhǔn)差來(lái)衡量其精度，則估計(jì) 的精度的階為 。 注注3 3 這里，定積分的解，就對(duì)應(yīng)我們選定的隨機(jī)變量的概率值。 pnbn,0abMp1121n2022-5-27例例 的計(jì)算的計(jì)算1.單位圓的面積等于 2.3.0.20.40.60.810.20.40.60.8141102dxx4)1 (1102dxx用隨機(jī)投點(diǎn)法求用隨機(jī)投點(diǎn)法求 的值的值 sjtdf_pi1(1000) piguji = 3.0520 sjtdf_pi1(10000) piguji = 3.1204 sjtdf_pi1(100000) piguji = 3.1296function piguji=sjtdf_pi1(mm)%m

12、m 是隨機(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 sjtdf_pi2(1000) piguji = 3.2120 sjtdf_pi2(10000) piguji = 3.1260 sjtdf_pi2(100000) piguji = 3.1373function piguji=sjtdf_pi2(mm)%mm 是隨機(jī)實(shí)驗(yàn)次數(shù)frq=0;xrandnum = unifrnd(0,1,1,mm);yrandnu

13、m = 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樣本平均值法樣本平均值法 基本原理基本原理：對(duì)積分 ，設(shè) 是 上的一個(gè)密度函數(shù)，改寫 可見(jiàn)，任一積分均可以表示為某個(gè)隨機(jī)變量（函數(shù)）的期望。由矩法，若有 個(gè)來(lái)自 的觀測(cè)值，則可給出 的一個(gè)矩估計(jì)。 最簡(jiǎn)單的，若 ，有限，可取 。 dxxfba xgba, XgXfEdxxgxgxfban xgab abxg /1 設(shè) 是來(lái)自 的隨機(jī)數(shù)，則 的一個(gè)估計(jì)為 具體步驟為 注注 可證

14、 是 的無(wú)偏估計(jì)。一般而言，樣本均值法要比隨機(jī)投點(diǎn)法更有效。 nxx ,.,1baU, niiniiixfnabxgxfn1121 用上式來(lái)估計(jì)，和計(jì)算，隨機(jī)數(shù)個(gè)獨(dú)立地產(chǎn)生nixfuabaxuuUniiin,.,1,.,1 , 012求解定積分的算例求解定積分的算例例例 計(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,100000) result = 7.2578dxx400 .

15、2cos7.2432 = 0 . 4sin0 . 8function 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)/mm2022-5-27例例 的計(jì)算的計(jì)算1.單位圓的面積等于 2.3.0.20.40.60.810.20.40.60.8141102dxx4)1 (1102dxx用樣本平均值法求用樣本平均值法求 的值的值func

16、tion result=ybjzf1(a,b,m,mm)%a是積分的下限 %b是積分的上限%積分函數(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*4ybjzf1(0,1,1,100) result = 3.08745746887753ybjzf1(0,1,1,1000) result = 3.15500646827616ybjzf1(0,1,1,10000) result = 3.11911714237706ybjzf1(0,1,1,10

17、0000) result = 3.14014654862983ybjzf1(0,1,1,1000000) result = 3.14093979612119function 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*4ybjzf2(0,1,1,100) result = 3.04500162146

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

19、隨機(jī)數(shù)出現(xiàn)概率小，從而提高抽樣的效。 xgixix2 ixf幾種降低估計(jì)方差的幾種降低估計(jì)方差的MCMC方法方法重要抽樣法重要抽樣法 關(guān)鍵因素在于 的選取，使得估計(jì)的方差較小。 重要抽樣法的基本思想，就是通過(guò)選取與 形狀接近的密度函數(shù) 來(lái)降低估計(jì)的方差。 xg xf xg幾種降低估計(jì)方差的幾種降低估計(jì)方差的MCMC方法方法分層抽樣法分層抽樣法 同樣是利用貢獻(xiàn)率大小來(lái)降低估計(jì)方差的方法。它首先是把樣本空間 分成一些小區(qū)間 ，且諸 不交， ，然后在各個(gè)小區(qū)間內(nèi)的抽樣數(shù)由其貢獻(xiàn)大小決定。對(duì) 貢獻(xiàn)大的 抽樣多，可提高抽樣效率。如果能夠提出較好抽樣區(qū)間的分配和各子區(qū)間內(nèi)抽樣次數(shù)的分配方案，分層抽樣法估計(jì)

20、積分可以達(dá)到非常令人滿意的效果。 DmDD ,.,1iDDUDiiD幾種降低估計(jì)方差的幾種降低估計(jì)方差的MCMC方法方法關(guān)聯(lián)抽樣法關(guān)聯(lián)抽樣法 將需要估計(jì)的積分分解成兩個(gè)積分之差， 對(duì) 的估計(jì)轉(zhuǎn)化為對(duì) ， 的估計(jì)的差。則相應(yīng)的，其估計(jì)的方差的大小則與 ， 的估計(jì)的正相關(guān)度有關(guān)，若兩者的相關(guān)程度越高，則 的估計(jì)方差越小。這便是關(guān)聯(lián)抽樣法的基本出發(fā)點(diǎn)。 2121IIdxxfdxxfdxxfbababa1I2I1I2I一個(gè)元件(或系統(tǒng))能正常工作的概率稱為元件(或系統(tǒng))的可靠性系統(tǒng)由元件組成,常見(jiàn)的元件連接方式：串聯(lián)并聯(lián)12212 2、系統(tǒng)的可靠性計(jì)算問(wèn)題、系統(tǒng)的可靠性計(jì)算問(wèn)題 例例 設(shè)兩系統(tǒng)都是由

21、4 個(gè)元件組成,每個(gè)元件正常工作的概率為 p=0.5 , 每個(gè)元件是否正常工作相互獨(dú)立.兩系統(tǒng)的連接方式如下圖所示，比較兩系統(tǒng)的可靠性.A1A2B2B1S1:)()()()()()(2121212121211BBAAPBBPAAPBBAAPSP)2 (22242ppppfunction 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

22、,0.5,1,mm); Rguji=frq/mmfunction 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(1,ii)=1)&(randnuma2(1,ii)=1) pass1=1; else pass1=0; end if (randnumb1

23、(1,ii)=1)&(randnumb2(1,ii)=1) pass2=1; else pass2=0; end if (pass1+pass2)=1 frq=frq+1; end End,Rguji=frq/mm例例 設(shè)兩系統(tǒng)都是由 4 個(gè)元件組成,每個(gè)元件正常工作的概率為 p , 每個(gè)元件是否正常工作相互獨(dú)立.兩系統(tǒng)的連接方式如下圖所示，比較兩系統(tǒng)的可靠性.A1A2B2B1S1:)()()()()()(2121212121211BBAAPBBPAAPBBAAPSP)2 (22242ppppA1A2B2B1S2:212)()(iiiBAPSP)()(12SPSP22)2(pp)2 (

24、22pp222pp例例 設(shè)兩系統(tǒng)都是由 4 個(gè)元件組成,每個(gè)元件的壽命服從參數(shù)為的指數(shù)分布，每個(gè)元件是否正常工作相互獨(dú)立.兩系統(tǒng)的連接方式如下圖所示，求兩系統(tǒng)壽命大于T=100的概率.A1A2B2B1S1:)()()()()()(2121212121211BBAAPBBPAAPBBAAPSP例例 設(shè)兩系統(tǒng)都是由 4 個(gè)元件組成,每個(gè)元件的壽命服從參數(shù)為的指數(shù)分布，每個(gè)元件是否正常工作相互獨(dú)立.兩系統(tǒng)的連接方式如下圖所示，求兩系統(tǒng)壽命大于T=100的概率.A1A2B2B1S1:)()()()()()(2121212121211BBAAPBBPAAPBBAAPSPA1A2B2B1S2:212)()

25、(iiiBAPSPfunction 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

26、,1,mm);for ii=1:mm if (randnuma1(1,ii)t)&(randnuma2(1,ii)t) pass1=1; else pass1=0; end if (randnumb1(1,ii)t)&(randnumb2(1,ii)t) pass2=1; else pass2=0; end if (pass1+pass2)=1 frq=frq+1; end End,Rguji=frq/mmfunction Rguji=litiR2(t,thetaa1,thetaa2,thetab1,thetab2,mm)%t 是要求系統(tǒng)生存的壽命%thetaa1 是元件A1的數(shù)學(xué)期望%thetaa2 是元件A2的數(shù)學(xué)期望%thetab1 是元件B1的數(shù)學(xué)期望 %thetab2