蒙特卡羅方法課件(PPT 56頁)

1、1目 錄 第一章 蒙特卡羅方法概述第二章 隨機(jī)數(shù)的產(chǎn)生第三章 EM算法和MCMC方法參考書 : 茆詩松等, 高等數(shù)理統(tǒng)計(jì)(第6章), 高等教育出版社,1998;2.徐鐘濟(jì)，蒙特卡羅方法，上海科學(xué)技術(shù)出版社第1頁，共56頁。2第一章 蒙特卡羅方法概述 蒙特卡羅方法又稱隨機(jī)抽樣技巧或統(tǒng)計(jì)試驗(yàn)方法。 蒙特卡羅方法是一種計(jì)算方法，但與一般數(shù)值計(jì)算方法有很大區(qū)別。它以概率統(tǒng)計(jì)理論為基礎(chǔ)。由于蒙特卡羅方法能夠比較逼真地描述事物的特點(diǎn)及物理實(shí)驗(yàn)過程，解決一些數(shù)值方法難以解決的問題，因而該方法的應(yīng)用領(lǐng)域日趨廣泛。第2頁，共56頁。31.蒙特卡羅方法的基本思想 理論基礎(chǔ)：大數(shù)定律；中心極限定理； F(X)U(0

2、,1)?；舅枷耄?.當(dāng)所求問題的解是某個(gè)事件的概率，或者是某個(gè)隨機(jī)變量的期望，或與概率、數(shù)學(xué)期望有關(guān)的量時(shí)，通過某種試驗(yàn)的方法，得出該事件發(fā)生的頻率，或該隨機(jī)變量若干個(gè)觀察值的算術(shù)平均值，根據(jù)大數(shù)定律得到問題的解；2. 要生成分布函數(shù)為F(x)的隨機(jī)數(shù)，可先生成U(0,1)隨機(jī)數(shù)F，則可得到隨機(jī)數(shù)X=F-1(F) 。第3頁，共56頁。4例（利用MC進(jìn)行歐式期權(quán)定價(jià)）設(shè)股票價(jià)格St服從風(fēng)險(xiǎn)中性測(cè)度下的幾何Brown運(yùn)動(dòng)：其離散化形式為根據(jù)金融工程理論，設(shè)現(xiàn)在股票價(jià)格為S0，T時(shí)刻到期（單位天），敲定價(jià)為K的歐式看漲期權(quán)的價(jià)格為MC方案：按照（1）遞推產(chǎn)生n條風(fēng)險(xiǎn)中性測(cè)度下的軌道，提取出ST (

3、n)；（2）第4頁，共56頁。52. 蒙特卡羅方法的誤差 根據(jù)中心極限定理如果隨機(jī)變量序列X1，X2，XN獨(dú)立同分布，且具有有限非零的方差2 ，即則當(dāng)N充分大時(shí)，有如下的近似式它表明，誤差收斂速度的階為 以概率1-成立。第5頁，共56頁。6通常，蒙特卡羅方法的誤差定義為關(guān)于蒙特卡羅方法的誤差需說明兩點(diǎn)：第一，蒙特卡羅方法的誤差為概率誤差，這與其他數(shù)值計(jì)算方法是有區(qū)別的。第二，誤差中的均方差是未知的，必須使用其估計(jì)值來代替，在計(jì)算所求量的同時(shí)，可計(jì)算出 。 第6頁，共56頁。7減小方差的各種技巧 顯然，當(dāng)給定置信度后，誤差由和N決定。要減小，或者是增大N，或者是減小方差2。在固定的情況下，要把精

4、度提高一個(gè)數(shù)量級(jí)，試驗(yàn)次數(shù)N需增加兩個(gè)數(shù)量級(jí)。因此，單純?cè)龃驨不是一個(gè)有效的辦法。降低方差的各種技巧，引起了人們的普遍注意。一般來說，降低方差的技巧，往往會(huì)使觀察一個(gè)子樣的時(shí)間增加。在固定時(shí)間內(nèi)，使觀察的樣本數(shù)減少。所以，一種方法的優(yōu)劣，需要由方差和觀察一個(gè)子樣的費(fèi)用（使用計(jì)算機(jī)的時(shí)間）兩者來衡量。這就是蒙特卡羅方法中效率的概念。它定義為 其中c是觀察一個(gè)子樣的平均費(fèi)用。第7頁，共56頁。8蒙特卡羅方法的特點(diǎn)優(yōu)點(diǎn)能夠比較逼真地描述具有隨機(jī)性質(zhì)的事物的特點(diǎn)及物理實(shí)驗(yàn)過程。受幾何條件限制小。收斂速度與問題的維數(shù)無關(guān)。誤差容易確定。程序結(jié)構(gòu)簡單，易于實(shí)現(xiàn)。 缺點(diǎn)收斂速度慢。誤差具有概率性。第8頁，共

