【原】Python用MCMC馬爾科夫鏈蒙特卡洛、拒絕抽樣和Metropolis-Hastings采樣算法

拓端數(shù)據(jù) 2022-06-20 發(fā)布于上海

展開(kāi)全文

原文鏈接：http:///?p=27267

我們將研究?jī)煞N對(duì)分布進(jìn)行抽樣的方法：拒絕抽樣和使用 Metropolis Hastings 算法的馬爾可夫鏈蒙特卡洛方法 (MCMC)。像往常一樣，我將提供直觀的解釋、理論和一些帶有代碼的示例。

相關(guān)視頻

背景

在我們進(jìn)入主題之前，讓我們將馬爾可夫鏈蒙特卡羅（MCMC）這個(gè)術(shù)語(yǔ)分解為它的基本組成部分：蒙特卡羅方法和馬爾可夫鏈。

馬爾可夫鏈

Markov Chain 是“在狀態(tài)空間上經(jīng)歷從一種狀態(tài)到另一種狀態(tài)的轉(zhuǎn)換的隨機(jī)過(guò)程”。

正如你所看到的，它看起來(lái)就像一個(gè)有限狀態(tài)機(jī)，只是我們用概率注釋了狀態(tài)轉(zhuǎn)換。例如，我們可以查看今天是否晴天，明天晴天的概率為 0.9，下雨的概率為 0.1。同樣在雨天開(kāi)始。應(yīng)該清楚的是，從給定的狀態(tài)，所有傳出的轉(zhuǎn)換應(yīng)該總計(jì) 1.0，因?yàn)樗且粋€(gè)適當(dāng)?shù)姆植肌?/p>

表示此信息的另一種方法是通過(guò)轉(zhuǎn)移矩陣 P：

將其表示為矩陣的有趣之處在于，我們可以通過(guò)矩陣乘法來(lái)模擬馬爾可夫鏈。例如，假設(shè)我們從陽(yáng)光明媚的狀態(tài)開(kāi)始，我們可以將其表示為行向量：x(0)=[10]。這隱含地表示我們處于晴天狀態(tài)的概率為 1，因此處于下雨?duì)顟B(tài)的概率為 0?，F(xiàn)在，如果我們執(zhí)行矩陣乘法，我們可以在一步之后找出處于每個(gè)狀態(tài)的概率：

我們可以看到明天有 0.9 的機(jī)會(huì)晴天（根據(jù)我們的簡(jiǎn)單模型），有 0.1 的機(jī)會(huì)下雨。實(shí)際上，我們可以繼續(xù)將轉(zhuǎn)換矩陣相乘，以在 k 步之后找到太陽(yáng)/雨的機(jī)會(huì)：

我們可以很容易地計(jì)算 x(k) 的各種 k 值，使用 numpy：

import numpy as npPa = nps.ardraasy(\[\[0.9, 0.1\], \[0.5, 0.5\]\])istsdatea = np.arasdray(\[1, 0\])siasdmulatase_markasdov(istaasdte, P, 10)

我們可以在這里看到一個(gè)有趣的現(xiàn)象，當(dāng)我們?cè)跔顟B(tài)機(jī)中采取更多步驟時(shí)，晴天或下雨的概率似乎會(huì)收斂。您可能認(rèn)為這與我們所處的初始狀態(tài)有關(guān)，但實(shí)際上并非如此。如果我們將初始狀態(tài)初始化為隨機(jī)值，我們將得到相同的結(jié)果

siasmdulasteds_marksov(nap.arsdray(\[r, 1 - r\]), P, 10)

這種穩(wěn)態(tài)分布稱(chēng)為 stationary distribution 通常用 π 表示。可以通過(guò)多種方式找到該穩(wěn)態(tài)向量π。最直接的方法是在 n 接近無(wú)窮大時(shí)取極限。

下一種方法就是求解方程。由于根據(jù)定義 q是穩(wěn)態(tài)，因此乘以 P 應(yīng)該返回相同的值：

