蒙特卡洛方法的應(yīng)用原稿PPT課件

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

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

3、數(shù)學(xué)模擬數(shù)學(xué)模擬 計(jì)算機(jī)模擬可以反復(fù)進(jìn)行，改變系統(tǒng)的結(jié)構(gòu)和系數(shù)都比較容易。第4頁(yè)/共114頁(yè) 蒙特卡洛（蒙特卡洛（Monte CarloMonte Carlo）方法）方法是一種應(yīng)用隨機(jī)數(shù)來(lái)進(jìn)行計(jì)算機(jī)模擬的方法此方法對(duì)研究的系統(tǒng)進(jìn)行隨機(jī)觀察抽樣，通過(guò)對(duì)樣本值的觀察統(tǒng)計(jì)，求得所研究系統(tǒng)的某些參數(shù) 蒙特卡洛方法也稱(chēng)為蒙特卡洛方法也稱(chēng)為隨機(jī)模擬方法,其起源最早可以追溯到18世紀(jì)下半葉的Buffon試驗(yàn).第5頁(yè)/共114頁(yè)例 在1777年,法國(guó)學(xué)者Buffon提出用試驗(yàn)方法求圓周率鸕鬧.其原理如下：假設(shè)平面上有元數(shù)條距離為1的等矩平行線,現(xiàn)向該平面隨機(jī)地投擲一根長(zhǎng)度為KI1的針,則我們可以計(jì)算該針與任一

