房愛(ài)東，謝士春

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

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

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

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

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

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

(1)

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

(2)

(3)

2 獨(dú)立采樣法

利用采樣技術(shù)的目的是采集指定概率分布p(x)(目標(biāo)分布)的樣本點(diǎn)，以用于近似計(jì)算.針對(duì)不同形式的概率分布，有不同的解決方案：

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

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

圖2 拒絕采樣法

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

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

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

3 馬氏鏈?zhǔn)矫商乜_(MCMC)類采樣法

MCMC類采樣法通過(guò)設(shè)計(jì)轉(zhuǎn)移矩陣，生成一條馬爾可夫鏈?zhǔn)蛊錁O限分布等于目標(biāo)分布.主要包括MH算法和Gibbs算法.MH算法適用于一維隨機(jī)變量的采樣，對(duì)于高維分布需要使用Gibbs算法，該方法是前者的特例.

3.1 馬爾可夫鏈

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

我們關(guān)心一類特殊的馬氏鏈.如果馬氏鏈可以從每個(gè)其他狀態(tài)到達(dá)每個(gè)狀態(tài)，稱為不可約馬爾可夫鏈.一個(gè)具有所有正返回狀態(tài)的不可約馬氏鏈稱為可遍歷的馬氏鏈[8].我們將只使用可遍歷馬氏鏈，因?yàn)樗鼈兙哂形ㄒ坏姆€(wěn)態(tài)分布.

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

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

(4)

公式(4)推導(dǎo)過(guò)程中的最后一個(gè)等式表明此時(shí)系統(tǒng)處在平穩(wěn)狀態(tài)，即分布p(x)對(duì)馬氏鏈?zhǔn)遣蛔兊?，這意味著轉(zhuǎn)移概率不會(huì)改變p(·)分布.為方便實(shí)現(xiàn)，我們?cè)囍褩l件加強(qiáng)，由此導(dǎo)出公式(5)的細(xì)節(jié)平衡條件(detailed balance)：

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

(5)

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

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

3.2 Metropolis算法

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

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

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

3.3 Metropolis-Hastings算法

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

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

(6)

下面證明構(gòu)造的馬氏鏈中狀態(tài)x(t)→x*的概率等于狀態(tài)x*→x(t)的概率，即滿足細(xì)節(jié)平衡條件.為行文方便，引入以下符號(hào)：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)，即滿足細(xì)節(jié)平衡條件.

(1)當(dāng)p(x*)q(x*,x)=p(x)q(x,x*)時(shí)，α(x,x*)=α(x*,x)=1，

(7)

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

(2)當(dāng)p(x*)q(x*,x)>p(x)q(x,x*)時(shí)，α(x,x*)=1，α(x*,x)<1，

(8)

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

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

注意到Metropolis算法是Metropolis-Hastings算法的特例，因?yàn)樘嶙h分布是對(duì)稱的，所以q(x(t)|x*)=q(x*|x(t))，校正因子c=1.

算法描述如下：

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

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

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

3.4 Gibbs采樣法(Gibbs Sampler)

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

令提議分布q(·)為

算法描述從略[10].

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

盡管我們采取上述技術(shù)可以減少樣本依賴性，但不能完全消除.自相關(guān)性的存在將降低采樣精度，因?yàn)楦咦韵嚓P(guān)同樣意味著我們采樣得到的樣本平均值將高于從后驗(yàn)分布中抽取獨(dú)立樣本的均值.

4 實(shí)驗(yàn)

假設(shè)葉斯后驗(yàn)分布有如下形式：

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

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

(9)

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

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

圖4 Bibbs采樣點(diǎn)軌跡圖

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

5 結(jié)語(yǔ)

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