其中 I 是單位矩陣。如果您擴(kuò)展我們的向量/矩陣符號(hào)，您會(huì)發(fā)現(xiàn)這只是一個(gè)方程組以及 π1,π2,...,πn 總和為 1 (即 π 形成概率分布）。在我們的例子中只有兩個(gè)狀態(tài)：π1+π2=1。

然而，并不是每個(gè)馬爾可夫鏈都有一個(gè)平穩(wěn)的分布，甚至是唯一的但是，如果我們向馬爾可夫鏈添加兩個(gè)額外的約束，我們可以保證這些屬性：

_不可約_：我們必須能夠最終從任何其他狀態(tài)到達(dá)任何一種狀態(tài)（即期望步數(shù)是有限的）。
_非周期性_：系統(tǒng)永遠(yuǎn)不會(huì)返回到具有固定周期的相同狀態(tài)（例如，不會(huì)每 5 步確定性地返回開(kāi)始“晴天”）。

一個(gè)重要的定理說(shuō)，如果馬爾可夫鏈?zhǔn)潜闅v的，那么它有一個(gè)唯一的穩(wěn)態(tài)概率向量 π。在 MCMC 的上下文中，我們可以從任何狀態(tài)跳轉(zhuǎn)到任何其他狀態(tài)（以一定的有限概率），輕松滿足不可約性。

我們將使用的另一個(gè)有用的定義是 detailed balance and reversible Markov Chains. 如果存在滿足此條件的概率分布 π，則稱(chēng)馬爾可夫鏈?zhǔn)强赡娴模ㄒ卜Q(chēng)為詳細(xì)平衡條件）：

換句話說(shuō)，從長(zhǎng)遠(yuǎn)來(lái)看，你從狀態(tài) i 轉(zhuǎn)移到狀態(tài) j 的次數(shù)比例，與你從狀態(tài) j 轉(zhuǎn)移到狀態(tài) i 的次數(shù)比例相同。事實(shí)上，如果馬爾可夫鏈?zhǔn)强赡娴模敲次覀兙椭浪哂衅椒€(wěn)分布（這就是我們使用相同符號(hào) π 的原因）。

馬爾可夫鏈蒙特卡羅方法

馬爾可夫鏈蒙特卡羅 (MCMC) 方法只是一類(lèi)使用馬爾可夫鏈從特定概率分布（蒙特卡羅部分）中采樣的算法。他們通過(guò)創(chuàng)建一個(gè)馬爾可夫鏈來(lái)工作，其中限制分布（或平穩(wěn)分布）只是我們想要采樣的分布。

這是一張可能有助于描述該過(guò)程的圖片. 想象一下，我們正在嘗試制作一個(gè) MCMC 來(lái)嘗試使用 PDF f(x)對(duì)任意一維分布進(jìn)行采樣。在這種情況下，我們的狀態(tài)將是沿 x 軸的點(diǎn)，而我們的轉(zhuǎn)換概率將是從一種狀態(tài)到另一種狀態(tài)的機(jī)會(huì)。這是情況的簡(jiǎn)化圖：

該圖顯示了我們?cè)噲D用粗黑線逼近的密度函數(shù)，以及使用從橙色狀態(tài)過(guò)渡的藍(lán)線的馬爾可夫鏈的一部分的可視化。特別是，對(duì)于 i=-3,-2,-1,1,2,3，只是從狀態(tài) X0 到 Xi 的轉(zhuǎn)換。但是，x 軸線上的每個(gè)點(diǎn)實(shí)際上都是這個(gè)馬爾可夫鏈中的一個(gè)勢(shì)態(tài)。請(qǐng)注意，這意味著我們有一個(gè)無(wú)限的狀態(tài)空間，因此我們不能再將轉(zhuǎn)換很好地表示為矩陣。MCMC 方法的真正“訣竅”是我們想要設(shè)計(jì)狀態(tài)（或 x 軸上的點(diǎn)）之間的轉(zhuǎn)換概率，以便我們將大部分時(shí)間花在 f(x) 很大的區(qū)域中，并且在它較小的區(qū)域中的時(shí)間相對(duì)較少（即與我們的密度函數(shù)成精確比例）。

