蒙特卡洛方法基本思想

1、蒙特卡洛方法基本思想,實(shí)驗(yàn)?zāi)康?實(shí)驗(yàn)內(nèi)容,學(xué)習(xí)計算機(jī)模擬的基本過程與方法。,1、模擬的概念。,4、實(shí)驗(yàn)作業(yè)。,3、計算機(jī)模擬實(shí)例。,2、產(chǎn)生隨機(jī)數(shù)的計算機(jī)命令。,模擬的概念,模擬就是利用物理的、數(shù)學(xué)的模型來類比、模仿現(xiàn)實(shí)系統(tǒng)及其演變過程，以尋求過程規(guī)律的一種方法。,模擬的基本思想是建立一個試驗(yàn)?zāi)Ｐ?，這個模型包含所研究系統(tǒng)的主要特點(diǎn)．通過對這個實(shí)驗(yàn)?zāi)Ｐ偷倪\(yùn)行，獲得所要研究系統(tǒng)的必要信息,模擬的方法,1、物理模擬： 對實(shí)際系統(tǒng)及其過程用功

2、能相似的實(shí)物系統(tǒng)去模仿。例如，軍事演習(xí)、船艇實(shí)驗(yàn)、沙盤作業(yè)等。,物理模擬通常花費(fèi)較大、周期較長，且在物理模型上改變系統(tǒng)結(jié)構(gòu)和系數(shù)都較困難。而且，許多系統(tǒng)無法進(jìn)行物理模擬，如社會經(jīng)濟(jì)系統(tǒng)、生態(tài)系統(tǒng)等。,在實(shí)際問題中，面對一些帶隨機(jī)因素的復(fù)雜系統(tǒng)，用分析方法建模常常需要作許多簡化假設(shè)，與面臨的實(shí)際問題可能相差甚遠(yuǎn)，以致解答根本無法應(yīng)用。這時，計算機(jī)模擬幾乎成為唯一的選擇。,在一定的假設(shè)條件下，運(yùn)用數(shù)學(xué)運(yùn)算模擬系統(tǒng)的運(yùn)行，稱為數(shù)學(xué)模擬。現(xiàn)代的

3、數(shù)學(xué)模擬都是在計算機(jī)上進(jìn)行的，稱為計算機(jī)模擬。,2、數(shù)學(xué)模擬,計算機(jī)模擬可以反復(fù)進(jìn)行，改變系統(tǒng)的結(jié)構(gòu)和系數(shù)都比較容易。,蒙特卡洛（Monte Carlo）方法是一種應(yīng)用隨機(jī)數(shù)來進(jìn)行計算機(jī)模擬的方法．此方法對研究的系統(tǒng)進(jìn)行隨機(jī)觀察抽樣，通過對樣本值的觀察統(tǒng)計，求得所研究系統(tǒng)的某些參數(shù)．,蒙特卡洛方法也稱為隨機(jī)模擬方法,其起源最早可以追溯到18世紀(jì)下半葉的Buffon試驗(yàn).,用蒙特卡洛方法進(jìn)行計算機(jī)模擬的步驟：,[1] 設(shè)計一個邏輯框圖，即

4、模擬模型．這個框圖要正確反映系統(tǒng)各部分運(yùn)行時的邏輯關(guān)系。[2] 模擬隨機(jī)現(xiàn)象．可通過具有各種概率分布的模擬隨機(jī)數(shù)來模擬隨機(jī)現(xiàn)象．,產(chǎn)生模擬隨機(jī)數(shù)的計算機(jī)命令,,在Matlab軟件中，可以直接產(chǎn)生滿足各種分布的隨機(jī)數(shù)，命令如下：,2．產(chǎn)生m*n階離散均勻分布的隨機(jī)數(shù)矩陣：R = unidrnd(N)R = unidrnd(N,mm,nn),當(dāng)只知道一個隨機(jī)變量取值在（a，b）內(nèi)，但不知道（也沒理由假設(shè)）它在何處取值的概率大，在何處取

5、值的概率小，就只好用U（a，b）來模擬它。,1．產(chǎn)生m*n階［a，b］均勻分布U（a，b）的隨機(jī)數(shù)矩陣： unifrnd (a,b,m, n)產(chǎn)生一個[a，b]均勻分布的隨機(jī)數(shù)：unifrnd (a,b),當(dāng)研究對象視為大量相互獨(dú)立的隨機(jī)變量之和，且其中每一種變量對總和的影響都很小時，可以認(rèn)為該對象服從正態(tài)分布。,機(jī)械加工得到的零件尺寸的偏差、射擊命中點(diǎn)與目標(biāo)的偏差、各種測量誤差、人的身高、體重等，都可近似看成服從正態(tài)分布。,,若連續(xù)型

