1、蒙特卡洛方法初探實驗目的實驗目的實驗內容實驗內容學習計算機模擬的基本過程與方法。學習計算機模擬的基本過程與方法。1 1、模擬的概念。、模擬的概念。4 4、實驗作業。、實驗作業。3 3、計算機模擬實例。、計算機模擬實例。2 2、產生隨機數的計算機命令。、產生隨機數的計算機命令。模擬的概念模擬的概念 模擬就是利用物理的、數學的模型來類比、模仿現實系統及其演變過程，以尋求過程規律的一種方法。 模擬的基本思想是建立一個試驗模型，這個模型包含所研究系統的主要特點通過對這個實驗模型的運行，獲得所要研究系統的必要信息模擬的方法模擬的方法1、物理模擬物理模擬： 對實際系統及其過程用功能相似的實物系統去模仿。例

2、如，軍事演習、船艇實驗、沙盤作業等。 物理模擬通常花費較大、周期較長，且在物理模型上改變系統結構和系數都較困難。而且，許多系統無法進行物理模擬，如社會經濟系統、生態系統等。 在實際問題中，面對一些帶隨機因素的復雜系統，用分析方法建模常常需要作許多簡化假設，與面臨的實際問題可能相差甚遠，以致解答根本無法應用。這時，計算機模擬幾乎成為唯一的選擇。 在一定的假設條件下，運用數學運算模擬系統在一定的假設條件下，運用數學運算模擬系統的運行，稱為數學模擬。現代的數學模擬都是在的運行，稱為數學模擬。現代的數學模擬都是在計算機上進行的，稱為計算機模擬。計算機上進行的，稱為計算機模擬。2、數學模擬數學模擬 計算

3、機模擬可以反復進行，改變系統的結構和系數都比較容易。 蒙特卡洛（蒙特卡洛（Monte CarloMonte Carlo）方法）方法是一種應用隨機數來進行計算機模擬的方法此方法對研究的系統進行隨機觀察抽樣，通過對樣本值的統計分析，求得所研究系統的某些參數用蒙特卡洛方法進行計算機模擬的步驟用蒙特卡洛方法進行計算機模擬的步驟：1 設計一個邏輯框圖，即模擬模型這個框圖要正確反映系統各部分運行時的邏輯關系。2 模擬隨機現象可通過具有各種概率分布的模擬隨機數來模擬隨機現象產生模擬隨機數的計算機命令產生模擬隨機數的計算機命令 在Matlab軟件中，可以直接產生滿足各種分布的隨機數，命令如下：1產生m*n階(

4、a，b)均勻分布U（a，b）的隨機數矩陣： unifrnd (a,b,m, n); 產生一個a，b均勻分布的隨機數：unifrnd (a,b) 當只知道一個隨機變量取值在（a，b）內，但不知道（也沒理由假設）它在何處取值的概率大，在何處取值的概率小，就只好用U（a，b）來模擬它。2產生mm*nn階離散均勻分布的隨機數矩陣：R = unidrnd(N)R = unidrnd(N,mm,nn)當研究對象視為大量相互獨立的隨機變量之和，且其中每一種變量對總和的影響都很小時，可以認為該對象服從正態分布。若連續型隨機變量X的概率密度函數為 其中 0為常數，則稱X服從參數為 的指數指數分布分布。 0001