就我們的人物而言，我們希望將大部分時(shí)間花在中心周?chē)?，而較少時(shí)間花在外面。事實(shí)上，如果我們模擬我們的馬爾可夫鏈足夠長(zhǎng)的時(shí)間，狀態(tài)的限制分布應(yīng)該近似于我們?cè)噲D采樣的 PDF。因此，使用 MCMC 方法進(jìn)行采樣的基本算法為：

從任意點(diǎn) x 開(kāi)始。
以一定的轉(zhuǎn)移概率跳轉(zhuǎn)到點(diǎn) x'（這可能意味著保持相同的狀態(tài)）。
轉(zhuǎn)到第 2 步，直到我們轉(zhuǎn)換了 T 時(shí)間。
記錄當(dāng)前狀態(tài) x′，進(jìn)行步驟 2。

現(xiàn)在，我們?cè)诿總€(gè)點(diǎn) x 軸上花費(fèi)的比例次數(shù)應(yīng)該是我們?cè)噲D模擬的 PDF 的近似值，即如果我們繪制 x 值的直方圖，我們應(yīng)該得到相同的形狀。

拒絕抽樣

現(xiàn)在，在我們進(jìn)入 MCMC 方法的具體算法之前，我想介紹另一種對(duì)概率分布進(jìn)行采樣的方法，我們稍后將使用它，稱(chēng)為 rejection sampling. 主要思想是，如果我們?cè)噲D從分布 f(x) 中采樣，我們將使用另一個(gè)工具分布 g(x) 來(lái)幫助從 f(x) 中采樣。唯一的限制是對(duì)于某些 M>1，f(x)<Mg(x)。它的主要用途是當(dāng) f(x) 的形式難以直接采樣時(shí)（但仍然可以在任何點(diǎn) xx 對(duì)其進(jìn)行評(píng)估）。

以下是算法的細(xì)分：

從 g(x) 中采樣 x。
從 U(0,Mg(x)) 中采樣 y（均勻分布）。
如果 y<f(x)，則接受 x 作為 f(x) 的樣本，否則轉(zhuǎn)到步驟 1。

另一種看待它的方法是我們采樣點(diǎn) x0 的概率。這與從 g 中采樣 x0 的概率乘以我們接受的次數(shù)的比例成正比，它簡(jiǎn)單地由 f(x0) 和 Mg(x0) 之間的比率給出：

等式告訴我們對(duì)任意點(diǎn)進(jìn)行采樣的概率與 f(x0) 成正比。在對(duì)許多點(diǎn)進(jìn)行采樣并找到我們看到 x0 的次數(shù)比例之后，常數(shù) M 被歸一化，我們得到了 PDF f(x) 的正確結(jié)果。

讓我們通過(guò)一個(gè)例子更直觀地看一下它。我們要從中采樣的目標(biāo)分布 f(x) 是 double gamma 分布，基本上是一個(gè)雙邊伽馬分布。我們將使用正態(tài)分布 g(x) 作為我們的包絡(luò)分布。下面的代碼向我們展示了如何找到縮放常數(shù) M，并為我們描繪了拒絕抽樣在概念上是如何工作的。

# 目標(biāo) = 雙伽馬分布dsg = stats.dgamma(a=1)# 生成PDF的樣本x = np.linspace# 繪圖ax = df.plot(style=\['--', '-'\]

從圖中，一旦我們找到 g(x)的樣本（在這種情況下 x=2），我們從范圍等于 Mg(x) 高度的均勻分布中繪制. 如果它在目標(biāo) PDF 的高度內(nèi)，我們接受它（綠色），否則拒絕（拒絕）。

點(diǎn)擊標(biāo)題查閱往期內(nèi)容

R語(yǔ)言實(shí)現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣

左右滑動(dòng)查看更多

實(shí)施拒絕抽樣

下面的代碼為我們的目標(biāo)雙伽馬分布實(shí)現(xiàn)拒絕采樣。它繪制標(biāo)準(zhǔn)化直方圖并將其與我們應(yīng)該得到的理論 PDF 匹配。

# 從拒絕采樣算法生成樣本sdampales = \[rejeasdctioan_samplaing for x in range(10000)\]# 繪制直方圖與目標(biāo) PDFdf\['Target'\].plot