6、隨機(jī)變量X的概率密度函數(shù)為 其中 >0為常數(shù)，則稱X服從參數(shù)為 的指數(shù)分布。,指數(shù)分布的期望值為,,排隊服務(wù)系統(tǒng)中顧客到達(dá)率為常數(shù)時的到達(dá)間隔、故障率為常數(shù)時零件的壽命都服從指數(shù)分布。,指數(shù)分布在排隊論、可靠性分析中有廣泛應(yīng)用。,注意：Matlab中，產(chǎn)生參數(shù)為 的指數(shù)分布的命令為exprnd( ),例 顧客到達(dá)某商店的間隔時間服從參數(shù)為10的指數(shù)分布,指數(shù)分布的均值為10。 指兩個顧客到達(dá)商店的平均間隔時

7、間是10個單位時間.即平均10個單位時間到達(dá)1個顧客. 顧客到達(dá)的間隔時間可用exprnd(10)模擬。,設(shè)離散型隨機(jī)變量X的所有可能取值為0,1,2,…,且取各個值的概率為其中 >0為常數(shù)，則稱X服從參數(shù)為 的帕松分布。,帕松分布在排隊系統(tǒng)、產(chǎn)品檢驗(yàn)、天文、物理等領(lǐng)域有廣泛應(yīng)用。,帕松分布的期望值為,1 事件的頻率 在一組不變的條件下，重復(fù)作n次試驗(yàn)，記m是n次試驗(yàn)中事件A發(fā)生的次數(shù)。 頻率

8、 f=m/n,2.頻率的穩(wěn)定性,擲一枚均勻硬幣，記錄擲硬幣試驗(yàn)中頻率P*的波動情況。 R = binornd(N,P,mm,nn),例1 頻率的穩(wěn)定性,3 概率的頻率定義,在一組不變的條件下，重復(fù)作n次試驗(yàn)，記m是n次試驗(yàn)中事件A發(fā)生的次數(shù)。當(dāng)試驗(yàn)次數(shù)n很大時，如果頻率m/n穩(wěn)定地在某數(shù)值p附近擺動，而且一般地說，隨著試驗(yàn)次數(shù)的增加，這種擺動的幅度越來越小，稱數(shù)值p為事件A在這一組不變的條件下發(fā)生的概率，記作P(A)=p.,4 頻率的

9、基本性質(zhì),（1） 對任意事件A，有,,（2）,,,（3）若A1，A2，…，An是互不相容的，則,,頻率定義的意義:,(1) 提供了估計概率的方法;(2)提供了一種檢驗(yàn)理論正確與否的準(zhǔn)則.,理論依據(jù)： 大數(shù)定律,大量的隨機(jī)現(xiàn)象中平均結(jié)果的穩(wěn)定性,大數(shù)定律的客觀背景,大量拋擲硬幣正面出現(xiàn)頻率,字母使用頻率,生產(chǎn)過程中的廢品率,……,大數(shù)定律,貝努里（Bernoulli） 大數(shù)定律,設(shè) nA 是 n 次獨(dú)立重復(fù)試驗(yàn)中事件 A 發(fā)生的

