3.5　馬爾科夫鏈蒙特卡洛

在以貝葉斯方法為基礎(chǔ)的各種機(jī)器學(xué)習(xí)技術(shù)中，常常需要計(jì)算后驗(yàn)概率，再通過最大后驗(yàn)概率（MAP）的方法進(jìn)行參數(shù)推斷和決策。然而，在很多時(shí)候，后驗(yàn)分布的形式可能非常復(fù)雜，這個(gè)時(shí)候?qū)ふ移渲械淖畲蠛篁?yàn)估計(jì)或者對后驗(yàn)概率進(jìn)行積分等計(jì)算往往非常困難，此時(shí)可以通過采樣的方法來求解。這其中非常重要的一個(gè)基礎(chǔ)就是馬爾科夫鏈蒙特卡洛（Markov chain Monte Carlo，MCMC）。

3.5.1　重要性采樣

前面已經(jīng)詳細(xì)介紹了基于隨機(jī)算法進(jìn)行定積分求解的技術(shù)。這里主要用到其中的平均值法。這里僅做簡單回顧。

在計(jì)算定積分時(shí)，會(huì)把g（x）拆解成兩項(xiàng)的乘積，即g（x）=f（x）p（x），其f（x）是某種形式的函數(shù)，而p（x）是關(guān)于隨機(jī)變量X的概率分布（也就是PDF或PMF）。如此一來，上述定積分就可以表示為求f（x）的期望，即

當(dāng)然，g（x）的分布函數(shù)可能具有很復(fù)雜的形式，仍然無法直接求解，這時(shí)就可以用采樣的方法去近似。這時(shí)積分可以表示為

在貝葉斯分析中，蒙特卡洛積分運(yùn)算常常被用來在對后驗(yàn)概率進(jìn)行積分時(shí)做近似估計(jì)。比如要計(jì)算

便可以使用下面這個(gè)近似形式

不難發(fā)現(xiàn)，在利用蒙特卡洛法進(jìn)行積分求解時(shí)，非常重要的一個(gè)環(huán)節(jié)就是從特定的分布中進(jìn)行采樣。這里的“采樣”意思也就是生成滿足某種分布的觀測值。之前已經(jīng)介紹過“逆采樣”和“拒絕采樣”等方法。

采樣的目的很多時(shí)候都是為了做近似積分運(yùn)算，前面的采樣方法（逆采樣和拒絕采樣）都是先對分布進(jìn)行采樣，然后再用采樣的結(jié)果近似計(jì)算積分。下面要介紹的另外一種方法“重要性采樣”（Importance sampling）則兩步并做一步，直接做近似計(jì)算積分。

現(xiàn)在的目標(biāo)是計(jì)算下面的積分

E［f（x）］=∫f（x）p（x）dx

按照蒙特卡洛求定積分的方法，將從滿足p（x）的概率分布中獨(dú)立地采樣出一系列隨機(jī)變量x_i，然后便有

但是現(xiàn)在的困難是對滿足p（x）的概率分布進(jìn)行采樣非常困難，畢竟實(shí)際中很多p（x）的形式都相當(dāng)復(fù)雜。這時(shí)該怎么做呢？于是想到做等量變換，于是有

如果把其中的看成是一個(gè)新的函數(shù)h（x），則有

其中被稱為是x_i的權(quán)重或重要性權(quán)重（Importance Weights）。所以這種采樣的方法就被稱為是重要性采樣（Importance Sampling）。

如圖3-19所示，在使用重要性采樣時(shí)，并不會(huì)拒絕掉某些采樣點(diǎn)，這與在使用拒絕采樣時(shí)不同。此時(shí)，所有的采樣點(diǎn)都將為我們所用，但是它們的權(quán)重是不同的。因?yàn)闄?quán)重為，所以在圖3-19中，當(dāng)p（x_i）＞q（x_i）時(shí)，采樣點(diǎn)x_i的權(quán)重就大（深色線在淺色線上方時(shí)），反之就小。重要性采樣就是通過這種方式來從一個(gè)“參考分布”q（x）中獲得“目標(biāo)分布”p（x）的。

3.5.2　馬爾科夫鏈蒙特卡洛的基本概念