總的來(lái)說(shuō)，我們的拒絕采樣器非常適合。與理論分布相比，抽取更多樣本會(huì)改善擬合。

拒絕抽樣的很大一部分是它很容易實(shí)現(xiàn)（在 Python 中只需幾行代碼），但有一個(gè)主要缺點(diǎn)：它很慢。

Metropolis-Hastings 算法

這 Metropolis-Hastings Algorithm (MH) 是一種 MCMC 技術(shù)，它從難以直接采樣的概率分布中抽取樣本。與拒絕抽樣相比，對(duì) MH 的限制實(shí)際上更加寬松：對(duì)于給定的概率密度函數(shù) p(x)，我們只要求我們有一個(gè)與 p_成正比_的函數(shù) f(x)。這在對(duì)后驗(yàn)分布進(jìn)行采樣時(shí)非常有用。

Metropolis-Hastings 算法的推導(dǎo)

為了推導(dǎo)出 Metropolis-Hastings 算法，我們首先從最終目標(biāo)開(kāi)始：創(chuàng)建一個(gè)馬爾可夫鏈，其中穩(wěn)態(tài)分布等于我們的目標(biāo)分布 p(x)。就馬爾可夫鏈而言，我們已經(jīng)知道狀態(tài)空間將是什么：概率分布的支持，即 x 值。因此（假設(shè)馬爾可夫鏈的構(gòu)造正確）我們最終得到的穩(wěn)態(tài)分布將只是 p(x)。剩下的是確定這些 x 值之間的轉(zhuǎn)換概率，以便我們可以實(shí)現(xiàn)這種穩(wěn)態(tài)行為。

馬爾可夫鏈的詳細(xì)平衡條件，這里用另一種方式寫(xiě)成：

這里 p(x)是我們的目標(biāo)分布，P(x→x′) 是從點(diǎn) x到點(diǎn) x′ 的轉(zhuǎn)移概率。所以我們的目標(biāo)是確定P(x→x′)的形式。既然我們要構(gòu)建馬爾可夫鏈，讓我們從使用等式 5 作為該構(gòu)建的基礎(chǔ)開(kāi)始。請(qǐng)記住，詳細(xì)的平衡條件保證我們的馬爾可夫鏈將具有平穩(wěn)分布（它存在）。此外，如果我們也包括遍歷性（不以固定間隔重復(fù)狀態(tài)，并且每個(gè)狀態(tài)最終都能夠達(dá)到任何其他狀態(tài)），我們將建立一個(gè)具有唯一平穩(wěn)分布 p(x)的馬爾可夫鏈.

我們可以將方程重新排列為：

這里我們使用 f(x)來(lái)表示一個(gè) 與 p(x)_成正比的函數(shù)。_這是為了強(qiáng)調(diào)我們并不明確需要 p(x)，只是需要與它成比例的東西，這樣比率才能達(dá)到相同的效果?，F(xiàn)在這里的“技巧”是我們將把 P(x→x′)分解為兩個(gè)獨(dú)立的步驟：一個(gè)提議分布 g(x→x′) 和接受分布 A(x→x′)（類(lèi)似于拒絕抽樣的工作原理）。由于它們是獨(dú)立的，我們的轉(zhuǎn)移概率只是兩者的乘積：

此時(shí)，我們必須弄清楚 g(x)和 A(x) 的合適選擇是什么。由于 g(x) 是“建議分布”，它決定了我們可能采樣的下一個(gè)點(diǎn)。因此，重要的是它具有與我們的目標(biāo)分布 p(x)（遍歷性條件）相同的支持。這里的一個(gè)典型選擇是以當(dāng)前狀態(tài)為中心的正態(tài)分布。現(xiàn)在給定一個(gè)固定的提議分布 g(x)，我們希望找到一個(gè)匹配的 A(x)。

雖然不明顯，但滿足公式的 A(x) 的典型選擇是：

我們可以通過(guò)考慮 f(x′)g(x′→x)小于等于 1 和大于 1 的情況。當(dāng)小于等于 1 時(shí)，它的倒數(shù)大于 1，因此 LHS 的分母，A(x′→ x), 等式 8 為 1，而分子等于 RHS?；蛘撸?dāng)f(x′)g(x′→x)是大于 1 LHS 的分子是 1，而分母正好是 RHS 的倒數(shù)，導(dǎo)致 LHS 等于 RHS。

