[科普中國]-馬爾科夫鏈蒙特卡洛方法

從理論上說，貝葉斯推斷和分析是容易實(shí)施的，即對于任何先驗(yàn)分布，只需要計(jì)算所需后驗(yàn)分布的性質(zhì)，如后驗(yàn)分布的矩(如后驗(yàn)均值、后驗(yàn)方差)、后驗(yàn)概率密度函數(shù)等，而這些計(jì)算木質(zhì)上就是計(jì)算后驗(yàn)分布某一函數(shù)的高維積分。但在實(shí)踐中，鑒于未知參數(shù)的后驗(yàn)分布多為高維、復(fù)雜的非常見分布，對這些高維積分進(jìn)行計(jì)算十分困難，這一困難使得貝葉斯推斷方法在實(shí)踐中的應(yīng)用受到很大的限制，在很長一段時(shí)間，貝葉斯推斷主要用于處理簡單低維的問題，以避免計(jì)算上的困難。MCMC方法突破了這一原木極為困難的計(jì)算問題，它通過模擬的方式對高維積分進(jìn)行計(jì)算，進(jìn)而使原本異常復(fù)雜的高維積分計(jì)算問題迎刃而解，使貝葉斯方法僅適用于解決簡單低維問題的狀況大有改觀為貝葉斯方法的應(yīng)用開辟了新的道路。1

MCMC方法的結(jié)構(gòu)如圖1所示。2

MCMC方法是使用馬爾科夫鏈的蒙特卡洛積分，其基木思想是：構(gòu)造一條Markov鏈，使其平穩(wěn)分布為待估參數(shù)的后驗(yàn)分布，通過這條馬爾科夫鏈產(chǎn)生后驗(yàn)分布的樣木，并基于馬爾科夫鏈達(dá)到平穩(wěn)分布時(shí)的樣木(有效樣本)進(jìn)行蒙特卡洛積分。

設(shè) 為某一空間，n為產(chǎn)生的總樣木數(shù)，m為鏈條達(dá)到平穩(wěn)時(shí)的樣本數(shù)，則MCMC方法的基木思路可概括為：

(1)構(gòu)造Markov鏈。構(gòu)造一條Markov鏈，使其收斂到平穩(wěn)分布 ；

(2)產(chǎn)生樣本：由 中的某一點(diǎn) 出發(fā)，用(1)中的Markov鏈進(jìn)行抽樣模擬，產(chǎn)生點(diǎn)序列： ；

(3)蒙特卡洛積分：任一函數(shù)f(x)的期望估計(jì)為： 。1

MCMC中采樣算法如圖2所示。

在采用MCMC方法時(shí)，馬爾科夫鏈轉(zhuǎn)移核的構(gòu)造至關(guān)重要，不同的轉(zhuǎn)移核構(gòu)造方法，將產(chǎn)生不同的MCMC方法，目前常用的MCMC方法主要有兩種：Gibbs抽樣和Metropolis-Hastings算法。1

Metropolis 算法是馬爾科夫鏈蒙特卡羅的基石。它是由 Metropolos 等人在 1953年的一篇僅 4 頁的文章中提出。Metropolis 算法用一個(gè)對稱的建議分布 T(x,y)來產(chǎn)生一個(gè)潛在的轉(zhuǎn)移點(diǎn)，然后根據(jù)特定的接受拒絕方法來決定是否轉(zhuǎn)移到該潛在點(diǎn)。最初的 Metropolis 算法將 T 取為對稱的函數(shù)，而 Metropolis-Hasting 方法將之推廣到非對稱的 T，其中 T 滿足 T(x,y)>0 的充要條件是 T(y,x)>0。已知當(dāng)前狀態(tài) ，該算法的迭代步驟如下：

第 1 步：從建議分布 抽取潛在轉(zhuǎn)移點(diǎn) y。

第 2 步：在[0,1]均勻分布中抽取一個(gè)隨機(jī)數(shù) ，并且更新 ，若滿足 ，否則令 。對于 r(x,y)，Metropolis 和 Hastings 建議采用如下形式的：

Gibbs 算法是一種特殊的 Metropolis 算法。對于一個(gè) d 維隨機(jī)變量 ，Gibbs 算法在生成 的第 i 個(gè)分量時(shí)選擇建議分布 T 為第 i 個(gè)分量基于其他所有分量的條件分布。即已知當(dāng)前狀態(tài) ， Gibbs 算法通過如下的方法更新 d 個(gè)分量來更新 ：

第 1 步 ： 對 于 固 定 的 i（初值為1），從條件分件=(|)中抽樣出。

第 2 步：令 i=i+1，重復(fù)第 1 步。3

MCMC方法依賴于模擬的收斂性，下面介紹三種常用的收斂性診斷方法。

這種方法的思路是選取多個(gè)不同的初值，同時(shí)產(chǎn)生多條馬爾科夫鏈，經(jīng)過一段時(shí)間后，若這幾條鏈穩(wěn)定下來，則說明算法收斂了。在實(shí)際操作中，可在同一個(gè)二維圖中畫出這些不同的馬爾科夫鏈產(chǎn)生的后驗(yàn)樣本值對迭代次數(shù)的散點(diǎn)圖，如果經(jīng)過若干次迭代后，這些散點(diǎn)圖基木穩(wěn)定，重合在一起，則可判斷其算法收斂。

這種方法的思路是選取多個(gè)不同的初值，同時(shí)產(chǎn)生T條馬爾科夫鏈，記第j條鏈的方差估計(jì)為，鏈內(nèi)方差的均值為W，鏈間方差為B，則：。R的值趨近1，則表明MCMC模擬收斂，且比較穩(wěn)定，通常R