Monte Carlo（蒙特卡洛方法）

1、Monte Carlo Simulation Methods(蒙特卡羅模擬方法) 一. 概述與思想 二. 隨機(jī)數(shù)的生成. 三. 實(shí)例-港口模型 四. 作業(yè)引例：電梯系統(tǒng) 然而這種做法可能是難以接受的，因?yàn)樵谑占y(tǒng)計(jì)數(shù)據(jù)時(shí)要再三驚擾乘客，并且電梯運(yùn)行模式的不斷變化也會(huì)使乘客感到迷惑 我們可以提出若干供選擇的電梯運(yùn)行模式，如設(shè)定停偶數(shù)層、奇數(shù)層的電梯或直達(dá)電梯理論上，對(duì)每種供選擇的模式都能夠做若干次試驗(yàn)，以確定哪一種模式能為那些要到達(dá)特定樓層的乘客提供最好的服務(wù)。一 概述與思想 所以在某些情況下，對(duì)對(duì)象的行為進(jìn)行直接觀測(cè)或重復(fù)試驗(yàn)可能是不可行的。 與此有關(guān)的另一個(gè)問(wèn)題是大城市交通控制系統(tǒng)可供選擇的

2、運(yùn)行模式的檢驗(yàn)，為了做試驗(yàn)而不停地改變單行道的交通方向和配置交通信號(hào)將是不現(xiàn)實(shí)的 在對(duì)象的行為不能做分析性的解釋，或數(shù)據(jù)無(wú)法直接收集的情況下，建模者可以用某種方式間接地模擬其行為，試驗(yàn)所研究的供選擇的各種方案，以估計(jì)它們?cè)鯓佑绊憣?duì)象的行為，然后收集數(shù)據(jù)來(lái)確定哪種方案是最好的 例如，為了得到一艘擬建造的潛艇受到的阻力，造一個(gè)原型是不可行的，我們可以按比例建一個(gè)模型，去模擬實(shí)際的潛艇的行為 這里將研究另外一種形式的模擬蒙特卡羅(MonteCarlo)模擬，一般是借助于計(jì)算機(jī)完成的 蒙特卡洛(Monte Carlo)方法，或稱計(jì)算機(jī)隨機(jī)模擬方法，是一種基于“隨機(jī)數(shù)”的計(jì)算方法。 這一方法源于美國(guó)在第

3、二次世界大戰(zhàn)中研制原子彈的“曼哈頓計(jì)劃”。該計(jì)劃的主持人之一、數(shù)學(xué)家馮諾伊曼用馳名世界的賭城摩納哥的Monte Carlo來(lái)命名這種方法，為它蒙上了一層神秘色彩。從Buffon(蒲豐)投針問(wèn)題談起 試驗(yàn)者時(shí)間(年)針長(zhǎng)投針次數(shù)相交次數(shù)的估計(jì)值Wolf18500.80500025323.15956Smith18550.60320412183.15665Fox18840.7510304893.15951Lazzarini19250.83340818083.141592921 確定行為的模擬例：曲線下的面積 本節(jié)以曲線下的面積為例說(shuō)明蒙特卡羅模擬在確定行為建模中的應(yīng)用 下面的算法給出了用蒙特卡羅方法