這樣，我們已經(jīng)證明，我們創(chuàng)建的馬爾可夫鏈的穩(wěn)定狀態(tài)將等于我們的目標(biāo)分布 (p(x))，因?yàn)樵敿?xì)的平衡條件通過(guò)構(gòu)造得到滿足。所以整體算法將是（與上面的 MCMC 算法非常匹配）：

通過(guò)選擇一個(gè)隨機(jī) x 來(lái)初始化初始狀態(tài)。
根據(jù)g(x→x′)找到新的x′。
根據(jù) A(x→x′) 以均勻概率接受 x′。如果接受到 x' 的轉(zhuǎn)換，否則保持狀態(tài) x。
進(jìn)行第 2 步，T 次。
將狀態(tài) x 保存為樣本，轉(zhuǎn)至步驟 2 對(duì)另一個(gè)點(diǎn)進(jìn)行采樣。

預(yù)燒和相關(guān)樣本

在我們繼續(xù)實(shí)現(xiàn)之前，我們需要討論關(guān)于 MCMC 方法的兩個(gè)非常重要的話題。第一個(gè)主題與我們選擇的初始狀態(tài)有關(guān)。由于我們隨機(jī)選擇 x 的值，它很可能位于 p(x) 非常小的區(qū)域（想想我們的雙伽馬分布的尾部）。如果從這里開(kāi)始，它可能會(huì)花費(fèi)不成比例的時(shí)間來(lái)遍歷低密度的 x 值，從而錯(cuò)誤地給我們一種感覺(jué)，即這些 x 值應(yīng)該比它們更頻繁地出現(xiàn)。所以解決這個(gè)問(wèn)題的方法是“預(yù)燒”采樣器通過(guò)生成一堆樣本并將它們?nèi)拥?。樣本的?shù)量將取決于我們?cè)噲D模擬的分布的細(xì)節(jié)。

我們上面提到的第二個(gè)問(wèn)題是兩個(gè)相鄰樣本之間的相關(guān)性。由于根據(jù)我們的轉(zhuǎn)換函數(shù) P(x→x′)的定義，繪制 x′ 取決于當(dāng)前狀態(tài) x。因此，我們失去了樣本的一項(xiàng)重要屬性：獨(dú)立性。為了糾正這一點(diǎn)，我們抽取 Tth 個(gè)樣本，并且只記錄最后抽取的樣本。假設(shè) T 足夠大，樣本應(yīng)該是相對(duì)獨(dú)立的。與預(yù)燒一樣，T 的值取決于目標(biāo)和提議分布。

實(shí)現(xiàn) Metropolis-Hastings 算法

讓我們使用上面的雙伽馬分布示例。讓我們將我們的提議分布定義為以 x 為中心、標(biāo)準(zhǔn)差為 2、N(x, 2) 的正態(tài)分布（記住 x 是當(dāng)前狀態(tài)）：

給定 f(x) 與我們的基礎(chǔ)分布 p(x) 成比例，我們的接受分布如下所示：

由于正態(tài)分布是對(duì)稱(chēng)的，因此正態(tài)分布的 PDF 在其各自點(diǎn)進(jìn)行評(píng)估時(shí)會(huì)相互抵消?，F(xiàn)在讓我們看一些代碼：

# 模擬與雙伽馬分布成比例的 f(x)f = ambd x: g.df(x* mat.i

采樣器 = mhspler()

＃ 樣本

sames = \[nex(saper) for x in range(10000)\]# 繪制直方圖與目標(biāo) PDFdf\['Target'\].plot

來(lái)自我們的 MH 采樣器的樣本很好地近似于我們的雙伽馬分布。此外，查看自相關(guān)圖，我們可以看到它在整個(gè)樣本中非常小，表明它們是相對(duì)獨(dú)立的。如果我們沒(méi)有為 T 選擇一個(gè)好的值或沒(méi)有預(yù)燒期，我們可能會(huì)在圖中看到較大的值。