4、平行線相交的概率.此處隨機(jī)投針可以這樣理解z針的中心點(diǎn)與最近的平行線間的距離Z均勻地分布在區(qū)間0.1/2上,針與平行線的夾角以不管相交與否)均勻地分布在區(qū)間0,而上(見(jiàn)圖6。.于是,針與線相交的充要條件是本寸,從而針線相交概率為1第6頁(yè)/共114頁(yè)用蒙特卡洛方法進(jìn)行計(jì)算機(jī)模擬的步驟用蒙特卡洛方法進(jìn)行計(jì)算機(jī)模擬的步驟：1 設(shè)計(jì)一個(gè)邏輯框圖，即模擬模型這個(gè)框圖要正確反映系統(tǒng)各部分運(yùn)行時(shí)的邏輯關(guān)系。2 模擬隨機(jī)現(xiàn)象可通過(guò)具有各種概率分布的模擬隨機(jī)數(shù)來(lái)模擬隨機(jī)現(xiàn)象第7頁(yè)/共114頁(yè)產(chǎn)生模擬隨機(jī)數(shù)的計(jì)算機(jī)命令 在Matlab軟件中，可以直接產(chǎn)生滿足各種分布的隨機(jī)數(shù)，命令如下：2產(chǎn)生m*n階離散均勻分布的

5、隨機(jī)數(shù)矩陣：R = unidrnd(N)R = unidrnd(N,mm,nn) 當(dāng)只知道一個(gè)隨機(jī)變量取值在（a，b）內(nèi)，但不知道（也沒(méi)理由假設(shè)）它在何處取值的概率大，在何處取值的概率小，就只好用U（a，b）來(lái)模擬它。1產(chǎn)生m*n階a，b均勻分布U（a，b）的隨機(jī)數(shù)矩陣： unifrnd (a,b,m, n)產(chǎn)生一個(gè)a，b均勻分布的隨機(jī)數(shù)：unifrnd (a,b)第8頁(yè)/共114頁(yè)當(dāng)研究對(duì)象視為大量相互獨(dú)立的隨機(jī)變量之和，且其中每一種變量對(duì)總和的影響都很小時(shí)，可以認(rèn)為該對(duì)象服從正態(tài)分布。機(jī)械加工得到的零件尺寸的偏差、射擊命中點(diǎn)與目標(biāo)的偏差、各種測(cè)量誤差、人的身高、體重等，都可近似看成服從正態(tài)

6、分布。第9頁(yè)/共114頁(yè)若連續(xù)型隨機(jī)變量X的概率密度函數(shù)為 其中 0為常數(shù)，則稱(chēng)X服從參數(shù)為 的指數(shù)分指數(shù)分布布。0001)(/xxexft指數(shù)分布的期望值為 第10頁(yè)/共114頁(yè)排隊(duì)服務(wù)系統(tǒng)中顧客到達(dá)率為常數(shù)時(shí)的到達(dá)間隔、故障率為常數(shù)時(shí)零件的壽命都服從指數(shù)分布。指數(shù)分布在排隊(duì)論、可靠性分析中有廣泛應(yīng)用。注意：注意：Matlab中，產(chǎn)生參數(shù)為 的指數(shù)分布的命令為exprnd( )例例 顧客到達(dá)某商店的間隔時(shí)間服從參數(shù)為顧客到達(dá)某商店的間隔時(shí)間服從參數(shù)為1010的指的指數(shù)分布數(shù)分布 指數(shù)分布的均值為指數(shù)分布的均值為1010。 指兩個(gè)顧客到達(dá)商店的平均間隔時(shí)間是指兩個(gè)顧客到達(dá)商店的平均間隔時(shí)間是1

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

8、 n次試驗(yàn)中事件次試驗(yàn)中事件A A發(fā)生的次數(shù)。發(fā)生的次數(shù)。 頻率頻率 f=m/n f=m/n 2.2.頻率的穩(wěn)定性頻率的穩(wěn)定性 擲一枚均勻硬幣，記錄擲硬幣試驗(yàn)中頻率擲一枚均勻硬幣，記錄擲硬幣試驗(yàn)中頻率P P* *的波動(dòng)情況。的波動(dòng)情況。 R = binornd(N,P,mm,nn) 例例1 1 頻率的穩(wěn)定性頻率的穩(wěn)定性第13頁(yè)/共114頁(yè)3 3 概率的頻率定義概率的頻率定義 在一組不變的條件下，重復(fù)作在一組不變的條件下，重復(fù)作n n次試驗(yàn)，記次試驗(yàn)，記m m是是n n次次試驗(yàn)中事件試驗(yàn)中事件A A發(fā)生的次數(shù)。當(dāng)試驗(yàn)次數(shù)發(fā)生的次數(shù)。當(dāng)試驗(yàn)次數(shù)n n很大時(shí)，很大時(shí)，如果頻率如果頻率m/nm/n穩(wěn)定

9、地在某數(shù)值穩(wěn)定地在某數(shù)值p p附近擺動(dòng)，而且一附近擺動(dòng)，而且一般地說(shuō)，隨著試驗(yàn)次數(shù)的增加，這種擺動(dòng)的幅度般地說(shuō)，隨著試驗(yàn)次數(shù)的增加，這種擺動(dòng)的幅度越來(lái)越小，稱(chēng)數(shù)值越來(lái)越小，稱(chēng)數(shù)值p p為事件為事件A A在這一組不變的條件在這一組不變的條件下發(fā)生的概率，記作下發(fā)生的概率，記作P(A)=p.P(A)=p.第14頁(yè)/共114頁(yè)4 頻率的基本性質(zhì)頻率的基本性質(zhì) （1） 對(duì)任意事件對(duì)任意事件A，有，有 1)(0AP（2）1)(SP0)(P（3）若）若A1，A2，An是互不相容的，則是互不相容的，則 )()(11nkknkkAPAP 頻率定義的意義:(1) 提供了估計(jì)概率的方法;(2)提供了一種檢驗(yàn)理論正

10、確與否的準(zhǔn)則.理論依據(jù)： 大數(shù)定律 第15頁(yè)/共114頁(yè) 大量的隨機(jī)現(xiàn)象中平均結(jié)果的穩(wěn)定性 大數(shù)定律的客觀背景大量拋擲硬幣正面出現(xiàn)頻率字母使用頻率生產(chǎn)過(guò)程中的廢品率第16頁(yè)/共114頁(yè)大數(shù)定律大數(shù)定律貝努里（Bernoulli） 大數(shù)定律設(shè) nA 是 n 次獨(dú)立重復(fù)試驗(yàn)中事件 A 發(fā)生的次數(shù), p 是每次試驗(yàn)中 A 發(fā)生的概率，則0有0limpnnPAn或1limpnnPAn第17頁(yè)/共114頁(yè)在概率的統(tǒng)計(jì)定義中，事件 A 發(fā)生的頻率“ 穩(wěn)定于”事件 A 在一次試驗(yàn)中發(fā)生的概率是指：nnAnnA頻率與 p 有較大偏差pnnA是小概率事件, 因而在 n 足夠大時(shí), 可以用頻率近似代替 p . 這

11、種穩(wěn)定稱(chēng)為依概率穩(wěn)定.貝努里（貝努里（Bernoulli） 大數(shù)定律的意義大數(shù)定律的意義：第18頁(yè)/共114頁(yè)定義a 是一常數(shù)，0limaYPnn(或)1limaYPnn則稱(chēng)隨機(jī)變量序列,21nYYY依概率收斂于常數(shù) a , 記作aYnPn故pnnnPA,21nYYY是一系列隨機(jī)變量，設(shè)0有若第19頁(yè)/共114頁(yè)在 Bernoulli 定理的證明過(guò)程中， Y n 是相互獨(dú)立的服從 0-1分布的隨機(jī)變量序列 Xk 的算術(shù)平均值, Y n 依概率收斂于其數(shù)學(xué)期望 p . 結(jié)果同樣適用于服從其它分布的獨(dú)立隨機(jī)變量序列第20頁(yè)/共114頁(yè)Chebyshev 大數(shù)定律,21nXXX相互獨(dú)立，設(shè)隨機(jī)變量序

12、列(指任意給定 n 1, 相互獨(dú)立)，且具有相同的數(shù)學(xué)期望和方差nXXX,21, 2 , 1,)(,)(2kXDXEkk則0有01lim1nkknXnP或11lim1nkknXnP第21頁(yè)/共114頁(yè)定理的意義定理的意義：當(dāng) n 足夠大時(shí)，算術(shù)平均值幾乎就是一個(gè)常數(shù),可以用算術(shù)平均值近似地代替數(shù)學(xué)期望.具有相同數(shù)學(xué)期望和方差的獨(dú)立隨機(jī)變量序列的算術(shù)平均值依概率收斂于數(shù)學(xué)期望.第22頁(yè)/共114頁(yè) 例如要估計(jì)某地區(qū)的平均畝產(chǎn)量，例如要估計(jì)某地區(qū)的平均畝產(chǎn)量，要收割某些有代表性的地塊，例如要收割某些有代表性的地塊，例如n n 塊塊. . 計(jì)算其平均畝產(chǎn)量，則當(dāng)計(jì)算其平均畝產(chǎn)量，則當(dāng)n n 較大時(shí)，

13、可用較大時(shí)，可用它作為整個(gè)地區(qū)平均畝產(chǎn)量的一個(gè)估計(jì)它作為整個(gè)地區(qū)平均畝產(chǎn)量的一個(gè)估計(jì). .第23頁(yè)/共114頁(yè)辛欽大數(shù)定律辛欽大數(shù)定律 設(shè),21nXXX相互獨(dú)立，服從同一分布，且具有數(shù)學(xué)期望 E(X k) = , k= 1,2,則對(duì)任意正數(shù) 001lim1nkknXnP第24頁(yè)/共114頁(yè),21nXXX相互獨(dú)立，注3: 設(shè)隨機(jī)變量序列, 2 , 1,)(iXEkki則0有01lim1knikinXnP具有相同的分布，且記knikiMXn11第25頁(yè)/共114頁(yè)11nPA),(21kAAAgnP),(21kg則則22nPAknPkA),(21kxxxg連續(xù)，若第26頁(yè)/共114頁(yè) 大數(shù)定律以嚴(yán)格

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

15、114頁(yè)function 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)在在MatlabMatlab中編輯中編輯.m.m文件輸入以下命令：文件輸入以下命令：第29頁(yè)/共114頁(yè)在在MatlabMatlab命令行中輸入以下命令：命令行中輸入以下命令：liti1(1,0.5,1000)第30頁(yè)/共114頁(yè)在在MatlabMatlab命令行中輸入以下命令：命令行中輸入以下命令：lit

