MCMC采樣技術及其在貝葉斯推斷上的應用

房愛東，謝士春

(宿州學院信息工程學院，安徽 宿州 234000)

從理論上講，貝葉斯統(tǒng)計是非常簡單的——后驗分布正比例于樣本似然分布與先驗分布.比例因子往往很難得到，所以不能直接用于精確推斷.大多數(shù)情況下，比例因子要求求解一個數(shù)值積分，然而當有多個參數(shù)時可能是困難的.盡管統(tǒng)計學家們從他們的決策理論研究中發(fā)現(xiàn)貝葉斯統(tǒng)計在理論上提供了真正優(yōu)勢，但在實踐中這些優(yōu)勢并不是真的可用，計算邊緣后驗分布密度函數(shù)的困難成為阻礙貝葉斯方法廣泛應用的最大障礙.隨著計算機算力的增強，計算貝葉斯統(tǒng)計技術的出現(xiàn)改變了這一切——統(tǒng)計推斷可以基于后驗分布抽取的隨機樣本，而且從采樣樣本中計算出的估計量可以通過設置足夠大的樣本大小來實現(xiàn)任何所需的精度[1].馬氏鏈式蒙特卡羅采樣法的基本思想是：構造一個馬氏鏈，并將后驗分布作為其穩(wěn)態(tài)分布[2-4].只要我們讓馬氏鏈運行足夠長，從鏈中提取的樣本可以被認為是從目標(后驗)分布中的樣本抽取，然后基于這些樣本就可以做出各種貝葉斯統(tǒng)計推斷.

1 蒙特卡羅類方法背后的思想

引入采樣技術的目的是獲取指定概率分布(密度)的采樣點.大數(shù)定理指出，如果能從未知的概率分布p(x)中獲取獨立同分布的采樣點xi，那么隨著采集樣本數(shù)的增加，樣本分布將收斂到真正分布[5].圖1表明，隨著采樣點的增多，每個樣本點出現(xiàn)的頻率漸近于該點的概率.

圖1 采樣點出現(xiàn)的頻率漸近于該點的實際概率(右側(cè)圖)

也就是說，樣本點x出現(xiàn)的頻率PN(x)隨著采樣次數(shù)的增多近似于其概率PN(x)[6].

(1)

在貝葉斯推斷中，常常需要計算各種統(tǒng)計量，例如求解函數(shù)f(x)關于隨機變量x的期望，且p(x):

(2)

(3)

2 獨立采樣法

利用采樣技術的目的是采集指定概率分布p(x)(目標分布)的樣本點，以用于近似計算.針對不同形式的概率分布，有不同的解決方案：

(1)如果能夠直接從目標分布抽樣，例如離散型隨機模型，或者連續(xù)型密度函數(shù)p(x)的累積分布函數(shù)存在解析形式的逆函數(shù)，借助于類似輪盤賭法或逆采樣法即可獲得獨立同分布樣本點.

(2)如果不能，只能求數(shù)值解，可以考慮用拒絕/接受或者重要采樣法.事先選定一個已知的且易抽樣(概率統(tǒng)計教科書中介紹的常見分布)的提議分布q(x)，從中抽取樣本點x(i)作為候選樣本點，并應用某種判優(yōu)規(guī)則： 若判定結果略優(yōu)于均勻采樣，那么此樣本就可以接受為p(x)的樣本點，否則拒絕.這種抽樣的準確率很大程度上取決于p(x)和q(x)的相似度，若相差很大，拒絕率會很高，特別是高維分布.

圖2 拒絕采樣法

算法思想：尋找易采樣的提議分布q(x)，確定一個常數(shù)M，使得Mq(x)完全“裹”住p(x).重復如下步驟，直至獲取N個樣本點：

(1)采樣x(i)～q(x),u～Unif(0,1)；(2)如果u接受x(i)；否則拒絕，轉(zhuǎn)(1).

圖2闡釋了拒絕采樣的基本思想[7]：對于q(x)產(chǎn)生的采樣點x(i)，由于p(x(i))≥uMq(x(i))，所以x(i)作為目標分布p(x)的高密度點(High Density Area，高密區(qū))被選中，作為低密點的x(j)卻被拒絕.事實上，低密點代表稀有事件，高密點才是較具代表性的樣本點.但由于u是均勻分布，一定概率的低密度點同樣也可能被選中，以滿足充分采樣的需要，否則我們采樣得到的樣本平均值將高于從后驗分布中抽取獨立樣本的均值.這種采樣法的缺點一是M較難確定，二是要求提議分布q(x)的目標分布p(x)形狀相似(兩者的高密區(qū)低密區(qū)相似)，否則拒絕率較高.

3 馬氏鏈式蒙特卡羅(MCMC)類采樣法

MCMC類采樣法通過設計轉(zhuǎn)移矩陣，生成一條馬爾可夫鏈使其極限分布等于目標分布.主要包括MH算法和Gibbs算法.MH算法適用于一維隨機變量的采樣，對于高維分布需要使用Gibbs算法，該方法是前者的特例.

