蒙特卡羅與計算機模擬有代碼

第六講蒙特卡羅與計算機模擬內(nèi)容：計算機模擬(或稱仿真)是一種廣義數(shù)值計算方法，適合解決一些規(guī)模大、難以解析化以及受隨機因素影響的不確定數(shù)學(xué)模型目的：了解蒙特卡羅方法的基本思想，掌握利用Matlab對離散/連續(xù)系統(tǒng)進行模擬的方法要求：掌握Matlab隨機數(shù)函數(shù)，處理應(yīng)用問題了解蒙特卡羅方法的起源和基本思想掌握離散系統(tǒng)和連續(xù)系統(tǒng)計算機模擬實例掌握隨機函數(shù)randunifrndnormrndexprnd了解Matlab仿真模塊Simulink蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第1頁!蒙特卡羅方法的起源和基本思想

蒙特卡羅方法(MonteCarlomethod)，或稱計算機隨機模擬方法，是一種基于“隨機數(shù)”的計算方法。源于美國在第二次世界大戰(zhàn)研制原子彈的“曼哈頓計劃”，該計劃的主持人之一數(shù)學(xué)家馮·諾伊曼用馳名世界的賭城—摩納哥的MonteCarlo—來命名這種方法，為它蒙上了一層神秘色彩。蒙特卡羅方法的基本思想很早以前就被人們所發(fā)現(xiàn)和利用。早在17世紀，人們就知道用事件發(fā)生的“頻率”來決定事件的“概率”。19世紀人們用蒲豐投針的方法來計算圓周率π，上世紀40年代電子計算機的出現(xiàn)，特別是近年來高速電子計算機的出現(xiàn)，使得用數(shù)學(xué)方法在計算機上大量、快速地模擬這樣的試驗成為可能。蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第2頁!蒲豐投針實驗近似計算圓周率π蒲豐投針實驗：法國科學(xué)家蒲豐(Buffon)在1777年提出的蒲豐投針實驗是早期幾何概率一個非常著名的例子。蒲豐投針實驗的重要性并非是為了求得比其它方法更精確的π值，而是它開創(chuàng)了使用隨機數(shù)處理確定性數(shù)學(xué)問題的先河，是用偶然性方法去解決確定性計算的前導(dǎo)，由此可以領(lǐng)略到從“概率土壤”上開出的一朵瑰麗的鮮花——蒙特卡羅方法(MC)蒲豐投針實驗可歸結(jié)為下面的數(shù)學(xué)問題：平面上畫有距離為a的一些平行線，向平面上任意投一根長為l(l<a)的針，假設(shè)針落在任意位置的可能性相同，試求針與平行線相交的概率P(從而求π)蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第3頁!蒲豐投針實驗近似計算圓周率π蒲豐投針實驗：如右圖所示，以M表示針落下后的中點，以x表示M到最近一條平行線的距離，以φ表示針與此線的交角：針落地的所有可能結(jié)果滿足：其樣本空間視作矩形區(qū)域Ω,面積是:針與平行線相交的條件：它是樣本空間Ω子集A，面積是：symslphi;int('l/2*sin(phi)',phi,0,pi);%ans=l因此，針與平行線相交的概率為：從而有：特別當(dāng)時蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第4頁!蒲豐投針實驗近似計算圓周率π蒲豐投針實驗的計算機模擬：意大利數(shù)學(xué)家拉澤里尼得到了準確到6位小數(shù)的π值，不過他的實驗因為太準確而受到了質(zhì)疑蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第5頁!常見分布的隨機數(shù)產(chǎn)生語句蒙特卡羅方法的關(guān)鍵步驟在于隨機數(shù)的產(chǎn)生，計算機產(chǎn)生的隨機數(shù)都不是真正的隨機數(shù)(由算法確定的緣故)，如果偽隨機數(shù)能夠通過一系列統(tǒng)計檢驗，我們也可以將其當(dāng)作真正的隨機數(shù)使用：蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第6頁!常見分布的隨機數(shù)產(chǎn)生語句④產(chǎn)生m×n階均值為mu方差為sigma的正態(tài)分布的隨機數(shù)矩陣