4、求曲線下面積的計(jì)算機(jī)模擬的計(jì)算格式 在給定區(qū)間上曲線y=cosx下面積的真值是2注意到即使對(duì)于產(chǎn)生的相當(dāng)多的點(diǎn)數(shù)，誤差也是可觀的對(duì)單變量函數(shù)，一般說(shuō)來(lái)，蒙特卡羅方法無(wú)法與在數(shù)值分析中學(xué)到的積分方法相比，沒有誤差界以及難以求出函數(shù)的上界M也是它的缺點(diǎn)然而，蒙特卡羅方法可以推廣到多變量函數(shù)，在那里它變得更加實(shí)用Monte Carlo數(shù)值積分的優(yōu)點(diǎn)與一般的數(shù)值積分方法比較，Monte Carlo方法具有以下優(yōu)點(diǎn)： 蒙特卡羅方法是隨機(jī)的，為使預(yù)測(cè)值與真值之差變小需要做大量的試驗(yàn)在最終的估計(jì)中，要討論保證預(yù)先給定的置信水平所要求的試驗(yàn)次數(shù)，需要統(tǒng)計(jì)學(xué)的背景知識(shí)，然而作為一般準(zhǔn)則，結(jié)果的精度提高一倍(即誤

5、差減少一半)，試驗(yàn)次數(shù)大約需要增至4倍做任何的蒙特卡羅模擬，都要用到隨機(jī)數(shù)2.同樣拋100次的結(jié)果也不近相同。蒙特卡羅模擬是一種隨機(jī)模型1.拋100次硬幣得到51個(gè)正面，并且接下來(lái)的10次(即使不太可能剛巧10次)全為正面的情況是可能出現(xiàn)的，這樣，用110次的結(jié)果進(jìn)行估計(jì)實(shí)際上比用100次要差 要記住，對(duì)于根據(jù)模擬結(jié)果的預(yù)測(cè)寄予太多的信任是有危險(xiǎn)的，特別是在模擬中包含的假設(shè)沒有清楚表明的時(shí)候還有，由于用了大量的數(shù)據(jù)和龐大的計(jì)算，再加上非專業(yè)人員理解模擬模型和計(jì)算機(jī)輸出相對(duì)容易，所以常會(huì)導(dǎo)致對(duì)模擬結(jié)果的過(guò)分相信二.隨機(jī)數(shù)的生成1.蒙特卡羅模擬的關(guān)鍵是生成優(yōu)良的隨機(jī)數(shù)。2.在計(jì)算機(jī)實(shí)現(xiàn)中,我們是通

6、過(guò)確定性的算法生成 隨機(jī)數(shù),所以這樣生成的序列在本質(zhì)上不是隨機(jī) 的,只是很好的模仿了隨機(jī)數(shù)的性質(zhì)(如可以通過(guò) 統(tǒng)計(jì)檢驗(yàn))。我們通常稱之為偽隨機(jī)數(shù)(pseudo-random numbers)。3.在模擬中，我們需要產(chǎn)生各種概率分布的隨機(jī)數(shù)，而大多數(shù)概率分布的隨機(jī)數(shù)產(chǎn)生均基于均勻分布U(0,1)的隨機(jī)數(shù)。U(0,1)隨機(jī)數(shù)的生成乘同余法： 稱為種子,a 是乘因子,m是模數(shù)一個(gè)簡(jiǎn)單的例子一個(gè)簡(jiǎn)單的例子(續(xù))上面的例子中，第一個(gè)隨機(jī)數(shù)生成器的周期長(zhǎng)度是 10，而后兩個(gè)生成器的周期長(zhǎng)度只有它的一半。我們自然希望生成器的周期越長(zhǎng)越好，這樣我們得到的分布就更接近于真實(shí)的均勻分布。線性同余生成器(混合同余法

7、)(Linear Congruential Generator )c是非負(fù)整數(shù).通過(guò)適當(dāng)選取參數(shù)c可以改善隨機(jī)數(shù)的統(tǒng)計(jì)性質(zhì)(獨(dú)立性,均勻性).常用的線性同余生成器Modulus mMultiplier aReference231-1=214748364716807Lewis, Goodman, and Miller39373LEcuyer742938285Fishman and Moore950706376Fishman and Moore1226874159Fishman and Moore214748339940692LEcuyer214748356340014LEcuyer復(fù)雜一些的生成

