數(shù)學(xué)建模之Monte_Carlo模擬

1、內(nèi)容提綱1.引言2.Monte Carlo模擬基本思想3.隨機(jī)數(shù)生成函數(shù)4.應(yīng)用實(shí)例舉例5.排隊(duì)論模擬6.Monte Carlo模擬求解規(guī)劃問題Monte Carlo方法：方法：蒙特卡羅方法，又稱隨機(jī)模擬方法，屬于計(jì)算數(shù)學(xué)的一個(gè)分支，它是在上世紀(jì)四蒙特卡羅方法，又稱隨機(jī)模擬方法，屬于計(jì)算數(shù)學(xué)的一個(gè)分支，它是在上世紀(jì)四十年代中期為了適應(yīng)當(dāng)時(shí)原子能事業(yè)的發(fā)展而發(fā)展起來的。十年代中期為了適應(yīng)當(dāng)時(shí)原子能事業(yè)的發(fā)展而發(fā)展起來的。亦稱統(tǒng)計(jì)模擬方法，亦稱統(tǒng)計(jì)模擬方法，statistical simulation method 利用隨機(jī)數(shù)進(jìn)行數(shù)值模擬的方法利用隨機(jī)數(shù)進(jìn)行數(shù)值模擬的方法Monte Carlo名字

2、的由來：名字的由來：Monte Carlo是摩納哥（是摩納哥（monaco)的首都，該城以賭博聞名的首都，該城以賭博聞名Nicholas Metropolis (1915-1999)Monte-Carlo, MonacoMonte CarloMonte Carlo方法的基本思想方法的基本思想 蒙特卡羅方法，或稱計(jì)算機(jī)隨機(jī)模擬方法，是一種基蒙特卡羅方法，或稱計(jì)算機(jī)隨機(jī)模擬方法，是一種基于于“隨機(jī)數(shù)隨機(jī)數(shù)”的計(jì)算方法。源于美國(guó)在第的計(jì)算方法。源于美國(guó)在第二二次世界大戰(zhàn)研次世界大戰(zhàn)研制原子彈的制原子彈的“曼哈頓計(jì)劃曼哈頓計(jì)劃”，該計(jì)劃的主持人之一數(shù)學(xué)家該計(jì)劃的主持人之一數(shù)學(xué)家馮馮諾伊曼用馳名世界的賭

3、城諾伊曼用馳名世界的賭城摩納哥的摩納哥的Monte CarloMonte Carlo來來命名這種方法，為它蒙上了一層神秘色彩。命名這種方法，為它蒙上了一層神秘色彩。 蒙特卡羅蒙特卡羅方法的基本思想很早以前就被人們所發(fā)現(xiàn)和方法的基本思想很早以前就被人們所發(fā)現(xiàn)和利用。早在利用。早在1717世紀(jì)，人們就知道用事件發(fā)生的世紀(jì)，人們就知道用事件發(fā)生的“頻率頻率”來來決定事件的決定事件的“概率概率”。1 19 9世紀(jì)人們用世紀(jì)人們用蒲豐蒲豐投針的方法來投針的方法來計(jì)計(jì)算算圓周率圓周率，上，上世紀(jì)世紀(jì)4040年代電子計(jì)算機(jī)的出現(xiàn)，特別是近年代電子計(jì)算機(jī)的出現(xiàn)，特別是近年來高速電子計(jì)算機(jī)的出現(xiàn)，使得用數(shù)學(xué)方法

4、在計(jì)算機(jī)上年來高速電子計(jì)算機(jī)的出現(xiàn)，使得用數(shù)學(xué)方法在計(jì)算機(jī)上大量、快速地模擬這樣的試驗(yàn)成為可能。大量、快速地模擬這樣的試驗(yàn)成為可能。蒲豐投針實(shí)驗(yàn)蒲豐投針實(shí)驗(yàn)： 法國(guó)科學(xué)家蒲豐法國(guó)科學(xué)家蒲豐(Buffon)(Buffon)在在17771777年提出的蒲豐年提出的蒲豐投針實(shí)驗(yàn)是早期幾何概率一個(gè)非常著名的例子。投針實(shí)驗(yàn)是早期幾何概率一個(gè)非常著名的例子。蒲豐蒲豐投針實(shí)驗(yàn)投針實(shí)驗(yàn)的重要性并非是為了求得比其它方法更精確的重要性并非是為了求得比其它方法更精確的的值，而是它開創(chuàng)了使用隨機(jī)數(shù)處理確定性數(shù)學(xué)問值，而是它開創(chuàng)了使用隨機(jī)數(shù)處理確定性數(shù)學(xué)問題的先河，是用偶然性方法去解決確定性計(jì)算的前導(dǎo)。題的先河，是用偶