normrnd(mu,sigma,m,n)產(chǎn)生一個均值為mu方差為sigma的正態(tài)分布的隨機數(shù)normrnd(mu,sigma)⑤產(chǎn)生m×n階期望值為mu(mu=1/λ)的指數(shù)分布的隨機數(shù)矩陣

exprnd(mu,m,n)產(chǎn)生一個期望值為mu的指數(shù)分布的隨機數(shù)

exprnd(mu)注意:

產(chǎn)生一個參數(shù)為λ的指數(shù)分布的隨機數(shù)應(yīng)輸入

exprnd(1/λ)蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第7頁!MATLAB隨機數(shù)的“重置”問題

Matlab的隨機數(shù)是偽隨機數(shù)，但在一定的信度之下可以看作真正的隨機數(shù)。問題是rand函數(shù)產(chǎn)生的隨機數(shù)從一個隨機數(shù)序列中取出來，而每次啟動Matlab時，rand的狀態(tài)都會被重置(相當(dāng)于把序列的指針移到了隨機數(shù)序列的開始)，換言之次啟動Matlab調(diào)用的第n次rand函數(shù)與下一次啟動調(diào)用的第n個rand函數(shù)產(chǎn)生相同的數(shù)值。如果想打亂這種狀態(tài)，可以為rand指定一個與當(dāng)前時間相關(guān)的初始狀態(tài)，而不用默認狀態(tài)：rand('state',sum(100*clock));或者rand('state',sum(100*clock)*rand);蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第8頁!非常見分布的隨機數(shù)的產(chǎn)生[2]離散型隨機變量(以p117離散分布為例)：x=[2,4,6,8];px=[0.1,0.4,0.3,0.2];%以下為程序片段Fx=0;forn=1:length(px),Fx=[Fx,sum(px(1:n))];endr=rand;index=find(r<Fx);x(index(1)-1)%已編寫通用離散分布隨機數(shù)產(chǎn)生程序scatrnd(x,px,n)蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第9頁!離散系統(tǒng)的計算機模擬實例范例海港系統(tǒng)的卸載貨物問題1

(p110-111,p119-129)程序片段(船只到港時間均勻分布,船只卸貨時間均勻分布)ShipBetweenTime(1)=unifrnd(15,145,1,1);%船只到港間隔時間隨機化(均勻分布)ShipUnloadTime(1)=unifrnd(45,90,1,1);%船只卸貨時間隨機化(均勻分布)通用程序haibor.m可實現(xiàn)多次模擬，并將結(jié)果保存到H.txtdeleteH.txt%清除歷史數(shù)據(jù)harbor(100,15,145,45,90)loadH.txt;Hmean=mean(H);

%導(dǎo)入H并按列取平均值蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第10頁!離散系統(tǒng)的計算機模擬實例范例海港系統(tǒng)的卸載貨物問題3

(p110-111,p119-129)程序片段(船只到港時間離散分布,船只卸貨時間離散分布)[1]

編寫船只到港間隔離散累積分布函數(shù)并作階梯圖：xs=15:10:145;fori=1:length(xs)-1,x(i)=(xs(i)+xs(i+1))/2;endpx=[0.009,0.029,0.035,0.051,0.090,0.161,0.200,0.172,0.125,0.071,0.037,0.017,0.003];Fx=0;fori=1:length(px),Fx=[Fx,sum(px(1:i))];endplot([10,x],Fx,'-rs');holdon;stairs([0,x-5,145],[Fx,1]);set(gca,'xtick',0:5:145);set(gca,'xgrid','on');axistight;蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第11頁!離散系統(tǒng)的計算機模擬實例范例海港系統(tǒng)的卸載貨物問題3

(p110-111,p119-129)程序片段(船只到港時間離散分布,船只卸貨時間離散分布)[3]

