Monte Carlo（蒙特卡洛方法）

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

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

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

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

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

6、過確定性的算法生成 隨機(jī)數(shù),所以這樣生成的序列在本質(zhì)上不是隨機(jī) 的,只是很好的模仿了隨機(jī)數(shù)的性質(zhì)(如可以通過 統(tǒng)計檢驗)。我們通常稱之為偽隨機(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ù)一個簡單的例子一個簡單的例子(續(xù))上面的例子中，第一個隨機(jī)數(shù)生成器的周期長度是 10，而后兩個生成器的周期長度只有它的一半。我們自然希望生成器的周期越長越好，這樣我們得到的分布就更接近于真實的均勻分布。線性同余生成器(混合同余法

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

8、器Multiple recursive generator算法實現(xiàn)許多程序語言中都自帶生成隨機(jī)數(shù)的方法，如 c 中的 random() 函數(shù)，Matlab中的rand()函數(shù)等。但這些生成器生成的隨機(jī)數(shù)效果很不一樣，比如 c 中的函數(shù)生成的隨機(jī)數(shù)性質(zhì)就比較差，如果用 c ，最好自己再編一個程序。Matlab 中的 rand() 函數(shù)，經(jīng)過了很多優(yōu)化?？梢援a(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)的分點.若隨機(jī)變量 ,有令則有據(jù)此,可得

9、產(chǎn)生 的隨機(jī)數(shù)的具體過程為:每產(chǎn)生一個(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é)的目的說明如何對簡單的隨機(jī)行為建模，以形成直覺和理解考察3個簡單的隨機(jī)模型： 1拋一枚正規(guī)的硬幣 2擲一個或一對正規(guī)的骰子 3擲一個或一對不正規(guī)的骰子 拋一枚硬幣時得到正面或反面的機(jī)會是1/2，于是拋很多次時出正面次數(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)的一個數(shù)，隨機(jī)賦值時我們可利用這個函數(shù)的累積性質(zhì)拋很多次時能夠得到如下的出現(xiàn)百分比：2.1 一枚正規(guī)的硬幣2.3 擲一個不正規(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、基本思路：通過一個容易生成的概率分布 g 和一個取舍準(zhǔn)則生成另一個與 g 相近的概率分布 f 。具體步驟：下面我們驗證由上述步驟生成的隨機(jī)數(shù) Y 確實具有密度函數(shù) f(x) 所以為了提高舍取法的效率，我們應(yīng)該使 c 的取值盡可能的小，也就是使 f 和 g 的分布更為相近。3.標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)的生成正態(tài)分布是概率統(tǒng)計中最重要的分布，在此我們著重討論如何生成標(biāo)準(zhǔn)正態(tài)分布隨機(jī)數(shù)。引理：Box-Muller 算法利用中心極限定理設(shè) 是n個相互獨立的在(0,1)上均勻分布的隨機(jī)變量,由中心極限定理知 漸近服從正態(tài)分布N(0,1).一般取n=10即可,若取n=12,則上式簡化為 再由公式 即可得到正態(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實例一 例（港口系統(tǒng)）考察一個帶有船只卸貨設(shè)備的小港口，任何時間僅能為一艘船只卸

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

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

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

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

17、or i=1:n if rand2+rand2 =0,則可以通過模擬估算.構(gòu)造一個矩形包含曲邊梯形,dmax f(x).產(chǎn)生n(足夠大)個在矩形區(qū)域內(nèi)的點,如果落在由函數(shù)f(x)構(gòu)成的曲邊梯形內(nèi)的點為m個,則所求定積分為 .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è) 三 趕上火車的概率問題描述:如圖,一列火車從A站開往B站,某人每天趕往B站上這趟火車.他已了解到:(

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