蒙特卡羅與計(jì)算機(jī)模擬有代碼

1、第六講 蒙特卡羅與計(jì)算機(jī)模擬內(nèi)容：計(jì)算機(jī)模擬(或稱仿真)是一種廣義數(shù)值計(jì)算 方法，適合解決一些規(guī)模大、難以解析化 以及受隨機(jī)因素影響的不確定數(shù)學(xué)模型目的：了解蒙特卡羅方法的基本思想，掌握利用 Matlab對(duì)離散/連續(xù)系統(tǒng)進(jìn)行模擬的方法要求：掌握Matlab隨機(jī)數(shù)函數(shù)，處理應(yīng)用問(wèn)題了解蒙特卡羅方法的起源和基本思想掌握離散系統(tǒng)和連續(xù)系統(tǒng)計(jì)算機(jī)模擬實(shí)例掌握隨機(jī)函數(shù) rand unifrnd normrnd exprnd了解Matlab仿真模塊 Simulink 蒙特卡羅方法的起源和基本思想 蒙特卡羅方法(Monte Carlo method)，或稱計(jì)算機(jī)隨機(jī)模擬方法，是一種基于“隨機(jī)數(shù)”的計(jì)算方法。

2、源于美國(guó)在第二次世界大戰(zhàn)研制原子彈的“曼哈頓計(jì)劃”，該計(jì)劃的主持人之一數(shù)學(xué)家馮諾伊曼用馳名世界的賭城摩納哥的Monte Carlo來(lái)命名這種方法，為它蒙上了一層神秘色彩。 蒙特卡羅方法的基本思想很早以前就被人們所發(fā)現(xiàn)和利用。早在17世紀(jì)，人們就知道用事件發(fā)生的“頻率”來(lái)決定事件的“概率”。19世紀(jì)人們用蒲豐投針的方法來(lái)計(jì)算圓周率，上世紀(jì)40年代電子計(jì)算機(jī)的出現(xiàn)，特別是近年來(lái)高速電子計(jì)算機(jī)的出現(xiàn)，使得用數(shù)學(xué)方法在計(jì)算機(jī)上大量、快速地模擬這樣的試驗(yàn)成為可能。蒲豐投針實(shí)驗(yàn)近似計(jì)算圓周率蒲豐投針實(shí)驗(yàn)： 法國(guó)科學(xué)家蒲豐(Buffon)在1777年提出的蒲豐投針實(shí)驗(yàn)是早期幾何概率一個(gè)非常著名的例子。蒲豐投

3、針實(shí)驗(yàn)的重要性并非是為了求得比其它方法更精確的值，而是它開(kāi)創(chuàng)了使用隨機(jī)數(shù)處理確定性數(shù)學(xué)問(wèn)題的先河，是用偶然性方法去解決確定性計(jì)算的前導(dǎo)，由此可以領(lǐng)略到從“概率土壤”上開(kāi)出的一朵瑰麗的鮮花蒙特卡羅方法(MC) 蒲豐投針實(shí)驗(yàn)可歸結(jié)為下面的數(shù)學(xué)問(wèn)題：平面上畫(huà)有距離為a的一些平行線，向平面上任意投一根長(zhǎng)為l(la)的針，假設(shè)針落在任意位置的可能性相同，試求針與平行線相交的概率P(從而求)蒲豐投針實(shí)驗(yàn)近似計(jì)算圓周率蒲豐投針實(shí)驗(yàn)： 如右圖所示，以M表示針落下后的中點(diǎn)，以x表示M到最近一條平行線的距離，以表示針與此線的交角：針落地的所有可能結(jié)果滿足：其樣本空間視作矩形區(qū)域, 面積是:針與平行線相交的條件：它

4、是樣本空間子集A，面積是：syms l phi; int(l/2*sin(phi),phi,0,pi); %ans=l因此，針與平行線相交的概率為：從而有： 特別當(dāng) 時(shí)蒲豐投針實(shí)驗(yàn)近似計(jì)算圓周率蒲豐投針實(shí)驗(yàn)的計(jì)算機(jī)模擬：format long; %設(shè)置15位顯示精度a=1; l=0.6; %兩平行線間的寬度和針長(zhǎng)figure; axis(0,pi,0,a/2); %初始化繪圖板set(gca,nextplot,add); %初始化繪圖方式為疊加counter=0; n=2010; %初始化計(jì)數(shù)器和設(shè)定投針次數(shù)x=unifrnd(0,a/2,1,n); phi=unifrnd(0,pi,1,n)