5、56頁。9第二章 隨機(jī)數(shù)的產(chǎn)生2.1 逆變換法設(shè)隨機(jī)變量X的分布函數(shù)為F(x),定義定理2.1 設(shè)隨機(jī)變量U服從U(0,1)分布，則的分布函數(shù)為F(x).由定理2.1,要生成分布函數(shù)為F(x)的隨機(jī)數(shù)，可先生成U(0,1)隨機(jī)數(shù)U，則可得到隨機(jī)數(shù)X=F-1(U) 第9頁，共56頁。102.2 合成法如果X的密度函數(shù)p(x)難于抽樣,而X關(guān)于Y的條件密度函數(shù)p(x|y)以及Y的密度函數(shù)g(y)均易于抽樣，則X 的隨機(jī)數(shù)可如下產(chǎn)生：Step1 由Y的分布g(y)抽取y;Step2 由X關(guān)于Y的條件密度函數(shù)p(x|y)抽取x.例2.1 設(shè)X的密度函數(shù)為由合成法，X的隨機(jī)數(shù)可如下抽?。?）取uU(0,

6、1); 2)取 ,確定i，使3) 由pi(x)抽取x.第10頁，共56頁。112.3 篩選抽樣 當(dāng)p(x)難以直接抽樣時(shí)，如果可以將p(x) 表示成p(x)=ch(x)g(x),其中h(.)是一密度函數(shù)且易于抽樣，而0g(y),回到1）上述方法就是篩選抽樣法，它是一種非常重要的抽樣方法，可解決許多難以直接抽樣的分布的抽樣問題。第11頁，共56頁。12h(x)的的選取有多種方法。一種直觀的方法是：如果存在一個(gè)函數(shù)M(x)，滿足p(x)M(x),且令h(x)=M(x)/c, 若h(x)易于抽樣，則篩選抽樣變?yōu)?）由U(0,1)抽取u，由h(y)抽取y；2）如果up(y)/M(y),則x=y停止；3

7、）如果u p(y)/M(y),回到1）。篩選抽樣的理論依據(jù)如下：定理 設(shè)X的密度函數(shù)為p(x),且p(x)=ch(x)g(x)，其中01時(shí)，如果 ，則x=y， 否則轉(zhuǎn)到1）；第14頁，共56頁。152.4 隨機(jī)向量的抽樣法設(shè)X1,Xk的聯(lián)合概率密度為定理2.4 設(shè)U1,Uk是獨(dú)立同分布的U(0,1)變量， X1,Xk是方程的解，其中 是對(duì)應(yīng)于 的分布函數(shù)，則X1,Xk的分布為(2.4).(2.4)(2.5)第15頁，共56頁。16隨機(jī)向量的逆變換抽樣法：由U(0,1)分布獨(dú)立地抽取u1,uk;用方程（2.5)解x1,xk例2.3 設(shè)X1,X2的聯(lián)合密度函數(shù)為試生成X1,X2的隨機(jī)數(shù)。解：第16

8、頁，共56頁。17相應(yīng)的邊際分布函數(shù)和條件分布函數(shù)分別為方程(2.5)變?yōu)榇朔匠滩灰捉?，不妨交換兩自變量的次序第17頁，共56頁。18相應(yīng)的邊際分布函數(shù)和條件分布函數(shù)分別為方程(2.5)變?yōu)閷?duì)服從特定分布的隨機(jī)向量有一些特殊的抽樣方法。第18頁，共56頁。19例2.6 試生成k維正態(tài)分布 的隨機(jī)數(shù)。解：注意到若 ，則存在下三角陣使其中C可由迭代實(shí)現(xiàn)：首先，由 ，有從而。因于是得依此類推，第19頁，共56頁。20一般迭代公式為至此，我們可以給出k維正態(tài)分布的抽樣步驟：1)迭代計(jì)算 ；2)由N(0,1)分布獨(dú)立抽取k個(gè)隨機(jī)數(shù) ；3)計(jì)算第20頁，共56頁。212.5 隨機(jī)模擬計(jì)算2.5.1 隨機(jī)投