編寫船只卸貨時間離散累積分布函數(shù)并作階梯圖：xs=45:5:90;fori=1:length(xs)-1,x(i)=(xs(i)+xs(i+1))/2;endpx=[0.017,0.045,0.095,0.086,0.130,0.185,0.208,0.143,0.091];Fx=0;fori=1:length(px),Fx=[Fx,sum(px(1:i))];endplot([40,x],Fx,'-rs');holdon;stairs([40,x-2.5,90],[Fx,1]);set(gca,'xtick',40:2.5:90);set(gca,'xgrid','on');axistight;蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第12頁!離散系統(tǒng)的計算機模擬實例范例海港系統(tǒng)的卸載貨物問題3

(p110-111,p119-129)程序片段(船只到港時間指數(shù)分布,船只卸貨時間均勻分布)[5]模擬船只到港間隔/卸貨時間均為離散分布的海港系統(tǒng)ShipBetweenTime(1)=harborrnd(sbtxs,sbtpx,1);%船只到港間隔時間隨機化(離散分布)ShipUnloadTime(1)=harborrnd(sutxs,sutpx,1);%船只卸貨時間隨機化(離散分布)通用程序haibor3.m可實現(xiàn)多次模擬，結(jié)果保存到H3.txtdeleteH3.txt%清除歷史數(shù)據(jù)loadharbor.mat%載入數(shù)據(jù)harbor3(100,sbtxs,sbtpx,sutxs,sutpx)loadH3.txt;Hmean3=mean(H3);

%導(dǎo)入H3并按列取平均值蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第13頁!連續(xù)系統(tǒng)的計算機模擬實例[1]

將隨等距時間τ連續(xù)變化的狀態(tài)變量軌跡x(t),y(t)用歐拉法離散化：蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第14頁!蒲豐投針實驗近似計算圓周率π蒲豐投針實驗的計算機模擬：formatlong;%設(shè)置15位顯示精度a=1;

l=0.6;

%兩平行線間的寬度和針長figure;axis([0,pi,0,a/2]);%初始化繪圖板set(gca,'nextplot','add');

%初始化繪圖方式為疊加counter=0;

n=2010;

%初始化計數(shù)器和設(shè)定投針次數(shù)x=unifrnd(0,a/2,1,n);phi=unifrnd(0,pi,1,n);

%樣本空間Ωfori=1:nifx(i)<l*sin(phi(i))/2

%滿足此條件表示針與線的相交

plot(phi(i),x(i),'r.');frame(i)=getframe;%描點并取幀

counter=counter+1;%統(tǒng)計針與線相交的次數(shù)endendfren=counter/n;

pihat=2*l/(a*fren)

%用頻率近似計算π%movie(frame,1)%播放幀動畫1次蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第15頁!蒲豐投針實驗計算圓周率π蒙特卡羅投點法是蒲豐投針實驗的推廣：在一個邊長為a的正方形內(nèi)隨機投點，該點落在此正方形的內(nèi)切圓中的概率應(yīng)為該內(nèi)切圓與正方形的面積比值，即n=10000;a=2;m=0;fori=1:n

x=rand(1)*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)蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第16頁!常見分布的隨機數(shù)產(chǎn)生語句MATLAB可以直接產(chǎn)生滿足各種分布的隨機數(shù)具體命令如下:①產(chǎn)生m×n階[0,1]上均勻分布的隨機數(shù)矩陣rand(m,n)產(chǎn)生一個[0,1]上均勻分布的隨機數(shù)

rand②產(chǎn)生m×n階[a,b]上均勻分布的隨機數(shù)矩陣

unifrnd(a,b,m,n)

產(chǎn)生一個[a,b]上均勻分布的隨機數(shù)

unifrnd(a,b)③產(chǎn)生一個1:n的隨機排列(元素均出現(xiàn)且不重復(fù))

p=randperm(n)注意:randperm(6)與unifrnd(1,6,1,6)的區(qū)別蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第17頁!常見分布的隨機數(shù)產(chǎn)生語句⑥

