機器學習主題模型之LDA引數求解——Gibbs取樣

機器學習主題模型之LDA引數求解——Gibbs取樣

 

LDA引數推導的Gibbs取樣方法基於馬爾科夫鏈蒙特卡洛方法,因此首先學習MCMC方法。

一、馬爾科夫鏈蒙特卡洛方法

MCMC(Markov Chain Monte Carlo)方法是構造適合的馬爾科夫鏈,使其平穩分佈為待估引數的後驗分佈,抽樣並使用蒙特卡洛方法進行積分計算,實現了抽樣分佈隨模擬的進行而改變的動態模擬,彌補了傳統蒙特卡洛積分只能靜態模擬的缺陷。

1、蒙特卡洛方法

蒙特卡洛方法是以概率統計的理論方法為基礎的一種數值計算方法,將所要求解的問題同一定的概率模型相聯絡,用計算機實現統計模擬或抽樣,以獲得問題的近似解,故又稱隨機抽樣法或統計試驗法。這種方法求得的是近似解,對於一些複雜問題往往是唯一可行的方法,模擬的樣本數越大越接近真實值,但同時計算量也大幅上升。

蒙特卡洛方法的三個步驟:

(1)構造/描述概率過程

對於本身就具有隨機性質的問題,如粒子運輸問題,只要正確描述和模擬這個概率過程即可;對於本來不是隨機性質的確定性問題,如定積分計算,就需要構造一個人為的概率過程,將不具有隨機性質的問題轉化為隨機性質的問題。

例如要求f(x)在[a,b]上的積分,但原函式F(x)難以求解;如果可以得到x在[a,b]上的概率分佈函式p(x),則定積分可以近似求解為:

        \int_{a}^{b}f(x)dx=\int_{a}^{b}\frac{f(x)}{p(x)}p(x)dx\approx \frac{1}{n}\sum_{i=0}^{n-1}\frac{f(x_{i})}{p(x_{i})}

(2)實現從已知概率分佈中抽樣

構造概率模型後,就可以基於概率分佈去取樣所需的n個x的樣本集,帶入求和式子即可求解。最基本的一個概率分佈是(0,1)上的均勻分佈,可以使用線性同餘發生器LCG生成偽隨機數樣本:

        x:=(a*x c)%m,  0<a<m, 0≤c<m, m>0

對於一些常見的連續分佈,如正態分佈、t分佈、F分佈、Beta分佈、\Gamma分佈等都可以通過uniform(0,1)的樣本轉換而得。

對於不是常見的分佈,可以使用接受-拒絕取樣來得到該分佈的樣本,即使用一個可以取樣的分佈q(x)與一個常數c相乘後,使得真實分佈總在cq(x)曲線的下方。

       

重複以下步驟直到獲得n個被接受的取樣點:

        ①從cq(x)中獲得一個隨機取樣點x_{0},值為cq(x_{0})

        ②計算接受概率acceptance probability:\alpha =\frac{p(x_{0})}{cq(x_{0})}

        ③從uniform(0,1)中隨機生成一個值u,如果ɑ≥u,則接受取樣x_{0},否則拒絕x_{0}並回到第一步。

從圖上可以看出,cq(x)與p(x)曲線接近的地方接受概率較高,所以這樣的地方採到的點就會比較多;而在接受概率較低的、差距較大的地方採到的點就會比較少,這也就保證了這個方法的有效性。

(3)建立各種估計量。

在構造了概率模型並從中抽樣後,就可以確定一個隨機變數,作為所要求問題的解——無偏估計。建立各種估計量,相當於對模擬實驗的結果進行考察和登記,從中得到問題的解。

 

2、馬爾科夫鏈

當後驗分佈的樣本相互獨立時,可以使用大數定律求解;當樣本不獨立時,就要使用馬爾科夫鏈的平穩分佈進行抽樣。