9、點(diǎn)法考慮積分 ，設(shè)a,b有限，0f(x)M,令=(x,y):axb,0yM,并設(shè)(X,Y)是在上均勻分布的二維隨機(jī)向量，其聯(lián)合密度函數(shù)為則易見， 是中曲線f(x)下方面積。 假設(shè)我們向中投點(diǎn)，若點(diǎn)落在y=f(x)下方稱為中的，則點(diǎn)中的概率為第21頁，共56頁。22若我們進(jìn)行了n次投點(diǎn)，其中n0次中的，則可以得到一個(gè)估計(jì)不難看出， 是的無偏估計(jì)，且其方差為(2.5.1)第22頁，共56頁。232.5.2 樣本均值法于是，積分注意到，若XU(a,b),則由大數(shù)定律，若 ，則MC方法為：1) 獨(dú)立產(chǎn)生n個(gè)U(a,b)隨機(jī)數(shù)2)按(2.5.2)估計(jì)。(2.5.2)第23頁，共56頁。24可證，在0f(

10、x)M條件下，2.5.3 降低方差的技術(shù)Monte Carlo 方法中一類重要的研究課題是考慮一些降低估計(jì)方差的技術(shù)。常用的方法有：重要抽樣法，分層抽樣法，關(guān)聯(lián)抽樣法等。一 重要抽樣法由上節(jié)，樣本平均法比投點(diǎn)法有效，將樣本平均法做更一般的推廣，設(shè)g(x)是(a,b)上的密度函數(shù)，改寫第24頁，共56頁。25由大數(shù)定律，若 ，則MC方法為：1)選擇適當(dāng)?shù)膅(x),獨(dú)立產(chǎn)生n個(gè)g(x)隨機(jī)數(shù)2)由(2.5.3)估計(jì)。顯然(2.5.3)第25頁，共56頁。26從理論上看，因，若f(x)0,取則有因?yàn)槲粗?，這是作不到的，但它提示我們?nèi)(x)與f(x)形狀接近，應(yīng)能降低方差。這就是重要抽樣法的基本思想

11、。其方差與g(x)有關(guān)。問題變?yōu)椋绾芜x擇g(.)使估計(jì)的方差最小。第26頁，共56頁。27例2.5.1 分別用投點(diǎn)法，均值法，重要抽樣法，求積分 ,比較各種方法的有效性。解 i)投點(diǎn)法1)產(chǎn)生隨機(jī)數(shù) 2) 對(duì)每對(duì) ，記 的次數(shù)為n0.則Gii)均值法1)產(chǎn)生隨機(jī)數(shù) 2)第27頁，共56頁。28iii)重要抽樣法由重要抽樣法的思想，需選擇一個(gè)與 相似的密度函數(shù)。由Taylor展開式 取1)產(chǎn)生隨機(jī)數(shù) 2) 取則（數(shù)值計(jì)算）真值投點(diǎn)法均值法重要抽樣法1.718281.87561.81451.7219模擬結(jié)果第28頁，共56頁。29二、分層抽樣法另一種利用貢獻(xiàn)率大小來降低估計(jì)方差的方法是分層抽樣法

12、。它首先把樣本空間D分成一些不交的小區(qū)間 ，然后在各小區(qū)間內(nèi)的抽樣數(shù)由其貢獻(xiàn)大小決定。即，定義 ，則Di內(nèi)的抽樣數(shù)ni應(yīng)與pi成正比?？紤]積分將0,1分成m個(gè)小區(qū)間：則記 為第i個(gè)小區(qū)間的長度，i=1,m.在每個(gè)小區(qū)間上的積分值可用均值法估計(jì)出來，然后將其相加即可給出的一個(gè)估計(jì)。具體步驟為：第29頁，共56頁。301) 獨(dú)立產(chǎn)生U(0,1)隨機(jī)數(shù)2)計(jì)算3) 計(jì)算于是可得的估計(jì)為(2.5.4)易見， 是的無偏估計(jì)，其方差為(2.5.5)(2.5.6)第30頁，共56頁。31續(xù)例2.5.1 考察分層抽樣法求積分 的方差。解：先將區(qū)間0,1劃分成兩個(gè)小區(qū)間0,0.5,0.5,1,則設(shè)一共抽n個(gè)隨機(jī)

