貝葉斯統計：Gibbs sampling 原理到實踐

多項分佈公式怎麼來的

經過前面兩篇文章，我們終於來到了 Gibbs samping，為什麼我這麼興奮呢？因為我當初看貝葉斯推斷就是因為 LDA 模型，裡面用到了 Gibbs sampling 的方法來解引數問題。

由於頭條對 Markdown 支援不是太好，可以去原文連結檢視：https：//www。zybuluo。com/zhuanxu/note/1027270。

好了下面開始介紹 Gibbs 取樣。

Gibbs 取樣

前面一篇文章我介紹了細緻平穩條件，即：

Gibbs 取樣的演算法如下：

我們來證明 Gibbs 取樣演算法也滿足細緻平穩條件。

假設 x = x 1 ， 。 。 。 ， x D ，當我們取樣第 k 個數據的時候，

此時我們的接受率為：

上面公式一個關鍵的部分是：

帶入就可以得到 1，即 gibbs 取樣是一定接受的取樣。

下面我們照慣例還是來一個例子。

例子

假設有我們有資料 x1，x2。。。xN，其中 1-n 的數服從一個泊松分佈，n+1-N 的數服從另一個泊松分佈。

而泊松分佈的引數服從 Gamma 分佈，我們總結下目前的先驗假設：

此時後驗機率是：

我們下一步開始我們的取樣， 先計算下 logP：

接著來取樣

我們此處只取了跟當前取樣引數有關的項，因為其他都是一些常數，作用只是將機率分佈歸一化，不影響取樣。

此處都服從 Gamma 分佈。

最後是 n：

n 沒有顯示的分佈資訊，但是我們簡單將其認為是一個多項分佈。

下面是關鍵 Python 程式碼：

完整程式碼見 Gibbs 。

取樣後輸出如下圖：

LDA

下面我們來最激動人心的 LDA 模型，看怎麼用 Gibbs 取樣來解。

先看 LDA 的模型：

整個過程可以描述為

對於文件 d，從分佈中取樣出，即文件 d 的主題分佈

對於文件的每個詞從多項分佈中取樣出主題 k，即

對於主題 k，從中取樣出，即主題 k 的詞分佈

對於每個詞，從主題分佈中取樣出具體的單詞 t

上面整個模型中，模型引數有：

為了做 Gibbs 取樣，我們先看這幾個引數的聯合分佈。

上面公式中是多項分佈，而是 Dir 分佈，我們知道 Dir 分佈和多項分佈是共軛分佈，因此後驗分佈比較容易寫出來。

下面我們開始計算每個引數的機率，主要是計算下面 3 個機率：

先來計算第一個主題的機率分佈：

其中

表示第 m 篇文件中，屬於每個主題的詞個數。

接著計算每個主題的詞分佈

其中

第 m 篇文件中屬於 k 個主題的，每個單詞的數量，共 （1-V） 個單詞。最後來估算每個詞的主題。

我們可以看到 p（zmj=k） 是一個多項分佈，每一項的機率都是 beta*theta，而他們本身也是需要從 Dir 分佈中取樣出來的，一個自然的想法就是我們可以用估計值來代替，根據 Dir 分佈我們能夠很方便的計算出機率來。

此處需要的數學基礎可以看頭條號中的主題模型：LDA 數學基礎。

裡面一個點是 Dir 分佈和其數學期望

我們上面在 Gibbs 取樣中計算分別取樣都可以用 Dir 分佈的的期望來作為新的引數值。

編碼

介紹完數學基礎後，我麼就能來看如何實現了，下面是一個虛擬碼。

完整的程式碼可以見玩點高階的 —— 帶你入門 Topic 模型 LDA（小改進 + 附原始碼）

pymc3 實現

下面我們再來用 pymc3 來實現下。

可以說 pymc3 寫出來的程式碼真是簡潔。but。就是太慢了，完整的程式碼可以看 gibbs-lda。

總結

本文介紹了 mh 演算法的特例 Gibbs 取樣，並且給出了證明為什麼 Gibbs 取樣 work，最後我們用 Gibbs 取樣來解決了 LDA 的引數估計問題。

你的鼓勵是我繼續寫下去的動力，期待我們共同進步。