5、; %樣本空間for i=1:n if x(i)l*sin(phi(i)/2 %滿足此條件表示針與線的相交 plot(phi(i),x(i),r.); frame(i)=getframe; %描點(diǎn)并取幀 counter=counter+1; %統(tǒng)計(jì)針與線相交的次數(shù) endendfren=counter/n; pihat=2*l/(a*fren) %用頻率近似計(jì)算%movie(frame,1) %播放幀動(dòng)畫(huà)1次蒲豐投針實(shí)驗(yàn)近似計(jì)算圓周率蒲豐投針實(shí)驗(yàn)的計(jì)算機(jī)模擬： 意大利數(shù)學(xué)家拉澤里尼得到了準(zhǔn)確到6位小數(shù)的值，不過(guò)他的實(shí)驗(yàn)因?yàn)樘珳?zhǔn)確而受到了質(zhì)疑蒲豐投針實(shí)驗(yàn)計(jì)算圓周率蒙特卡羅投點(diǎn)法是蒲豐投針實(shí)驗(yàn)的推

6、廣： 在一個(gè)邊長(zhǎng)為a的正方形內(nèi)隨機(jī)投點(diǎn)，該點(diǎn)落在此正方形的內(nèi)切圓中的概率應(yīng)為該內(nèi)切圓與正方形的面積比值，即n=10000; a=2; m=0; for i=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(投點(diǎn)法近似計(jì)算的為: ,num2str(4*m/n);xyo(a/2,a/2)常見(jiàn)分布的隨機(jī)數(shù)產(chǎn)生語(yǔ)句 蒙特卡羅方法的關(guān)鍵步驟在于隨機(jī)數(shù)的產(chǎn)生，計(jì)算機(jī)產(chǎn)生的隨機(jī)數(shù)都不是真正的隨機(jī)數(shù)(由算法確定的緣故)，如果偽隨機(jī)數(shù)能夠通過(guò)一系列統(tǒng)計(jì)檢驗(yàn)，我們也可以將其當(dāng)作真正的隨機(jī)數(shù)使用： 常見(jiàn)

7、分布的隨機(jī)數(shù)產(chǎn)生語(yǔ)句MATLAB可以直接產(chǎn)生滿足各種分布的隨機(jī)數(shù)具體命令如下: 產(chǎn)生mn階0,1上均勻分布的隨機(jī)數(shù)矩陣 rand(m,n)產(chǎn)生一個(gè)0,1上均勻分布的隨機(jī)數(shù) rand 產(chǎn)生mn階a,b上均勻分布的隨機(jī)數(shù)矩陣 unifrnd (a,b,m, n) 產(chǎn)生一個(gè)a,b上均勻分布的隨機(jī)數(shù) unifrnd(a,b) 產(chǎn)生一個(gè)1:n的隨機(jī)排列(元素均出現(xiàn)且不重復(fù)) p=randperm(n)注意: randperm(6)與unifrnd (1,6,1, 6)的區(qū)別常見(jiàn)分布的隨機(jī)數(shù)產(chǎn)生語(yǔ)句 產(chǎn)生mn階均值為mu方差為sigma的正態(tài)分布的隨機(jī)數(shù)矩陣 normrnd(mu,sigma,m,n)產(chǎn)生

8、一個(gè)均值為mu方差為sigma的正態(tài)分布的隨機(jī)數(shù) normrnd(mu,sigma) 產(chǎn)生mn階期望值為mu (mu=1/)的指數(shù)分布的隨機(jī)數(shù)矩陣 exprnd(mu,m,n)產(chǎn)生一個(gè)期望值為mu的指數(shù)分布的隨機(jī)數(shù) exprnd(mu)注意: 產(chǎn)生一個(gè)參數(shù)為的指數(shù)分布的隨機(jī)數(shù)應(yīng)輸入 exprnd(1/)常見(jiàn)分布的隨機(jī)數(shù)產(chǎn)生語(yǔ)句 產(chǎn)生mn階參數(shù)為A1,A2,A3的指定分布name的隨機(jī)數(shù)矩陣 random(name,A1,A2,A3,m,n)產(chǎn)生一個(gè)參數(shù)為為A1,A2,A3的指定分布name的隨機(jī)數(shù) random(name,A1,A2,A3)舉例: 產(chǎn)生24階的均值為0方差為1的正態(tài)分布的隨機(jī)數(shù)

9、矩陣 random(Normal,0,1,2,4)name的取值可以是(詳情參見(jiàn)help random)：norm or Normal / unif or Uniformpoiss or Poisson / beta or Betaexp or Exponential / gam or Gammageo or Geometric / unid or Discrete UniformMATLAB隨機(jī)數(shù)的“重置”問(wèn)題 Matlab的隨機(jī)數(shù)是偽隨機(jī)數(shù)，但在一定的信度之下可以看作真正的隨機(jī)數(shù)。問(wèn)題是rand函數(shù)產(chǎn)生的隨機(jī)數(shù)從一個(gè)隨機(jī)數(shù)序列中取出來(lái)，而每次啟動(dòng)Matlab時(shí)，rand的狀態(tài)都會(huì)被重置(相

10、當(dāng)于把序列的指針移到了隨機(jī)數(shù)序列的開(kāi)始)，換言之第一次啟動(dòng)Matlab調(diào)用的第n次rand函數(shù)與下一次啟動(dòng)調(diào)用的第n個(gè)rand函數(shù)產(chǎn)生相同的數(shù)值。 如果想打亂這種狀態(tài)，可以為rand指定一個(gè)與當(dāng)前時(shí)間相關(guān)的初始狀態(tài)，而不用默認(rèn)狀態(tài)：rand(state,sum(100*clock);或者rand(state,sum(100*clock)*rand);非常見(jiàn)分布的隨機(jī)數(shù)的產(chǎn)生 對(duì)于常見(jiàn)分布隨機(jī)數(shù)，可由相應(yīng)Matlab函數(shù)直接產(chǎn)生，對(duì)于非常見(jiàn)分布隨機(jī)數(shù)可如下處理：1 連續(xù)型隨機(jī)變量(以p116指數(shù)分布為例)：syms t x lambda;Fx=int(lambda*exp(-lambda*t),

11、t,0,x) %分布函數(shù)syms r;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ù)分布的隨機(jī)數(shù)%指數(shù)分布隨機(jī)數(shù)產(chǎn)生函數(shù)已經(jīng)提供 exprnd(1/3,1,1)非常見(jiàn)分布的隨機(jī)數(shù)的產(chǎn)生2 離散型隨機(jī)變量(以p117離散分布為例)：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);endr=rand; ind

