1、1Monte Carlo Simulation Methods(蒙特卡羅模擬方法)主要內(nèi)容: 1.各種隨機(jī)數(shù)的生成方法. 2.MCMC方法.2從Buffon 投針問(wèn)題談起 220, /2 0, sin,:sin.llaXAX隨機(jī)投針可以理解成針的中心點(diǎn)與最近的平行線的距離X是均勻地分布在區(qū)間 上的r.v.，針與平行線的夾角 是均勻地分布在區(qū)間 上的r.v.，且X與 相互獨(dú)立，于是針與平行線相交的充要條件為 即相交 3Buffon 投針問(wèn)題2sin0022(sin)2lllpP Xdxdaa 于是有：2lap4試驗(yàn)者時(shí)間(年)針長(zhǎng)投針次數(shù)相交次數(shù)的估計(jì)值Wolf18500.80500025323

2、.15956Smith18550.60320412183.15665Fox18840.7510304893.15951Lazzarini19250.83340818083.141592925數(shù)值積分問(wèn)題1( )()( ) 0,1 . . .1() nniinf x dxEf Xf xXUknii df Unwith probability as n 10kk 計(jì)算積分 我們可以將此積分看成 的數(shù)學(xué)期望。其中 (均勻分布)。于是可以將上式積分看作是f(X)的數(shù)學(xué)期望.若U ,1為U U0,1.則可以取作為 的估計(jì)，由大數(shù)定律，可以保證收斂性，即：這表明可以用隨機(jī)模擬的方法計(jì)算積分。6Monte

3、Carlo數(shù)值積分的優(yōu)點(diǎn)與一般的數(shù)值積分方法比較，Monte Carlo方法具有以下優(yōu)點(diǎn)：1. Monte Carlo一般的數(shù)值方法很難推廣到高維積分的情形，而方法很容易推廣到高維情形2/1/22. ()() dO nO n一般的數(shù)值積分方法的收斂階為 ，而由中心極限定理可以保證 Monte Carlo 方法的收斂階為 。此收斂階與維數(shù)無(wú)關(guān)，且在高維時(shí)明顯優(yōu)于一般的數(shù)值方法。7隨機(jī)模擬計(jì)算的基本思路1.針對(duì)實(shí)際問(wèn)題建立一個(gè)簡(jiǎn)單且便于實(shí)現(xiàn)的概率統(tǒng)計(jì)模型，使所求的量（或解）恰好是該模型某個(gè)指標(biāo)的概率分布或者數(shù)字特征。2.對(duì)模型中的隨機(jī)變量建立抽樣方法，在計(jì)算機(jī)上進(jìn)行模擬測(cè)試，抽取足夠多的隨機(jī)數(shù)，對(duì)

4、有關(guān)事件進(jìn)行統(tǒng)計(jì)3.對(duì)模擬試驗(yàn)結(jié)果加以分析，給出所求解的估計(jì)及其精度(方差)的估計(jì)4.必要時(shí)，還應(yīng)改進(jìn)模型以降低估計(jì)方差和減少試驗(yàn)費(fèi)用，提高模擬計(jì)算的效率8隨機(jī)數(shù)的生成1.蒙特卡羅模擬的關(guān)鍵是生成優(yōu)良的隨機(jī)數(shù)。2.在計(jì)算機(jī)實(shí)現(xiàn)中,我們是通過(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ù)。9U(0,1)隨機(jī)數(shù)的生成一個(gè)簡(jiǎn)單的隨機(jī)數(shù)生成器：1101

5、 mod , , /iiiiixaxmuaxxmxm其中 均為整數(shù), 可以任意選取。111 ( ) , () iiiixf xug x隨機(jī)數(shù)生成器的一般形式：10一個(gè)簡(jiǎn)單的例子1i+116 mod11, u/11 6,11iiixxxam（）0 1 ,x 1,6,3,7,9,10,5,8,4當(dāng)時(shí) 得到序列:,1,6,3.,2.003,1 ,1,3,9.3,2,2,1,3,9,5,42,6,7,10,8,6.axax如果令 得到序列:如果令 得到序列：11一個(gè)簡(jiǎn)單的例子(續(xù))上面的例子中，第一個(gè)隨機(jī)數(shù)生成器的周期長(zhǎng)度是 10，而后兩個(gè)生成器的周期長(zhǎng)度只有它的一半。我們自然希望生成器的周期越長(zhǎng)越好