16、i1(1,0.5,10000)第31頁(yè)/共114頁(yè)練習(xí)練習(xí) 頻率的穩(wěn)定性頻率的穩(wěn)定性1 1 事件的頻率事件的頻率 R = binornd(N,P,mm,nn) 在一組不變的條件下，重復(fù)作在一組不變的條件下，重復(fù)作n n次試驗(yàn)，次試驗(yàn)，記記m m是是n n次試驗(yàn)中事件次試驗(yàn)中事件A A發(fā)生的次數(shù)。發(fā)生的次數(shù)。 頻率頻率 f=m/n f=m/n 2.2.頻率的穩(wěn)定性頻率的穩(wěn)定性 練習(xí)練習(xí)擲一枚不均勻硬幣，正面出現(xiàn)概率為擲一枚不均勻硬幣，正面出現(xiàn)概率為0.30.3，記錄前記錄前10001000次擲硬幣試驗(yàn)中正面頻率的波動(dòng)次擲硬幣試驗(yàn)中正面頻率的波動(dòng)情況，并畫(huà)圖。情況，并畫(huà)圖。 第32頁(yè)/共114頁(yè)在

17、在MatlabMatlab命令行中輸入以下命令：命令行中輸入以下命令：liti1(1,0.3,1000)第33頁(yè)/共114頁(yè)例例2 2 擲擲兩枚不均勻硬幣，每枚正面出現(xiàn)概率兩枚不均勻硬幣，每枚正面出現(xiàn)概率為為0.40.4，記錄前，記錄前10001000次擲硬幣試驗(yàn)中兩枚都為次擲硬幣試驗(yàn)中兩枚都為正面頻率的波動(dòng)情況，并畫(huà)圖。正面頻率的波動(dòng)情況，并畫(huà)圖。 在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)*ran

18、dnum(2,i); pro(i)=a/i;end pro=pro，num=1:mm;plot(num,pro)第34頁(yè)/共114頁(yè) 熊宇樂(lè) 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; end y=y; num=1:1000; plot(num,y)第35頁(yè)/共114頁(yè)孟亞function bino (n,p,m)x=binornd(n,p,1,m);y=binornd(n,p,1,m);for i=1:

19、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)第36頁(yè)/共114頁(yè)liti2(1,0.4,100)第37頁(yè)/共114頁(yè)liti2(1,0.4,10000)第38頁(yè)/共114頁(yè)在一袋中有在一袋中有10 10 個(gè)相同的球，分別標(biāo)有號(hào)碼個(gè)相同的球，分別標(biāo)有號(hào)碼1,2,1,2,10,10。每次任取一個(gè)球，記錄其號(hào)碼。每次任取一個(gè)球，記錄其號(hào)碼后放回袋中，再任取下一個(gè)。這種取法叫做后放回袋中，再任取下一個(gè)。這種取法叫做“有放回抽取有放回抽取”。今有放回抽取。今有放回抽取

20、3 3個(gè)球，求個(gè)球，求這這3 3個(gè)球的號(hào)碼均為偶數(shù)的概率。個(gè)球的號(hào)碼均為偶數(shù)的概率。（用頻率估計(jì)概率） 例3:解：令A(yù)=有放回抽取3個(gè)球，求這3個(gè)球的號(hào)碼均為偶數(shù)=(2,2,2)，(2,2,4)，.，(10,10,10)81)(AP第39頁(yè)/共114頁(yè)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=f

21、rq/mm第40頁(yè)/共114頁(yè)例例4 4 兩盒火柴，每盒兩盒火柴，每盒2020根。每次隨機(jī)在任一根。每次隨機(jī)在任一盒中取出一根火柴。問(wèn)其中一盒中火柴被取盒中取出一根火柴。問(wèn)其中一盒中火柴被取完而另一盒中至少還有完而另一盒中至少還有5 5根火柴的概率有多大？根火柴的概率有多大？（用頻率估計(jì)概率）（用頻率估計(jì)概率） 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第41頁(yè)/共114頁(yè)

22、function proguji=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 (a120)&(a2=5 frq=frq+1; end% a1=a1,a2=a2,frq % pauseend proguji=frq/mm第42頁(yè)/共114頁(yè)二二. 幾何概率幾何概率1.1.定義 向任一可度量區(qū)域向任一可度量區(qū)域G G內(nèi)投一點(diǎn)，如果所投內(nèi)投一點(diǎn)，如果所投的點(diǎn)落在的點(diǎn)落在G G

23、中任意可度量區(qū)域中任意可度量區(qū)域g g內(nèi)的可能內(nèi)的可能性與性與g g的度量成正比，而與的度量成正比，而與g g的位置和形的位置和形狀無(wú)關(guān)，則稱(chēng)這個(gè)隨機(jī)試驗(yàn)為幾何型隨狀無(wú)關(guān)，則稱(chēng)這個(gè)隨機(jī)試驗(yàn)為幾何型隨機(jī)試驗(yàn)。或簡(jiǎn)稱(chēng)為幾何概型。機(jī)試驗(yàn)?；蚝?jiǎn)稱(chēng)為幾何概型。第43頁(yè)/共114頁(yè)2. 概率計(jì)算概率計(jì)算 1. P(A)=A的度量的度量/S的度量的度量?jī)扇思s定于兩人約定于12點(diǎn)到點(diǎn)到1點(diǎn)到某地會(huì)面，先點(diǎn)到某地會(huì)面，先到者等到者等20分鐘后離去，試求兩人能會(huì)面分鐘后離去，試求兩人能會(huì)面的概率？的概率？ 例例5解：設(shè)x,y分別為甲、乙到達(dá)時(shí)刻(分鐘)令A(yù)=兩人能會(huì)面兩人能會(huì)面=(x,y)|x-y|20，x60,