8、器Multiple recursive generator算法實(shí)現(xiàn)許多程序語(yǔ)言中都自帶生成隨機(jī)數(shù)的方法，如 c 中的 random() 函數(shù)，Matlab中的rand()函數(shù)等。但這些生成器生成的隨機(jī)數(shù)效果很不一樣，比如 c 中的函數(shù)生成的隨機(jī)數(shù)性質(zhì)就比較差，如果用 c ，最好自己再編一個(gè)程序。Matlab 中的 rand() 函數(shù)，經(jīng)過(guò)了很多優(yōu)化。可以產(chǎn)生性質(zhì)很好的隨機(jī)數(shù)，可以直接利用。從U(0,1)到其它概率分布的隨機(jī)數(shù)1.離散型隨機(jī)數(shù)的模擬2.連續(xù)型隨機(jī)數(shù)的模擬3.正態(tài)隨機(jī)數(shù)的模擬1.離散型隨機(jī)數(shù)的模擬設(shè)隨機(jī)變量 的分布律為令將 作為區(qū)間(0,1)的分點(diǎn).若隨機(jī)變量 ,有令則有據(jù)此,可得

9、產(chǎn)生 的隨機(jī)數(shù)的具體過(guò)程為:每產(chǎn)生一個(gè)(0,1)區(qū)間上均勻分布隨機(jī)數(shù) ,若 則令 取值 .例1:離散型隨機(jī)變量X有如下分布律: X 0 1 2P(x) 0.3 0.3 0.4設(shè) 是(0,1)上均勻分布的隨機(jī)數(shù),令則 是具有X分布律的隨機(jī)數(shù).2 隨機(jī)行為的模擬 本節(jié)的目的說(shuō)明如何對(duì)簡(jiǎn)單的隨機(jī)行為建模，以形成直覺和理解考察3個(gè)簡(jiǎn)單的隨機(jī)模型： 1拋一枚正規(guī)的硬幣 2擲一個(gè)或一對(duì)正規(guī)的骰子 3擲一個(gè)或一對(duì)不正規(guī)的骰子 拋一枚硬幣時(shí)得到正面或反面的機(jī)會(huì)是1/2，于是拋很多次時(shí)出正面次數(shù)的比例接近0.5設(shè)x為0,1內(nèi)的隨機(jī)數(shù)，f(x)定義如下：當(dāng)0 x0.5, f(x)=正面，否則 f(x)=反面， f

10、(x)將結(jié)果是正面或反面賦值到0，1內(nèi)的一個(gè)數(shù)，隨機(jī)賦值時(shí)我們可利用這個(gè)函數(shù)的累積性質(zhì)拋很多次時(shí)能夠得到如下的出現(xiàn)百分比：2.1 一枚正規(guī)的硬幣2.3 擲一個(gè)不正規(guī)的骰子2.連續(xù)型隨機(jī)數(shù)的模擬a.逆變換方法(常用) (Inverse Transform Method)b.舍取方法 (Acceptance-Rejection Method)定理: 設(shè)隨機(jī)變量Y的分布函數(shù)F(y)是連續(xù)函數(shù), 而U是在(0,1)上均勻分布的隨機(jī)變量, 令, 則Y與X有相同的分布.舍取方法舍取方法(Acceptance-Rejection )方法最早由 Von Neumann提出，現(xiàn)在已經(jīng)廣泛應(yīng)用于各種隨機(jī)數(shù)的生成。

11、基本思路：通過(guò)一個(gè)容易生成的概率分布 g 和一個(gè)取舍準(zhǔn)則生成另一個(gè)與 g 相近的概率分布 f 。具體步驟：下面我們驗(yàn)證由上述步驟生成的隨機(jī)數(shù) Y 確實(shí)具有密度函數(shù) f(x) 所以為了提高舍取法的效率，我們應(yīng)該使 c 的取值盡可能的小，也就是使 f 和 g 的分布更為相近。3.標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)的生成正態(tài)分布是概率統(tǒng)計(jì)中最重要的分布，在此我們著重討論如何生成標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)。引理：Box-Muller 算法利用中心極限定理設(shè) 是n個(gè)相互獨(dú)立的在(0,1)上均勻分布的隨機(jī)變量,由中心極限定理知 漸近服從正態(tài)分布N(0,1).一般取n=10即可,若取n=12,則上式簡(jiǎn)化為 再由公式 即可得到正態(tài)分