5、然性方法去解決確定性計(jì)算的前導(dǎo)。由此由此可以可以領(lǐng)略到從領(lǐng)略到從“概率土壤概率土壤”上開出的一朵瑰麗的鮮上開出的一朵瑰麗的鮮花花- -蒙特卡羅蒙特卡羅方法方法( (MC)MC) 蒲豐投針實(shí)驗(yàn)可歸結(jié)為下面的數(shù)學(xué)問題：平面上畫蒲豐投針實(shí)驗(yàn)可歸結(jié)為下面的數(shù)學(xué)問題：平面上畫有距離為有距離為a a的一些平行線，向平面上任意投一根長(zhǎng)為的一些平行線，向平面上任意投一根長(zhǎng)為l l (la)(la)的針，假設(shè)針落在任意位置的可能性相同，試求針的針，假設(shè)針落在任意位置的可能性相同，試求針與平行線相交的概率與平行線相交的概率P(P(從而求從而求) )蒲豐投針實(shí)驗(yàn)蒲豐投針實(shí)驗(yàn)： 如右圖所示，以如右圖所示，以M M表示

6、針落下表示針落下后的中點(diǎn)，以后的中點(diǎn)，以x x表示表示M M到最近一條平行到最近一條平行線的距離，以線的距離，以表示針與此線的交角：表示針與此線的交角：針落地的所有可能結(jié)果滿足：針落地的所有可能結(jié)果滿足：其樣本空間視作矩形區(qū)域其樣本空間視作矩形區(qū)域, , 面積是面積是: :針與平行線相交的條件：針與平行線相交的條件：它是樣本空間它是樣本空間子集子集A A，面積是：，面積是：syms l phi; int(l/2*sin(phi),phi,0,pi) %ans=l因此，針與平行線相交的概率為：因此，針與平行線相交的概率為： 從而有：從而有： 0, 2/0ax 2/aS0, 2/sin0lx 02