12、ex=find(rFx); x(index(1)-1)%已編寫(xiě)通用離散分布隨機(jī)數(shù)產(chǎn)生程序 scatrnd(x,px,n)離散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例范例 海港系統(tǒng)的卸載貨物問(wèn)題(p110-111,p119-129)離散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例范例 海港系統(tǒng)的卸載貨物問(wèn)題1 (p110-111,p119-129)程序片段(船只到港時(shí)間均勻分布,船只卸貨時(shí)間均勻分布)ShipBetweenTime(1)=unifrnd(15,145,1,1);%船只到港間隔時(shí)間隨機(jī)化(均勻分布)ShipUnloadTime(1)=unifrnd(45,90,1,1);%船只卸貨時(shí)間隨機(jī)化(均勻分布)通用程序haibor.

13、m可實(shí)現(xiàn)多次模擬，并將結(jié)果保存到H.txtdelete H.txt %清除歷史數(shù)據(jù)harbor(100,15,145,45,90)load H.txt;Hmean=mean(H); %導(dǎo)入H并按列取平均值離散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例范例 海港系統(tǒng)的卸載貨物問(wèn)題2 (p110-111,p119-129)程序片段(船只到港時(shí)間指數(shù)分布,船只卸貨時(shí)間均勻分布)ShipBetweenTime(1)=exprnd(60,1,1);%船只到港間隔時(shí)間隨機(jī)化(指數(shù)分布)ShipUnloadTime(1)=unifrnd(45,90,1,1);%船只卸貨時(shí)間隨機(jī)化(均勻分布)通用程序haibor2.m可實(shí)現(xiàn)多次模