產(chǎn)生m×n階參數(shù)為A1,A2,A3的指定分布'name'的隨機數(shù)矩陣random('name',A1,A2,A3,m,n)產(chǎn)生一個參數(shù)為為A1,A2,A3的指定分布'name'的隨機數(shù)random('name',A1,A2,A3)舉例:產(chǎn)生2×4階的均值為0方差為1的正態(tài)分布的隨機數(shù)矩陣random('Normal',0,1,2,4)'name'的取值可以是(詳情參見helprandom)：'norm'or'Normal'/'unif'or'Uniform''poiss'or'Poisson'/'beta'or'Beta''exp'or'Exponential'/'gam'or'Gamma''geo'or'Geometric'/'unid'or'DiscreteUniform'……蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第18頁!非常見分布的隨機數(shù)的產(chǎn)生

對于常見分布隨機數(shù)，可由相應(yīng)Matlab函數(shù)直接產(chǎn)生，對于非常見分布隨機數(shù)可如下處理：[1]連續(xù)型隨機變量(以p116指數(shù)分布為例)：symstxlambda;Fx=int('lambda*exp(-lambda*t)',t,0,x)%分布函數(shù)symsr;Fxinv=finverse(Fx,x);%求反函數(shù)Fxinv=subs(Fxinv,x,r)%替換反函數(shù)變量x為rFxinv=inline(Fxinv)x=Fxinv(3,rand)%產(chǎn)生參數(shù)lambda=3指數(shù)分布的隨機數(shù)%指數(shù)分布隨機數(shù)產(chǎn)生函數(shù)已經(jīng)提供exprnd(1/3,1,1)蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第19頁!離散系統(tǒng)的計算機模擬實例范例海港系統(tǒng)的卸載貨物問題(p110-111,p119-129)蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第20頁!離散系統(tǒng)的計算機模擬實例范例海港系統(tǒng)的卸載貨物問題2

(p110-111,p119-129)程序片段(船只到港時間指數(shù)分布,船只卸貨時間均勻分布)ShipBetweenTime(1)=exprnd(60,1,1);%船只到港間隔時間隨機化(指數(shù)分布)ShipUnloadTime(1)=unifrnd(45,90,1,1);%船只卸貨時間隨機化(均勻分布)通用程序haibor2.m可實現(xiàn)多次模擬，結(jié)果保存到H2.txtdeleteH2.txt%清除歷史數(shù)據(jù)harbor2(100,60,45,90)loadH2.txt;Hmean2=mean(H2);

%導(dǎo)入H2并按列取平均值蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第21頁!離散系統(tǒng)的計算機模擬實例范例海港系統(tǒng)的卸載貨物問題3

(p110-111,p119-129)程序片段(船只到港時間離散分布,船只卸貨時間離散分布)[2]

編寫船只到港間隔離散累積分布反函數(shù)并作線性插值：Fxi=0:0.001:1-eps;xi=interp1(Fx,[0,x],Fxi,'linear');r=rand(1,n);rnd=[];fori=1:nindex=find(r(i)<=Fxi);ifxs(1)<=xi(index(1)-1)<=xs(length(xs))rnd=[rnd,xi(index(1)-1)];endend%以上程序已編寫通用M函數(shù)文件harborrnd(xs,px,n)%即給出n個滿足離散分布(x,px)的船只到港間隔隨機數(shù)蒙特卡羅與計算機模擬有代碼共25頁，您現(xiàn)在瀏覽的是第22頁!離散系統(tǒng)的計算機模擬實例范例海港系統(tǒng)的卸載貨物問題3

(p110-111,p119-129)程序片段(船只到港時間離散分布,船只卸貨時間離散分布)[4]

編寫船只卸貨時間離散累積分布反函數(shù)并作線性插值：Fxi=0:0.001:1-eps;xi=interp1(Fx,[0,x],Fxi,'linear');r=rand(1,n);rnd=[];fori=1:nindex=find(r(i)<