24、60,y y6060P(A)=A的面積/S的面積=（602-402）/602=5/9=0.5556第44頁(yè)/共114頁(yè)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;proguji=0;for ii=1:mm if abs(randnum(ii,1)=20 frq=frq+1; endendproguji=frq/mmliti5(10000)proguji = 0.5557第45頁(yè)/共114頁(yè)例

25、6 6在我方某前沿防守地域，敵人以一個(gè)炮排（含兩門(mén)火炮）為單位對(duì)我方進(jìn)行干擾和破壞為躲避我方打擊，敵方對(duì)其陣地進(jìn)行了偽裝并經(jīng)常變換射擊地點(diǎn) 經(jīng)過(guò)長(zhǎng)期觀察發(fā)現(xiàn)，我方指揮所對(duì)敵方目標(biāo)的指示有50是準(zhǔn)確的，而我方火力單位，在指示正確時(shí)，有1/3的射擊效果能毀傷敵人一門(mén)火炮，有1/6的射擊效果能全部消滅敵人 現(xiàn)在希望能用某種方式把我方將要對(duì)敵人實(shí)施的1打擊結(jié)果顯現(xiàn)出來(lái)，利用頻率穩(wěn)定性，確定有效射擊的概率第46頁(yè)/共114頁(yè)分析分析： 這是一個(gè)概率問(wèn)題，可以通過(guò)理論計(jì)算得到相應(yīng)的概率和期望值. 為了能顯示我方射擊的過(guò)程，現(xiàn)采用模擬的方式。第47頁(yè)/共114頁(yè) 需要模擬出以下兩件事： 1. 1. 問(wèn)題分析

26、問(wèn)題分析1 1 觀察所對(duì)目標(biāo)的指示正確與否觀察所對(duì)目標(biāo)的指示正確與否模擬試驗(yàn)有兩種結(jié)果，每一種結(jié)果出現(xiàn)的概率都是1/2 因此，可用投擲一枚硬幣的方式予以確可用投擲一枚硬幣的方式予以確定定，當(dāng)硬幣出現(xiàn)正面時(shí)為指示正確，反之為不正確第48頁(yè)/共114頁(yè)2 2 當(dāng)指示正確時(shí)，我方火力單位的射擊結(jié)當(dāng)指示正確時(shí)，我方火力單位的射擊結(jié)果情況果情況 模擬試驗(yàn)有三種結(jié)果：毀傷一門(mén)火炮的可能性為1/3(即2/6)，毀傷兩門(mén)的可能性為1/6，沒(méi)能毀傷敵火炮的可能性為1/2(即3/6) 這時(shí)可用投擲骰子的方法來(lái)確定可用投擲骰子的方法來(lái)確定：如果出現(xiàn)的是、三個(gè)點(diǎn)：則認(rèn)為沒(méi)能擊中敵人；如果出現(xiàn)的是、點(diǎn)：則認(rèn)為毀傷敵人一

27、門(mén)火炮；若出現(xiàn)的是點(diǎn)：則認(rèn)為毀傷敵人兩門(mén)火炮第49頁(yè)/共114頁(yè)2. 2. 符號(hào)假設(shè)符號(hào)假設(shè)i：要模擬的打擊次數(shù)； k1：沒(méi)擊中敵人火炮的射擊總數(shù)； k2：擊中敵人一門(mén)火炮的射擊總數(shù)；k3：擊中敵人兩門(mén)火炮的射擊總數(shù)E：有效射擊比率； 第50頁(yè)/共114頁(yè)3. 3. 模擬框圖模擬框圖初始化:i=0,k1=0,k2=0,k3=0i=i+1骰子點(diǎn)數(shù)?k1=k1+1k2=k2+1k3=k3+1k1=k1+1imm?E=(k2+k3)/mm 停止硬幣正面?YNNY1，2，34，56第51頁(yè)/共114頁(yè)function liti6(p,mm)efreq=zeros(1,mm);randnum1 = bi

28、nornd(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)在在MatlabMatlab中編輯中編輯.m.m文件輸入以下命令：文件輸入以下命令：第52頁(yè)/共114頁(yè)在在MatlabMatlab命令行中輸入以下

29、命令：命令行中輸入以下命令：liti6(0.5,2000)第53頁(yè)/共114頁(yè)在在MatlabMatlab命令行中輸入以下命令：命令行中輸入以下命令：liti6(0.5,20000)第54頁(yè)/共114頁(yè)5. 5. 理論計(jì)算理論計(jì)算第55頁(yè)/共114頁(yè)6. 6. 結(jié)果比較結(jié)果比較 模擬結(jié)果與理論計(jì)算近似一致，能更模擬結(jié)果與理論計(jì)算近似一致，能更加真實(shí)地表達(dá)實(shí)際戰(zhàn)斗動(dòng)態(tài)過(guò)程加真實(shí)地表達(dá)實(shí)際戰(zhàn)斗動(dòng)態(tài)過(guò)程 第56頁(yè)/共114頁(yè)3. 某廠生產(chǎn)的燈泡能用1000小時(shí)的概率為0.8, 能用1500小時(shí)的概率為0.4 , 求已用1000小時(shí)的燈泡能用到1500小時(shí)的概率2. 在一袋中有10 個(gè)相同的球，分別標(biāo)

30、有號(hào)碼1,2,10。今有放回任取兩個(gè)球，求取得的第一個(gè)球號(hào)碼為奇數(shù)，第二個(gè)球的號(hào)碼為偶數(shù)的概率。1. 擲三枚不均勻硬幣，每枚正面出現(xiàn)概率為0.3，記錄前1000次擲硬幣試驗(yàn)中至少兩枚都為正面頻率的波動(dòng)情況，并畫(huà)圖。 作業(yè):第57頁(yè)/共114頁(yè) 某廠生產(chǎn)的燈泡能用1000小時(shí)的概率為0.8, 能用1500小時(shí)的概率為0.4 , 求已用1000小時(shí)的燈泡能用到1500小時(shí)的概率解解 令 A 燈泡能用到1000小時(shí) B 燈泡能用到1500小時(shí)所求概率為)()(APABPABPAB 218 . 04 . 0)()(APBP例例 第58頁(yè)/共114頁(yè)例:在一袋中有在一袋中有10 個(gè)相同的球，分別標(biāo)有號(hào)個(gè)

31、相同的球，分別標(biāo)有號(hào)碼碼1,2,10。今不放回任取兩個(gè)球，求。今不放回任取兩個(gè)球，求取得的第一個(gè)球號(hào)碼為奇數(shù)，第二個(gè)球取得的第一個(gè)球號(hào)碼為奇數(shù)，第二個(gè)球的號(hào)碼為偶數(shù)的概率。的號(hào)碼為偶數(shù)的概率。解：令解：令A(yù)=A=抽取抽取2 2個(gè)球，第一個(gè)球號(hào)碼為奇數(shù)，第二個(gè)球個(gè)球，第一個(gè)球號(hào)碼為奇數(shù)，第二個(gè)球的號(hào)碼為偶數(shù)的號(hào)碼為偶數(shù)=(1,2)=(1,2)，(1,4)(1,4)，. .，(9,10)(9,10)185)(AP第59頁(yè)/共114頁(yè)2021-11-19圓的面積單位圓的面積等于 計(jì)算第一象限內(nèi)的單位圓的面積，方法為：把它分成n個(gè)窄的曲邊梯形，計(jì)算S大，S小，其中n可以為1000, 10000, .實(shí)

32、驗(yàn)二: 的計(jì)算第60頁(yè)/共114頁(yè)2021-11-19數(shù)值積分法4)1(141102102dxxdxx第61頁(yè)/共114頁(yè)2021-11-19梯形公式inabaxnabxfbfafdxxfiniiba,)(2/)()()(11第62頁(yè)/共114頁(yè)2021-11-19辛普森公式nabxfxfbfafdxxfniiniiba6)(2)(4)()()(11105 . 0第63頁(yè)/共114頁(yè)2021-11-19無(wú)窮級(jí)數(shù)法 arctg x =x-x3/3+x5/5-x7/7+x9/9-. /4 = arctg 1 =1-1/3+1/5-1/7+1/9- 收斂太慢！ |x| 應(yīng)當(dāng)比 1 小很多，級(jí)數(shù)收斂才

33、快。 /4 = arctg 1/2 + arctg 1/3 /4 = 4 arctg 1/5 - arctg 1/239第64頁(yè)/共114頁(yè)2021-11-19 某次隨機(jī)投點(diǎn)的結(jié)果某次隨機(jī)投點(diǎn)的結(jié)果0.20.40.60.810.20.40.60.81第65頁(yè)/共114頁(yè) Monte-Carlo 方法謝文賢西北工業(yè)大學(xué)理學(xué)院 高性能計(jì)算中心2004-05-142004-05-14第66頁(yè)/共114頁(yè)Monte-Carlo 方法 一、MC 的起源和發(fā)展 Buffon試驗(yàn) 排隊(duì)系統(tǒng)模擬：M/G/1M/G/1排隊(duì)系統(tǒng) 二、MC 的原理 三、隨機(jī)數(shù)的產(chǎn)生原理與IMSL庫(kù) 均勻分布U(0,1)的隨機(jī)數(shù)的產(chǎn)

34、生 其他各種分布的隨機(jī)數(shù)的產(chǎn)生 隨機(jī)過(guò)程模擬 四、MC的應(yīng)用舉例 定積分的MC計(jì)算 隨機(jī)微分方程的數(shù)值模擬 五、EM算法及其MCMC方法第67頁(yè)/共114頁(yè)一、MC 的起源和發(fā)展 隨機(jī)模擬方法，也稱(chēng)為Monte Carlo方法，是一種基于“隨機(jī)數(shù)”的計(jì)算方法。這一方法源于美國(guó)在第一次世界大戰(zhàn)進(jìn)行的研制原子彈的“曼哈頓計(jì)劃”。該計(jì)劃的主持人之一、數(shù)學(xué)家馮諾伊曼用馳名世界的賭城摩納哥的Monte Carlo來(lái)命名這種方法，為它蒙上了一層神秘色彩。馮諾伊曼是公理化方法和計(jì)算機(jī)體系的領(lǐng)袖人物，Monte Carlo方法也是他的功勞。 第68頁(yè)/共114頁(yè) 事實(shí)上，Monte Carlo方法的基本思想很

35、早以前就被人們所發(fā)現(xiàn)和利用。早在17世紀(jì)，人們就知道用事件發(fā)生的“頻率”來(lái)決定事件的“概率”。18世紀(jì)下半葉的法國(guó)學(xué)者Buffon提出用投針試驗(yàn)的方法來(lái)確定圓周率的值。這個(gè)著名的Buffon試驗(yàn)是Monte Carlo方法的最早的嘗試！第69頁(yè)/共114頁(yè) 歷史上曾有幾位學(xué)者相繼做過(guò)這樣的試驗(yàn)。不過(guò)呢，他們的試驗(yàn)是費(fèi)時(shí)費(fèi)力的，同時(shí)精度不夠高，實(shí)施起來(lái)也很困難。然而，隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展，人們不需要具體實(shí)施這些試驗(yàn)，而只要在計(jì)算機(jī)上進(jìn)行大量的、快速的模擬試驗(yàn)就可以了。第70頁(yè)/共114頁(yè) 在大眾的心目中，科學(xué)的代言人是心不在焉的牛頓或者爆炸式發(fā)型的愛(ài)因斯坦，但這只是傳統(tǒng)形象，比他們更了解現(xiàn)代

36、計(jì)算技術(shù)的馮諾伊曼是個(gè)衣著考究，風(fēng)度翩翩的人物，他說(shuō)：純粹數(shù)學(xué)和應(yīng)用數(shù)學(xué)的許多分支非常需要計(jì)算工具，用以打破目前由于純粹分析的研究方法不能解決非線性問(wèn)題而形成的停滯狀態(tài)。 Monte Carlo方法是現(xiàn)代計(jì)算技術(shù)的最為杰出的成果之一，它在工程領(lǐng)域的作用是不可比擬的。 第71頁(yè)/共114頁(yè)Buffon試驗(yàn) 假設(shè)平面上有無(wú)數(shù)條距離為1的等距平行線，現(xiàn)向該平面隨機(jī)投擲一根長(zhǎng)度為 的針( )， 則我們可計(jì)算該針與任一平行線相交的概率。這里，隨機(jī)投針指的是：針的中心點(diǎn)與最近的平行線間的距離 均勻的分布在區(qū)間 。上，針與平行線的夾角 （不管相交與否）均勻的分布在區(qū)間 上。 因此，針與線相交的充要條件是 l

37、1l21 , 0 x, 021sinx第72頁(yè)/共114頁(yè)Buffon試驗(yàn) 從而針線相交的概率為 根據(jù)上式，若我們做大量的投針試驗(yàn)并記錄針與線相交的次數(shù)，則由大數(shù)定理可以估計(jì)出針線相交的概率 ，從而得到 的估計(jì)值。 針與線的位置關(guān)系：ldxdlXPpl22sin20sin20 xlp第73頁(yè)/共114頁(yè)排隊(duì)系統(tǒng)模擬：M/G/1M/G/1排隊(duì)系統(tǒng)服務(wù)規(guī)則：先到先服務(wù)假設(shè)：（）顧客到達(dá)遵循Poisson分布；（）服務(wù)時(shí)間服從一般分布；（）到達(dá)間隔與服務(wù)時(shí)間相互獨(dú)立1第74頁(yè)/共114頁(yè)排隊(duì)系統(tǒng)模擬：M/G/1M/G/1排隊(duì)系統(tǒng)關(guān)心的指標(biāo)：（）時(shí)刻t時(shí)，系統(tǒng)中的顧客數(shù)；即隊(duì)長(zhǎng)的分布；（）顧客的等待時(shí)

38、間；（）服務(wù)的忙碌程度；（） 用最樸素的Monte-Carlo方法可以得到這些指標(biāo)的估計(jì)第75頁(yè)/共114頁(yè)二、MC 的原理 應(yīng)用Monte Carlo方法求解工程技術(shù)問(wèn)題可以分為兩類(lèi)： 確定性問(wèn)題 隨機(jī)性問(wèn)題第76頁(yè)/共114頁(yè)確定性系統(tǒng)隨機(jī)性系統(tǒng)模擬自然界Monte-Carlo模擬，即隨機(jī)模擬（重復(fù)“試驗(yàn)”）重復(fù)試驗(yàn)計(jì)算機(jī)模擬第77頁(yè)/共114頁(yè)思路： 1、 針對(duì)實(shí)際問(wèn)題建立一個(gè)簡(jiǎn)單且便于實(shí)現(xiàn)的概率統(tǒng)計(jì)模型，使問(wèn)題的解對(duì)應(yīng)于該模型中隨機(jī)變量的概率分布或其某些數(shù)字特征，比如，均值和方差等。所構(gòu)造的模型在主要特征參量方面要與實(shí)際問(wèn)題或系統(tǒng)相一致的。 2、 根據(jù)模型中各個(gè)隨機(jī)變量的分布，在計(jì)算機(jī)

39、上產(chǎn)生隨機(jī)數(shù)，實(shí)現(xiàn)一次模擬過(guò)程所需的足夠數(shù)量的隨機(jī)數(shù)。通常先產(chǎn)生均勻分布的隨機(jī)數(shù)，然后生成服從某一分布的隨機(jī)數(shù)，再進(jìn)行隨機(jī)模擬試驗(yàn)。 第78頁(yè)/共114頁(yè)思路： 3、 根據(jù)概率模型的特點(diǎn)和隨機(jī)變量的分布特性，設(shè)計(jì)和選取合適的抽樣方法，并對(duì)每個(gè)隨機(jī)變量進(jìn)行抽樣（包括直接抽樣、分層抽樣、相關(guān)抽樣、重要抽樣等）。 4、 按照所建立的模型進(jìn)行仿真試驗(yàn)、計(jì)算，求出問(wèn)題的隨機(jī)解。 5、 統(tǒng)計(jì)分析模擬試驗(yàn)結(jié)果，給出問(wèn)題的估計(jì)以及其精度估計(jì)。 6、 必要時(shí)，還應(yīng)改進(jìn)模型以降低估計(jì)方差和減少試驗(yàn)費(fèi)用，提高模擬計(jì)算的效率。 第79頁(yè)/共114頁(yè) 收斂性: : 由大數(shù)定律, Monte-Carlo模擬的收斂是以概率

40、而言的 誤差: : 用頻率估計(jì)概率時(shí)誤差的估計(jì)，可由中心極限定理，給定置信水平 的條件下，有： 模擬次數(shù): :由誤差公式得NU2/1|)(XgVar22/1)(UN第80頁(yè)/共114頁(yè)三、隨機(jī)數(shù)的產(chǎn)生原理與IMSL庫(kù) 隨機(jī)數(shù)的產(chǎn)生是實(shí)現(xiàn)MC計(jì)算的先決條件。而大多數(shù)概率分布的隨機(jī)數(shù)的產(chǎn)生都是基于均勻分布U(0,1)的隨機(jī)數(shù)。 首先，介紹服從均勻分布U(0,1)的隨機(jī)數(shù)的產(chǎn)生方法。 其次，介紹服從其他各種分布的隨機(jī)數(shù)的產(chǎn)生方法。以及服從正態(tài)分布的隨機(jī)數(shù)的產(chǎn)生方法。 接下來(lái)，介紹隨機(jī)過(guò)程模擬。 最后，關(guān)于隨機(jī)數(shù)的幾點(diǎn)注。第81頁(yè)/共114頁(yè)均勻分布U(0,1)的隨機(jī)數(shù)的產(chǎn)生 產(chǎn)生均勻分布的標(biāo)準(zhǔn)算法在

41、很多高級(jí)計(jì)算機(jī)語(yǔ)言的書(shū)都可以看到。算法簡(jiǎn)單，容易實(shí)現(xiàn)。使用者可以自己手動(dòng)編程實(shí)現(xiàn)。IMSL統(tǒng)計(jì)庫(kù)中也提供給我們用于產(chǎn)生均勻分布的各種函數(shù)。我們的重點(diǎn)是怎樣通過(guò)均勻分布產(chǎn)生服從其他分布的隨機(jī)數(shù)。因此，直接使用IMSL提供的可靠安全的標(biāo)準(zhǔn)函數(shù)，當(dāng)然不用費(fèi)事了。 第82頁(yè)/共114頁(yè)IMSL庫(kù)中的函數(shù)使用 RNSET: 種子的設(shè)定 CALL RNSET (ISEED) RNOPT: 產(chǎn)生器的類(lèi)型的設(shè)定 CALL RNOPT (IOPT) RNUN/DRNUN: 產(chǎn)生均勻分布的隨機(jī)數(shù) CALL RNUN (NR, R) 第83頁(yè)/共114頁(yè)其他各種分布的隨機(jī)數(shù)的產(chǎn)生 基本方法有如下三種： 逆變換法 合

42、成法 篩選法 第84頁(yè)/共114頁(yè)逆變換法 設(shè)隨機(jī)變量 的分布函數(shù)為 ，定義 定理 設(shè)隨機(jī)變量 服從 上的均勻分布，則 的分布函數(shù)為 。 因此，要產(chǎn)生來(lái)自 的隨機(jī)數(shù)，只要先產(chǎn)生來(lái)自 的隨機(jī)數(shù)，然后計(jì)算 即可。 其步驟為 X xF 10 ,:inf1yyxFxyFU) 1 , 0( UFX1 xF xF1 , 0U uF1 uFxuU11 , 0計(jì)算，抽取由第85頁(yè)/共114頁(yè)合成法 合成法的應(yīng)用最早見(jiàn)于Butlter 的書(shū)中。 構(gòu)思如下： 如果 的密度函數(shù) 難于抽樣，而 關(guān)于 的條件密度函數(shù) 以及 的密度函數(shù) 均易于抽樣，則 的隨機(jī)數(shù)可如下產(chǎn)生： 可以證明由此得到 的服從 。X xpXYyxp

43、Y ygX xyxpyygY抽取由條件分布，抽取的分布由X xp第86頁(yè)/共114頁(yè)篩選抽樣 假設(shè)我們要從 抽樣，如果可以將 表示成 ，其中 是一個(gè)密度函數(shù)且易于抽樣，而 ， 是常數(shù)，則 的抽樣可如下進(jìn)行： 定理 設(shè) 的密度函數(shù) ，且 ，其中 ， ， 是一個(gè)密度函數(shù)。令 和 分別服從 和 ，則在 的條件 下， 的條件密度為 xp xgxhcxp xph 10 xg1cX 。，回到如果，停止，則如果，抽取，由抽取由1321 , 01yguyxyguyyhuUX xp xgxhcxph 10 xg1cUY1 , 0U yh YgU xpYgUxpYY第87頁(yè)/共114頁(yè)標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)數(shù) 的隨機(jī)