13、數(shù)，其中在0,0.5)上抽n1個(gè)，則使用分層抽樣法求得 的方差為第31頁，共56頁。32對(duì)n1求導(dǎo)易知，在n固定下，當(dāng) 時(shí)的方差最小，為如果我們將區(qū)間進(jìn)行10等份，并確定出最優(yōu)的抽樣次數(shù)分配： ，則可得到分層抽樣法估計(jì)的方差為.一般地，若諸 已知，在n固定下，當(dāng)時(shí)，估計(jì)的方差最小，為第32頁，共56頁。33分層抽樣法在實(shí)施上有兩個(gè)主要問題，其一是怎樣劃分區(qū)間，簡單而常用的方法是將區(qū)間等分；另一個(gè)問題是在區(qū)間劃分好后如何確定抽樣次數(shù)的分配。由于在實(shí)際中 總是未知的，因而前面最優(yōu)分配的結(jié)論無法應(yīng)用。即使如此，分層抽樣法還是有其作用的?？梢宰C明，即使取簡單的分配也有事實(shí)上，取 ，代入(2.5.5)得

14、由Cauchy-Schwarz不等式，有據(jù)此，在(2.5.6)式兩端各乘以 并相加得 于是第33頁，共56頁。34三、關(guān)聯(lián)抽樣法考慮積分差若用 估計(jì)，則其方差為顯然，在 確定后， 正相關(guān)度越高，則 的方差越小。這便是關(guān)聯(lián)抽樣法的基本出發(fā)點(diǎn)。考慮用重要抽樣法來估計(jì)I1,I2,即改寫為產(chǎn)生n個(gè)U(0,1)隨機(jī)數(shù) ；令則第34頁，共56頁。35第三章 數(shù)據(jù)添加算法 在Bayes統(tǒng)計(jì)或極大似然估計(jì)的計(jì)算中，經(jīng)常會(huì)遇到這樣一類問題：設(shè)我們能觀測(cè)到的數(shù)據(jù)是Y，關(guān)于Y的后驗(yàn)分布p(|Y)很復(fù)雜，難以直接進(jìn)行各種統(tǒng)計(jì)計(jì)算.假如我們能假定一些沒有能觀察到的潛在數(shù)據(jù)Z為已知（譬如，Y為某變量的截尾觀測(cè)值，Z為該變

15、量的真值），則可能得到一個(gè)關(guān)于的簡單的添加后驗(yàn)分布p(|Y,Z),利用p(|Y,Z)的簡單性我們可以進(jìn)行各種計(jì)算，如極大化，抽樣等，然后回過頭來，又可以對(duì)Z的假定做檢查或改進(jìn)。如此進(jìn)行，我們就將一個(gè)復(fù)雜的極大化問題轉(zhuǎn)變?yōu)橐幌盗泻唵蔚臉O大化或抽樣。在統(tǒng)計(jì)上，這種處理問題的方法稱為“數(shù)據(jù)添加算法”。 常用的“數(shù)據(jù)添加算法”有EM算法和Markov Chain Monte Carlo方法。第35頁，共56頁。363.1 EM算法 先考慮一種簡單情形。設(shè)某元件的失效時(shí)間Y關(guān)于變量x有直線回歸關(guān)系，假設(shè)在一次試驗(yàn)中得到一批數(shù)據(jù)，如圖， “”表示該元件失效時(shí)間坐標(biāo)， ”“表示對(duì)應(yīng)元件的截尾時(shí)間（小于失效時(shí)