對於一個非週期的馬爾科夫鏈,任何兩個狀態i,j∈E是相通的(不可約),那麼它的狀態轉移概率的極限存在且與i無關,則稱此馬爾科夫鏈具有遍歷性:\lim_{n\rightarrow \infty }p^{(n)}(i,j)=p_{j},此時若滿足p_{j}\geqslant 0,\; \sum p_{j}=1,則稱p_{j} 為轉移概率的極限分佈。

平穩分佈:若有限或無限數列\left \{ q_{j},j=1,2,... \right \} 滿足q_{j}\geq 0,\; \sum q_{j}=1,則稱它是概率分佈;如果此概率分佈滿足q_{j}=\sum_{i}^{ }q_{i}p(i,j),則稱它是平穩分佈。

有限馬爾科夫鏈轉移概率的極限分佈一定是平穩分佈;無限馬爾科夫鏈轉移概率的極限分佈不一定是平穩分佈。

在得到平穩分佈的馬爾科夫鏈狀態轉移矩陣後,就可以從這個平穩分佈中取樣蒙特卡洛方法所需的樣本集了。

 

       

 

3、Metropolis-Hastings演算法

對不滿足細緻平穩條件的馬爾科夫鏈狀態轉移矩陣p,引入ɑ(i,j)使其可以取等號,即得到滿足條件的目標矩陣P:

        q_{i}p(i,j)\neq q_{j}p(j,i)\\\\q_{i}p(i,j)\alpha (i,j)= q_{j}p(j,i)\alpha (j,i)\\\\\Rightarrow \left\{\begin{matrix} \alpha (i,j)=q_{j}p(j,i)\\ \alpha (j,i)=q_{i}p(i,j) \end{matrix}\right.\\\\\therefore P(i,j)=\alpha (i,j)\cdot p(i,j)

ɑ(i,j)∈[0,1]稱為接受率,P可以通過任意一個馬爾科夫鏈狀態轉移矩陣p乘以ɑ(i,j)得到,注意每次轉移後需要歸一化矩陣使得每一行之和為1。但由於ɑ(i,j)可能非常小,導致大部分取樣都被拒絕轉移,馬爾科夫鏈需要取樣很多次才能收斂,取樣效率很低。

MH取樣將接受率做了如下改進,解決了接受率過低的問題:

        \alpha (i,j)=\min \left \{ \frac{q_{j}p(j,i)}{q_{i}p(i,j)},1 \right \}

MH演算法在特徵很多的情況下,接受率的計算量很大;且在特徵維度很高,難以求出特徵的聯合概率分佈,只能得到特徵之間的條件概率分佈時,需要採取Gibbs取樣的方法解決這類問題。

       

 

4、Gibbs Sampling

馬爾科夫鏈通過轉移概率矩陣可以收斂到穩定的概率分佈,這意味著MCMC可以藉助馬爾科夫鏈的平穩分佈特性模擬高維概率分佈p(x);當馬爾科夫鏈經過burn-in階段,消除初始引數的影響、達到平穩狀態後,每一次狀態轉移都可以生成待模擬分佈的一個樣本。

Gibbs取樣是MCMC演算法的一個特例,只對2維以上的情況有效,它交替隨機固定某一維度,通過其他維度的值來抽樣該維度的值。設p(x,y)是一個二維的聯合概率分佈,由下圖可知:

       

        p(x_{1},y_{1})p(y_{2}|x_{1})=p(x_{1})p(y_{1}|x_{1})p(y_{2}|x_{1})=p(x_{1},y_{2})p(y_{1}|x_{1})

        p(x_{1},y_{1})p(x_{2}|y_{1})=p(y_{1})p(x_{1}|y_{1})p(x_{2}|y_{1})=p(x_{2},y_{1})p(x_{1}|y_{1})

        \Rightarrow \left\{\begin{matrix} p(A)p(y_{2}|x_{1})=p(B)p(y_{1}|x_{1})\\ p(A)p(x_{2}|y_{1})=p(C)p(x_{1}|y_{1}) \end{matrix}\right.   

所以在 x=x_{1} 這條直線上,如果用條件概率分佈 p(y|x_{1}) 作為馬爾科夫鏈的狀態轉移概率,則任意兩個點之間的轉移滿足細緻平穩條件;在y=y_{1} 這條直線上也是如此。Gibbs取樣和MH的不同之處在於,MH只是單獨考慮某兩個狀態的取樣(相當於直接從A到D),而Gibbs取樣是沿著不同維度進行輪換取樣的,且接受率ɑ=1 是MH演算法的特例情況。最終的近似平穩分佈可以通過取樣結果\left \{(x_{n},y_{n}),(x_{n 1},y_{n 1}),... \right \} 近似表達。

              

 

二、Gibbs取樣演算法求解LDA思路

在Gibbs取樣演算法求解LDA的方法中,假設α、η是已知的先驗輸入,目的是推斷隱變數的分佈,即估計p(θ,z,β,w|ɑ,η),可以通過對θ、β積分將它們邊緣化去掉,這種方法稱為collapse Gibbs sampling。

       

        p(\beta ,z,\theta ,w|\alpha ,\eta )=\prod_{k=1}^{K}p(\beta _{k}|\eta )\prod_{m=1}^{M}p(\theta _{m}|\alpha )\prod_{n=1}^{N}p(z_{mn}|\theta _{m})p(w_{mn}|z_{mn},\beta _{k})

        p(z,w|\alpha ,\eta )\\\\=\int \int p(\theta ,z,w,\beta |\alpha ,\eta )d\theta d\beta \\\\=\int \prod_{m=1}^{M}p(\theta _{m}|\alpha )\prod_{n=1}^{N}p(z_{mn}|\theta_{m})d\mathbf{\theta}\cdot \int \prod_{k=1}^{K}p(\beta _{k}|\eta)\prod_{m=1}^{M}\prod_{n=1}^{N}p(w_{mn}|z_{mn},\beta_{k})d\beta

假設每篇文件長度N相同,分別考慮θ、β這兩部分的積分,令n_{n}^{k} 代表第m篇文件中分配給第k個主題的詞的個數,對應的多項分佈的計數為:n_{m}=\left ( n_{m}^{1}, n_{m}^{2},..., n_{m}^{K} \right )

        \sum_{n=1}^{N}p(z_{mn}|\theta _{m})=\prod_{k=1}^{K}\theta _{mk}^{n_{m}^{k}}

先對θ部分進行計算,根據α→θ→z的Dir-Mult共軛分佈關係,省略M的連乘可得:

        \int p(\theta _{m}|\alpha )\prod_{n=1}^{N}p(z_{mn}|\theta _{m})d\theta _{m}\\\\=\int \frac{\Gamma (\sum_{k=1}^{K}\alpha _{k})}{\prod_{k=1}^{K}\Gamma (\alpha _{k})}\prod_{k=1}^{K}\theta _{mk}^{\alpha _{k}-1}\prod_{k=1}^{K}\theta _{mk}^{n_{m}^{k}}d\theta _{m}\\\\=\int \frac{\Gamma (\sum_{k=1}^{K}\alpha _{k})}{\prod_{k=1}^{K}\Gamma (\alpha _{k})}\prod_{k=1}^{K}\theta _{mk}^{n_{m}^{k} \alpha _{k}-1}d\theta _{m}

根據Dirichlet分佈的積分性質可知:

        \int \frac{\Gamma (\sum_{k=1}^{K}n_{m}^{k} \alpha _{k})}{\prod_{k=1}^{K}\Gamma (n_{m}^{k} \alpha _{k})}\prod_{k=1}^{K}\theta _{mk}^{n_{m}^{k} \alpha _{k}-1}d\theta _{m}=1

因此θ部分積分可化為:

        \int p(\theta _{m}|\alpha )\prod_{n=1}^{N}p(z_{mn}|\theta _{m})d\theta _{m}= \frac{\Gamma (\sum_{k=1}^{K}\alpha _{k})}{\prod_{k=1}^{K}\Gamma (\alpha _{k})}\cdot \frac{\prod_{k=1}^{K}\Gamma (n_{m}^{k} \alpha _{k})}{\Gamma (\sum_{k=1}^{K}n_{m}^{k} \alpha _{k})}

 

同樣地,令n_{k}^{v} 代表第k個主題對應詞典中第v個詞的詞個數,對應的多項分佈的計數為:n_{k}=\left ( n_{k}^{1}, n_{k}^{2},..., n_{k}^{V} \right ),可將β部分的積分寫為:

        \int \prod_{k=1}^{K}p(\beta _{k}|\eta)\prod_{m=1}^{M}\prod_{n=1}^{N}p(w_{mn}|z_{mn},\beta_{k})d\beta \\\\=\prod_{k=1}^{K}\int p(\beta _{k}|\eta)\prod_{m=1}^{M}\prod_{n=1}^{N}p(w_{mn}|z_{mn},\beta_{k})d\beta_{k}\\\\=\prod_{k=1}^{K}\int \frac{\Gamma (\sum_{v=1}^{V}\eta_{v})}{\prod_{v=1}^{V}\Gamma (\eta_{v})}\prod_{v=1}^{V}\beta _{mk}^{n_{k}^{v} \eta _{v}-1}d\beta _{k}\\\\=\prod_{k=1}^{K}\frac{\Gamma (\sum_{v=1}^{V}\eta_{v})}{\prod_{v=1}^{V}\Gamma (\eta_{v})}\cdot \frac{\prod_{v=1}^{V}\Gamma (n_{k}^{v} \eta _{v})}{\Gamma (\sum_{v=1}^{V}n_{k}^{v} \eta _{v})}

最終得到主題和詞的聯合分佈如下:

        p(z,w|\alpha ,\eta )=\prod_{m=1}^{M}\frac{\Gamma (\sum_{k=1}^{K}\alpha _{k})}{\prod_{k=1}^{K}\Gamma (\alpha _{k})}\cdot \frac{\prod_{k=1}^{K}\Gamma (n_{m}^{k} \alpha _{k})}{\Gamma (\sum_{k=1}^{K}n_{m}^{k} \alpha _{k})} \times \prod_{k=1}^{K}\frac{\Gamma (\sum_{v=1}^{V}\eta_{v})}{\prod_{v=1}^{V}\Gamma (\eta_{v})}\cdot \frac{\prod_{v=1}^{V}\Gamma (n_{k}^{v} \eta _{v})}{\Gamma (\sum_{v=1}^{V}n_{k}^{v} \eta _{v})}

 

Gibbs取樣的目標是接近分佈p(z|w,ɑ,η),由於w,ɑ,η已知,令z_{mn} 表示第m篇文件的第n個詞對應的主題,z_{-mn} 代表z去掉z_{mn} 後的主題分佈。

        p(z_{mn}|z_{-mn},w,\alpha ,\eta )=\frac{p(z_{mn},z_{-mn},w|\alpha ,\eta)}{p(z_{-mn},w|\alpha ,\eta )}

注意到Gibbs取樣只需要為z_{mn} 採一個樣,因此不用得到上式準確的值,只需得到z_{mn} 取值的比值即可,將上式簡化為: 

        p(z_{mn}=k|z_{-mn},w,\alpha ,\eta )\\\\\propto p(z_{mn}=k,z_{-mn},w|\alpha ,\eta )\\\\=\left ( \frac{\Gamma (\sum_{i=1}^{K}\alpha _{i})}{\prod_{i=1}^{K}\Gamma (\alpha _{i})} \right )^{M}\prod_{j\neq m}^{ }\frac{\prod_{i=1}^{K}\Gamma (n_{j}^{i} \alpha _{i})}{\Gamma (\sum_{i=1}^{K}n_{j}^{i} \alpha _{i})}\left ( \frac{\Gamma (\sum_{r=1}^{V}\eta _{r})}{\prod_{r=1}^{V}\Gamma (\eta _{r})} \right )^{K}\\\\\cdot \prod_{i=1}^{K}\prod_{r\neq v}^{ }\Gamma (n_{i}^{r} \eta _{r})\frac{\prod_{i=1}^{K}\Gamma (n_{m}^{i} \alpha _{i})}{\Gamma (\sum_{i=1}^{K}n_{m}^{i} \alpha _{i})}\prod_{i=1}^{K}\frac{\Gamma (n_{i}^{v} \eta _{v})}{\Gamma (\sum_{r=1}^{V}n_{i}^{r} \eta _{r})}

        \propto \prod_{i=1}^{K}\Gamma (n_{m}^{i} \alpha _{i})\frac{\Gamma (n_{i}^{v} \eta _{v})}{\Gamma (\sum_{r=1}^{V}n_{i}^{r} \eta _{r})}

        \propto \prod_{i\neq k}^{ }\Gamma (n_{m}^{i,-mn} \alpha _{i})\frac{\Gamma (n_{i,-mn}^{v} \eta _{v})}{\Gamma (\sum_{r=1}^{V}n_{i,-mn}^{r} \eta _{r})}\cdot \Gamma (n_{m}^{k,-mn} \alpha _{k} 1)\frac{\Gamma (n_{k,-mn}^{v} \eta _{v} 1)}{\Gamma (\sum_{r=1}^{V}n_{k,-mn}^{r} \eta _{r} 1)}

        = \prod_{i}^{ }\Gamma (n_{m}^{i,-mn} \alpha _{i})\frac{\Gamma (n_{i,-mn}^{v} \eta _{v})}{\Gamma (\sum_{r=1}^{V}n_{i,-mn}^{r} \eta _{r})}\cdot (n_{m}^{k,-mn} \alpha _{k})\frac{ (n_{k,-mn}^{v} \eta _{v})}{(\sum_{r=1}^{V}n_{k,-mn}^{r} \eta _{r})}

        \propto (n_{m}^{k,-mn} \alpha _{k})\frac{ (n_{k,-mn}^{v} \eta _{v})}{(\sum_{r=1}^{V}n_{k,-mn}^{r} \eta _{r})}

Gibbs取樣公式的物理意義就是在K條doc→topic→word的路徑中進行取樣,當取樣收斂後可得所有詞最終的主題z,統計每篇文件中的主題概率p(topic|doc)可得文件主題分佈\theta _{m},統計每個主題中所有的詞概率p(word|topic)可得主題詞分佈\beta _{k}

       

 

三、LDA Gibbs取樣演算法流程總結

(1)選擇適當的主題總數K,及LDA模型引數ɑ、η;

(2)對應文集中的每一篇文件的每個詞w,隨機賦予一個主題編號z;

(3)重新掃描文集,對每個詞w使用Gibbs取樣公式更新它的主題編號,並更新文集中的詞;

(4)重複取樣過程直至收斂;

(5)統計文集中topic-word共現頻率矩陣,即可得到LDA模型。

通常訓練LDA模型時,取Gibbs取樣收斂之後的n個迭代的結果進行平均來做引數估計,這樣模型質量更高。

 

當LDA模型訓練好後,要確定一篇新文件的主題時,由於各個主題的詞分佈\beta _{k} 已經確定,只需要對\theta _{m} 進行取樣計算即可。

選擇適合的超引數K非常關鍵,對於簡單的語義區分可選擇較小的K,對於複雜的語義區分則選擇較大的K,且需要足夠多的語料。

Gibbs取樣容易並行化,可以很方便的使用大資料平臺來分散式地訓練海量文件的LDA模型。

 

 

參考資料

https://en.wikipedia.org/wiki/Latent_Dirichlet_allocation

lda數學八卦

https://blog.csdn.net/pipisorry/article/details/51373090

http://www.cnblogs.com/pinard/p/6831308.html