7、/sinldlASalSASp/2)(/ )(apl /2蒲豐投針實(shí)驗(yàn)蒲豐投針實(shí)驗(yàn)的計(jì)算機(jī)模擬：的計(jì)算機(jī)模擬：format long; %format long; %設(shè)置設(shè)置1515位顯示精度位顯示精度a=1;a=1; l=0.6;l=0.6; % %兩平行線間的寬度兩平行線間的寬度和和針長(zhǎng)針長(zhǎng)figure; axis(0,pi,0,a/2); %figure; axis(0,pi,0,a/2); %初始化繪圖板初始化繪圖板set(gca,nextplot,set(gca,nextplot,a adddd);); % %初始化繪圖方式為疊加初始化繪圖方式為疊加counter=0;counter

8、=0; n=n=20102010; ; % %初始化計(jì)數(shù)器和設(shè)定初始化計(jì)數(shù)器和設(shè)定投針次數(shù)投針次數(shù)x=unifrnd(0,a/2,1,n);x=unifrnd(0,a/2,1,n); p phihi=unifrnd(0,pi,1,n);=unifrnd(0,pi,1,n); % %樣本空間樣本空間for i=1:nfor i=1:n if x(i)l if x(i)l* *sin(sin(p phihi(i)/2(i)/2 % %滿足此條件表示針與線的相交滿足此條件表示針與線的相交 plot(phi(i),x(i),plot(phi(i),x(i),r.r.) ); ; counter=cou

9、nter+1; counter=counter+1; % %統(tǒng)計(jì)針與線相交的次數(shù)統(tǒng)計(jì)針與線相交的次數(shù)frame(frame(countercounter)=getframe; %)=getframe; %描點(diǎn)并取幀描點(diǎn)并取幀 end endendendfren=counter/n;fren=counter/n; pihat=2pihat=2* *l/(al/(a* *fren)fren) % %用用頻率頻率近似計(jì)算近似計(jì)算figure(2)figure(2)movie(frame,1) %movie(frame,1) %播放幀動(dòng)畫播放幀動(dòng)畫1 1次次一些人進(jìn)行了實(shí)驗(yàn)，其結(jié)果列于下表一些人進(jìn)行了

10、實(shí)驗(yàn)，其結(jié)果列于下表 ：實(shí)驗(yàn)者年份投計(jì)次數(shù)的實(shí)驗(yàn)值沃爾弗(Wolf)185050003.1596斯密思(Smith)185532043.1553?？怂?Fox)189411203.1419拉查里尼(Lazzarini)190134083.1415929蒙特卡羅投點(diǎn)法是蒙特卡羅投點(diǎn)法是蒲豐投針實(shí)驗(yàn)蒲豐投針實(shí)驗(yàn)的推廣：的推廣： 在一個(gè)邊長(zhǎng)為在一個(gè)邊長(zhǎng)為a的正方形內(nèi)隨機(jī)投點(diǎn)，的正方形內(nèi)隨機(jī)投點(diǎn)，該點(diǎn)落在此正方形的內(nèi)切圓中的概率該點(diǎn)落在此正方形的內(nèi)切圓中的概率應(yīng)為該內(nèi)切圓與正方形的面積比值，應(yīng)為該內(nèi)切圓與正方形的面積比值，即即n=10000; a=2; m=0; for i=1:n x=rand(1)

11、*a; y=rand(1)*a; if ( (x-a/2)2+(y-a/2)2 = (a/2)2 ) m=m+1; endenddisp(投點(diǎn)法近似計(jì)算的投點(diǎn)法近似計(jì)算的為: ,num2str(4*m/n);xyo(a/2,a/2)/4a:a/222基本思想基本思想 由上面的例子可以看出，當(dāng)所求問題的解由上面的例子可以看出，當(dāng)所求問題的解是某個(gè)事件的概率，或者是某個(gè)隨機(jī)變量的數(shù)是某個(gè)事件的概率，或者是某個(gè)隨機(jī)變量的數(shù)學(xué)期望，或者是與之有關(guān)的量時(shí)，通過某種試學(xué)期望，或者是與之有關(guān)的量時(shí)，通過某種試驗(yàn)的方法，得出該事件發(fā)生的頻率，再通過它驗(yàn)的方法，得出該事件發(fā)生的頻率，再通過它得到問題的解。這就是

12、蒙特卡羅方法的基本思得到問題的解。這就是蒙特卡羅方法的基本思想。想。 蒙特卡羅方法的關(guān)鍵步驟在于隨機(jī)數(shù)的產(chǎn)蒙特卡羅方法的關(guān)鍵步驟在于隨機(jī)數(shù)的產(chǎn)生，計(jì)算機(jī)產(chǎn)生的隨機(jī)數(shù)都不是真正的隨機(jī)數(shù)生，計(jì)算機(jī)產(chǎn)生的隨機(jī)數(shù)都不是真正的隨機(jī)數(shù)( (由算法確定的緣故由算法確定的緣故) )，如果偽隨機(jī)數(shù)能夠通過，如果偽隨機(jī)數(shù)能夠通過一系列統(tǒng)計(jì)檢驗(yàn)，我們也可以將其當(dāng)作真正的一系列統(tǒng)計(jì)檢驗(yàn)，我們也可以將其當(dāng)作真正的隨機(jī)數(shù)使用。隨機(jī)數(shù)使用。rand(seed,0.1);rand(1)每次運(yùn)行程序產(chǎn)生的值是相同的每次運(yùn)行程序產(chǎn)生的值是相同的rand(state,sum(100*clock)*rand);rand(1) %每次

13、重新啟動(dòng)每次重新啟動(dòng)matlab時(shí)，輸出的隨機(jī)數(shù)不一樣時(shí)，輸出的隨機(jī)數(shù)不一樣注意注意: : 產(chǎn)生一個(gè)參數(shù)為產(chǎn)生一個(gè)參數(shù)為的指數(shù)分布的隨機(jī)數(shù)應(yīng)輸入的指數(shù)分布的隨機(jī)數(shù)應(yīng)輸入 exprnd(1/ exprnd(1/) )e0( )00txf xx產(chǎn)生產(chǎn)生m mn n階階參數(shù)為參數(shù)為A1,A2,A3A1,A2,A3的指定分布的指定分布namename的的隨機(jī)數(shù)矩陣隨機(jī)數(shù)矩陣 random(name,A1,A2,A3,m,n) random(name,A1,A2,A3,m,n)舉例舉例: : 產(chǎn)生產(chǎn)生2 24 4階的均值為階的均值為0 0方差為方差為1 1的正態(tài)分布的隨機(jī)數(shù)矩陣的正態(tài)分布的隨機(jī)數(shù)矩陣 r

14、andom(Normal,0,1,2,4)random(Normal,0,1,2,4) namename 的取值可以是的取值可以是( (詳情參見詳情參見help random)help random)：norm or Normal / unif or Uniformnorm or Normal / unif or Uniformpoiss or Poisson / beta or Betapoiss or Poisson / beta or Betaexp or Exponential / gam or Gammaexp or Exponential / gam or Gammageo or

15、Geometric / unid or Discrete Uniformgeo or Geometric / unid or Discrete Uniform非常見分布的隨機(jī)數(shù)的產(chǎn)生非常見分布的隨機(jī)數(shù)的產(chǎn)生 逆變換方法逆變換方法1 1 ( ) (0,1) ( ) F xUuFu由定理，要產(chǎn)生來自的隨機(jī)數(shù)，只要先產(chǎn)生來自隨機(jī)數(shù)，然后計(jì)算 即可。具體步驟如下：(1) (0,1)U上生成均勻分布的隨機(jī)數(shù) 。-1(2) , ( ) ( ) XXUFFx計(jì)算則為來自 分布的隨機(jī)數(shù).Acceptance-Rejection 方法方法最早由 Von Neumann提出，現(xiàn)在已經(jīng)廣泛應(yīng)用于各種隨機(jī)數(shù)的生成?；?/p>

16、本思路：通過一個(gè)容易生成的概率分布 g 和一個(gè)取舍準(zhǔn)則生成另一個(gè)與 g 相近的概率分布 f 。 為要模擬服從給定分布的隨機(jī)變量,用生成一個(gè)易于實(shí)現(xiàn)的不可約遍歷鏈 作為隨機(jī)樣本，使其平穩(wěn)分布為給定分布的方法，稱為馬氏鏈蒙特卡羅方法馬氏鏈蒙特卡羅方法.0,nXXn馬氏鏈蒙特卡羅方法馬氏鏈蒙特卡羅方法1 連續(xù)型隨機(jī)變量連續(xù)型隨機(jī)變量(以指數(shù)分布為例以指數(shù)分布為例)：e01-e0( )( )( )0000txxtxf tf t dt F xtx syms t x lambda;Fx=int(lambda*exp(-lambda*t),t,0,x) %分布函數(shù)分布函數(shù)syms r;Fxinv=finve

17、rse(Fx,x); %求反函數(shù)求反函數(shù)Fxinv=subs(Fxinv,x,r) %替換反函數(shù)變量替換反函數(shù)變量x為為rFxinv=inline(Fxinv)x=Fxinv(3,rand) %產(chǎn)生參數(shù)產(chǎn)生參數(shù) lambda=3 指數(shù)分布的隨機(jī)數(shù)指數(shù)分布的隨機(jī)數(shù)%指數(shù)分布隨機(jī)數(shù)產(chǎn)生函數(shù)已經(jīng)提供指數(shù)分布隨機(jī)數(shù)產(chǎn)生函數(shù)已經(jīng)提供 exprnd(1/3,1,1)2 離散型隨機(jī)變量離散型隨機(jī)變量(以離散分布為例以離散分布為例)：x=2,4,6,8; px=0.1,0.4,0.3,0.2; %以下為程序片段以下為程序片段Fx=0; for n=1:length(px),Fx=Fx,sum(px(1:n);

18、endr=rand; index=find(rMAXK或PMAXP時(shí)停止迭代框框 圖圖初始化:給定MAXK,MAXP;k=0,p=0,Q:大整數(shù)xj=aj+R(bj-aj) j=1,2,nj=0j=j+1,p=p+1PMAXP?YNxj=aj+R(bj-aj)gi(X)0?i=1,2nYNjMAXK?YN輸出X,Q,停止YN 例例 max 21212221382xxxxxxz s.t 10321 xx 01x 02x 在Matlab軟件包中編程，共需三個(gè)文件:randlp.m, mylp.m, lpconst.m.主程序?yàn)閞andlp.m.% mylp.m% mylp.mfunction z=mylp(x) %目標(biāo)函數(shù)z=2*x(1)2+x(2)2-x(1)*x(2)-8*x(1)-3*x(2)