馬爾科夫鏈蒙特卡洛方法是一類用于從一個(gè)概率分布中進(jìn)行采樣的算法，該類方法以構(gòu)造馬爾科夫鏈為基礎(chǔ)，而被構(gòu)造的馬爾科夫鏈分布應(yīng)該同需要之分布等價(jià)。

MCMC構(gòu)造馬爾科夫鏈，使其穩(wěn)態(tài)分布等于要采樣的分布，這樣就可以通過馬爾科夫鏈來采樣。這種等價(jià)如何理解是深入探討具體操作方法之前，需要先攻克的一個(gè)問題。在此之前，希望讀者對馬爾科夫鏈已經(jīng)有了一個(gè)比較清晰的認(rèn)識(shí)。

現(xiàn)在，用下面的式子來表示每一步（時(shí)刻推進(jìn)）中從狀態(tài)s_i到狀態(tài)s_j的轉(zhuǎn)移概率

p（i，j）=p（i→j）=P（X_t+1=s_j|X_t=s_i）

這里的一步是指時(shí)間從時(shí)刻t過渡到下一時(shí)刻t+1。

馬爾科夫鏈在時(shí)刻t處于狀態(tài)j的概率可以表示為π_j（t）=P（X_t=s_j）。這里用向量π（t）來表示在時(shí)刻t各個(gè)狀態(tài)的概率向量。在初始時(shí)刻，需要選取一個(gè)特定的π（0），通常情況下可以使向量中一個(gè)元素為1，其他元素均為0，即從某個(gè)特定的狀態(tài)開始。隨著時(shí)間的推移，狀態(tài)會(huì)通過轉(zhuǎn)移矩陣，向各個(gè)狀態(tài)擴(kuò)散。

馬爾科夫鏈在t+1時(shí)刻處于狀態(tài)s_i的概率可以通過t時(shí)刻的狀態(tài)概率和轉(zhuǎn)移概率來求得，并可通過查普曼-柯爾莫哥洛夫等式得到

用轉(zhuǎn)移矩陣寫成矩陣的形式如下

π（t+1）=π（t）P

其中，轉(zhuǎn)移矩中的元素p（i，j）表示p（i→j）。因此，π（t）=π（0）P^t。此外，表示矩陣Pⁿ中第ij個(gè)元素，即。

一條馬爾科夫鏈有一個(gè)平穩(wěn)分布π^?，是指給定任意一個(gè)初始分布π（0），經(jīng)過有限步之后，最終都會(huì)收斂到平穩(wěn)分布π^?。平穩(wěn)分布具有性質(zhì)π^?=π^?P。可以結(jié)合本章前面談到的矩陣極限的概念來理解馬爾科夫鏈的平穩(wěn)分布。

一條馬爾科夫鏈擁有平穩(wěn)分布的一個(gè)充分條件是對于任意兩個(gè)狀態(tài)i和j，其滿足細(xì)致平衡（Detailed Balance）：p（j→i）π_j=p（i→j）π_i。可以看到，此時(shí)

所以π=πP，明顯滿足平穩(wěn)分布的條件。如果一條馬爾科夫鏈滿足細(xì)致平衡，就說它是可逆的。

最后來總結(jié)一下MCMC的基本思想。在拒絕采樣和重要性采樣中，當(dāng)前生成的樣本點(diǎn)與之前生成的樣本點(diǎn)之間是沒有關(guān)系的，它的采樣都是獨(dú)立進(jìn)行的。然而，MCMC是基于馬爾科夫鏈進(jìn)行的采樣。這就表明，當(dāng)前的樣本點(diǎn)生成與上一時(shí)刻的樣本點(diǎn)是有關(guān)的。如圖3-20所示，假設(shè)當(dāng)前時(shí)刻生成的樣本點(diǎn)是x，下一次時(shí)刻采樣到x′的（轉(zhuǎn)移）概率就是p（x|x′），或者寫成p（x→x′）。我們希望的是這個(gè)過程如此繼續(xù)下去，最終的概率分布收斂到π（x），也就是要采樣的分布。

易見，轉(zhuǎn)移概率p（x|x′）與采樣分布π（x）之間必然存在某種聯(lián)系，因?yàn)橄Ｍ勒?i>p（x|x′）來進(jìn)行轉(zhuǎn)移采樣最終能收斂到π（x）。那這個(gè)關(guān)系到底是怎么樣的呢？首先，給定所有可能的x，然后求從x轉(zhuǎn)移到x′的概率，所得之結(jié)果就是x′出現(xiàn)的概率，如果用公式表示即