6、，這樣我們得到的分布就更接近于真實(shí)的均勻分布。0 (max在給定 的情況下，生成器的周期與和初值種子)選擇有關(guān)。12線性同余生成器(Linear Congruential Generator )111 () mod /iiiixaxcmuxm一般形式：0 cc 1.可以證明 的選擇對(duì)生成的隨機(jī)數(shù)的均勻性影響不大，所以為了提高計(jì)算速度，一般都令 。02. 1 mmaxm線性同余器可以達(dá)到的最長(zhǎng)周期為 ，我們可以通過(guò)適當(dāng)?shù)倪x擇 和，使無(wú)論選取怎樣的初值 都可以達(dá)到最大周期(一般選取 為質(zhì)數(shù))13常用的線性同余生成器Modulus mMultiplier aReference231-1=214748

7、364716807Lewis, Goodman, and Miller39373LEcuyer742938285Fishman and Moore950706376Fishman and Moore1226874159Fishman and Moore214748339940692LEcuyer214748356340014LEcuyer14復(fù)雜一些的生成器（一）1.Combining Generators:,1,1,111,11111111111 mod, / 1,2.( 1) mod(1)/ 0 (1)/ 0( )j ijj ijj ij ijJjij ijiiiijJxa xmuxmjJ

8、xxmuxmxmmxmm考慮個(gè)簡(jiǎn)單線性同余生成器：其中為所有中最大的一個(gè)15復(fù)雜一些的生成器（二）2.Multiple recursive generator1122102(.) mo,.)d/iiiki kikkixa xa xa xmuxxxxm需(要選取種子12-1-1 (,.) 1ijiii kkkxmxxxmm每個(gè)有種選擇，所以向量 可以取個(gè)不同的值，所以這樣的隨機(jī)數(shù)生成器的最大周期可以達(dá)到 ，大大提高了簡(jiǎn)單同余生成器的周期。16算法實(shí)現(xiàn)許多程序語(yǔ)言中都自帶生成隨機(jī)數(shù)的方法，如 c 中的 random() 函數(shù)，Matlab中的rand()函數(shù)等。但這些生成器生成的隨機(jī)數(shù)效果很不一樣

9、，比如 c 中的函數(shù)生成的隨機(jī)數(shù)性質(zhì)就比較差，如果用 c ，最好自己再編一個(gè)程序。Matlab 中的 rand() 函數(shù)，經(jīng)過(guò)了很多優(yōu)化。可以產(chǎn)生性質(zhì)很好的隨機(jī)數(shù)，可以直接利用。17由rand()函數(shù)生成的U0,1隨機(jī)數(shù)18由rand函數(shù)生成的2維隨機(jī)點(diǎn)19從U(0,1)到其它概率分布的隨機(jī)數(shù)U(0,1)的均勻分布的隨機(jī)數(shù)，是生成其他概率分布隨機(jī)數(shù)的基礎(chǔ)，下面我們主要介紹兩種將U(0,1)隨機(jī)數(shù)轉(zhuǎn)換為其他分布的隨機(jī)數(shù)的方法。1.逆變換方法 (Inverse Transform Method)2.舍取方法 (Acceptance-Rejection Method)20Inverse Transf

10、orm Method11 ( ) ( )inf :( ), 01( )( )XF xFyx F xyyFyF x設(shè)隨機(jī)變量 的分布函數(shù)為，定義：即是的反函數(shù)。11 (0,1) ( ) ( ) UXFUF x定理 ：設(shè)隨機(jī)變量服從 上的均勻分布，則的分布函數(shù)為 。11( ) ()( )( )( )FUP XxP FUxP UF xF x 由 的定義和均勻分布的分布函明數(shù)可證：得：21Inverse Transform Method1 1 ( ) (0,1) ( ) F xUuFu由定理，要產(chǎn)生來(lái)自的隨機(jī)數(shù)，只要先產(chǎn)生來(lái)自隨機(jī)數(shù)，然后計(jì)算 即可。具體步驟如下：(1) (0,1)U上生成均勻分布的隨

11、機(jī)數(shù) 。-1(2) , ( ) ( ) XXUFFx計(jì)算則為來(lái)自 分布的隨機(jī)數(shù).22幾個(gè)具體例子（一）-11 ( , ) 0 ( ) 1, ( )() 01 (0,1) ( - ) ( , ) XU a bxaxaF xaxbbaxbFyaba yyUUab a UU a b例 ： 設(shè)，則其分布函數(shù)為 ，生成隨機(jī)數(shù)，則是來(lái)自的隨機(jī)數(shù)。23幾個(gè)具體例子（二）/12 exp( ) ( )1, 0 ( )log(1) log(1) ( ) 1 xXXF xexFyyXUUUU 例 ：設(shè)服從指數(shù)分布，則的分布函數(shù)為：通過(guò)計(jì)算得，則：服從指數(shù)分布 其中服從均勻分布又因?yàn)?和有著同樣的分布，所以也可以取：