44、數(shù)產(chǎn)生方法很多。簡(jiǎn)要介紹三種。 法1、 變換法（Box 和Muller 1958） 設(shè) ， 是獨(dú)立同分布的 變量，令 則 與 獨(dú)立，均服從標(biāo)準(zhǔn)正態(tài)分布。 法2、 結(jié)合合成法與篩選法。（略） 法3 3、 近似方法（利用中心極限定理） 即用 個(gè) 變量產(chǎn)生一個(gè) 變量。 其中 是抽自 的隨機(jī)數(shù)， 可近似為一個(gè) 變量。1 , 0N1U2U1 , 0U2122112sinln22cosln22121UUXUUX1X2Xn1 , 0U1 , 0N2112unxiu1 , 0Uniiunu111 , 0N第88頁(yè)/共114頁(yè)IMSL庫(kù)中的函數(shù)使用 RNSET: 種子的設(shè)定 CALL RNSET (ISEED)

45、 RNOPT: 產(chǎn)生器的類(lèi)型的設(shè)定 CALL RNOPT (IOPT) RNNOA: 產(chǎn)生服從標(biāo)準(zhǔn)正態(tài)分布的偽隨機(jī)數(shù) CALL RNNOA (NR, R) 第89頁(yè)/共114頁(yè)隨機(jī)過(guò)程模擬 高斯白噪聲的產(chǎn)生： 一種簡(jiǎn)單有效的模擬高斯白噪聲過(guò)程的方法是由獨(dú)立的單位正態(tài)隨機(jī)序列 ，按照如圖所示的方式連接，其中 式中 為白噪聲的強(qiáng)度。kktDkD第90頁(yè)/共114頁(yè) 由于正態(tài)隨機(jī)數(shù)的獨(dú)立性，當(dāng) 很小時(shí)，按(a)、(b)所構(gòu)成的過(guò)程的相關(guān)時(shí)間很短，從而具有很高的截止頻率。 當(dāng)然， 過(guò)小，樣本的計(jì)算量過(guò)大。因此，選取 適當(dāng)小即可。ttt第91頁(yè)/共114頁(yè)關(guān)于隨機(jī)數(shù)的幾點(diǎn)注 注1 1 由于均勻分布的隨機(jī)