3.1 馬爾可夫鏈

馬爾可夫鏈(簡稱馬氏鏈)是一種描述隨機過程的數(shù)學模型.從t0時刻起，以一定的概率轉(zhuǎn)移到狀態(tài)空間S={s1,s2,…,sn}的某一個，以后各個時刻均如此.馬氏鏈具有“無記憶”特點，即給定過程的過去和現(xiàn)在狀態(tài)，將來狀態(tài)只取決于現(xiàn)在的狀態(tài)，又稱之為馬爾可夫?qū)傩?由于馬氏鏈的轉(zhuǎn)移概率僅取決于當前的狀態(tài)，可為之建立一步狀態(tài)轉(zhuǎn)移概率矩陣T(si→sj).我們把t0時刻初始狀態(tài)和以后各時刻的狀態(tài)拉成一條鏈，稱為馬氏鏈.

我們關心一類特殊的馬氏鏈.如果馬氏鏈可以從每個其他狀態(tài)到達每個狀態(tài)，稱為不可約馬爾可夫鏈.一個具有所有正返回狀態(tài)的不可約馬氏鏈稱為可遍歷的馬氏鏈[8].我們將只使用可遍歷馬氏鏈，因為它們具有唯一的穩(wěn)態(tài)分布.

假設p(x)是我們的目標分布函數(shù)，我們擬獲取一系列服從該分布的采樣點.MCMC類算法將采樣點的產(chǎn)生過程構成一個可遍歷馬爾可夫鏈，它呈平穩(wěn)分布且收斂于穩(wěn)態(tài)p(x)，所以采樣點序列服從p(x).p(x)是已知的，問題的關鍵是如何構造可遍歷馬爾可夫鏈.

平穩(wěn)狀態(tài)及細節(jié)平衡條件.假如我們模擬大量物理粒子在馬氏鏈中的隨機行為[9]，粒子不斷地從一個狀態(tài)躍遷到另一狀態(tài)，最后這些移動將達到一個動態(tài)平衡，稱為平穩(wěn)狀態(tài).在達到平穩(wěn)狀態(tài)后，從某個狀態(tài)x流出去的粒子數(shù)，將會等于流回該狀態(tài)的粒子數(shù)，這時系統(tǒng)滿足 “一般平衡條件”(公式(4)的第一個等式)：

(4)

公式(4)推導過程中的最后一個等式表明此時系統(tǒng)處在平穩(wěn)狀態(tài)，即分布p(x)對馬氏鏈是不變的，這意味著轉(zhuǎn)移概率不會改變p(·)分布.為方便實現(xiàn)，我們試著把條件加強，由此導出公式(5)的細節(jié)平衡條件(detailed balance)：

P(x)T(x→y)=P(y)T(y→x)

(5)

可以證明，滿足細節(jié)平衡條件必然滿足一般平衡條件.

MCMC類算法有兩種主要方法：Metropolis-Hastings算法和吉布斯抽樣算法.

3.2 Metropolis算法

算法思想：引進一個提議分布q(x)，抽取候選狀態(tài)(樣本點)x*～q(x*丨x(t)).評估目標分布是否在候選狀態(tài)x*附近具有足夠大的密度，若滿足，則接受x*，并將其設置為下一個狀態(tài).如果候選狀態(tài)密度低于當前狀態(tài)x(t)，那么很可能(但不保證)它將被拒絕.接受或拒絕候選狀態(tài)的標準由以下啟發(fā)式定義：

(1)如果p(x*)≥p(x(t))，接受x*作為目標分布的樣本點，并作為提議分布的新狀態(tài)，意味著上述操作將鼓勵跳轉(zhuǎn)到目標分布的高密區(qū)狀態(tài).

需要注意的是，Metropolis算法中唯一限制是提議分布必須是對稱的，滿足條件的這類分布有正態(tài)分布、柯西分布、t分布和均勻分布.但是如果目標分布的支持集是不對稱的，例如x∈(0，+∞)，我們將考慮使用Metropolis-Hastings算法.

3.3 Metropolis-Hastings算法

為了能夠使用非對稱的提議分布，Metropolis-Hastings算法定義了附加校正因子c，

校正因子c調(diào)整轉(zhuǎn)移算子，以確保x(t)→x*的概率等于x*→x(t)的概率，而不管提議分布是否對稱.此時候選狀態(tài)的x*接受率為

(6)

下面證明構造的馬氏鏈中狀態(tài)x(t)→x*的概率等于狀態(tài)x*→x(t)的概率，即滿足細節(jié)平衡條件.為行文方便，引入以下符號：x?x(t)，提議分布q(x,x*)?q(x*|x(t))，狀態(tài)轉(zhuǎn)移算子T(x,x*)?T(x(t)→x*)，則接受率