π（x′）=∫π（x）p（x′|x）dx

這其實(shí)也就是查普曼-柯爾莫哥洛夫等式所揭示的關(guān)系。如果馬爾科夫鏈?zhǔn)瞧椒€(wěn)的，那么基于該馬爾科夫鏈的采樣過程最終將收斂到平穩(wěn)狀態(tài)

π^?（x′）=∫π^?（x′）p（x|x′）dx′

而平穩(wěn)狀態(tài)的充分條件又是滿足細(xì)致平衡

π（x）p（x′|x）dx=π（x）p（x|x′）

可見細(xì)致平衡的意思是說從x轉(zhuǎn)移到x′的概率應(yīng)該等于從x′轉(zhuǎn)移到x的概率，在這樣的情況下，采樣過程將最終收斂到目標(biāo)分布。也就是說，設(shè)計(jì)的MCMC算法，只要驗(yàn)證其滿足細(xì)致平衡的條件，那么用這樣的方法做基于馬爾科夫鏈的采樣，最終就能收斂到一個(gè)穩(wěn)定狀態(tài)，即目標(biāo)分布。

實(shí)際應(yīng)用中，有兩個(gè)馬爾科夫鏈蒙特卡洛采樣算法十分常用，它們是Metropolis-Hastings（梅特羅波利斯-黑斯廷斯）算法與Gibbs Sampling采樣（吉布斯），可以證明吉布斯采樣是梅特羅波利斯-黑斯廷斯算法的一種特殊情況，而且兩者都滿足細(xì)致平衡的條件。

3.5.3　Metropolis-Hastings算法

對給定的概率分布，如果從中直接采樣比較困難，那么Metropolis-Hastings算法就是一個(gè)從中獲取一系列隨機(jī)樣本的MCMC方法。這里所說的Metropolis-Hastings算法和模擬退火方法中所使用的Metropolis-Hastings算法本質(zhì)上是一回事。用于模擬退火和MCMC的原始算法最初是由Metropolis提出的，后來Hastings對其進(jìn)行了推廣。Hastings提出沒有必要使用一個(gè)對稱的參考分布（proposal distribution）。當(dāng)然達(dá)到均衡分布的速度與參考分布的選擇有關(guān)。

Metropolis-Hastings算法的執(zhí)行步驟如下。

Metropolis-Hastings算法的執(zhí)行步驟是先隨便指定一個(gè)初始的樣本點(diǎn)x⁰，u是來自均勻分布的一個(gè)隨機(jī)數(shù)，x^?是來自p分布的樣本點(diǎn)，樣本點(diǎn)是否從前一個(gè)位置跳轉(zhuǎn)到x^?，由α（x^?）和u的大小關(guān)系決定。如果接受，則跳轉(zhuǎn)，即新的采樣點(diǎn)為x^?，否則就拒絕，即新的采用點(diǎn)仍為上一輪的采樣點(diǎn)。這個(gè)算法看起來非常簡單，難免讓人疑問照此步驟來執(zhí)行，最終采樣分布能收斂到目標(biāo)分布嗎？根據(jù)之前的講解，可知要想驗(yàn)證這一點(diǎn)，就需要檢查細(xì)致平衡是否滿足。下面是具體的證明過程，不再贅言。

既然已經(jīng)從理論上證明Metropolis-Hastings算法的確實(shí)可行，下面舉一個(gè)簡單的例子來看看實(shí)際中Metropolis-Hastings算法的效果。稍微有點(diǎn)不一樣的地方是，這里并未出現(xiàn)π（x）p（x′|x）這樣的后驗(yàn)概率形式，作為一個(gè)簡單的示例，下面只演示從柯西分布中進(jìn)行采樣的做法。

柯西分布也叫作柯西-洛倫茲分布，它是以奧古斯丁·路易·柯西與亨德里克·洛倫茲名字命名的連續(xù)概率分布，其概率密度函數(shù)為

其中，x₀是定義分布峰值位置的位置參數(shù)，γ是最大值一半處的一半寬度的尺度參數(shù)。x₀=0且γ=1的特例稱為標(biāo)準(zhǔn)柯西分布，其概率密度函數(shù)為

