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

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

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

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

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

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

6、1); 2)取 ,確定i，使3) 由pi(x)抽取x.第10頁，共56頁。112.3 篩選抽樣 當p(x)難以直接抽樣時，如果可以將p(x) 表示成p(x)=ch(x)g(x),其中h(.)是一密度函數(shù)且易于抽樣，而0g(y),回到1）上述方法就是篩選抽樣法，它是一種非常重要的抽樣方法，可解決許多難以直接抽樣的分布的抽樣問題。第11頁，共56頁。12h(x)的的選取有多種方法。一種直觀的方法是：如果存在一個函數(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ù)如下：定理 設X的密度函數(shù)為p(x),且p(x)=ch(x)g(x)，其中01時，如果 ，則x=y， 否則轉(zhuǎn)到1）；第14頁，共56頁。152.4 隨機向量的抽樣法設X1,Xk的聯(lián)合概率密度為定理2.4 設U1,Uk是獨立同分布的U(0,1)變量， X1,Xk是方程的解，其中 是對應于 的分布函數(shù)，則X1,Xk的分布為(2.4).(2.4)(2.5)第15頁，共56頁。16隨機向量的逆變換抽樣法：由U(0,1)分布獨立地抽取u1,uk;用方程（2.5)解x1,xk例2.3 設X1,X2的聯(lián)合密度函數(shù)為試生成X1,X2的隨機數(shù)。解：第16

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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