46、數(shù)的產(chǎn)生總是采用某個(gè)確定的模型進(jìn)行的，從理論上講，總會(huì)有周期現(xiàn)象出現(xiàn)的。初值確定后，所有隨機(jī)數(shù)也隨之確定，并不滿足真正隨機(jī)數(shù)的要求。因此通常把由數(shù)學(xué)方法產(chǎn)生的隨機(jī)數(shù)成為偽隨機(jī)數(shù)。 但其周期又相當(dāng)長(zhǎng)，在實(shí)際應(yīng)用中幾乎不可能出現(xiàn)。因此，這種由計(jì)算機(jī)產(chǎn)生的偽隨機(jī)數(shù)可以當(dāng)作真正的隨機(jī)數(shù)來(lái)處理。 注2 2 應(yīng)對(duì)所產(chǎn)生的偽隨機(jī)數(shù)作各種統(tǒng)計(jì)檢驗(yàn)，如獨(dú)立性檢驗(yàn)，分布檢驗(yàn)，功率譜檢驗(yàn)等等。 第92頁(yè)/共114頁(yè)四、MC的應(yīng)用舉例 例一、定積分的MC計(jì)算 隨機(jī)投點(diǎn)法 樣本平均值法 幾種降低估計(jì)方差的MC方法 例二、 隨機(jī)微分方程的數(shù)值模擬 系統(tǒng)的可靠性數(shù)值模擬計(jì)算問(wèn)題第93頁(yè)/共114頁(yè)定積分的MC計(jì)算 事實(shí)上，