10、次數(shù), p 是每次試驗(yàn)中 A 發(fā)生的概率，則,有,或,,貝努里（Bernoulli） 大數(shù)定律的意義：,定義,a 是一常數(shù)，,(或,故,在 Bernoulli 定理的證明過程中， Y n 是相互獨(dú)立的服從 0-1分布的隨機(jī)變量序列 {Xk} 的算術(shù)平均值, Y n 依概率收斂于其數(shù)學(xué)期望 p .,結(jié)果同樣適用于服從其它分布的獨(dú)立隨機(jī)變量序列,Chebyshev 大數(shù)定律,(指任意給定 n > 1,

11、 相互獨(dú)立)，且具有相同的數(shù)學(xué)期望和方差,或,定理的意義：,當(dāng) n 足夠大時，算術(shù)平均值幾乎就是一個常數(shù),可以用算術(shù)平均值近似地代替數(shù)學(xué)期望.,具有相同數(shù)學(xué)期望和方差的獨(dú)立隨機(jī)變量序列的算術(shù)平均值依概率收斂于數(shù)學(xué)期望.,例如要估計某地區(qū)的平均畝產(chǎn)量，要收割某些有代表性的地塊，例如n 塊. 計算其平均畝產(chǎn)量，則當(dāng)n 較大時，可用它作為整個地區(qū)平均畝產(chǎn)量的一個估計.,辛欽大數(shù)定律,具有相同的分布，且,記,則,則,大

12、數(shù)定律以嚴(yán)格的數(shù)學(xué)形式表達(dá)了隨機(jī)現(xiàn)象最根本的性質(zhì)之一：,它是隨機(jī)現(xiàn)象統(tǒng)計規(guī)律的具體表現(xiàn).,大數(shù)定律在理論和實(shí)際中都有廣泛的應(yīng)用.,平均結(jié)果的穩(wěn)定性,例1 頻率的穩(wěn)定性,1 事件的頻率 在一組不變的條件下，重復(fù)作n次試驗(yàn)，記m是n次試驗(yàn)中事件A發(fā)生的次數(shù)。 頻率 f=m/n,2.頻率的穩(wěn)定性,擲一枚均勻硬幣，記錄擲硬幣試驗(yàn)中頻率P*的波動情況。 R = binornd(N,P,mm,nn),functi

13、on liti1(n,p,mm)pro=zeros(1,mm);randnum = binornd(n,p,1,mm)a=0;for i=1:mm a=a+randnum(1,i); pro(i)=a/i;end pro=pronum=1:mm;plot(num,pro),在Matlab中編輯.m文件輸入以下命令：,在Matlab命令行中輸入以下命令：,liti1(1,0.5,1000),在Matl

14、ab命令行中輸入以下命令：,liti1(1,0.5,10000),練習(xí) 頻率的穩(wěn)定性,1 事件的頻率 R = binornd(N,P,mm,nn) 在一組不變的條件下，重復(fù)作n次試驗(yàn)，記m是n次試驗(yàn)中事件A發(fā)生的次數(shù)。 頻率 f=m/n,2.頻率的穩(wěn)定性,練習(xí)擲一枚不均勻硬幣，正面出現(xiàn)概率為0.3，記錄前1000次擲硬幣試驗(yàn)中正面頻率的波動情況，并畫圖。,在Matlab命令行中輸入以下命令

15、：,liti1(1,0.3,1000),例2 擲兩枚不均勻硬幣，每枚正面出現(xiàn)概率為0.4，記錄前1000次擲硬幣試驗(yàn)中兩枚都為正面頻率的波動情況，并畫圖。,在Matlab中編輯.m文件輸入以下命令：,function liti2(n,p,mm)pro=zeros(1,mm);randnum = binornd(n,p,2,mm)；a=0; for i=1:mm a=a+randnum(1,i)*randnum(2,i); p

16、ro(i)=a/i;end pro=pro，num=1:mm;plot(num,pro),熊宇樂y=zeros(1,1000);a=binornd(1,0.4,1,1000);b=binornd(1,0.4,1,1000);c=0;d=0;for i=1:1000 c=c+a(1,i).*b(1,i); y(i)=c/i;endy=y;num=1:1000;plot(num,y),孟亞func

17、tion bino (n,p,m)x=binornd(n,p,1,m);y=binornd(n,p,1,m);for i=1:m if x(i)==1&y(i)==1 s(i)=1; else s(i)=0; endend for i=1:m y(i)=sum(s(1,1:i))/i;endplot(y),liti2(1,0.4,100),liti

18、2(1,0.4,10000),在一袋中有10 個相同的球，分別標(biāo)有號碼1,2,…,10。每次任取一個球，記錄其號碼后放回袋中，再任取下一個。這種取法叫做“有放回抽取”。今有放回抽取3個球，求這3個球的號碼均為偶數(shù)的概率。（用頻率估計概率）,例3:,解：令A(yù)={有放回抽取3個球，求這3個球的號碼均為偶數(shù)}={(2,2,2)，(2,2,4)，….，(10,10,10)},function proguji=liti3(n,mm)frq=0;

19、randnum=unidrnd(n,mm,3);proguji=0;for i=1:mm a=(randnum(i,1)+1)*(randnum(i,2)+1)*(randnum(i,3)+1)； if mod(a,2)==1 frq=frq+1 endend; proguji=frq/mm,例4 兩盒火柴，每盒20根。每次隨機(jī)在任一盒中取出一根火柴。問其中一盒中火柴被取完而另一盒中至少還有5根

20、火柴的概率有多大？（用頻率估計概率）,>> liti4(20,5,100)proguji = 0.4800>> liti4(20,5,1000)proguji = 0.4970>> liti4(20,5,10000)proguji = 0.4910>> liti4(20,5,100000)proguji = 0.4984,function proguji=

21、liti4(nn,num,mm)%nn 是每盒中的火柴數(shù) %num 是剩余的火柴數(shù)%mm 是隨機(jī)實(shí)驗(yàn)次數(shù)frq=0; randnum=binornd(1,0.5,mm,2*nn);proguji=0;for i=1:mm a1=0;a2=0;j=1; while (a1=5 frq=frq+1; end% a1=a1,a2=a2,frq % pauseend; progu

22、ji=frq/mm,二. 幾何概率,,1.定義,向任一可度量區(qū)域G內(nèi)投一點(diǎn)，如果所投的點(diǎn)落在G中任意可度量區(qū)域g內(nèi)的可能性與g的度量成正比，而與g的位置和形狀無關(guān)，則稱這個隨機(jī)試驗(yàn)為幾何型隨機(jī)試驗(yàn)。或簡稱為幾何概型。,2. 概率計算,1.     P(A)=[A的度量]/[S的度量],兩人約定于12點(diǎn)到1點(diǎn)到某地會面，先到者等20分鐘后離去，試求兩人能會面的概率？,例5,解：設(shè)x,y分別為甲、

23、乙到達(dá)時刻(分鐘),令A(yù)={兩人能會面}={(x,y)||x-y|≤20，x≤60,y≤60},P(A)=A的面積/S的面積=（602-402）/602=5/9=0.5556,function proguji=liti5(mm)%mm 是隨機(jī)實(shí)驗(yàn)次數(shù)frq=0;randnum1=unifrnd(0,60,mm,1);randnum2=unifrnd(0,60,mm,1);randnum=randnum1-randnum2;

24、proguji=0;for ii=1:mm if abs(randnum(ii,1))<=20 frq=frq+1; endendproguji=frq/mm,liti5(10000)proguji = 0.5557,例6　在我方某前沿防守地域，敵人以一個炮排（含兩門火炮）為單位對我方進(jìn)行干擾和破壞．為躲避我方打擊，敵方對其陣地進(jìn)行了偽裝并經(jīng)常變換射擊地點(diǎn)．,經(jīng)過長期觀察發(fā)現(xiàn)，我方指

25、揮所對敵方目標(biāo)的指示有50％是準(zhǔn)確的，而我方火力單位，在指示正確時，有1/3的射擊效果能毀傷敵人一門火炮，有1/6的射擊效果能全部消滅敵人．,現(xiàn)在希望能用某種方式把我方將要對敵人實(shí)施的1打擊結(jié)果顯現(xiàn)出來，利用頻率穩(wěn)定性，確定有效射擊的概率,分析： 這是一個概率問題，可以通過理論計算得到相應(yīng)的概率和期望值.,為了能顯示我方射擊的過程，現(xiàn)采用模擬的方式。,需要模擬出以下兩件事：,1. 問題分析,[1] 觀察所對目標(biāo)的指示正確與否,模擬

26、試驗(yàn)有兩種結(jié)果，每一種結(jié)果出現(xiàn)的概率都是1/2．,因此，可用投擲一枚硬幣的方式予以確定，當(dāng)硬幣出現(xiàn)正面時為指示正確，反之為不正確．,[2] 當(dāng)指示正確時，我方火力單位的射擊結(jié)果情況,模擬試驗(yàn)有三種結(jié)果：毀傷一門火炮的可能性為1/3(即2/6)，毀傷兩門的可能性為1/6，沒能毀傷敵火炮的可能性為1/2(即3/6)．,這時可用投擲骰子的方法來確定：如果出現(xiàn)的是１、２、３三個點(diǎn)：則認(rèn)為沒能擊中敵人；如果出現(xiàn)的是４、５點(diǎn)：則認(rèn)為毀傷敵人一門

27、火炮；若出現(xiàn)的是６點(diǎn)：則認(rèn)為毀傷敵人兩門火炮．,2. 符號假設(shè),i：要模擬的打擊次數(shù)； k1：沒擊中敵人火炮的射擊總數(shù)； k2：擊中敵人一門火炮的射擊總數(shù)；k3：擊中敵人兩門火炮的射擊總數(shù)．E：有效射擊比率；,3. 模擬框圖,function liti6(p,mm)efreq=zeros(1,mm);randnum1 = binornd(1,p,1,mm);randnum2 = unidrnd(6,1,

28、mm);k1=0;k2=0;k3=0;for i=1:mm if randnum1(i)==0 k1=k1+1; else if randnum2(i)<=3 k1=k1+1; elseif randnum2(i)==6 k3=k3+1; else k2=k2+1;

29、 end end efreq(i)=(k2+k3)/i; end num=1:mm;plot(num,efreq),在Matlab中編輯.m文件輸入以下命令：,在Matlab命令行中輸入以下命令：,liti6(0.5,2000),在Matlab命令行中輸入以下命令：,liti6(0.5,20000),5. 理論計算,6. 結(jié)果比較,模擬結(jié)果與理論計算近似一致，能更加真實(shí)地表達(dá)實(shí)際戰(zhàn)斗動態(tài)過程．

30、,3. 某廠生產(chǎn)的燈泡能用1000小時的概率為0.8, 能用1500小時的概率為0.4 , 求已用1000小時的燈泡能用到1500小時的概率（頻率估計概率）。,2. 在一袋中有10 個相同的球，分別標(biāo)有號碼1,2,…,10。今有放回任取兩個球，求取得的第一個球號碼為奇數(shù)，第二個球的號碼為偶數(shù)的概率（頻率估計概率）,擲三枚不均勻硬幣，每枚正面出現(xiàn)概率為0.3，記錄前1000次擲硬幣試驗(yàn)中至少兩枚都為正面頻率的波動情況，并畫圖。,

31、作業(yè):,4 : 兩船欲停靠同一個碼頭, 設(shè)兩船到達(dá)碼頭的時間各不相干，而且到達(dá)碼頭的時間在一晝夜內(nèi)是等可能的. 如果兩船到達(dá)碼頭后需在碼頭停留的時間分別是1 小時與2 小 時，試求在一晝夜內(nèi)，任一船到達(dá)時，需要等待空出碼頭的概率. （頻率估計概率）,5 :在0,1,2,3,…..,9中不重復(fù)地任取4個數(shù)，求它們能排成首位非零的四位偶數(shù)的概率. （頻率估計概率）,6:從1,2,….. ,10十個數(shù)字中有放回地任取5個數(shù)字, 求取出的5

32、個數(shù)字中按由小到大排列, 中間的那個數(shù)等于 4 的概率. （頻率估計概率）,某廠生產(chǎn)的燈泡能用1000小時的概率為0.8, 能用1500小時的概率為0.4 , 求已用1000小時的燈泡能用到1500小時的概率,解 令 A 燈泡能用到1000小時 B 燈泡能用到1500小時,所求概率為,,例:,在一袋中有10 個相同的球，分別標(biāo)有號碼1,2,…,10。今有放回任取兩個球，求取得的第一個球號碼為奇數(shù)，第二個球的

33、號碼為偶數(shù)的概率。,解：令A(yù)={抽取2個球，第一個球號碼為奇數(shù)，第二個球的號碼為偶數(shù)}={(1,2)，(1,4)，….，(9,10)},例4 兩船欲?？客粋€碼頭, 設(shè)兩船到達(dá)碼頭的時間各不相干，而且到達(dá)碼頭的時間在一晝夜內(nèi)是等可能的. 如果兩船到達(dá)碼頭后需在碼頭停留的時間分別是1 小時與2 小 時，試求在一晝夜內(nèi)，任一船到達(dá)時，需 要等待空出碼頭的概率.,解 設(shè)船1 到達(dá)碼頭的瞬時為 x ，0 ? x < 24

34、 船2 到達(dá)碼頭的瞬時為 y ，0 ? y < 24,設(shè)事件 A 表示任一船到達(dá)碼頭時需要等待 空出碼頭,解,例5 在0,1,2,3, ,9中不重復(fù)地任取4個數(shù)，求它們能排成首位非零的四位偶數(shù)的概率.,設(shè) A為“能排成首位非零的四位偶數(shù)”,解 設(shè)取出的5個數(shù)按由小到大排列為,1,1,2,3,3;,1,1,2,3,4;,所取的5個數(shù)字中至少有3個數(shù)字不大于4,例6 從1,2, ,1