首先有T(x,x*)=q(x,x*)α(x,x*)，T(x*,x)=q(x*,x)α(x*,x)，只需證明p(x)T(x,x*)=p(x*)T(x*,x)，即滿足細節(jié)平衡條件.

(1)當p(x*)q(x*,x)=p(x)q(x,x*)時，α(x,x*)=α(x*,x)=1，

(7)

比較公式(7)兩組等式的最后項，可知p(x)T(x,x*)=p(x*)T(x*,x)，滿足細節(jié)平衡條件.

(2)當p(x*)q(x*,x)>p(x)q(x,x*)時，α(x,x*)=1，α(x*,x)<1，

(8)

比較公式(8)兩組等式的最后項，可知p(x)T(x,x*)=p(x*)T(x*,x)，同樣滿足細節(jié)平衡條件.

(3)證明過程同(2).

注意到Metropolis算法是Metropolis-Hastings算法的特例，因為提議分布是對稱的，所以q(x(t)|x*)=q(x*|x(t))，校正因子c=1.

算法描述如下：

(a)從提議分布抽樣候選樣本點x*～q(x*|x(t))；

(d)從均勻分布抽取u～Unif(0,1)，如果u≤α接受x*作為新狀態(tài)，否則x*=x(t)，轉(zhuǎn)(a).

經(jīng)過充分一段時間的迭代(也稱burned in，煅燒成型時間)，馬氏鏈最終進入并收斂到平穩(wěn)狀態(tài)，只是不同的初始狀態(tài)，煅燒成型經(jīng)歷的時間不同.

3.4 Gibbs采樣法(Gibbs Sampler)

生成一維隨機變量并不困難，對于高維分布p(x)，x=(x1,x2,…,xn)，生成各分量非獨立的隨機向量就非常困難.Gibbs采樣法把一次一個n維樣本變?yōu)椤皀次一維”樣本，但前提是我們事先知道完全條件概率分布p(xj|x1,…,xj-1,xj+1…,xn)(數(shù)學上稱為半共軛，為方便常簡寫為p(xj|p(x-j)).不同于MH算法，該算法不引入新的提議分布q(·)，僅通過條件分布得到以給定分布p(x)為不變分布的馬氏鏈的轉(zhuǎn)移概率.

令提議分布q(·)為

算法描述從略[10].

需要說明的是，MCMC類采樣技術獲取的樣本點違反樣本獨立性假設，因為樣本點間存在自相關(autocorrelation)，這是由馬氏鏈的固有性質(zhì)決定的——采樣點依賴于前一個.為了減少這種輕微依賴性，第一種方法采用刪減(thinning)技術，對于生成的馬氏鏈每隔K個保留一個，其余的全丟棄；第二種塊采樣(Blocked Gibbs Sampler)技術，若xi,xj在給定xk的前提下并不條件獨立，可將xi,xj兩個變量組合在一起，并根據(jù)它們的聯(lián)合分布對所有其他變量進行采樣，而不單獨采樣每個變量；第三種折疊采樣(CollapsedGibbsSampler)技術[11]，若在給定xk的前提下，xi,xj不條件獨立，但在xk未知的情況下，xi,xj條件獨立，此時可積分掉xk，使得xi,xj條件獨立，相關內(nèi)容請參照概率圖模型[4].

盡管我們采取上述技術可以減少樣本依賴性，但不能完全消除.自相關性的存在將降低采樣精度，因為高自相關同樣意味著我們采樣得到的樣本平均值將高于從后驗分布中抽取獨立樣本的均值.

4 實驗

假設葉斯后驗分布有如下形式：

聯(lián)合分布密度函數(shù)等高線見圖3.

采用Gibbs 采樣技術，先求條件概率分布：

(9)

這里A(y),B(x)是關于y,x函數(shù)，當指定y,x時，它們是常數(shù).很明顯，公式(9)中隨機變量x,y都符合正態(tài)分布，即

圖3 聯(lián)合分布密度函數(shù)的等高線圖

圖4 Bibbs采樣點軌跡圖

編寫并運行Gibbs采樣法，初值設定在(1，6)，采樣4*105次.為消除自相關，采用刪減技術，K=20.觀察采樣點軌跡圖(圖4)，采樣點大多集中在高密區(qū)，它們是較具代表性的樣本點.

5 結語

采樣技術提供給統(tǒng)計學家一個不同的貝葉斯推斷方法，把他們從因數(shù)學計算的方便性而被限制到那些可以從理論分析中找到后驗分布的模型(例如共軛分布)中解放出來，從而可以選擇使用基于觀測模型得到的更現(xiàn)實的先驗分布.MCMC類采樣法是了解其參數(shù)的后驗分布的、對復雜模型進行推理的有效方法，這使得統(tǒng)計學家只需關注模型的設計而不必擔心計算能力，該方法的靈活運用大大提高了應用統(tǒng)計學、模式識別、機器學習理論領域處理來自現(xiàn)實世界復雜模型的能力.