16、間）。如果直線斜率和截矩的估計(jì)值已知，則我們可以在真實(shí)數(shù)據(jù)不小于截尾數(shù)據(jù)的前提下將各個(gè)被截尾的失效時(shí)間估計(jì)出來，從而得到所謂的”完全數(shù)據(jù)“，由此完全數(shù)據(jù)，重新對(duì)直線的斜率及截矩進(jìn)行估計(jì)，再依據(jù)新的估計(jì)量，得到新的”完全數(shù)據(jù)“。如此循環(huán)往復(fù)，則將一個(gè)復(fù)雜的估計(jì)問題替換成一系列簡單的估計(jì)問題。將之一般化，就給出EM算法。第36頁，共56頁。37EM算法是一種迭代方法，主要用來求后驗(yàn)分布的眾數(shù)（即極大似然估計(jì)）。它的每一步迭代由兩步組成：E步（求期望）和M步（極大化）。一般地，以p(|Y)表示基于Y的的后驗(yàn)密度，稱為觀測(cè)后驗(yàn)分布； p(|Y,Z)表示添加數(shù)據(jù)Z后得到的的后驗(yàn)密度，稱為添加后驗(yàn)分布；

17、p(Z|,Y)表示在給定觀測(cè)數(shù)據(jù)Y和參數(shù)條件下Z的條件密度。我們的目的是計(jì)算p(|Y)的眾數(shù)。于是EM算法如下進(jìn)行。記 為第i+1次迭代開始時(shí)后驗(yàn)眾數(shù)的估計(jì)值，則第i+1次迭代的兩步為E步：將p(|Y,Z)或log p(|Y,Z)關(guān)于Z的條件分布求期望，從而把Z積掉，即第37頁，共56頁。38M步：將 極大化，即找到一個(gè)點(diǎn) ，使將上述E，M步循環(huán)進(jìn)行，直至 充分小為止。例3.1 設(shè)總體X的分布律為其中(0,1),現(xiàn)進(jìn)行了X1234pk197次試驗(yàn)，觀察到1,2,3,4的頻數(shù)為取的先驗(yàn)分布()為U(0,1)分布,則的觀察后驗(yàn)分布為第38頁，共56頁。39現(xiàn)假設(shè)X=1可以分解為兩部分，其發(fā)生概率分

18、別為1/2和/4,令和y1-Z分別表示試驗(yàn)結(jié)果中落入這兩部分的次數(shù)（是不能觀測(cè)到的潛在數(shù)據(jù)），則的添加后驗(yàn)分布為(3.1.1)(3.1.2)顯然，用(3.1.2)式求極值比(3.1.1)式簡單。迭代如下：第39頁，共56頁。40E步：在給定下，M步：將 關(guān)于極大化得可以證明，在關(guān)于logp(|Y)的很一般的條件下，由算法得到的估計(jì)序列 收斂到的穩(wěn)定點(diǎn)。（不能保證是極大值點(diǎn)）。較為可行的辦法是選幾個(gè)不同的初值迭代，然后在諸估計(jì)值中加以選擇，這可減輕初值選取對(duì)結(jié)果的影響）第40頁，共56頁。41估計(jì)的精度假設(shè)EM算法最后的結(jié)果是 ，則根據(jù)似然估計(jì)的漸近正態(tài)性，其漸近方差可用 Fisher觀測(cè)信息的

19、倒數(shù)近似。（證明見高等數(shù)理統(tǒng)計(jì)p126定理2.5.4)第41頁，共56頁。423.2 Markov Chain Monte Carlo方法對(duì)于較簡單的后驗(yàn)分布，可直接計(jì)算或靜態(tài)MC等近似計(jì)算方法。但在實(shí)際中，觀測(cè)后驗(yàn)分布往往是復(fù)雜的，高維的，非標(biāo)準(zhǔn)形式的分布，上述方法都難以實(shí)施。對(duì)于這類問題，一種簡單且行之有效的Bayes計(jì)算方法就是MCMC。EM算法得到的是后驗(yàn)分布的眾數(shù)，有時(shí)我們希望得到其它一些后驗(yàn)量如后驗(yàn)均值，方差，后驗(yàn)分布的分位數(shù)等。計(jì)算這些后驗(yàn)量都可歸結(jié)為關(guān)于后驗(yàn)分布積分的計(jì)算。具體地，設(shè) 為后驗(yàn)密度，我們要計(jì)算的后驗(yàn)量可寫成某函數(shù)f(x)關(guān)于的期望(3.2.1)第42頁，共56頁。

