Monte_Carlo模擬

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

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

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

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

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

6、針落下表示針落下后的中點，以后的中點，以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蒲豐投針實驗蒲豐投針實驗的計算機模擬：的計算機模擬：format long; %format long; %設(shè)置設(shè)置1515位顯示精度位顯示精度a=1;a=1; l=0.6;l=0.6; % %兩平行線間的寬度兩平行線間的寬度和和針長針長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; ; % %初始化計數(shù)器和設(shè)定初始化計數(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),r.)plot(phi(i),x(i),r.); ; counter=coun

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

10、驗，其結(jié)果列于下表 ：實驗者年份投計次數(shù)的實驗值沃爾弗(Wolf)185050003.1596斯密思(Smith)185532043.1553?？怂?Fox)189411203.1419拉查里尼(Lazzarini)190134083.1415929蒙特卡羅投點法是蒙特卡羅投點法是蒲豐投針實驗蒲豐投針實驗的推廣：的推廣： 在一個邊長為在一個邊長為a的正方形內(nèi)隨機投點，的正方形內(nèi)隨機投點，該點落在此正方形的內(nèi)切圓中的概率該點落在此正方形的內(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(投點法近似計算的投點法近似計算的為: ,num2str(4*m/n);xyo(a/2,a/2)/4a:a/222基本思想基本思想 由上面的例子可以看出，當所求問題的解由上面的例子可以看出，當所求問題的解是某個事件的概率，或者是某個隨機變量的數(shù)是某個事件的概率，或者是某個隨機變量的數(shù)學期望，或者是與之有關(guān)的量時，通過某種試學期望，或者是與之有關(guān)的量時，通過某種試驗的方法，得出該事件發(fā)生的頻率，再通過它驗的方法，得出該事件發(fā)生的頻率，再通過它得到問題的解。這就是蒙

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

13、新啟動每次重新啟動matlab時，輸出的隨機數(shù)不一樣時，輸出的隨機數(shù)不一樣注意注意: : 產(chǎn)生一個參數(shù)為產(chǎn)生一個參數(shù)為的指數(shù)分布的隨機數(shù)應(yīng)輸入的指數(shù)分布的隨機數(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的的隨機數(shù)矩陣隨機數(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)分布的隨機數(shù)矩陣的正態(tài)分布的隨機數(shù)矩陣 ra

14、ndom(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 G

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

16、思路：通過一個容易生成的概率分布 g 和一個取舍準則生成另一個與 g 相近的概率分布 f 。 為要模擬服從給定分布的隨機變量,用生成一個易于實現(xiàn)的不可約遍歷鏈 作為隨機樣本，使其平穩(wěn)分布為給定分布的方法，稱為馬氏鏈蒙特卡羅方法馬氏鏈蒙特卡羅方法.0,nXXn馬氏鏈蒙特卡羅方法馬氏鏈蒙特卡羅方法1 連續(xù)型隨機變量連續(xù)型隨機變量(以指數(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=finver

17、se(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ù)分布的隨機數(shù)指數(shù)分布的隨機數(shù)%指數(shù)分布隨機數(shù)產(chǎn)生函數(shù)已經(jīng)提供指數(shù)分布隨機數(shù)產(chǎn)生函數(shù)已經(jīng)提供 exprnd(1/3,1,1)2 離散型隨機變量離散型隨機變量(以離散分布為例以離散分布為例)：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);e

18、ndr=rand; index=find(rMAXK或PMAXP時停止迭代框框 圖圖初始化:給定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軟件包中編程，共需三個文件:randlp.m, mylp.m, lpconst.m.主程序為randlp.m.% mylp.m% mylp.mfunction z=mylp(x) %目標函數(shù)z=2*x(1)2+x(2)2-x(1)*x(2)-8*x(1)-3*x(2)