下面是在R中進(jìn)行基于Metropolis-Hastings算法對柯西分布進(jìn)行采樣的示例代碼。

圖3-21所示為每次采樣點(diǎn)的數(shù)值分布情況。

為了更清晰明確地看到采樣結(jié)果符合預(yù)期分布，這里也繪制了分布的密度圖并與實(shí)際的分布密度進(jìn)行對照，如圖3-22所示，深色實(shí)線是采樣分布的密度圖，而淺色虛線則是實(shí)際柯西分布的密度圖，可見兩者吻合地相當(dāng)好。

當(dāng)然，作為一個(gè)簡單例子的開始，上面的例子中并沒有涉及p（y|x），接下來的這個(gè)例子則會(huì)演示涉及后驗(yàn)概率的基于Metropolis-Hastings算法的MCMC。這里將用Metropolis-Hastings算法從瑞利分布（Rayleigh Distribution）中采樣，瑞利分布的密度為

然后取自由度為X_t的卡方分布為參考分布。實(shí)現(xiàn)代碼如下。

然后要驗(yàn)證一下，生產(chǎn)的采樣數(shù)據(jù)是否真的符合瑞利分布。注意在R中使用瑞利分布的相關(guān)函數(shù)，需要加裝VGAM包。下面的代碼可以讓讀者直觀地感受到采樣結(jié)果的分布情況。

從圖3-23中不難看出，采樣點(diǎn)分布確實(shí)符合預(yù)期。當(dāng)然，這僅僅是一個(gè)演示用的小例子。顯然，它并不高效。因?yàn)椴捎米杂啥葹?i>X_t的卡方分布來作為參考函數(shù)，大約有40%的采樣點(diǎn)都被拒絕掉了。如果換做其他更加合適的參考函數(shù)，可以大大提升采樣的效率。

3.5.4　Gibbs采樣

Gibbs采樣是一種用以獲取一系列觀察值的MCMC算法，這些觀察值近似于來自一個(gè)指定的多變量概率分布，但從該分布中直接采樣較為困難。Gibbs采樣可以被看成是Metropolis-Hastings算法的一個(gè)特例。而這個(gè)特殊之處就在于，使用Gibbs采樣時(shí)，通常是對多變量分布進(jìn)行采樣的。比如說，現(xiàn)在有一個(gè)分布p（z₁，z₂，z₃），可見它是定義在三個(gè)變量上的分布。這種分布在實(shí)際中還有很多，比如一維的正態(tài)分布也是關(guān)于其均值和方差這兩個(gè)變量的分布。在算法的第t步，選出了三個(gè)值，，。這時(shí)，既然和是已知的，那么可以由此獲得一個(gè)新的，并用這個(gè)新的值替代原始的，而新的可以由下面這個(gè)條件分布來獲得

接下來同樣的道理再用一個(gè)新的來替代原始的，而這個(gè)新的可以由下面這個(gè)條件分布來獲得

可見剛剛獲得的新值已經(jīng)被直接用上了。同理，最后用來更新，而新的值由下面這個(gè)條件分布來獲得

按照這個(gè)順序如此往復(fù)即可。下面給出更為完整的吉布斯采樣的算法描述。

當(dāng)然如果要從理論上證明Gibbs采樣確實(shí)可以得到預(yù)期的分布，就應(yīng)該考察它是否滿足細(xì)致平衡。但是前面也講過，Gibbs采樣是Metropolis-Hastings算法的一個(gè)特例。所以其實(shí)甚至無需費(fèi)力去考察其是否滿足細(xì)致平衡，如果能夠證明吉布斯采樣就是Metropolis-Hastings算法，而Metropolis-Hastings算法是滿足細(xì)致平衡，其實(shí)也就得到了我們想要的。下面是證明的過程，其中x_-k表示一個(gè)向量（x₁，…，x_k-1，x_k+1，…，x_D），也就是從中剔除了x_k元素。

令x=x₁，…，x_D。

當(dāng)采樣第k個(gè)元素，

當(dāng)采樣第k個(gè)元素，

以Gibbs采樣為基礎(chǔ)實(shí)現(xiàn)的MCMC在自然語言處理中的LDA里有重要應(yīng)用。如果讀者正在或者后續(xù)準(zhǔn)備研究LDA的話，這些內(nèi)容將是非常重要的基礎(chǔ)。