12、 log( )XU 24幾個(gè)具體例子（三）12013 () ( ) ,. () 0, ( ) niiiijjiiXF xc ccP XcpqqpF cq例 離散分布隨機(jī)數(shù)的生成設(shè)為取值的離散隨機(jī)變量且，令 ，則有，我們按照如下步驟生成隨機(jī)數(shù)：(0,1) ;U(1) 生成上均勻分布隨機(jī)數(shù)1 , 1;kkkqUqkn(2) 尋找滿足 ( )KXcXF xkk-1kk(3) 令，顯然P(X=c )=P(q 0.5,Y=-X,1;exp( 0.5(1) )U U UXUUXUUXUUX 由下述步驟生成Y：生成 上的均勻分布且相互獨(dú)立的樣本 ；令若 ，（生成指數(shù)分布的樣且則取返回步驟若，且則取返回步驟若

13、，則舍去X,返本）回步驟139隨機(jī)向量的抽樣方法（一）12(,.) nXXX若隨機(jī)向量 獨(dú)立，則可用生成一維隨機(jī)數(shù)的方法分別生成向量的每一維即可。121212121121121121(,.) (,.) ( ,.)() (|). (|,.)(),(|). (|,.)nnnnnnnXXXXXXp x xxp x p xxp xx xxF xF xxF xx xx若隨機(jī)向量 不獨(dú)立。設(shè)概率密度函數(shù)為：為對(duì)應(yīng)的分布函數(shù)40隨機(jī)向量的抽樣方法（二）1212111211212,. (0,1) ,. () (|,.) 2,3.,.( ,.)nniiinnU UUUXXXF xUF xx xxUinXXXp

14、x xx定理2:設(shè) 是獨(dú)立同分布的變量，是方程組：的解，則 的分布為 1121221. (0,1)(,.),. 2. , . nnnUu uuxXXXxx從分布中獨(dú)立抽取 解上面的方程組得到樣本 生成的基本步驟：41生成多維正態(tài)隨機(jī)數(shù)的方法（一）122121(,.) ( , ),11 ( )exp()()2(2 ) |()|nnnTijXXXXnNf xxx 若 服從 維正態(tài)分布 則密度函數(shù)為:其中 為已知均值向量，為已知協(xié)方差陣，為其行列式。12(,.) (0,), ( , )TnnnnZZ ZZNIAAXAZN 若 則可證明：42生成多維正態(tài)隨機(jī)數(shù)的方法（二）生成多維正態(tài)隨機(jī)數(shù)的具體步驟：

15、1. TAAA 找到一個(gè)矩陣 使其滿足 122. (0,1) ,.nNnz zz由 分布獨(dú)立抽取 個(gè)隨機(jī)數(shù)123. ( ,.)nxAzzz zz計(jì)算 ，其中43v用Monte Carlo方法求解Laplace方程v參見(jiàn)書(shū)上5.8節(jié) P213P21544馬氏鏈在馬氏鏈在Monte Carlo隨機(jī)模擬中的應(yīng)用隨機(jī)模擬中的應(yīng)用定義定義 為要模擬服從給定分布的隨機(jī)變量,用生成一個(gè)易于實(shí)現(xiàn)的不可約遍歷鏈 作為隨機(jī)樣本,使其平穩(wěn)分布為 的方法，稱為馬氏鏈蒙特卡羅方法馬氏鏈蒙特卡羅方法.0,nXXn蒙特卡羅方法的一個(gè)首要步驟是產(chǎn)生服從給定的概率分布函數(shù) 的隨機(jī)變量(或稱為隨機(jī)樣本),由概率論知識(shí),熟知下面的

16、結(jié)論.45引理引理 生成隨機(jī)變量U,使其分布滿足U0,1,記為UU0,1, F(x)是給定的一個(gè)分布函數(shù),記 為F(x)的反函數(shù),則X=F-1(U)分布函數(shù)為F(x).)1 , 0(,)(:sup)(1yyxFxyF46).().1 ();1 (,)()()2(;1 , 0),() 1 (:,),()(, 0)(,)(:xZYZYMgYUUUygYZxxMgxMxgx則返回步驟否則直接返回步驟取若抽取量由下列步驟生成隨機(jī)變滿足和常數(shù)的分布函數(shù)選取一個(gè)易于生成是給定的分布函數(shù)假設(shè)命題命題47).),()()()(/ )(,)(,)(,)(, )()()(:11nZEhyZEhxdxhnZhynZ

17、EhZZnxZxZEhxdxhPnnkknn大數(shù)定律知由的近似值作為積分取則對(duì)充分大的存在若樣本個(gè)獨(dú)立的只需生成服從的概率分布函數(shù)為其中計(jì)算積分應(yīng)用舉例48米特羅波利斯(Metropolis)等人在1953年最早給出了通過(guò)生成一馬氏鏈實(shí)現(xiàn)從分布 中采樣(生成相關(guān)的樣本)這一重要基本思想.隨后,哈斯汀(Hastings)將其推廣到更一般的形式.下面僅敘述狀態(tài)空間S為至多可數(shù)的情形:49.,0,).1 (,);1 (,),(,1 , 0)2(;),()(),()(, 1min),(,) ,().0,() 1 (:0,).(),),(,),(11平穩(wěn)分布為是不可約遍歷馬氏鏈則返回步驟取否則舍去返回步

18、驟則令若抽取并計(jì)算抽取由給定由下列步驟生成為參照矩陣稱陣不可約遍歷概率轉(zhuǎn)移矩的為任選的易于實(shí)現(xiàn)條件而為將要生成的概率分布記nXXXXYYXYXUUUYXTXXYTYYXYXTnSXXnXXSjijiTSiinnnnnnnnnnnnnTT50.)(),()(),()(min),()(),()(, 1min),()(),(),()()().,()()(,),)(,(),(,:jiijjiijijPjijTjjiTijiTiijTjjiTijijiTiPiSjiPjPiSjijijiTPX事實(shí)上即性有對(duì)稱轉(zhuǎn)移概率矩陣其一步是不可約遍歷馬氏鏈易證證明梗概.0,),(.)()(,的平穩(wěn)分布是由上式即證得

19、求和上式兩邊對(duì)nXXSiiPjiSjnSjji51Markov chain Monte Carlo (MCMC)問(wèn)題提出：( ) ( ) ( ) ( ) ( ) nnRxRf xx dxxx設(shè) 為 上的密度函數(shù)，在很多問(wèn)題中我們需要計(jì)算高維積分 ，其中常常是非常復(fù)雜的多維概率分布，我們無(wú)法用前面介紹的方法生成 分布的隨機(jī)數(shù)。因此，我們必須探討一些新的生成方法。52MCMC方法的基本思路MCMC 是一種簡(jiǎn)單有效的計(jì)算方法，在統(tǒng)計(jì)物理，Bayes 統(tǒng)計(jì)計(jì)算，顯著性檢驗(yàn)，極大似然估計(jì)等領(lǐng)域都有著廣泛的應(yīng)用。基本思路：( ) Markov ( ) xx 通過(guò)建立一個(gè)平穩(wěn)分布為的鏈來(lái)得到的樣本，基于這些

20、樣就可以做各種統(tǒng)計(jì)推斷。53概率轉(zhuǎn)移核的構(gòu)造（一）MCMC的方法有很多，在此我們只介紹Metropolis-Hastings方法。基本思路： ( , )() (0()1)( ,) ( ,)( ,) ( ,) ( , )1( ,) ( ,) ( ,) x xqx xp x xq x xx xxxp x xq x xx x dxxxp x x 任意選擇一個(gè)不可約的轉(zhuǎn)移概率以及一個(gè)函數(shù)，。對(duì)任一組合，定義：則易見(jiàn)構(gòu)成一個(gè)概率轉(zhuǎn)移核。54概率轉(zhuǎn)移核的構(gòu)造（二） ( ) ( | )( ,)1-( ,)ttxXxqxxxx xxx xxx 此方法的實(shí)施比較直觀：如果鏈在時(shí)刻處于狀態(tài)，即，則首先由產(chǎn)生一個(gè)潛

21、在的轉(zhuǎn)移，然后根據(jù)概率接受作為鏈下一時(shí)刻的狀態(tài)值，而以概率拒絕轉(zhuǎn)移到，從而鏈在下一時(shí)刻仍然處于狀態(tài)。 ( )( , )( , )xq 我們的目標(biāo)是使成為馬氏鏈的平穩(wěn)分布，下面我們就介紹在給定后，如何選擇55概率轉(zhuǎn)移核的構(gòu)造（三）Metropolis-Hastings 算法：( ) ( , )( ,)min1, ( ) ( ,)x q x xx xx q x x 一個(gè)常用的選擇：( ,) ( ) ( , )( ) ( ,)( ,)( )( , ) ( ) ( , )( ) ( ,)( )q x xx q x xx q x xp x xxq x xx q x xx q x xx 此時(shí)有：56 Markov ( ) ( ,)( ) ( , ) (*)( )Markov x p x xx p x xx 定理：由上述過(guò)程產(chǎn)生的鏈?zhǔn)强赡娴模矗呵沂擎湹钠椒€(wěn)分布。:,*( ) ( , ) ( ) ( ,)( ) ( ,)min 1,( ) ( ,) min ( ) ( ,), ( ) ( , )( ) ( ( ) ( , )minproofxxxxx q x xx p x xx q x