14、擬，結(jié)果保存到H2.txtdelete H2.txt %清除歷史數(shù)據(jù)harbor2(100,60,45,90)load H2.txt;Hmean2=mean(H2); %導(dǎo)入H2并按列取平均值離散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例范例 海港系統(tǒng)的卸載貨物問(wèn)題3 (p110-111,p119-129)程序片段(船只到港時(shí)間離散分布,船只卸貨時(shí)間離散分布)1 編寫(xiě)船只到港間隔離散累積分布函數(shù)并作階梯圖：xs=15:10:145; for i=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

15、.172,0.125,0.071,0.037,0.017,0.003; Fx=0; for i=1:length(px),Fx=Fx,sum(px(1:i);endplot(10,x,Fx,-rs); hold on; stairs(0,x-5,145,Fx,1);set(gca,xtick,0:5:145); set(gca,xgrid,on); axis tight;離散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例范例 海港系統(tǒng)的卸載貨物問(wèn)題3 (p110-111,p119-129)程序片段(船只到港時(shí)間離散分布,船只卸貨時(shí)間離散分布)2 編寫(xiě)船只到港間隔離散累積分布反函數(shù)并作線性插值：Fxi=0:0.001:1

16、-eps;xi=interp1(Fx,0,x,Fxi,linear); r=rand(1,n);rnd=;for i=1:n index=find(r(i)=Fxi); if xs(1)=xi(index(1)-1)=xs(length(xs) rnd=rnd,xi(index(1)-1); endend%以上程序已編寫(xiě)通用M函數(shù)文件 harborrnd(xs,px,n)%即給出n個(gè)滿足離散分布(x,px)的船只到港間隔隨機(jī)數(shù)離散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例范例 海港系統(tǒng)的卸載貨物問(wèn)題3 (p110-111,p119-129)程序片段(船只到港時(shí)間離散分布,船只卸貨時(shí)間離散分布)3 編寫(xiě)船只卸貨時(shí)間離

17、散累積分布函數(shù)并作階梯圖：xs=45:5:90; for i=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; for i=1:length(px),Fx=Fx,sum(px(1:i);endplot(40,x,Fx,-rs); hold on; stairs(40,x-2.5,90,Fx,1);set(gca,xtick,40:2.5:90); set(gca,xgrid,on); axis tight;離散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例范例 海

18、港系統(tǒng)的卸載貨物問(wèn)題3 (p110-111,p119-129)程序片段(船只到港時(shí)間離散分布,船只卸貨時(shí)間離散分布)4 編寫(xiě)船只卸貨時(shí)間離散累積分布反函數(shù)并作線性插值：Fxi=0:0.001:1-eps;xi=interp1(Fx,0,x,Fxi,linear); r=rand(1,n);rnd=;for i=1:n index=find(r(i)=Fxi); if xs(1)=xi(index(1)-1)=xs(length(xs) rnd=rnd,xi(index(1)-1); endend%以上程序已編寫(xiě)通用M函數(shù)文件 harborrnd(xs,px,n)%即給出n個(gè)滿足離散分布(x,p

19、x)的船只卸貨時(shí)間隨機(jī)數(shù)離散系統(tǒng)的計(jì)算機(jī)模擬實(shí)例范例 海港系統(tǒng)的卸載貨物問(wèn)題3 (p110-111,p119-129)程序片段(船只到港時(shí)間指數(shù)分布,船只卸貨時(shí)間均勻分布)5 模擬船只到港間隔 / 卸貨時(shí)間均為離散分布的海港系統(tǒng)ShipBetweenTime(1)=harborrnd(sbtxs,sbtpx,1);%船只到港間隔時(shí)間隨機(jī)化(離散分布)ShipUnloadTime(1)=harborrnd(sutxs,sutpx,1);%船只卸貨時(shí)間隨機(jī)化(離散分布)通用程序haibor3.m可實(shí)現(xiàn)多次模擬，結(jié)果保存到H3.txtdelete H3.txt %清除歷史數(shù)據(jù)load harbor.mat %載入數(shù)據(jù)harbor3(100,sbtxs,sbtpx,sutx