5、)(/xxexfx 指數分布的期望值為 排隊服務系統中顧客到達間隔、質量與可靠性中電子元件的壽命通常服從指數分布。例例 顧客到達某商店的間隔時間服從參數為顧客到達某商店的間隔時間服從參數為10(10(分鐘分鐘) )的指數分布的指數分布( (指數分布的均值為指數分布的均值為10) - - 指兩個顧客到達商店的平均間隔時間指兩個顧客到達商店的平均間隔時間是是1010分鐘分鐘. .即平均即平均1010分鐘到達分鐘到達1 1個顧客個顧客. . 顧客顧客到達的間隔時間可用到達的間隔時間可用exprnd(10)exprnd(10)模擬。模擬。設離散型隨機變量X的所有可能取值為0,1,2,且取各個值的概率為

6、其中 0為常數，則稱X服從參數為 的帕松分布帕松分布。, 2 , 1 , 0,!)(kkekXPk帕松分布在排隊系統、產品檢驗、天文、物理等領域有廣泛應用。帕松分布的期望值為6產生1個參數為n,p的二項分布的隨機數binornd(n,p)，產生mn個參數為n,p的二項分布的隨機數binornd(n,p,m,n) 。 擲一枚均勻硬幣，正面朝上的次數擲一枚均勻硬幣，正面朝上的次數X X服從參數為，服從參數為，p p的二項分布的二項分布,XB(1,p),XB(1,p)頻率的穩定性模擬頻率的穩定性模擬 1. 1.事件的頻率事件的頻率 在一組不變的條件下，重復作在一組不變的條件下，重復作n n次試驗，次

7、試驗，記記m m是是n n次試驗中事件次試驗中事件A A發生的次數。發生的次數。 頻率頻率 f=m/n f=m/n 2.2.頻率的穩定性頻率的穩定性 擲一枚均勻硬幣，記錄擲硬幣試驗中擲一枚均勻硬幣，記錄擲硬幣試驗中頻率頻率P P* *的波動情況。的波動情況。 function liti1(p,mm)pro=zeros(1,mm);randnum = binornd(1,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)在在MatlabMatlab中編輯中編輯.m.m文件輸入以下

8、命令：文件輸入以下命令：在在MatlabMatlab命令行中輸入以下命令：命令行中輸入以下命令：liti1(0.5,1000)在在MatlabMatlab命令行中輸入以下命令：命令行中輸入以下命令：liti1(0.5,10000)練習練習擲一枚不均勻硬幣，正面出現概率擲一枚不均勻硬幣，正面出現概率為為0.30.3，記錄前，記錄前10001000次擲硬幣試驗中正面次擲硬幣試驗中正面頻率的波動情況，并畫圖。頻率的波動情況，并畫圖。 在在MatlabMatlab命令行中輸入以下命令：命令行中輸入以下命令：liti1(0.3,1000)例例2 2 擲兩枚不均勻硬幣，每枚正面出現概率擲兩枚不均勻硬幣，每

9、枚正面出現概率為為0.40.4，記錄前，記錄前10001000次擲硬幣試驗中兩枚都為次擲硬幣試驗中兩枚都為正面頻率的波動情況，并畫圖。正面頻率的波動情況，并畫圖。 在在Matlab中編輯中編輯.m文件輸入以下命令：文件輸入以下命令：function liti2(p,mm)pro=zeros(1,mm);randnum = binornd(1,p,2,mm)；a=0; for i=1:mm a=a+randnum(1,i)*randnum(2,i); pro(i)=a/i;end pro=pro，num=1:mm;plot(num,pro)liti2(0.4,100)liti2(0.4,1000

10、0)例例3: 在一袋中有在一袋中有10 10 個相同的球，分別標有號碼個相同的球，分別標有號碼1,2,1,2,10,10。每次任取一個球，記錄其號碼后放。每次任取一個球，記錄其號碼后放回袋中，再任取下一個。這種取法叫做回袋中，再任取下一個。這種取法叫做“有放有放回抽取回抽取”。今有放回抽取。今有放回抽取3 3個球，求這個球，求這3 3個球的個球的號碼均為偶數的概率。號碼均為偶數的概率。（用頻率估計概率）（用頻率估計概率） 解：有放回取解：有放回取3個球個球, 所有取法有所有取法有103種種; 有放回有放回取取3個偶數號碼的球個偶數號碼的球, 所有取法有所有取法有53種種. 所以所以81105)

11、(33AP古典概率模擬古典概率模擬 function proguji=liti3(n,mm)frq=0;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 4 兩盒火柴，每盒兩盒火柴，每盒2020根。每次隨機在任根。每次隨機在任一盒中取出一根火柴。問其中一盒中火柴一盒中取出一根火柴。問其中一盒中火柴被取完而另一盒中至少還有被取完而另一盒中至少還有5 5

12、根火柴的概根火柴的概率有多大？（用頻率估計概率）率有多大？（用頻率估計概率） function proguji=liti40(mm)%mm 是隨機實驗次數frq=0; randnum=binornd(1,0.5,mm,2*20);proguji=0;for i=1:mm a1=0;a2=0;j=1; while (a120)&(a2=5 frq=frq+1; endend; proguji=frq/mm liti40(100)proguji = 0.4800 liti40(1000)proguji = 0.4970 liti40(10000)proguji = 0.4910 liti4

13、0(100000)proguji = 0.4984例例4 4 兩盒火柴，每盒兩盒火柴，每盒n根。每次隨機在任根。每次隨機在任一盒中取出一根火柴。問其中一盒中火柴一盒中取出一根火柴。問其中一盒中火柴被取完而另一盒中至少還有被取完而另一盒中至少還有k根火柴的概根火柴的概率有多大？（用頻率估計概率）率有多大？（用頻率估計概率） function proguji=liti4(n,k,mm)%n是每盒中的火柴數 %k 是剩余的火柴數%mm 是隨機實驗次數frq=0; randnum=binornd(1,0.5,mm,2*n);proguji=0;for i=1:mm a1=0;a2=0;j=1; whi

14、le (a1n)&(a2=k , frq=frq+1; end% a1=a1,a2=a2,frq % pauseend; proguji=frq/mm 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 幾何概率模擬幾何概率模擬1.1.定義定義 向任一可度量區域向任一可度量區域G G內投一點，如果所投的點內投一點，如果所投的點落在落在G G中任意可度量區域中任意可度量區域g

15、 g內的可能性與內的可能性與g g的度的度量成正比，而與量成正比，而與g g的位置和形狀無關，則稱這的位置和形狀無關，則稱這個隨機試驗為幾何型隨機試驗。或簡稱為幾何個隨機試驗為幾何型隨機試驗。或簡稱為幾何概型。概型。2. 概率計算概率計算 P(A)=A的度量的度量/S的度量的度量例例5 兩人約定于兩人約定于12點到點到1點到某地會面，先到者點到某地會面，先到者等等20分鐘后離去，試求兩人能會面的概率？分鐘后離去，試求兩人能會面的概率？ 解：設解：設x, y分別為甲、乙到達時刻分別為甲、乙到達時刻(分鐘分鐘)令令A=兩人能會面兩人能會面=(x,y)|x-y|20，x60,60,y y6060P(

16、A)=A的面積的面積/S的面積的面積=（602-402）/602=5/9=0.5556function proguji=liti5(mm)%mm 是隨機實驗次數frq=0;randnum1=unifrnd(0,60,mm,1);randnum2=unifrnd(0,60,mm,1);randnum=randnum1-randnum2;proguji=0;for ii=1:mm if abs(randnum(ii,1)=20 frq=frq+1; endendproguji=frq/mmliti5(10000)proguji = 0.5557例例6 6在我方某前沿防守地域，敵人以一個炮排（含兩門

17、火炮）為單位對我方進行干擾和破壞為躲避我方打擊，敵方對其陣地進行了偽裝并經常變換射擊地點 經過長期觀察發現，我方指揮所對敵方目標的指示有50是準確的，而我方火力單位，在指示正確時，有1/3的概率能毀傷敵人一門火炮，有1/6的概率能全部消滅敵人 現在希望能用某種方式把我方將要對敵人實施的1次打擊結果顯現出來，利用頻率穩定性，確定有效射擊(毀傷一門炮或全部消滅)的概率.復雜概率模擬復雜概率模擬分析分析： 這是一個復雜概率問題，可以通過理論計算得到相應的概率. 為了直觀地顯示我方射擊的過程，現采用模擬的方式。 需要模擬出以下兩件事： 1. 1. 問題分析問題分析1 1 觀察所對目標的指示正確與否觀察

18、所對目標的指示正確與否 模擬試驗有兩種結果，每一種結果出現的概率都是1/2 因此，可用投擲一枚硬幣的方式予以確定可用投擲一枚硬幣的方式予以確定,當硬幣出現正面時為指示正確，反之為不正確 2 2 當指示正確時，我方火力單位的射擊結果當指示正確時，我方火力單位的射擊結果情況情況 模擬試驗有三種結果：毀傷一門火炮的可能性為1/3(即2/6)，毀傷兩門的可能性為1/6，沒能毀傷敵火炮的可能性為1/2(即3/6) 這時可用投擲骰子的方法來確定可用投擲骰子的方法來確定：如果出現的是、三個點：則認為沒能擊中敵人；如果出現的是、點：則認為毀傷敵人一門火炮；若出現的是點：則認為毀傷敵人兩門火炮2. 2. 符號假

19、設符號假設i：要模擬的打擊次數； k1：沒擊中敵人火炮的射擊總數； k2：擊中敵人一門火炮的射擊總數；k3：擊中敵人兩門火炮的射擊總數E：有效射擊(毀傷一門炮或兩門炮)的概率3. 3. 模擬框圖模擬框圖初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子點數?k1=k1+1k2=k2+1k3=k3+1k1=k1+1imm?E=(k2+k3)/mm 停止硬幣正面?YNNY1，2，34，56function liti6(p,mm)efreq=zeros(1,mm);randnum1 = binornd(1,p,1,mm);randnum2 = unidrnd(6,1,mm);k1=0;k2=

20、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; end end efreq(i)=(k2+k3)/i; end num=1:mm;plot(num,efreq)在在MatlabMatlab中編輯中編輯.m.m文件輸入以下命令：文件輸入以下命令：在在MatlabMatlab命令行中輸入以下命令：命令行中輸入以下命令：liti6(0.5,2000)在在MatlabMatlab命令行中輸入以下命令：命令行中輸入以

21、下命令：liti6(0.5,20000)5. 5. 理論計算理論計算6. 6. 結果比較結果比較 模擬結果與理論計算近似一致，能更加真模擬結果與理論計算近似一致，能更加真實地表達實際戰斗動態過程實地表達實際戰斗動態過程 三三. .蒙特卡洛模擬的理論基礎與模擬結果的誤差蒙特卡洛模擬的理論基礎與模擬結果的誤差大數定律大數定律中心極限定理中心極限定理大數定律大數定律貝努里（Bernoulli） 大數定律設 nA 是 n 次獨立重復試驗中事件 A 發生的次數, p 是每次試驗中 A 發生的概率，則0有0limpnnPAn或1limpnnPAn在概率的統計定義中，事件 A 發生的頻率“ 穩定于”事件 A

22、 在一次試驗中發生的概率是指：nnAnnA頻率與 p 有較大偏差pnnA是小概率事件, 因而在 n 足夠大時, 可以用頻率近似代替 p . 這種穩定稱為依概率穩定.貝努里（貝努里（Bernoulli） 大數定律的意義大數定律的意義： 設 X1，X2，XN是來自總體X (EX)的簡單隨機樣本，即X1，X2，XN 獨立同分布，則 EXXNXPNiiN111)(limXEXPNN辛欽 大數定律即中心極限定理中心極限定理 設隨機變量序列,21nXXX相互獨立，服從同一分布，且有期望和方差：, 2 , 1,0)(,)(2kXDXEkk則對于任意實數 x ,xtnkkndtexnnXP21221lim若令

23、) 1 , 0(1NnnXnkknkknXnX11則等價于) 1 , 0(/NnXn于是1)(unuXPN這表明，不等式近似地以概率1成立。上式也表明，nuXnnX收斂到的階為O(n -1/2)。通常，蒙特卡羅方法的誤差定義為nu 關于蒙特卡羅方法的誤差需說明兩點：（1），蒙特卡羅方法的誤差為概率誤差，也即蒙特卡羅方法的收斂是概率意義下的收斂，雖然不能斷言其誤差不超過某個值，但能指出其誤差以接近1的概率不超過某個界限。如=0.5，誤差./6745. 0n此時，誤差超過的概率與小于的概率1-相等，都等于0.5。 來代替，在計算所求量的同時，可計算出 。2112)1(1niiniiXnXn 關于蒙

24、特卡羅方法的誤差需說明兩點：（2）誤差中的均方差 是未知的，必須使用其估計值 顯然，當給定置信度后，誤差由 和n決定。要減小 ，或者是增大n ，或者是減小方差 2。在 固定的情況下，要把精度提高一個數量級，試驗次數n需增加兩個數量級。因此，單純增大n不是一個有效的辦法。四、減小方差的各種技巧四、減小方差的各種技巧 nu 另一方面，如能減小估計的均方差 ，比如降低一半，那誤差就減小一半，這相當于n增大四倍的效果(n=(u /)2)。因此降低方差的各種技巧，引起了人們的普遍注意。nu 一般來說，降低方差的技巧，往往會使觀察一個子樣的時間增加。在固定時間內，使觀察的樣本數減少。所以，一種方法的優劣，

25、需要由方差和觀察一個子樣的費用（使用計算機的時間）兩者來衡量。這就是蒙特卡羅方法中效率的概念。它定義為nc ，其中c是觀察一個子樣的平均費用。顯然 nc= (u /)2 2c它與 2c成正比。四、效率四、效率 總而言之，作為提高蒙特卡洛方法效率的重要方向，是在減小標準差的同時兼顧考慮費用大小，使 2c盡可能地小。nc= (u /)2 2c)五、蒙特卡洛方法的特點五、蒙特卡洛方法的特點1) Monte Carlo方法及其程序結構簡單 產生隨機數，計算相應的值。即通過大量的簡單重復抽樣和簡單計算實現該方法。2）收斂速度與問題維數無關）收斂速度與問題維數無關 Monte Carlo方法的收斂速度為O

26、(n -1/2)，與一般數值方法相比很慢。因此，用Monte Carlo方法不能解決精確度要求很高的問題 從nu可見，Monte Carlo方法的誤差只與標準差和樣本容量n有關，而與樣本所在空間無關，即Monte Carlo方法的收斂速度與問題維數無關。而其他數值方法則不然。因此，這就決定了Monte Carlo方法對多維問題的適用性。 在解題時受問題條線限制的影響較小，例如要計算s維空間中的任一區域Ds上的積分ssDdxdxdxxxxggs2121),( 3） Monte Carlo方法的適用性強時，無論區域Ds的形狀多么特殊，只要能給出描述Ds的幾何特征的條件，就可以從Ds中均勻產生n個點 , 得到積分的近似值),()()(2)(1isii