Python用MCMC馬爾科夫鏈蒙特卡洛、拒絕抽樣和Metropolis-Hastings采樣算法
原文鏈接：http://tecdat.cn/?p=27267
原文出處：拓端數(shù)據(jù)部落公眾號
相關(guān)視頻：馬爾可夫鏈原理可視化解釋與R語言區(qū)制轉(zhuǎn)換Markov regime switching實(shí)例
馬爾可夫鏈原理可視化解釋與R語言區(qū)制轉(zhuǎn)換Markov?regime?switching實(shí)例
，時長07:25
相關(guān)視頻
馬爾可夫鏈蒙特卡羅方法MCMC原理與R語言實(shí)現(xiàn)
，時長08:47
我們將研究兩種對分布進(jìn)行抽樣的方法：拒絕抽樣和使用 Metropolis Hastings 算法的馬爾可夫鏈蒙特卡洛方法 (MCMC)。像往常一樣，我將提供直觀的解釋、理論和一些帶有代碼的示例。
背景
在我們進(jìn)入主題之前，讓我們將馬爾可夫鏈蒙特卡羅（MCMC）這個術(shù)語分解為它的基本組成部分：蒙特卡羅方法和馬爾可夫鏈。
馬爾可夫鏈
?Markov Chain?是“在狀態(tài)空間上經(jīng)歷從一種狀態(tài)到另一種狀態(tài)的轉(zhuǎn)換的隨機(jī)過程”。
正如你所看到的，它看起來就像一個有限狀態(tài)機(jī)，只是我們用概率注釋了狀態(tài)轉(zhuǎn)換。例如，我們可以查看今天是否晴天，明天晴天的概率為 0.9，下雨的概率為 0.1。同樣在雨天開始。應(yīng)該清楚的是，從給定的狀態(tài)，所有傳出的轉(zhuǎn)換應(yīng)該總計(jì) 1.0，因?yàn)樗且粋€適當(dāng)?shù)姆植肌?/p>
表示此信息的另一種方法是通過轉(zhuǎn)移矩陣 P：
?
將其表示為矩陣的有趣之處在于，我們可以通過矩陣乘法來模擬馬爾可夫鏈。例如，假設(shè)我們從陽光明媚的狀態(tài)開始，我們可以將其表示為行向量：x(0)=[10]。這隱含地表示我們處于晴天狀態(tài)的概率為 1，因此處于下雨?duì)顟B(tài)的概率為 0?，F(xiàn)在，如果我們執(zhí)行矩陣乘法，我們可以在一步之后找出處于每個狀態(tài)的概率：?
我們可以看到明天有 0.9 的機(jī)會晴天（根據(jù)我們的簡單模型），有 0.1 的機(jī)會下雨。實(shí)際上，我們可以繼續(xù)將轉(zhuǎn)換矩陣相乘，以在 k?步之后找到太陽/雨的機(jī)會：?
我們可以很容易地計(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)
我們可以在這里看到一個有趣的現(xiàn)象，當(dāng)我們在狀態(tài)機(jī)中采取更多步驟時，晴天或下雨的概率似乎會收斂。您可能認(rèn)為這與我們所處的初始狀態(tài)有關(guān)，但實(shí)際上并非如此。如果我們將初始狀態(tài)初始化為隨機(jī)值，我們將得到相同的結(jié)果：
siasmdulasteds_marksov(nap.arsdray([r, 1 - r]), P, 10)x^(0) = [0.3653 0.6347]x^(1) = [0.6461 0.3539]x^(2) = [0.7584 0.2416]x^(3) = [0.8034 0.1966]x^(4) = [0.8214 0.1786]x^(5) = [0.8285 0.1715]x^(6) = [0.8314 0.1686]x^(7) = [0.8326 0.1674]x^(8) = [0.8330 0.1670]x^(9) = [0.8332 0.1668]
這種穩(wěn)態(tài)分布稱為?stationary distribution?通常用 π 表示?？梢酝ㄟ^多種方式找到該穩(wěn)態(tài)向量π。最直接的方法是在 nn 接近無窮大時取極限。
下一種方法就是求解方程。由于根據(jù)定義 q是穩(wěn)態(tài)，因此乘以 P 應(yīng)該返回相同的值：?
其中 I?是單位矩陣。如果您擴(kuò)展我們的向量/矩陣符號，您會發(fā)現(xiàn)這只是一個方程組以及 π1,π2,...,πn 總和為 1 (即 π 形成概率分布）。在我們的例子中只有兩個狀態(tài)：π1+π2=1。
然而，并不是每個馬爾可夫鏈都有一個平穩(wěn)的分布，甚至是唯一的?但是，如果我們向馬爾可夫鏈添加兩個額外的約束，我們可以保證這些屬性：
不可約：我們必須能夠最終從任何其他狀態(tài)到達(dá)任何一種狀態(tài)（即期望步數(shù)是有限的）。
非周期性：系統(tǒng)永遠(yuǎn)不會返回到具有固定周期的相同狀態(tài)（例如，不會每 5 步確定性地返回開始“晴天”）。
一個重要的定理說，如果馬爾可夫鏈?zhǔn)潜闅v的，那么它有一個唯一的穩(wěn)態(tài)概率向量 ππ。在 MCMC 的上下文中，我們可以從任何狀態(tài)跳轉(zhuǎn)到任何其他狀態(tài)（以一定的有限概率），輕松滿足不可約性。
我們將使用的另一個有用的定義是?detailed balance and reversible Markov Chains. 如果存在滿足此條件的概率分布 π，則稱馬爾可夫鏈?zhǔn)强赡娴模ㄒ卜Q為詳細(xì)平衡條件）：
換句話說，從長遠(yuǎn)來看，你從狀態(tài) i?轉(zhuǎn)移到狀態(tài) j?的次數(shù)比例，與你從狀態(tài) j 轉(zhuǎn)移到狀態(tài) i?的次數(shù)比例相同。事實(shí)上，如果馬爾可夫鏈?zhǔn)强赡娴?，那么我們就知道它具有平穩(wěn)分布（這就是我們使用相同符號 π 的原因）。
馬爾可夫鏈蒙特卡羅方法
馬爾可夫鏈蒙特卡羅 (MCMC) 方法只是一類使用馬爾可夫鏈從特定概率分布（蒙特卡羅部分）中采樣的算法。他們通過創(chuàng)建一個馬爾可夫鏈來工作，其中限制分布（或平穩(wěn)分布）只是我們想要采樣的分布。
這是一張可能有助于描述該過程的圖片. 想象一下，我們正在嘗試制作一個 MCMC 來嘗試使用 PDF f(x)對任意一維分布進(jìn)行采樣。在這種情況下，我們的狀態(tài)將是沿 x 軸的點(diǎn)，而我們的轉(zhuǎn)換概率將是從一種狀態(tài)到另一種狀態(tài)的機(jī)會。這是情況的簡化圖：
該圖顯示了我們試圖用粗黑線逼近的密度函數(shù)，以及使用從橙色狀態(tài)過渡的藍(lán)線的馬爾可夫鏈的一部分的可視化。特別是，對于 i=-3,-2,-1,1,2,3，只是從狀態(tài) X0 到 Xi 的轉(zhuǎn)換。但是，x 軸線上的每個點(diǎn)實(shí)際上都是這個馬爾可夫鏈中的一個勢態(tài)。請注意，這意味著我們有一個無限的狀態(tài)空間，因此我們不能再將轉(zhuǎn)換很好地表示為矩陣。MCMC 方法的真正“訣竅”是我們想要設(shè)計(jì)狀態(tài)（或 x 軸上的點(diǎn)）之間的轉(zhuǎn)換概率，以便我們將大部分時間花在 f(x) 很大的區(qū)域中，并且在它較小的區(qū)域中的時間相對較少（即與我們的密度函數(shù)成精確比例）。
就我們的人物而言，我們希望將大部分時間花在中心周圍，而較少時間花在外面。事實(shí)上，如果我們模擬我們的馬爾可夫鏈足夠長的時間，狀態(tài)的限制分布應(yīng)該近似于我們試圖采樣的 PDF。因此，使用 MCMC 方法進(jìn)行采樣的基本算法為：
從任意點(diǎn) x?開始。
以一定的轉(zhuǎn)移概率跳轉(zhuǎn)到點(diǎn) x'（這可能意味著保持相同的狀態(tài)）。
轉(zhuǎn)到第 2 步，直到我們轉(zhuǎn)換了 T?時間。
記錄當(dāng)前狀態(tài) x′，進(jìn)行步驟 2。
現(xiàn)在，我們在每個點(diǎn) x 軸上花費(fèi)的比例次數(shù)應(yīng)該是我們試圖模擬的 PDF 的近似值，即如果我們繪制 x 值的直方圖，我們應(yīng)該得到相同的形狀。
拒絕抽樣
現(xiàn)在，在我們進(jìn)入 MCMC 方法的具體算法之前，我想介紹另一種對概率分布進(jìn)行采樣的方法，我們稍后將使用它，稱為?rejection sampling. 主要思想是，如果我們試圖從分布 f(x) 中采樣，我們將使用另一個工具分布 g(x) 來幫助從 f(x) 中采樣。唯一的限制是對于某些 M>1，f(x)<Mg(x)。它的主要用途是當(dāng) f(x) 的形式難以直接采樣時（但仍然可以在任何點(diǎn) xx 對其進(jìn)行評估）。
以下是算法的細(xì)分：
從 g(x) 中采樣 x。
從 U(0,Mg(x)) 中采樣 y（均勻分布）。
如果 y<f(x)，則接受 x 作為 f(x) 的樣本，否則轉(zhuǎn)到步驟 1。
另一種看待它的方法是我們采樣點(diǎn) x0 的概率。這與從 g 中采樣 x0 的概率乘以我們接受的次數(shù)的比例成正比，它簡單地由 f(x0) 和 Mg(x0) 之間的比率給出：
等式告訴我們對任意點(diǎn)進(jìn)行采樣的概率與 f(x0) 成正比。在對許多點(diǎn)進(jìn)行采樣并找到我們看到 x0 的次數(shù)比例之后，常數(shù) M 被歸一化，我們得到了 PDF f(x) 的正確結(jié)果。
讓我們通過一個例子更直觀地看一下它。我們要從中采樣的目標(biāo)分布 f(x) 是?double gamma?分布，基本上是一個雙邊伽馬分布。我們將使用正態(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)，我們接受它（綠色），否則拒絕（拒絕）。
實(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
?
總的來說，我們的拒絕采樣器非常適合。與理論分布相比，抽取更多樣本會改善擬合。
拒絕抽樣的很大一部分是它很容易實(shí)現(xiàn)（在 Python 中只需幾行代碼），但有一個主要缺點(diǎn)：它很慢。
Metropolis-Hastings 算法
這?Metropolis-Hastings Algorithm?(MH) 是一種 MCMC 技術(shù)，它從難以直接采樣的概率分布中抽取樣本。與拒絕抽樣相比，對 MH 的限制實(shí)際上更加寬松：對于給定的概率密度函數(shù) p(x)，我們只要求我們有一個??與 p成正比的函數(shù) f(x)f(x) (x)p(x)！這在對后驗(yàn)分布進(jìn)行采樣時非常有用。
Metropolis-Hastings 算法的推導(dǎo)
為了推導(dǎo)出 Metropolis-Hastings 算法，我們首先從最終目標(biāo)開始：創(chuàng)建一個馬爾可夫鏈，其中穩(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ì)平衡條件，這里用另一種方式寫成：
這里 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ǔ)開始。請記住，詳細(xì)的平衡條件保證我們的馬爾可夫鏈將具有平穩(wěn)分布（它存在）。此外，如果我們也包括遍歷性（不以固定間隔重復(fù)狀態(tài)，并且每個狀態(tài)最終都能夠達(dá)到任何其他狀態(tài)），我們將建立一個具有唯一平穩(wěn)分布 p(x)的馬爾可夫鏈.
我們可以將方程重新排列為：
?
這里我們使用 f(x)來表示一個??與 p(x)成正比的函數(shù)。這是為了強(qiáng)調(diào)我們并不明確需要 p(x)，只是需要與它成比例的東西，這樣比率才能達(dá)到相同的效果?，F(xiàn)在這里的“技巧”是我們將把 P(x→x′)分解為兩個獨(dú)立的步驟：一個提議分布 g(x→x′) 和接受分布 A(x→x′)（類似于拒絕抽樣的工作原理）。由于它們是獨(dú)立的，我們的轉(zhuǎn)移概率只是兩者的乘積：?
?
此時，我們必須弄清楚 g(x)和 A(x) 的合適選擇是什么。由于 g(x) 是“建議分布”，它決定了我們可能采樣的下一個點(diǎn)。因此，重要的是它具有與我們的目標(biāo)分布 p(x)（遍歷性條件）相同的支持。這里的一個典型選擇是以當(dāng)前狀態(tài)為中心的正態(tài)分布。現(xiàn)在給定一個固定的提議分布 g(x)，我們希望找到一個匹配的 A(x)。
雖然不明顯，但滿足公式 的 A(x) 的典型選擇是：?
我們可以通過考慮 f(x′)g(x′→x)小于等于 1 和大于 1 的情況。當(dāng)小于等于 1 時，它的倒數(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ì)的平衡條件通過構(gòu)造得到滿足。所以整體算法將是（與上面的 MCMC 算法非常匹配）：
通過選擇一個隨機(jī) x 來初始化初始狀態(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 對另一個點(diǎn)進(jìn)行采樣。
預(yù)燒和相關(guān)樣本
在我們繼續(xù)實(shí)現(xiàn)之前，我們需要討論關(guān)于 MCMC 方法的兩個非常重要的話題。第一個主題與我們選擇的初始狀態(tài)有關(guān)。由于我們隨機(jī)選擇 xx 的值，它很可能位于 p(x)?非常小的區(qū)域（想想我們的雙伽馬分布的尾部）。如果從這里開始，它可能會花費(fèi)不成比例的時間來遍歷低密度的 x 值，從而錯誤地給我們一種感覺，即這些 x?值應(yīng)該比它們更頻繁地出現(xiàn)。所以解決這個問題的方法是“預(yù)燒”采樣器通過生成一堆樣本并將它們?nèi)拥?。樣本的?shù)量將取決于我們試圖模擬的分布的細(xì)節(jié)。
我們上面提到的第二個問題是兩個相鄰樣本之間的相關(guān)性。由于根據(jù)我們的轉(zhuǎn)換函數(shù) P(x→x′)的定義，繪制 x′ 取決于當(dāng)前狀態(tài) x。因此，我們失去了樣本的一項(xiàng)重要屬性：獨(dú)立性。為了糾正這一點(diǎn)，我們抽取 Tth 個樣本，并且只記錄最后抽取的樣本。假設(shè) T 足夠大，樣本應(yīng)該是相對獨(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)分布是對稱的，因此正態(tài)分布的 PDF 在其各自點(diǎn)進(jìn)行評估時會相互抵消。現(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
?
來自我們的 MH 采樣器的樣本很好地近似于我們的雙伽馬分布。此外，查看自相關(guān)圖，我們可以看到它在整個樣本中非常小，表明它們是相對獨(dú)立的。如果我們沒有為 T?選擇一個好的值或沒有預(yù)燒期，我們可能會在圖中看到較大的值。
結(jié)論
我希望你喜歡這篇關(guān)于使用拒絕抽樣和使用 Metropolis-Hastings 算法進(jìn)行 MCMC 抽樣的簡短文章。
最受歡迎的見解
1.使用R語言進(jìn)行METROPLIS-IN-GIBBS采樣和MCMC運(yùn)行
2.R語言中的Stan概率編程MCMC采樣的貝葉斯模型
3.R語言實(shí)現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣
4.R語言BUGS JAGS貝葉斯分析 馬爾科夫鏈蒙特卡洛方法（MCMC）采樣
5.R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
6.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析
7.R語言用Rcpp加速M(fèi)etropolis-Hastings抽樣估計(jì)貝葉斯邏輯回歸模型的參數(shù)
8.R語言使用Metropolis- Hasting抽樣算法進(jìn)行邏輯回歸
9.R語言中基于混合數(shù)據(jù)抽樣(MIDAS)回歸的HAR-RV模型預(yù)測GDP增長