12、布 的隨機(jī)數(shù).Matlab程序Function r=rnd-u(a,b)%產(chǎn)生在a,b間均勻分布的隨機(jī)數(shù)r=a+(b-a)*rand;returnMatlab程序Function r=rnd-beta(lamda)%模擬指數(shù)分布%lamda表示指數(shù)分布的參數(shù)r=-log(rand)/lamda;returnMatlab程序Function y= rnd(mean, segema)%模擬均值為mean,方差為segma的正態(tài)分布r=rand(1,12);x=sum(r)-6;y=segma*x+mean;return實(shí)例一 例（港口系統(tǒng)）考察一個(gè)帶有船只卸貨設(shè)備的小港口，任何時(shí)間僅能為一艘船只卸

13、貨船只進(jìn)港是為了卸貨，相鄰兩艘船到達(dá)的時(shí)間間隔在15分鐘到145分鐘之間變化一艘船只卸貨的時(shí)間由所卸貨物的類型決定，在45分鐘到90分鐘之間變化 需要回答以下問(wèn)題： 1每艘船只在港口的平均時(shí)間和最長(zhǎng)時(shí)間是多少? 2若一艘船只的等待時(shí)間是從到達(dá)到開始卸貨的時(shí)間，每艘船只的平均等待時(shí)間和最長(zhǎng)等待時(shí)間是多少? 3卸貨設(shè)備空閑時(shí)間的百分比是多少? 4船只排隊(duì)最長(zhǎng)的長(zhǎng)度是多少?3 實(shí)例二 為了得到一些合理的答案，利用計(jì)算機(jī)或可編程計(jì)算器來(lái)模擬港口的活動(dòng)假定相鄰兩艘船到達(dá)的時(shí)間間隔和每艘船只卸貨的時(shí)問(wèn)在它們各自的時(shí)間區(qū)問(wèn)內(nèi)均勻分布，例如兩艘船到達(dá)的時(shí)間間隔可以是15到145之間的任何整數(shù)，且這個(gè)區(qū)間內(nèi)的任

14、何整數(shù)等可能地出現(xiàn) 在給出模擬這個(gè)港口系統(tǒng)的一般算法之前，考慮有5艘船只的假想情況對(duì)每艘船只有以下數(shù)據(jù)： 設(shè)想碼頭設(shè)備的擁有者關(guān)心他們提供服務(wù)的質(zhì)量，并且要評(píng)價(jià)各種管理模式以確定為了改善服務(wù)是否值得增加費(fèi)用做一些統(tǒng)計(jì)可以幫助對(duì)服務(wù)質(zhì)量的評(píng)價(jià)。 例如，呆在港口時(shí)間最長(zhǎng)的船只是船5，呆了130分鐘，而平均是89分鐘(表5-14)通常顧客對(duì)等待時(shí)間的長(zhǎng)短非常在乎，例中最長(zhǎng)的等待時(shí)間為55分鐘，而平均是26分鐘如果排隊(duì)太長(zhǎng)有些顧客會(huì)改到別處去做生意，例中最長(zhǎng)的隊(duì)是2 用下面的蒙特卡羅模擬算法可以做這些統(tǒng)計(jì)，對(duì)各種管理模式進(jìn)行估價(jià) 現(xiàn)在假定你是碼頭設(shè)備擁有者的顧問(wèn)，如果能夠雇用更多的勞動(dòng)力，或者得到更好

