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

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

3、計算機模擬可以反復進行，改變系統的結構和系數都比較容易。 蒙特卡洛（蒙特卡洛（Monte CarloMonte Carlo）方法）方法是一種應用隨機數來進行計算機模擬的方法此方法對研究的系統進行隨機觀察抽樣，通過對樣本值的觀察統計，求得所研究系統的某些參數 蒙特卡洛方法也稱為蒙特卡洛方法也稱為隨機模擬方法,其起源最早可以追溯到18世紀下半葉的Buffon試驗.例 在1777年,法國學者Buffon提出用試驗方法求圓周率鸕鬧.其原理如下：假設平面上有元數條距離為1的等矩平行線,現向該平面隨機地投擲一根長度為KI1的針,則我們可以計算該針與任一平行線相交的概率.此處隨機投針可以這樣理解z針的中心點

4、與最近的平行線間的距離Z均勻地分布在區間0.1/2上,針與平行線的夾角以不管相交與否)均勻地分布在區間0,而上(見圖6。.于是,針與線相交的充要條件是本寸,從而針線相交概率為1用蒙特卡洛方法進行計算機模擬的步驟用蒙特卡洛方法進行計算機模擬的步驟：1 設計一個邏輯框圖，即模擬模型這個框圖要正確反映系統各部分運行時的邏輯關系。2 模擬隨機現象可通過具有各種概率分布的模擬隨機數來模擬隨機現象產生模擬隨機數的計算機命令產生模擬隨機數的計算機命令 在Matlab軟件中，可以直接產生滿足各種分布的隨機數，命令如下：2產生m*n階離散均勻分布的隨機數矩陣：R = unidrnd(N)R = unidrnd(

5、N,mm,nn) 當只知道一個隨機變量取值在（a，b）內，但不知道（也沒理由假設）它在何處取值的概率大，在何處取值的概率小，就只好用U（a，b）來模擬它。1產生m*n階a，b均勻分布U（a，b）的隨機數矩陣： unifrnd (a,b,m, n)產生一個a，b均勻分布的隨機數：unifrnd (a,b)當研究對象視為大量相互獨立的隨機變量之和，且其中每一種變量對總和的影響都很小時，可以認為該對象服從正態分布。機械加工得到的零件尺寸的偏差、射擊命中點與目標的偏差、各種測量誤差、人的身高、體重等，都可近似看成服從正態分布。若連續型隨機變量X的概率密度函數為 其中 0為常數，則稱X服從參數為 的指數

6、分指數分布布。0001)(/xxexft指數分布的期望值為 排隊服務系統中顧客到達率為常數時的到達間隔、故障率為常數時零件的壽命都服從指數分布。指數分布在排隊論、可靠性分析中有廣泛應用。注意：注意：Matlab中，產生參數為 的指數分布的命令為exprnd( )例例 顧客到達某商店的間隔時間服從參數為顧客到達某商店的間隔時間服從參數為1010的指的指數分布數分布 指數分布的均值為指數分布的均值為1010。 指兩個顧客到達商店的平均間隔時間是指兩個顧客到達商店的平均間隔時間是1010個單個單位時間位時間. .即平均即平均1010個單位時間到達個單位時間到達1 1個顧客個顧客. . 顧客顧客到達的

7、間隔時間可用到達的間隔時間可用exprnd(10)exprnd(10)模擬。模擬。設離散型隨機變量X的所有可能取值為0,1,2,且取各個值的概率為其中 0為常數，則稱X服從參數為 的帕松分布帕松分布。, 2 , 1 , 0,!)(kkekXPk帕松分布在排隊系統、產品檢驗、天文、物理等領域有廣泛應用。帕松分布的期望值為1 1 事件的頻率事件的頻率 在一組不變的條件下，重復作在一組不變的條件下，重復作n n次試驗，次試驗，記記m m是是n n次試驗中事件次試驗中事件A A發生的次數。發生的次數。 頻率頻率 f=m/n f=m/n 2.2.頻率的穩定性頻率的穩定性 擲一枚均勻硬幣，記錄擲硬幣試驗中

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

9、 p為事件為事件A A在這一組不變的條件在這一組不變的條件下發生的概率，記作下發生的概率，記作P(A)=p.P(A)=p.4 頻率的基本性質頻率的基本性質 （1） 對任意事件對任意事件A，有，有 1)(0AP（2）1)(SP0)(P（3）若）若A1，A2，An是互不相容的，則是互不相容的，則 )()(11nkknkkAPAP 頻率定義的意義頻率定義的意義:(1) 提供了估計概率的方法提供了估計概率的方法;(2)提供了一種檢驗理論正確與否的準則提供了一種檢驗理論正確與否的準則.理論依據：理論依據： 大數定律大數定律 大量的隨機現象中平均結果的穩定性大量的隨機現象中平均結果的穩定性 大數定律的客觀

10、背景大數定律的客觀背景大量拋擲硬幣大量拋擲硬幣正面出現頻率正面出現頻率字母使用頻率字母使用頻率生產過程中的生產過程中的廢品率廢品率大數定律大數定律貝努里（Bernoulli） 大數定律設 nA 是 n 次獨立重復試驗中事件 A 發生的次數, p 是每次試驗中 A 發生的概率，則0有0limpnnPAn或1limpnnPAn在概率的統計定義中，事件 A 發生的頻率“ 穩定于”事件 A 在一次試驗中發生的概率是指：nnAnnA頻率與 p 有較大偏差pnnA是小概率事件, 因而在 n 足夠大時, 可以用頻率近似代替 p . 這種穩定稱為依概率穩定.貝努里（貝努里（Bernoulli） 大數定律的意義

11、大數定律的意義：定義a 是一常數，0limaYPnn(或)1limaYPnn則稱隨機變量序列,21nYYY依概率收斂于常數 a , 記作aYnPn故pnnnPA,21nYYY是一系列隨機變量，設0有若在 Bernoulli 定理的證明過程中， Y n 是相互獨立的服從 0-1分布的隨機變量序列 Xk 的算術平均值, Y n 依概率收斂于其數學期望 p . 結果同樣適用于服從其它分布的獨立隨機變量序列Chebyshev 大數定律,21nXXX相互獨立，設隨機變量序列(指任意給定 n 1, 相互獨立)，且具有相同的數學期望和方差nXXX,21, 2 , 1,)(,)(2kXDXEkk則0有01li

12、m1nkknXnP或11lim1nkknXnP定理的意義定理的意義：當 n 足夠大時，算術平均值幾乎就是一個常數,可以用算術平均值近似地代替數學期望.具有相同數學期望和方差的獨立隨機變量序列的算術平均值依概率收斂于數學期望. 例如要估計某地區的平均畝產量，要例如要估計某地區的平均畝產量，要收割某些有代表性的地塊，例如收割某些有代表性的地塊，例如n n 塊塊. . 計算其平均畝產量，則當計算其平均畝產量，則當n n 較大時，可用較大時，可用它作為整個地區平均畝產量的一個估計它作為整個地區平均畝產量的一個估計. .辛欽大數定律辛欽大數定律 設,21nXXX相互獨立，服從同一分布，且具有數學期望 E

13、(X k) = , k= 1,2,則對任意正數 001lim1nkknXnP,21nXXX相互獨立，注3: 設隨機變量序列, 2 , 1,)(iXEkki則0有01lim1knikinXnP具有相同的分布，且記knikiMXn1111nPA),(21kAAAgnP),(21kg則則22nPAknPkA),(21kxxxg連續，若 大數定律以嚴格的數學形式表達了隨大數定律以嚴格的數學形式表達了隨機現象最根本的性質之一：機現象最根本的性質之一：它是隨機現象統計規律的具體表現它是隨機現象統計規律的具體表現.大數定律在理論和實際中都有廣泛的應用大數定律在理論和實際中都有廣泛的應用.平均結果的穩定性平均

14、結果的穩定性例例1 1 頻率的穩定性頻率的穩定性1 1 事件的頻率事件的頻率 在一組不變的條件下，重復作在一組不變的條件下，重復作n n次試驗，次試驗，記記m m是是n n次試驗中事件次試驗中事件A A發生的次數。發生的次數。 頻率頻率 f=m/n f=m/n 2.2.頻率的穩定性頻率的穩定性 擲一枚均勻硬幣，記錄擲硬幣試驗中頻率擲一枚均勻硬幣，記錄擲硬幣試驗中頻率P P* *的波動情況。的波動情況。 R = binornd(N,P,mm,nn) function liti1(n,p,mm)pro=zeros(1,mm);randnum = binornd(n,p,1,mm)a=0;for i

15、=1:mm a=a+randnum(1,i); pro(i)=a/i;end pro=pronum=1:mm;plot(num,pro)在在MatlabMatlab中編輯中編輯.m.m文件輸入以下命令：文件輸入以下命令：在在MatlabMatlab命令行中輸入以下命令：命令行中輸入以下命令：liti1(1,0.5,1000)在在MatlabMatlab命令行中輸入以下命令：命令行中輸入以下命令：liti1(1,0.5,10000)練習練習 頻率的穩定性頻率的穩定性1 1 事件的頻率事件的頻率 R = binornd(N,P,mm,nn) 在一組不變的條件下，重復作在一組不變的條件下，重復作n

16、n次試驗，次試驗，記記m m是是n n次試驗中事件次試驗中事件A A發生的次數。發生的次數。 頻率頻率 f=m/n f=m/n 2.2.頻率的穩定性頻率的穩定性 練習練習擲一枚不均勻硬幣，正面出現概率為擲一枚不均勻硬幣，正面出現概率為0.30.3，記錄前記錄前10001000次擲硬幣試驗中正面頻率的波動次擲硬幣試驗中正面頻率的波動情況，并畫圖。情況，并畫圖。 在在MatlabMatlab命令行中輸入以下命令：命令行中輸入以下命令：liti1(1,0.3,1000)例例2 2 擲兩枚不均勻硬幣，每枚正面出現概率擲兩枚不均勻硬幣，每枚正面出現概率為為0.40.4，記錄前，記錄前10001000次擲

17、硬幣試驗中兩枚都為次擲硬幣試驗中兩枚都為正面頻率的波動情況，并畫圖。正面頻率的波動情況，并畫圖。 在在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); pro(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

18、.4,1,1000); c=0;d=0; for i=1:1000 c=c+a(1,i).*b(1,i); y(i)=c/i; end y=y; num=1:1000; plot(num,y)孟亞function 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)liti2(1,0.4,10000)在一袋中有在一袋中有