47、不少的統(tǒng)計(jì)問(wèn)題，如計(jì)算概率、各階距等，最后都?xì)w結(jié)為定積分的近似計(jì)算問(wèn)題。 下面考慮一個(gè)簡(jiǎn)單的定積分 為了說(shuō)明問(wèn)題，我們首先介紹兩種求 的簡(jiǎn)單的MC方法，然后給出幾種較為復(fù)雜而更有效的MC方法。 dxxfba第94頁(yè)/共114頁(yè) 在計(jì)算積分上，MC的實(shí)用場(chǎng)合是計(jì)算重積分 其中 是 維空間的點(diǎn)，當(dāng) 較大時(shí)，用MC方法比一般的數(shù)值方法有優(yōu)點(diǎn)，主要是它的誤差與維數(shù) 無(wú)關(guān)。 dPPgIkPmmm第95頁(yè)/共114頁(yè)第96頁(yè)/共114頁(yè)隨機(jī)投點(diǎn)法 方法簡(jiǎn)述： 設(shè) ，有限， ， ，并設(shè) 是在 上均勻分布的二維隨機(jī)變量，其聯(lián)合密度函數(shù)為 。則易見(jiàn) 是 中 曲線下方的面積。假設(shè)我們向 中進(jìn)行隨機(jī)投點(diǎn)，若點(diǎn)落在