20、433.2.1 基本思路MCMC方法的基本思想是通過建立一個(gè)平穩(wěn)分布為(x)的Markov鏈來得到(x)的樣本，基于這些樣本可以作各種統(tǒng)計(jì)推斷。比如，若得到了平穩(wěn)分布為(x)的Markov鏈的樣本軌道 ，則(3.2.1)可估計(jì)為(3.2.2)注 由Markov鏈平穩(wěn)分布的概念可知，不論Markov鏈從什么初始狀態(tài)出發(fā)，經(jīng)過一段時(shí)間后，各個(gè)時(shí)間的邊際分布都是平穩(wěn)分布，因此可將經(jīng)過某個(gè)m時(shí)間之后的觀察值看作平穩(wěn)分布(x)的樣本。由遍歷性定理可知，MCMC的關(guān)鍵是如何構(gòu)造平穩(wěn)分布為的Markov鏈的轉(zhuǎn)移核p(x,y)第43頁，共56頁。44MCMC方法可概括為如下三步：(1) 在X上選一個(gè)“合適”的

21、Markov鏈，確定其轉(zhuǎn)移核p(x,y),使鏈的平穩(wěn)分布為。(2)由X中某一點(diǎn)X(0)出發(fā)，用(1)中的Markov鏈產(chǎn)生序列X1,Xn;(3) 對(duì)某個(gè)m和大的n，任一函數(shù)f(x)的期望估計(jì)如下MCMC有許多研究專題，如鏈的收斂性判斷（m大小的確定），鏈的長度（n的大小）的確定，估計(jì)誤差等等。以下主要討論轉(zhuǎn)移核的構(gòu)造。第44頁，共56頁。453.2.2 滿條件分布MCMC主要用于多變量，非標(biāo)準(zhǔn)形式，且各變量間相互不獨(dú)立時(shí)分布的模擬。令 ，我們總可以寫出其中 。如果(3.2.1)式中右端各個(gè)因子能夠直接模擬，則只需要進(jìn)行靜態(tài)模擬（抽樣過程中不改變抽樣分布）。實(shí)際中很難滿足上述條件，因此需進(jìn)行動(dòng)態(tài)

22、模擬（抽樣分布隨模擬的進(jìn)行而改變），如MCMC，此時(shí)滿條件分布扮演了一個(gè)重要角色。(3.2.1)在導(dǎo)出滿條件分布時(shí)，應(yīng)注意到這樣一個(gè)事實(shí)：記 ，(3.2.2)第45頁，共56頁。46等價(jià)地，若 ，且 ，則(3.2.3)一般地，用y表示觀測(cè)數(shù)據(jù)， ，其中 分別表示參數(shù)，超參數(shù)和缺損數(shù)據(jù)，則有其中， 表示完全數(shù)據(jù)的密度函數(shù)， 表示先驗(yàn)分布， 表示超參數(shù)的分布。有(3.2.2),各變量的滿條件分布如下：第46頁，共56頁。47例3.2.1 設(shè)(X1,X2)的聯(lián)合密度為且 ，則其滿條件分布為第47頁，共56頁。483.2.3 Gibbs 抽樣思想：設(shè) 的密度為 ，任意固定TN，在給定 條件下，如下定義

23、隨機(jī)變量 具有密度函數(shù) ，則對(duì)任一可測(cè)集B，因而X的密度也是 。上述過程定義了一個(gè)由X到X的轉(zhuǎn)移核，且其相應(yīng)的平穩(wěn)分布是。這樣構(gòu)造的MCMC稱為Gibbs抽樣。當(dāng)T只有一個(gè)元素時(shí)稱為單元素Gibbs抽樣。第48頁，共56頁。49單元素Gibbs抽樣具體步驟如下：在給定起始點(diǎn) 后，假定第t次迭代開始時(shí)的估計(jì)值為 ，則第t次迭代分為如下n步：（1）由滿條件分布 抽取 (i)由滿條件分布 抽取 ; (n)由滿條件分布 抽取記 ，則 是平穩(wěn)分布為的Markov鏈的實(shí)現(xiàn)值，其由x到x的轉(zhuǎn)移概率函數(shù)為第49頁，共56頁。503.2.3 Metroplis-Hastings方法在Gibbs抽樣中， 可能很難抽取，這時(shí)可采用更一般化的Me