15、的卸貨設(shè)備，使卸貨時(shí)間減少到每艘船3575分鐘，會(huì)有什么影響?表5-16給出了基于模擬算法的結(jié)果。從表5-16可以看到，每艘船的卸貨時(shí)間減少1520分鐘，使得船只呆在港口的時(shí)間特別是等待時(shí)間縮短了，而設(shè)備空閑時(shí)問(wèn)的百分比卻增加了近一倍船主對(duì)此是滿意的，因?yàn)檫@提高了長(zhǎng)期行駛時(shí)每艘船運(yùn)送貨物的效率 這樣，人港貿(mào)易好像會(huì)增加如果貿(mào)易量增加使得相鄰兩艘船到達(dá)的時(shí)間間隔縮短到10到120分鐘之間，模擬結(jié)果如表5-17從這個(gè)表可以看到，隨著貿(mào)易量的增加，船只又要在港口呆更長(zhǎng)的時(shí)間，但設(shè)備空閑時(shí)間少多了，于是船主和設(shè)備擁有者都隨著貿(mào)易量的增加而受益 假定我們對(duì)兩艘船到達(dá)的時(shí)間間隔和每艘船的卸貨時(shí)間分別在15

16、betweeni145和45unloadi90內(nèi)均勻分布不滿意，決定收集港口系統(tǒng)的經(jīng)驗(yàn)數(shù)據(jù)，并將結(jié)果并人我們的模型設(shè)想觀測(cè)了利用港口卸貨的1200艘船，收集的數(shù)據(jù)見表5-18 最后，按照表5-19和表5-20給出的規(guī)則生成betweeni和unloadi。(i=1,2,3,n)，將線性樣條子模型并入港口系統(tǒng)的模擬模型用這些子模型，表5-21給出了每次100艘船共6次獨(dú)立模擬的結(jié)果.作業(yè)一例 模擬求近似圓周率在邊長(zhǎng)為1的正方形內(nèi)有一半徑為0.5的內(nèi)切圓.現(xiàn)在模擬產(chǎn)生在正方形內(nèi)均勻分布的點(diǎn)n個(gè).如果有m個(gè)在圓內(nèi),則圓面積與正方形的面積比可近似為m/n.即/4m/n 4m/nn=10000m=0;F

17、or i=1:n if rand2+rand2 =0,則可以通過(guò)模擬估算.構(gòu)造一個(gè)矩形包含曲邊梯形,dmax f(x).產(chǎn)生n(足夠大)個(gè)在矩形區(qū)域內(nèi)的點(diǎn),如果落在由函數(shù)f(x)構(gòu)成的曲邊梯形內(nèi)的點(diǎn)為m個(gè),則所求定積分為 .n=106;a=0;b=1;d=1;m=0;for i=1:n x=a+rand*(b-a); y=d*rand; if y=x2 m=m+1; endends=m/n*(b-a)*dn=106;x=rand(1,n);y=x.2;s=sum(y)/n采用前面講的方法:作業(yè) 三 趕上火車的概率問(wèn)題描述:如圖,一列火車從A站開往B站,某人每天趕往B站上這趟火車.他已了解到:(

18、1)火車從A站到B站的運(yùn)行時(shí)間是均值為30分鐘,標(biāo)準(zhǔn)差為2分鐘的隨機(jī)變量;(2)火車在下午大約1點(diǎn)離開A站,離開時(shí)刻的頻率分布如下:出發(fā)時(shí)刻午后1:00午后1:05午后1:10頻率 0.7 0.2 0.1此人到達(dá)B站的時(shí)刻頻率分布為:時(shí)刻午后1:28午后1:30午后1:32午后1:34頻率0.30.40.20.1問(wèn)他能趕上火車的概率是多少?變量說(shuō)明 :火車從A站出發(fā)的時(shí)刻; :火車從A站到B站的運(yùn)行時(shí)間,單位:分鐘; :他到達(dá)B站的時(shí)刻問(wèn)題分析與假設(shè):此題包含多個(gè)隨機(jī)因數(shù).這里假設(shè) , , 都是隨機(jī)變量,其中 服從正態(tài)分布.模型建立很顯然,他能及時(shí)趕上火車的條件是: + .為了簡(jiǎn)化計(jì)算,將下午1點(diǎn)記為0時(shí)刻. 和 的分布律如下: /min 0 5 10 P(t) 0.7 0.2 0.1 /min 28 30 32 34 P(t) 0.3 0.4 0.2 0.1如果r為在(0,1)均勻分布的隨機(jī)數(shù),為了模擬隨機(jī)