19、10 10 個相同的球，分別標有號碼個相同的球，分別標有號碼1,2,101,2,10。每次任取一個球，記錄其號碼。每次任取一個球，記錄其號碼后放回袋中，再任取下一個。這種取法叫做后放回袋中，再任取下一個。這種取法叫做“有放回抽取有放回抽取”。今有放回抽取。今有放回抽取3 3個球，求個球，求這這3 3個球的號碼均為偶數的概率。個球的號碼均為偶數的概率。（用頻率估計（用頻率估計概率）概率） 例例3:解：令解：令A=有放回抽取有放回抽取3個球，求這個球，求這3個球的號碼個球的號碼均為偶數均為偶數=(2,2,2)，(2,2,4)，.，(10,10,10)81)(APfunction proguji=l

20、iti3(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根火柴的概率有多大？根火柴的概率有多大？（用頻率估計概率）（用頻率估計概

21、率） 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.4984function proguji=liti4(nn,num,mm)%nn 是每盒中的火柴數 %num 是剩余的火柴數%mm 是隨機實驗次數frq=0; randnum=binornd(1,0.5,mm,2*nn);proguji=0;for i=1:mm a1=0;a2=0;j=1; while (a120)&(a2=

22、5 frq=frq+1; end% a1=a1,a2=a2,frq % pauseend proguji=frq/mm二二. 幾何概率幾何概率1.1.定義定義 向任一可度量區域向任一可度量區域G G內投一點，如果所投內投一點，如果所投的點落在的點落在G G中任意可度量區域中任意可度量區域g g內的可能內的可能性與性與g g的度量成正比，而與的度量成正比，而與g g的位置和形的位置和形狀無關，則稱這個隨機試驗為幾何型隨狀無關，則稱這個隨機試驗為幾何型隨機試驗。或簡稱為幾何概型。機試驗。或簡稱為幾何概型。2. 概率計算概率計算 1. P(A)=A的度量的度量/S的度量的度量兩人約定于兩人約定于12

23、點到點到1點到某地會面，先點到某地會面，先到者等到者等20分鐘后離去，試求兩人能會面分鐘后離去，試求兩人能會面的概率？的概率？ 例例5解：設解：設x,y分別為甲、乙到達時刻分別為甲、乙到達時刻(分鐘分鐘)令令A=兩人能會面兩人能會面=(x,y)|x-y|20，x60,60,y y6060P(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=rand

24、num1-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在我方某前沿防守地域，敵人以一個炮排（含兩門火炮）為單位對我方進行干擾和破壞為躲避我方打擊，敵方對其陣地進行了偽裝并經常變換射擊地點 經過長期觀察發現，我方指揮所對敵方目標的指示有50是準確的，而我方火力單位，在指示正確時，有1/3的射擊效果能毀傷敵人一門火炮，有1/6的射擊效果能全部消滅敵人 現在希望能用某種方式把我方將要對敵人實施的1打擊結果顯現

25、出來，利用頻率穩定性，確定有效射擊的概率分析分析： 這是一個概率問題，可以通過理論計算得到相應的概率和期望值. 為了能顯示我方射擊的過程，現采用模擬的方式。 需要模擬出以下兩件事： 1. 1. 問題分析問題分析1 1 觀察所對目標的指示正確與否觀察所對目標的指示正確與否模擬試驗有兩種結果，每一種結果出現的概率都是1/2 因此，可用投擲一枚硬幣的方式予以確可用投擲一枚硬幣的方式予以確定定，當硬幣出現正面時為指示正確，反之為不正確2 2 當指示正確時，我方火力單位的射擊結當指示正確時，我方火力單位的射擊結果情況果情況 模擬試驗有三種結果：毀傷一門火炮的可能性為1/3(即2/6)，毀傷兩門的可能性為

26、1/6，沒能毀傷敵火炮的可能性為1/2(即3/6) 這時可用投擲骰子的方法來確定可用投擲骰子的方法來確定：如果出現的是、三個點：則認為沒能擊中敵人；如果出現的是、點：則認為毀傷敵人一門火炮；若出現的是點：則認為毀傷敵人兩門火炮2. 2. 符號假設符號假設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 停止硬幣正面?YNN

27、Y1，2，34，56function liti6(p,mm)efreq=zeros(1,mm);randnum1 = binornd(1,p,1,mm);randnum2 = unidrnd(6,1,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; end end efreq(i)=(k2+k3)/i; end num=1:mm;plot(num,efreq)在在MatlabMatla

28、b中編輯中編輯.m.m文件輸入以下命令：文件輸入以下命令：在在MatlabMatlab命令行中輸入以下命令：命令行中輸入以下命令：liti6(0.5,2000)在在MatlabMatlab命令行中輸入以下命令：命令行中輸入以下命令：liti6(0.5,20000)5. 5. 理論計算理論計算6. 6. 結果比較結果比較 模擬結果與理論計算近似一致，能更模擬結果與理論計算近似一致，能更加真實地表達實際戰斗動態過程加真實地表達實際戰斗動態過程 3. 某廠生產的燈泡能用某廠生產的燈泡能用1000小時的概率為小時的概率為0.8, 能能用用1500小時的概率為小時的概率為0.4 , 求已用求已用1000

29、小時的燈泡小時的燈泡能用到能用到1500小時的概率小時的概率2. 在一袋中有在一袋中有10 個相同的球，分別標有號碼個相同的球，分別標有號碼1,2,10。今有放回任取兩個球，求取得的第一個球號碼。今有放回任取兩個球，求取得的第一個球號碼為奇數，第二個球的號碼為偶數的概率。為奇數，第二個球的號碼為偶數的概率。 擲三枚不均勻硬幣，每枚正面出現概率為擲三枚不均勻硬幣，每枚正面出現概率為0.3，記，記 錄前錄前1000次擲硬幣試驗中至少兩枚都為正面頻率的次擲硬幣試驗中至少兩枚都為正面頻率的 波動情況，并畫圖。波動情況，并畫圖。 作業作業: 某廠生產的燈泡能用1000小時的概率為0.8, 能用1500小

30、時的概率為0.4 , 求已用1000小時的燈泡能用到1500小時的概率解解 令 A 燈泡能用到1000小時 B 燈泡能用到1500小時所求概率為)()(APABPABPAB218 . 04 . 0)()(APBP例例 例例:在一袋中有在一袋中有10 個相同的球，分別標有號個相同的球，分別標有號碼碼1,2,10。今不放回任取兩個球，求取。今不放回任取兩個球，求取得的第一個球號碼為奇數，第二個球的得的第一個球號碼為奇數，第二個球的號碼為偶數的概率。號碼為偶數的概率。解：令解：令A=A=抽取抽取2 2個球，第一個球號碼為奇數，第二個球個球，第一個球號碼為奇數，第二個球的號碼為偶數的號碼為偶數=(1,

31、2)=(1,2)，(1,4)(1,4)，.，(9,10)(9,10)185)(AP圓的面積 單位圓的面積等于 計算第一象限內的單位圓的面積，方法為：把它分成n個窄的曲邊梯形，計算S大，S小，其中n可以為1000, 10000, .實驗二實驗二: 的計算的計算數值積分法4)1(141102102dxxdxx梯形公式inabaxnabxfbfafdxxfiniiba,)(2/)()()(11辛普森公式nabxfxfbfafdxxfniiniiba6)(2)(4)()()(11105 . 0無窮級數法 arctg x =x-x3/3+x5/5-x7/7+x9/9-. /4 = arctg 1 =1-

32、1/3+1/5-1/7+1/9- 收斂太慢！ |x| 應當比 1 小很多，級數收斂才快。 /4 = arctg 1/2 + arctg 1/3 /4 = 4 arctg 1/5 - arctg 1/239 某次隨機投點的結果0.20.40.60.810.20.40.60.81 Monte-Carlo 方法方法謝文賢謝文賢西北工業大學理學院西北工業大學理學院 高性能計算中心高性能計算中心2004-05-142004-05-14Monte-Carlo 方法方法 一、一、MC 的起源和發展的起源和發展 Buffon試驗試驗 排隊系統模擬：排隊系統模擬：M/G/1M/G/1排隊系統排隊系統 二、二、M

33、C 的原理的原理 三、隨機數的產生原理與三、隨機數的產生原理與IMSL庫庫 均勻分布均勻分布U(0,1)的隨機數的產生的隨機數的產生 其他各種分布的隨機數的產生其他各種分布的隨機數的產生 隨機過程模擬隨機過程模擬 四、四、MC的應用舉例的應用舉例 定積分的定積分的MC計算計算 隨機微分方程的數值模擬隨機微分方程的數值模擬 五、五、EM算法及其算法及其MCMC方法方法一、一、MC 的起源和發展的起源和發展 隨機模擬方法，也稱為Monte Carlo方法，是一種基于“隨機數”的計算方法。這一方法源于美國在第一次世界大戰進行的研制原子彈的“曼哈頓計劃”。該計劃的主持人之一、數學家馮諾伊曼用馳名世界的

34、賭城摩納哥的Monte Carlo來命名這種方法，為它蒙上了一層神秘色彩。馮諾伊曼是公理化方法和計算機體系的領袖人物，Monte Carlo方法也是他的功勞。 事實上，Monte Carlo方法的基本思想很早以前就被人們所發現和利用。早在17世紀，人們就知道用事件發生的“頻率”來決定事件的“概率”。18世紀下半葉的法國學者Buffon提出用投針試驗的方法來確定圓周率的值。這個著名的Buffon試驗是Monte Carlo方法的最早的嘗試！ 歷史上曾有幾位學者相繼做過這樣的試驗。不過呢，他們的試驗是費時費力的，同時精度不夠高，實施起來也很困難。然而，隨著計算機技術的飛速發展，人們不需要具體實施這

35、些試驗，而只要在計算機上進行大量的、快速的模擬試驗就可以了。 在大眾的心目中，科學的代言人是心不在焉的牛頓或者爆炸式發型的愛因斯坦，但這只是傳統形象，比他們更了解現代計算技術的馮諾伊曼是個衣著考究，風度翩翩的人物，他說：純粹數學和應用數學的許多分支非常需要計算工具，用以打破目前由于純粹分析的研究方法不能解決非線性問題而形成的停滯狀態。 Monte Carlo方法是現代計算技術的最為杰出的成果之一，它在工程領域的作用是不可比擬的。 Buffon試驗試驗 假設平面上有無數條距離為1的等距平行線，現向該平面隨機投擲一根長度為 的針( )， 則我們可計算該針與任一平行線相交的概率。這里，隨機投針指的是

36、：針的中心點與最近的平行線間的距離 均勻的分布在區間 。上，針與平行線的夾角 （不管相交與否）均勻的分布在區間 上。 因此，針與線相交的充要條件是 l1l21 , 0 x, 021sinxBuffon試驗試驗 從而針線相交的概率為 根據上式，若我們做大量的投針試驗并記錄針與線相交的次數，則由大數定理可以估計出針線相交的概率 ，從而得到 的估計值。 針與線的位置關系：ldxdlXPpl22sin20sin20 xlp排隊系統模擬：排隊系統模擬：M/G/1M/G/1排隊系統排隊系統服務規則：服務規則：先到先服務假設假設：（）顧客到達遵循Poisson分布；（）服務時間服從一般分布；（）到達間隔與服

37、務時間相互獨立1排隊系統模擬：排隊系統模擬：M/G/1M/G/1排隊系統排隊系統關心的指標：關心的指標：（）時刻t時，系統中的顧客數；即隊長的分布；（）顧客的等待時間；（）服務的忙碌程度；（） 用最樸素的Monte-Carlo方法可以得到這些指標的估計二、二、MC 的原理的原理 應用Monte Carlo方法求解工程技術問題可以分為兩類： 確定性問題 隨機性問題確定性系統確定性系統隨機性系統隨機性系統模擬模擬自然界自然界Monte-Carlo模擬，即隨機模擬（重復模擬，即隨機模擬（重復“試驗試驗”）重復試驗重復試驗計算機模擬計算機模擬思路：思路： 1、 針對實際問題建立一個簡單且便于實現的概率

38、統計模型，使問題的解對應于該模型中隨機變量的概率分布或其某些數字特征，比如，均值和方差等。所構造的模型在主要特征參量方面要與實際問題或系統相一致的。 2、 根據模型中各個隨機變量的分布，在計算機上產生隨機數，實現一次模擬過程所需的足夠數量的隨機數。通常先產生均勻分布的隨機數，然后生成服從某一分布的隨機數，再進行隨機模擬試驗。 思路：思路： 3、 根據概率模型的特點和隨機變量的分布特性，設計和選取合適的抽樣方法，并對每個隨機變量進行抽樣（包括直接抽樣、分層抽樣、相關抽樣、重要抽樣等）。 4、 按照所建立的模型進行仿真試驗、計算，求出問題的隨機解。 5、 統計分析模擬試驗結果，給出問題的估計以及其

39、精度估計。 6、 必要時，還應改進模型以降低估計方差和減少試驗費用，提高模擬計算的效率。 收斂性收斂性: : 由大數定律, Monte-Carlo模擬的收斂是以概率而言的 誤差誤差: : 用頻率估計概率時誤差的估計，可由中心極限定理，給定置信水平 的條件下，有： 模擬次數模擬次數: :由誤差公式得NU2/1|)(XgVar22/1)(UN三、隨機數的產生原理與三、隨機數的產生原理與IMSL庫庫 隨機數的產生是實現MC計算的先決條件。而大多數概率分布的隨機數的產生都是基于均勻分布U(0,1)的隨機數。 首先，介紹服從均勻分布U(0,1)的隨機數的產生方法。 其次，介紹服從其他各種分布的隨機數的產

40、生方法。以及服從正態分布的隨機數的產生方法。 接下來，介紹隨機過程模擬。 最后，關于隨機數的幾點注。均勻分布均勻分布U(0,1)的隨機數的產生的隨機數的產生 產生均勻分布的標準算法在很多高級計算機語言的書都可以看到。算法簡單，容易實現。使用者可以自己手動編程實現。IMSL統計庫中也提供給我們用于產生均勻分布的各種函數。我們的重點是怎樣通過均勻分布產生服從其他分布的隨機數。因此，直接使用IMSL提供的可靠安全的標準函數，當然不用費事了。 IMSL庫中的函數使用庫中的函數使用 RNSET: 種子的設定 CALL RNSET (ISEED) RNOPT: 產生器的類型的設定 CALL RNOPT (

41、IOPT) RNUN/DRNUN: 產生均勻分布的隨機數 CALL RNUN (NR, R) 其他各種分布的隨機數的產生其他各種分布的隨機數的產生 基本方法有如下三種： 逆變換法 合成法 篩選法 逆變換法逆變換法 設隨機變量 的分布函數為 ，定義 定理定理 設隨機變量 服從 上的均勻分布，則 的分布函數為 。 因此，要產生來自 的隨機數，只要先產生來自 的隨機數，然后計算 即可。 其步驟為 X xF 10 ,:inf1yyxFxyFU) 1 , 0( UFX1 xF xF1 , 0U uF1 uFxuU11 , 0計算，抽取由合成法合成法 合成法的應用最早見于Butlter 的書中。 構思如下

42、： 如果 的密度函數 難于抽樣，而 關于 的條件密度函數 以及 的密度函數 均易于抽樣，則 的隨機數可如下產生： 可以證明由此得到 的服從 。X xpXYyxpY ygX xyxpyygY抽取由條件分布，抽取的分布由X xp篩選抽樣篩選抽樣 假設我們要從 抽樣，如果可以將 表示成 ，其中 是一個密度函數且易于抽樣，而 ， 是常數，則 的抽樣可如下進行： 定理定理 設 的密度函數 ，且 ，其中 ， ， 是一個密度函數。令 和 分別服從 和 ，則在 的條件 下， 的條件密度為 xp xgxhcxp xph 10 xg1cX 。，回到如果，停止，則如果，抽取，由抽取由1321 , 01yguyxyg

43、uyyhuUX xp xgxhcxph 10 xg1cUY1 , 0U yh YgU xpYgUxpYY標準正態分布的隨機數標準正態分布的隨機數 的隨機數產生方法很多。簡要介紹三種。 法法1、 變換法（Box 和Muller 1958） 設 ， 是獨立同分布的 變量，令 則 與 獨立，均服從標準正態分布。 法法2、 結合合成法與篩選法。（略） 法法3 3、 近似方法（利用中心極限定理） 即用 個 變量產生一個 變量。 其中 是抽自 的隨機數， 可近似為一個 變量。1 , 0N1U2U1 , 0U2122112sinln22cosln22121UUXUUX1X2Xn1 , 0U1 , 0N211

44、2unxiu1 , 0Uniiunu111 , 0NIMSL庫中的函數使用庫中的函數使用 RNSET: 種子的設定 CALL RNSET (ISEED) RNOPT: 產生器的類型的設定 CALL RNOPT (IOPT) RNNOA: 產生服從標準正態分布的偽隨機數 CALL RNNOA (NR, R) 隨機過程模擬隨機過程模擬 高斯白噪聲的產生： 一種簡單有效的模擬高斯白噪聲過程的方法是由獨立的單位正態隨機序列 ，按照如圖所示的方式連接，其中 式中 為白噪聲的強度。kktDkD 由于正態隨機數的獨立性，當 很小時，按(a)、(b)所構成的過程的相關時間很短，從而具有很高的截止頻率。 當然，

45、 過小，樣本的計算量過大。因此，選取 適當小即可。ttt關于隨機數的幾點注關于隨機數的幾點注 注注1 1 由于均勻分布的隨機數的產生總是采用某個確定的模型進行的，從理論上講，總會有周期現象出現的。初值確定后，所有隨機數也隨之確定，并不滿足真正隨機數的要求。因此通常把由數學方法產生的隨機數成為偽隨機數。 但其周期又相當長，在實際應用中幾乎不可能出現。因此，這種由計算機產生的偽隨機數可以當作真正的隨機數來處理。 注注2 2 應對所產生的偽隨機數作各種統計檢驗，如獨立性檢驗，分布檢驗，功率譜檢驗等等。 四、四、MC的應用舉例的應用舉例 例一、定積分的例一、定積分的MC計算計算 隨機投點法隨機投點法

46、樣本平均值法樣本平均值法 幾種降低估計方差的幾種降低估計方差的MC方法方法 例二、例二、 隨機微分方程的數值模擬隨機微分方程的數值模擬 系統的可靠性系統的可靠性數值模擬數值模擬計算問題計算問題定積分的定積分的MC計算計算 事實上，不少的統計問題，如計算概率、各階距等，最后都歸結為定積分的近似計算問題。 下面考慮一個簡單的定積分 為了說明問題，我們首先介紹兩種求 的簡單的MC方法，然后給出幾種較為復雜而更有效的MC方法。 dxxfba 在計算積分上，MC的實用場合是計算重積分 其中 是 維空間的點，當 較大時，用MC方法比一般的數值方法有優點，主要是它的誤差與維數 無關。 dPPgIkPmmm隨

47、機投點法隨機投點法 方法簡述方法簡述： 設 ，有限， ， ，并設 是在 上均勻分布的二維隨機變量，其聯合密度函數為 。則易見 是 中 曲線下方的面積。假設我們向 中進行隨機投點，若點落在 下方，(即 稱為中的，否則稱為不中，則點中的概率為 。若我們進行了 次投點，其中 次中的，則用頻率來估計概率 。即 。ab Mxf0Mybxayx0 ,:,YX,MybxaIabM0,1 dxxfba xfy xfy xfy abMpn0npnnp0 那么我們可以得到 的一個估計 具體試驗步驟為nnabM01 用上式來估計的個數統計，和，計算，隨機數，個獨立地產生0,.,11 , 02nyxfxfMvyabu

48、axnivuUniiiiiiiii 注注1 1 隨機投點法的思想簡單明了，且每次投點結果服從二項分布，故 ，其中 注注2 2 可證 是 的無偏估計。若用估計的標準差來衡量其精度，則估計 的精度的階為 。 注注3 3 這里，定積分的解，就對應我們選定的隨機變量的概率值。 pnbn,0abMp1121n樣本平均值法樣本平均值法 基本原理基本原理：對積分 ，設 是 上的一個密度函數，改寫 可見，任一積分均可以表示為某個隨機變量（函數）的期望。由矩法，若有 個來自 的觀測值，則可給出 的一個矩估計。 最簡單的，若 ，有限，可取 。 dxxfba xgba, XgXfEdxxgxgxfban xgab abxg /1 設 是來自 的隨機數，則 的一個估計為 具體步驟為 注注 可證 是 的無偏估計。一般而言，樣本均值法要比隨機投點法更有效。 nxx ,.,1baU, niiniiixfnabxgxfn1121 用上式來估計，和計算，隨機數個獨立地產生nixfuabaxuuUniiin,.,1,.,1 , 012求解定積分的算例求解定積分的算例 計算定積分 事實上，其精確解為 隨機投點法：sjtd_djf.f90，及動態圖 樣本平均值法：ybpjz_djf.f90 注注 所有計算中的隨機