48、下方，(即 稱(chēng)為中的，否則稱(chēng)為不中，則點(diǎn)中的概率為 。若我們進(jìn)行了 次投點(diǎn)，其中 次中的，則用頻率來(lái)估計(jì)概率 。即 。ab Mxf0Mybxayx0 ,:,YX,MybxaIabM0,1 dxxfba xfy xfy xfy abMpn0npnnp0第97頁(yè)/共114頁(yè) 那么我們可以得到 的一個(gè)估計(jì) 具體試驗(yàn)步驟為nnabM01 用上式來(lái)估計(jì)的個(gè)數(shù)統(tǒng)計(jì)，和，計(jì)算，隨機(jī)數(shù)，個(gè)獨(dú)立地產(chǎn)生0,.,11 , 02nyxfxfMvyabuaxnivuUniiiiiiiii第98頁(yè)/共114頁(yè) 注1 1 隨機(jī)投點(diǎn)法的思想簡(jiǎn)單明了，且每次投點(diǎn)結(jié)果服從二項(xiàng)分布，故 ，其中 注2 2 可證 是 的無(wú)偏估計(jì)。若用

49、估計(jì)的標(biāo)準(zhǔn)差來(lái)衡量其精度，則估計(jì) 的精度的階為 。 注3 3 這里，定積分的解，就對(duì)應(yīng)我們選定的隨機(jī)變量的概率值。 pnbn,0abMp1121n第99頁(yè)/共114頁(yè)樣本平均值法 基本原理：對(duì)積分 ，設(shè) 是 上的一個(gè)密度函數(shù)，改寫(xiě) 可見(jiàn)，任一積分均可以表示為某個(gè)隨機(jī)變量（函數(shù)）的期望。由矩法，若有 個(gè)來(lái)自 的觀測(cè)值，則可給出 的一個(gè)矩估計(jì)。 最簡(jiǎn)單的，若 ，有限，可取 。 dxxfba xgba, XgXfEdxxgxgxfban xgab abxg /1第100頁(yè)/共114頁(yè) 設(shè) 是來(lái)自 的隨機(jī)數(shù)，則 的一個(gè)估計(jì)為 具體步驟為 注 可證 是 的無(wú)偏估計(jì)。一般而言，樣本均值法要比隨機(jī)投點(diǎn)法更有效。 nxx ,.,1baU, niiniiixfnabxgxfn1121 用上式來(lái)估計(jì)，和計(jì)算，隨機(jī)數(shù)個(gè)獨(dú)立地產(chǎn)生nixfuabaxuuUniiin,.,1,.,1 , 012第101頁(yè)/共114頁(yè)求解定積分的算例 計(jì)算定積分 事實(shí)上，其精確解為 隨機(jī)投點(diǎn)法：sjtd_djf.f90，及動(dòng)態(tài)圖 樣本平均值法：ybpjz_djf.f90 注 所有計(jì)算中的隨機(jī)數(shù)的產(chǎn)生均來(lái)自IMSL庫(kù)。樣本數(shù)為10000。增加樣本數(shù)目，可提