狄利克雷分布公式_深入机器学习系列11－隐式狄利克雷分布
轉(zhuǎn)載請(qǐng)注明出處，該文章的官方來(lái)源：
LDA | Teaching ML
前言
LDA是一種概率主題模型：隱式狄利克雷分布（Latent Dirichlet Allocation，簡(jiǎn)稱(chēng)LDA）。LDA是2003年提出的一種主題模型，它可以將文檔集中每篇文檔的主題以概率分布的形式給出。 通過(guò)分析一些文檔，我們可以抽取出它們的主題（分布），根據(jù)主題（分布）進(jìn)行主題聚類(lèi)或文本分類(lèi)。同時(shí)，它是一種典型的詞袋模型，即一篇文檔是由一組詞構(gòu)成，詞與詞之間沒(méi)有先后順序的關(guān)系。一篇文檔可以包含多個(gè)主題，文檔中每一個(gè)詞都由其中的一個(gè)主題生成。
舉一個(gè)簡(jiǎn)單的例子，比如假設(shè)事先給定了這幾個(gè)主題：Arts、Budgets、Children、Education，然后通過(guò)學(xué)習(xí)的方式，獲取每個(gè)主題Topic對(duì)應(yīng)的詞語(yǔ)，如下圖所示：
然后以一定的概率選取上述某個(gè)主題，再以一定的概率選取那個(gè)主題下的某個(gè)單詞，不斷的重復(fù)這兩步，最終生成如下圖所示的一篇文章（不同顏色的詞語(yǔ)分別表示不同主題）。
我們看到一篇文章后，往往會(huì)推測(cè)這篇文章是如何生成的，我們通常認(rèn)為作者會(huì)先確定幾個(gè)主題，然后圍繞這幾個(gè)主題遣詞造句寫(xiě)成全文。LDA要干的事情就是根據(jù)給定的文檔，判斷它的主題分布。在LDA模型中，生成文檔的過(guò)程有如下幾步：
- 從狄利克雷分布 中生成文檔i的主題分布 ；
- 從主題的多項(xiàng)式分布 中取樣生成文檔i第j個(gè)詞的主題 ；
- 從狄利克雷分布 中取樣生成主題 對(duì)應(yīng)的詞語(yǔ)分布 ；
- 從詞語(yǔ)的多項(xiàng)式分布 中采樣最終生成詞語(yǔ) .
LDA的圖模型結(jié)構(gòu)如下圖所示：
LDA會(huì)涉及很多數(shù)學(xué)知識(shí)，后面的章節(jié)我會(huì)首先介紹LDA涉及的數(shù)學(xué)知識(shí)，然后在這些數(shù)學(xué)知識(shí)的基礎(chǔ)上詳細(xì)講解LDA的原理。
1 數(shù)學(xué)預(yù)備
1.1 Gamma函數(shù)
在高等數(shù)學(xué)中，有一個(gè)長(zhǎng)相奇特的Gamma函數(shù)
通過(guò)分部積分，可以推導(dǎo)gamma函數(shù)有如下遞歸性質(zhì)
通過(guò)該遞歸性質(zhì)，我們可以很容易證明，gamma函數(shù)可以被當(dāng)成階乘在實(shí)數(shù)集上的延拓，具有如下性質(zhì)
1.2 Digamma函數(shù)
如下函數(shù)被稱(chēng)為Digamma函數(shù)，它是Gamma函數(shù)對(duì)數(shù)的一階導(dǎo)數(shù)
這是一個(gè)很重要的函數(shù)，在涉及Dirichlet分布相關(guān)的參數(shù)的極大似然估計(jì)時(shí)，往往需要用到這個(gè)函數(shù)。Digamma函數(shù)具有如下一個(gè)漂亮的性質(zhì)
1.3 二項(xiàng)分布（Binomial distribution）
二項(xiàng)分布是由伯努利分布推出的。伯努利分布，又稱(chēng)兩點(diǎn)分布或0-1分布，是一個(gè)離散型的隨機(jī)分布，其中的隨機(jī)變量只有兩類(lèi)取值，即0或者1。二項(xiàng)分布是重復(fù)n次的伯努利試驗(yàn)。簡(jiǎn)言之，只做一次實(shí)驗(yàn)，是伯努利分布，重復(fù)做了n次，是二項(xiàng)分布。二項(xiàng)分布的概率密度函數(shù)為：
對(duì)于k=1,2，...,n，其中C(n,k)是二項(xiàng)式系數(shù)（這就是二項(xiàng)分布的名稱(chēng)的由來(lái)）
1.4 多項(xiàng)分布
多項(xiàng)分布是二項(xiàng)分布擴(kuò)展到多維的情況。多項(xiàng)分布是指單次試驗(yàn)中的隨機(jī)變量的取值不再是0-1，而是有多種離散值可能（1,2,3...,k）。比如投擲6個(gè)面的骰子實(shí)驗(yàn)，N次實(shí)驗(yàn)結(jié)果服從K=6的多項(xiàng)分布。其中：
多項(xiàng)分布的概率密度函數(shù)為：
1.5 Beta分布
1.5.1 Beta分布
首先看下面的問(wèn)題1（問(wèn)題1到問(wèn)題4都取自于文獻(xiàn)【1】）。
問(wèn)題1：
為解決這個(gè)問(wèn)題，可以嘗試計(jì)算
落在區(qū)間[x,x+delta x]的概率。首先，把[0,1]區(qū)間分成三段[0,x),[x,x+delta x]，(x+delta x,1]，然后考慮下簡(jiǎn)單的情形：即假設(shè)n個(gè)數(shù)中只有1個(gè)落在了區(qū)間[x,x+delta x]內(nèi)，由于這個(gè)區(qū)間內(nèi)的數(shù)X(k)是第k大的，所以[0,x)中應(yīng)該有k?1個(gè)數(shù)，(x+delta x,1]這個(gè)區(qū)間中應(yīng)該有n?k個(gè)數(shù)。 如下圖所示：上述問(wèn)題可以轉(zhuǎn)換為下述事件E：
對(duì)于上述事件E，有：
其中，o(delta x)表示delta x的高階無(wú)窮小。顯然，由于不同的排列組合，即n個(gè)數(shù)中有一個(gè)落在[x,x+delta x]區(qū)間的有n種取法，余下n?1個(gè)數(shù)中有k?1個(gè)落在[0,x)的有C(n-1,k-1)種組合。所以和事件E等價(jià)的事件一共有nC(n-1,k-1)個(gè)。
文獻(xiàn)【1】中證明，只要落在[x,x+delta x]內(nèi)的數(shù)字超過(guò)一個(gè)，則對(duì)應(yīng)的事件的概率就是o(delta x)。所以
的概率密度函數(shù)為：利用Gamma函數(shù)，我們可以將f(x)表示成如下形式：
在上式中，我們用alpha=k，beta=n-k+1替換，可以得到beta分布的概率密度函數(shù)
1.5.2 共軛先驗(yàn)分布
什么是共軛呢？軛的意思是束縛、控制。共軛從字面上理解，則是共同約束，或互相約束。在貝葉斯概率理論中，如果后驗(yàn)概率P(z|x)和先驗(yàn)概率p(z)滿(mǎn)足同樣的分布，那么，先驗(yàn)分布和后驗(yàn)分布被叫做共軛分布，同時(shí)，先驗(yàn)分布叫做似然函數(shù)的共軛先驗(yàn)分布。
1.5.3 Beta-Binomial 共軛
我們?cè)趩?wèn)題1的基礎(chǔ)上增加一些觀測(cè)數(shù)據(jù)，變成問(wèn)題2：
第2步的條件可以用另外一句話(huà)來(lái)表述，即“Yi中有m1個(gè)比X(k)小，m2個(gè)比X(k)大”，所以X(k)是
中k+m1大的數(shù)。根據(jù)1.5.1的介紹，我們知道事件p服從beta分布,它的概率密度函數(shù)為：
按照貝葉斯推理的邏輯，把以上過(guò)程整理如下：
- 1、p是我們要猜測(cè)的參數(shù)，我們推導(dǎo)出p的分布為f(p)=Beta(p|k,n-k+1),稱(chēng)為p的先驗(yàn)分布
- 2、根據(jù)Yi中有m1個(gè)比p小，有m2個(gè)比p大，Yi相當(dāng)是做了m次伯努利實(shí)驗(yàn)，所以m1服從二項(xiàng)分布B(m,p)
- 3、在給定了來(lái)自數(shù)據(jù)提供(m1,m2)知識(shí)后，p的后驗(yàn)分布變?yōu)閒(p|m1,m2)=Beta(p|k+m1,n-k+1+m2)
貝葉斯估計(jì)的基本過(guò)程是：
先驗(yàn)分布 + 數(shù)據(jù)的知識(shí) = 后驗(yàn)分布
以上貝葉斯分析過(guò)程的簡(jiǎn)單直觀的表示就是：
Beta(p|k,n-k+1) + BinomCount(m1,m2) = Beta(p|k+m1,n-k+1+m2)
更一般的，對(duì)于非負(fù)實(shí)數(shù)alpha和beta，我們有如下關(guān)系
Beta(p|alpha,beta) + BinomCount(m1,m2) = Beta(p|alpha+m1,beta+m2)
針對(duì)于這種觀測(cè)到的數(shù)據(jù)符合二項(xiàng)分布，參數(shù)的先驗(yàn)分布和后驗(yàn)分布都是Beta分布的情況，就是Beta-Binomial共軛。換言之，Beta分布是二項(xiàng)式分布的共軛先驗(yàn)概率分布。二項(xiàng)分布和Beta分布是共軛分布意味著，如果我們?yōu)槎?xiàng)分布的參數(shù)p選取的先驗(yàn)分布是Beta分布，那么以p為參數(shù)的二項(xiàng)分布用貝葉斯估計(jì)得到的后驗(yàn)分布仍然服從Beta分布。
1.6 Dirichlet 分布
1.6.1 Dirichlet 分布
Dirichlet分布，是beta分布在高維度上的推廣。Dirichlet分布的的密度函數(shù)形式跟beta分布的密度函數(shù)類(lèi)似：
其中
至此，我們可以看到二項(xiàng)分布和多項(xiàng)分布很相似，Beta分布和Dirichlet分布很相似。并且Beta分布是二項(xiàng)式分布的共軛先驗(yàn)概率分布。那么Dirichlet分布呢？Dirichlet分布是多項(xiàng)式分布的共軛先驗(yàn)概率分布。下文來(lái)論證這點(diǎn)。
1.6.2 Dirichlet-Multinomial 共軛
在1.5.3章問(wèn)題2的基礎(chǔ)上，我們更進(jìn)一步引入問(wèn)題3：
類(lèi)似于問(wèn)題1的推導(dǎo)，我們可以容易推導(dǎo)聯(lián)合分布。為了簡(jiǎn)化計(jì)算，我們?nèi)3滿(mǎn)足x1+x2+x3=1,x1和x2是變量。如下圖所示。
概率計(jì)算如下：
于是我們得到聯(lián)合分布為：
觀察上述式子的最終結(jié)果，可以看出上面這個(gè)分布其實(shí)就是3維形式的Dirichlet分布。令alpha1=k1,alpha2=k2,alpha3=n-k1-k2+1，分布密度函數(shù)可以寫(xiě)為：
為了論證Dirichlet分布是多項(xiàng)式分布的共軛先驗(yàn)概率分布，在上述問(wèn)題3的基礎(chǔ)上再進(jìn)一步，提出問(wèn)題4。
為了方便計(jì)算，我們記
根據(jù)問(wèn)題中的信息，我們可以推理得到p1,p2在X;Y這m+n個(gè)數(shù)中分別成為了第k1+m1,k1+k2+m1+m2大的數(shù)。后驗(yàn)分布p應(yīng)該為
同樣的，按照貝葉斯推理的邏輯，可將上述過(guò)程整理如下：
- 1 我們要猜測(cè)參數(shù)P=(p1,p2,p3)，其先驗(yàn)分布為Dir(p|k);
- 2 數(shù)據(jù)Yi落到三個(gè)區(qū)間[0,p1),[p1,p2],(p2,1]的個(gè)數(shù)分別是m1,m2,m3,所以m=(m1,m2,m3)服從多項(xiàng)分布Mult(m|p);
- 3 在給定了來(lái)自數(shù)據(jù)提供的知識(shí)m后，p的后驗(yàn)分布變?yōu)镈ir(P|k+m)
上述貝葉斯分析過(guò)程的直觀表述為：
Dir(p|k) + Multcount(m) = Dir(p|k+m)
針對(duì)于這種觀測(cè)到的數(shù)據(jù)符合多項(xiàng)分布，參數(shù)的先驗(yàn)分布和后驗(yàn)分布都是Dirichlet分布的情況，就是Dirichlet-Multinomial共軛。這意味著，如果我們?yōu)槎囗?xiàng)分布的參數(shù)p選取的先驗(yàn)分布是Dirichlet分布，那么以p為參數(shù)的多項(xiàng)分布用貝葉斯估計(jì)得到的后驗(yàn)分布仍然服從Dirichlet分布。
1.7 Beta和Dirichlet分布的一個(gè)性質(zhì)
如果p=Beta(t|alpha,beta)，那么
上式右邊的積分對(duì)應(yīng)到概率分布Beta(t|alpha+1,beta)，對(duì)于這個(gè)分布，我們有
把上式帶人E(p)的計(jì)算式，可以得到：
這說(shuō)明，對(duì)于Beta分布的隨機(jī)變量，其期望可以用上式來(lái)估計(jì)。Dirichlet分布也有類(lèi)似的結(jié)論。對(duì)于p=Dir(t|alpha)，有
這個(gè)結(jié)論在后文的推導(dǎo)中會(huì)用到。
1.8 總結(jié)
LDA涉及的數(shù)學(xué)知識(shí)較多，需要認(rèn)真體會(huì)，以上大部分的知識(shí)來(lái)源于文獻(xiàn)【1,2,3】,如有不清楚的地方，參見(jiàn)這些文獻(xiàn)以了解更多。
2 主題模型LDA
在介紹LDA之前，我們先介紹幾個(gè)基礎(chǔ)模型：Unigram model、mixture of unigrams model、pLSA model。為了方便描述，首先定義一些變量：
- 1 w表示詞，V表示所有詞的個(gè)數(shù)
- 2 z表示主題，k表示主題的個(gè)數(shù)
- 3 表示語(yǔ)料庫(kù)，M表示語(yǔ)料庫(kù)中的文檔數(shù)。
- 4 表示文檔，N表示文檔中詞的個(gè)數(shù)。
2.1 一元模型(Unigram model)
對(duì)于文檔
,用 表示 的先驗(yàn)概率，生成文檔W的概率為：其圖模型為（圖中被涂色的w表示可觀測(cè)變量，N表示一篇文檔中總共N個(gè)單詞，M表示M篇文檔）：
2.2 混合一元模型(Mixture of unigrams model)
該模型的生成過(guò)程是：給某個(gè)文檔先選擇一個(gè)主題Z，再根據(jù)該主題生成文檔，該文檔中的所有詞都來(lái)自一個(gè)主題。生成文檔的概率為：
其圖模型為（圖中被涂色的w表示可觀測(cè)變量，未被涂色的z表示未知的隱變量，N表示一篇文檔中總共N個(gè)單詞，M表示M篇文檔）：
2.3 pLSA模型
在混合一元模型中，假定一篇文檔只由一個(gè)主題生成，可實(shí)際中，一篇文章往往有多個(gè)主題，只是這多個(gè)主題各自在文檔中出現(xiàn)的概率大小不一樣。在pLSA中，假設(shè)文檔由多個(gè)主題生成。下面通過(guò)一個(gè)投色子的游戲（取自文獻(xiàn)【2】的例子）說(shuō)明pLSA生成文檔的過(guò)程。
首先，假定你一共有K個(gè)可選的主題，有V個(gè)可選的詞。假設(shè)你每寫(xiě)一篇文檔會(huì)制作一顆K面的“文檔-主題”骰子（扔此骰子能得到K個(gè)主題中的任意一個(gè)），和K個(gè)V面的“主題-詞項(xiàng)”骰子（每個(gè)骰子對(duì)應(yīng)一個(gè)主題，K個(gè)骰子對(duì)應(yīng)之前的K個(gè)主題，且骰子的每一面對(duì)應(yīng)要選擇的詞項(xiàng)，V個(gè)面對(duì)應(yīng)著V個(gè)可選的詞）。 比如可令K=3，即制作1個(gè)含有3個(gè)主題的“文檔-主題”骰子，這3個(gè)主題可以是：教育、經(jīng)濟(jì)、交通。然后令V = 3，制作3個(gè)有著3面的“主題-詞項(xiàng)”骰子，其中，教育主題骰子的3個(gè)面上的詞可以是：大學(xué)、老師、課程，經(jīng)濟(jì)主題骰子的3個(gè)面上的詞可以是：市場(chǎng)、企業(yè)、金融，交通主題骰子的3個(gè)面上的詞可以是：高鐵、汽車(chē)、飛機(jī)。
其次，每寫(xiě)一個(gè)詞，先扔該“文檔-主題”骰子選擇主題，得到主題的結(jié)果后，使用和主題結(jié)果對(duì)應(yīng)的那顆“主題-詞項(xiàng)”骰子，扔該骰子選擇要寫(xiě)的詞。先扔“文檔-主題”的骰子，假設(shè)以一定的概率得到的主題是：教育，所以下一步便是扔教育主題篩子，以一定的概率得到教育主題篩子對(duì)應(yīng)的某個(gè)詞大學(xué)。
- 上面這個(gè)投骰子產(chǎn)生詞的過(guò)程簡(jiǎn)化一下便是：“先以一定的概率選取主題，再以一定的概率選取詞”。事實(shí)上，一開(kāi)始可供選擇的主題有3個(gè)：教育、經(jīng)濟(jì)、交通，那為何偏偏選取教育這個(gè)主題呢？其實(shí)是隨機(jī)選取的，只是這個(gè)隨機(jī)遵循一定的概率分布。比如可能選取教育主題的概率是0.5，選取經(jīng)濟(jì)主題的概率是0.3，選取交通主題的概率是0.2，那么這3個(gè)主題的概率分布便是{教育：0.5，經(jīng)濟(jì)：0.3，交通：0.2}，我們把各個(gè)主題z在文檔d中出現(xiàn)的概率分布稱(chēng)之為主題分布，且是一個(gè)多項(xiàng)分布。
- 同樣的，從主題分布中隨機(jī)抽取出教育主題后，依然面對(duì)著3個(gè)詞：大學(xué)、老師、課程，這3個(gè)詞都可能被選中，但它們被選中的概率也是不一樣的。比如大學(xué)這個(gè)詞被選中的概率是0.5，老師這個(gè)詞被選中的概率是0.3，課程被選中的概率是0.2，那么這3個(gè)詞的概率分布便是{大學(xué)：0.5，老師：0.3，課程：0.2}，我們把各個(gè)詞語(yǔ)w在主題z下出現(xiàn)的概率分布稱(chēng)之為詞分布，這個(gè)詞分布也是一個(gè)多項(xiàng)分布。
- 所以，選主題和選詞都是兩個(gè)隨機(jī)的過(guò)程，先從主題分布{教育：0.5，經(jīng)濟(jì)：0.3，交通：0.2}中抽取出主題：教育，然后從該主題對(duì)應(yīng)的詞分布{大學(xué)：0.5，老師：0.3，課程：0.2}中抽取出詞：大學(xué)。
最后，你不停的重復(fù)扔“文檔-主題”骰子和”主題-詞項(xiàng)“骰子，重復(fù)N次（產(chǎn)生N個(gè)詞），完成一篇文檔，重復(fù)這產(chǎn)生一篇文檔的方法M次，則完成M篇文檔。
上述過(guò)程抽象出來(lái)即是pLSA的文檔生成模型。在這個(gè)過(guò)程中，我們并未關(guān)注詞和詞之間的出現(xiàn)順序，所以pLSA是一種詞袋模型。定義如下變量:
- 表示隱藏的主題；
- 表示海量文檔中某篇文檔被選中的概率；
- 表示詞 在文檔 中出現(xiàn)的概率；針對(duì)海量文檔，對(duì)所有文檔進(jìn)行分詞后，得到一個(gè)詞匯列表，這樣每篇文檔就是一個(gè)詞語(yǔ)的集合。對(duì)于每個(gè)詞語(yǔ)，用它在文檔中出現(xiàn)的次數(shù)除以文檔中詞語(yǔ)總的數(shù)目便是它在文檔中出現(xiàn)的概率；
- 表示主題 在文檔 中出現(xiàn)的概率；
- 表示詞 在主題 中出現(xiàn)的概率。與主題關(guān)系越密切的詞其條件概率越大。
我們可以按照如下的步驟得到“文檔-詞項(xiàng)”的生成模型：
- 1 按照 選擇一篇文檔 ；
- 2 選定文檔 之后，從主題分布中按照概率 選擇主題；
- 3 選定主題后，從詞分布中按照概率 P(w_{j}|z_{k})選擇一個(gè)詞。
利用看到的文檔推斷其隱藏的主題（分布）的過(guò)程，就是主題建模的目的：自動(dòng)地發(fā)現(xiàn)文檔集中的主題（分布）。文檔d和單詞w是可被觀察到的，但主題z卻是隱藏的。如下圖所示（圖中被涂色的d、w表示可觀測(cè)變量，未被涂色的z表示未知的隱變量，N表示一篇文檔中總共N個(gè)單詞，M表示M篇文檔）。
上圖中，文檔d和詞w是我們得到的樣本，可觀測(cè)得到，所以對(duì)于任意一篇文檔，其
是已知的。根據(jù)這個(gè)概率可以訓(xùn)練得到文檔-主題概率以及主題-詞項(xiàng)概率。即：故得到文檔中每個(gè)詞的生成概率為：
可以直接得出，而 和 未知，所以 就是我們要估計(jì)的參數(shù),我們要最大化這個(gè)參數(shù)。因?yàn)樵摯烙?jì)的參數(shù)中含有隱變量z，所以我們可以用EM算法來(lái)估計(jì)這個(gè)參數(shù)。2.4 LDA模型
LDA的不同之處在于，pLSA的主題的概率分布P(c|d)是一個(gè)確定的概率分布，也就是雖然主題c不確定，但是c符合的概率分布是確定的，比如符合某個(gè)多項(xiàng)分布，這個(gè)多項(xiàng)分布的各參數(shù)是確定的。 但是在LDA中，這個(gè)多項(xiàng)分布都是不確定的，高斯分布又服從一個(gè)狄利克雷先驗(yàn)分布(Dirichlet prior)。即LDA就是pLSA的貝葉斯版本,正因?yàn)長(zhǎng)DA被貝葉斯化了，所以才會(huì)加的兩個(gè)先驗(yàn)參數(shù)。
LDA模型中一篇文檔生成的方式如下所示:
- 1 按照 選擇一篇文檔 ；
- 2 從狄利克雷分布 中生成文檔i的主題分布 ；
- 3 從主題的多項(xiàng)式分布 中取樣生成文檔i第j個(gè)詞的主題 ；
- 4 從狄利克雷分布 中取樣生成主題 對(duì)應(yīng)的詞語(yǔ)分布 ；
- 5 從詞語(yǔ)的多項(xiàng)式分布 中采樣最終生成詞語(yǔ)
從上面的過(guò)程可以看出，LDA在pLSA的基礎(chǔ)上，為主題分布和詞分布分別加了兩個(gè)Dirichlet先驗(yàn)。
拿之前講解pLSA的例子進(jìn)行具體說(shuō)明。如前所述，在pLSA中，選主題和選詞都是兩個(gè)隨機(jī)的過(guò)程，先從主題分布{教育：0.5，經(jīng)濟(jì)：0.3，交通：0.2}中抽取出主題：教育，然后從該主題對(duì)應(yīng)的詞分布{大學(xué)：0.5，老師：0.3，課程：0.2}中抽取出詞：大學(xué)。 在LDA中，選主題和選詞依然都是兩個(gè)隨機(jī)的過(guò)程。但在LDA中，主題分布和詞分布不再唯一確定不變，即無(wú)法確切給出。例如主題分布可能是{教育：0.5，經(jīng)濟(jì)：0.3，交通：0.2}，也可能是{教育：0.6，經(jīng)濟(jì)：0.2，交通：0.2}，到底是哪個(gè)我們不能確定，因?yàn)樗请S機(jī)的可變化的。 但再怎么變化，也依然服從一定的分布，主題分布和詞分布由Dirichlet先驗(yàn)確定。
舉個(gè)文檔d產(chǎn)生主題z的例子。
在pLSA中，給定一篇文檔d，主題分布是一定的，比如{ P(zi|d), i = 1,2,3 }可能就是{0.4,0.5,0.1}，表示z1、z2、z3這3個(gè)主題被文檔d選中的概率都是個(gè)固定的值：P(z1|d) = 0.4、P(z2|d) = 0.5、P(z3|d) = 0.1，如下圖所示:
在LDA中，主題分布（各個(gè)主題在文檔中出現(xiàn)的概率分布）和詞分布（各個(gè)詞語(yǔ)在某個(gè)主題下出現(xiàn)的概率分布）是唯一確定的。LDA為提供了兩個(gè)Dirichlet先驗(yàn)參數(shù)，Dirichlet先驗(yàn)為某篇文檔隨機(jī)抽取出主題分布和詞分布。
給定一篇文檔d，現(xiàn)在有多個(gè)主題z1、z2、z3，它們的主題分布{ P(zi|d), i = 1,2,3 }可能是{0.4,0.5,0.1}，也可能是{0.2,0.2,0.6}，即這些主題被d選中的概率都不再是確定的值，可能是P(z1|d) = 0.4、P(z2|d) = 0.5、P(z3|d) = 0.1，也有可能是P(z1|d) = 0.2、P(z2|d) = 0.2、P(z3|d) = 0.6，而主題分布到底是哪個(gè)取值集合我們不確定，但其先驗(yàn)分布是dirichlet分布，所以可以從無(wú)窮多個(gè)主題分布中按照dirichlet先驗(yàn)隨機(jī)抽取出某個(gè)主題分布出來(lái)。如下圖所示
LDA在pLSA的基礎(chǔ)上給兩參數(shù)
和 加了兩個(gè)先驗(yàn)分布的參數(shù)。這兩個(gè)分布都是Dirichlet分布。 下面是LDA的圖模型結(jié)構(gòu)：3 LDA 參數(shù)估計(jì)
在spark中，提供了兩種方法來(lái)估計(jì)參數(shù)，分別是變分EM（期望最大）算法（見(jiàn)文獻(xiàn)【3】【4】）和在線(xiàn)學(xué)習(xí)算法（見(jiàn)文獻(xiàn)【5】）。下面將分別介紹這兩種算法以及其源碼實(shí)現(xiàn)。
3.1 變分EM算法
變分貝葉斯算法的詳細(xì)信息可以參考文獻(xiàn)【9】。
在上文中，我們知道LDA將變量theta和phi（為了方便起見(jiàn)，我們將上文LDA圖模型中的beta改為了phi）看做隨機(jī)變量，并且為theta添加一個(gè)超參數(shù)為alpha的Dirichlet先驗(yàn)，為phi添加一個(gè)超參數(shù)為eta的Dirichlet先驗(yàn)來(lái)估計(jì)theta和beta的最大后驗(yàn)（MAP）。 可以通過(guò)最優(yōu)化最大后驗(yàn)估計(jì)來(lái)估計(jì)參數(shù)。我們首先來(lái)定義幾個(gè)變量：
- 下式的gamma表示詞為w，文檔為j時(shí)，主題為k的概率，如公式**（3.1.1）**
- 表示詞w在文檔j中出現(xiàn)的次數(shù)；
- 表示詞w在主題k中出現(xiàn)的次數(shù)，如公式**（3.1.2）**
- 表示主題k在文檔j中出現(xiàn)的次數(shù)，如公式**（3.1.3）**
- 表示主題k中包含的詞出現(xiàn)的總次數(shù)，如公式**（3.1.4）**
- 表示文檔j中包含的主題出現(xiàn)的總次數(shù),如公式**（3.1.5）**
根據(jù)文獻(xiàn)【4】中2.2章節(jié)的介紹，我們可以推導(dǎo)出如下更新公式**(3.1.6)**，其中alpha和eta均大于1：
收斂之后，最大后驗(yàn)估計(jì)可以得到公式**(3.1.7)**：
變分EM算法的流程如下：
- 1 初始化狀態(tài)，即隨機(jī)初始化 和
- 2 E-步，對(duì)每一個(gè)（文檔，詞匯）對(duì)i，計(jì)算 ，更新gamma值
- 3 M-步，計(jì)算隱藏變量phi和theta。即計(jì)算 和
- 4 重復(fù)以上2、3兩步，直到滿(mǎn)足最大迭代數(shù)
第4.2章會(huì)從代碼層面說(shuō)明該算法的實(shí)現(xiàn)流程。
3.2 在線(xiàn)學(xué)習(xí)算法
3.2.1 批量變分貝葉斯
在變分貝葉斯推導(dǎo)(VB)中，根據(jù)文獻(xiàn)【3】，使用一種更簡(jiǎn)單的分布q(z,theta,beta)來(lái)估計(jì)真正的后驗(yàn)分布，這個(gè)簡(jiǎn)單的分布使用一組自由變量(free parameters)來(lái)定義。 通過(guò)最大化對(duì)數(shù)似然的一個(gè)下界（Evidence Lower Bound (ELBO)）來(lái)最優(yōu)化這些參數(shù)，如下公式**(3.2.1)**
最大化ELBO就是最小化q(z,theta,beta)和p(z,theta,beta|w,alpha,eta)的KL距離。根據(jù)文獻(xiàn)【3】，我們將q因式分解為如下**（3.2.2）**的形式：
后驗(yàn)z通過(guò)phi來(lái)參數(shù)化，后驗(yàn)theta通過(guò)gamma來(lái)參數(shù)化，后驗(yàn)beta通過(guò)lambda來(lái)參數(shù)化。為了簡(jiǎn)單描述，我們把lambda當(dāng)作“主題”來(lái)看待。公式**(3.2.2)分解為如下(3.2.3)**形式：
我們現(xiàn)在將上面的期望擴(kuò)展為變分參數(shù)的函數(shù)形式。這反映了變分目標(biāo)只依賴(lài)于
，即詞w出現(xiàn)在文檔d中的次數(shù)。當(dāng)使用VB算法時(shí)，文檔可以通過(guò)它們的詞頻來(lái)匯總（summarized），如公式**(3.2.4)**上面的公式中，W表示詞的數(shù)量，D表示文檔的數(shù)量。l表示文檔d對(duì)ELBO所做的貢獻(xiàn)。L可以通過(guò)坐標(biāo)上升法來(lái)最優(yōu)化，它的更新公式如**(3.2.5)**:
log(theta)和log(beta)的期望通過(guò)下面的公式**(3.2.6)**計(jì)算：
通過(guò)EM算法，我們可以將這些更新分解成E-步和M-步。E-步固定lambda來(lái)更新gamma和phi；M-步通過(guò)給定phi來(lái)更新lambda。批VB算法的過(guò)程如下**(算法1)**所示：
3.2.2 在線(xiàn)變分貝葉斯
批量變分貝葉斯算法需要固定的內(nèi)存，并且比吉布斯采樣更快。但是它仍然需要在每次迭代時(shí)處理所有的文檔，這在處理大規(guī)模文檔時(shí)，速度會(huì)很慢，并且也不適合流式數(shù)據(jù)的處理。 文獻(xiàn)【5】提出了一種在線(xiàn)變分推導(dǎo)算法。設(shè)定gamma(n_d,lambda)和phi(n_d,lambda)分別表示gamma_d和phi_d的值，我們的目的就是設(shè)定phi來(lái)最大化下面的公式**(3.2.7)**
我們?cè)谒惴?中介紹了在線(xiàn)VB算法。因?yàn)樵~頻的第t個(gè)向量
是可觀察的，我們?cè)贓-步通過(guò)固定lambda來(lái)找到gamma_t和phi_t的局部最優(yōu)解。 然后，我們計(jì)算lambda_cap。如果整個(gè)語(yǔ)料庫(kù)由單個(gè)文檔重復(fù)D次組成，那么這樣的lambda_cap設(shè)置是最優(yōu)的。之后，我們通過(guò)lambda之前的值以及l(fā)ambda_cap來(lái)更新lambda。我們給lambda_cap設(shè)置的權(quán)重如公式**(3.2.8)**所示：在線(xiàn)VB算法的實(shí)現(xiàn)流程如下算法2所示
那么在在線(xiàn)VB算法中，alpha和eta是如何更新的呢？參考文獻(xiàn)【8】提供了計(jì)算方法。給定數(shù)據(jù)集，dirichlet參數(shù)的可以通過(guò)最大化下面的對(duì)數(shù)似然來(lái)估計(jì)
其中，
有多種方法可以最大化這個(gè)目標(biāo)函數(shù)，如梯度上升，Newton-Raphson等。Spark使用Newton-Raphson方法估計(jì)參數(shù)，更新alpha。Newton-Raphson提供了一種參數(shù)二次收斂的方法， 它一般的更新規(guī)則如下公式**(3.3.3)**:
其中，H表示海森矩陣。對(duì)于這個(gè)特別的對(duì)數(shù)似然函數(shù)，可以應(yīng)用Newton-Raphson去解決高維數(shù)據(jù)，因?yàn)樗梢栽诰€(xiàn)性時(shí)間求出海森矩陣的逆矩陣。一般情況下，海森矩陣可以用一個(gè)對(duì)角矩陣和一個(gè)元素都一樣的矩陣的和來(lái)表示。 如下公式**(3.3.4)**，Q是對(duì)角矩陣，C11是元素相同的一個(gè)矩陣。
為了計(jì)算海森矩陣的逆矩陣，我們觀察到，對(duì)任意的可逆矩陣Q和非負(fù)標(biāo)量c，有下列式子**(3.3.5)**:
因?yàn)镼是對(duì)角矩陣，所以Q的逆矩陣可以很容易的計(jì)算出來(lái)。所以Newton-Raphson的更新規(guī)則可以重寫(xiě)為如下**(3.3.6)**的形式
其中b如下公式**(3.3.7)**，
4 LDA代碼實(shí)現(xiàn)
4.1 LDA使用實(shí)例
我們從官方文檔【6】給出的使用代碼為起始點(diǎn)來(lái)詳細(xì)分析LDA的實(shí)現(xiàn)。
import org.apache.spark.mllib.clustering.{LDA, DistributedLDAModel} import org.apache.spark.mllib.linalg.Vectors // 加載和處理數(shù)據(jù) val data = sc.textFile("data/mllib/sample_lda_data.txt") val parsedData = data.map(s => Vectors.dense(s.trim.split(' ').map(_.toDouble))) // 為文檔編號(hào)，編號(hào)唯一。List（（id，vector）....） val corpus = parsedData.zipWithIndex.map(_.swap).cache() // 將文檔聚類(lèi)為3類(lèi) val ldaModel = new LDA().setK(3).run(corpus) val topics = ldaModel.topicsMatrix for (topic <- Range(0, 3)) {print("Topic " + topic + ":")for (word <- Range(0, ldaModel.vocabSize)) { print(" " + topics(word, topic)); }println() }以上代碼主要做了兩件事：加載和切分?jǐn)?shù)據(jù)、訓(xùn)練模型。在樣本數(shù)據(jù)中，每一行代表一篇文檔，經(jīng)過(guò)處理后，corpus的類(lèi)型為L(zhǎng)ist((id,vector)*)，一個(gè)(id,vector)代表一篇文檔。將處理后的數(shù)據(jù)傳給org.apache.spark.mllib.clustering.LDA類(lèi)的run方法， 就可以開(kāi)始訓(xùn)練模型。run方法的代碼如下所示：
def run(documents: RDD[(Long, Vector)]): LDAModel = {val state = ldaOptimizer.initialize(documents, this)var iter = 0val iterationTimes = Array.fill[Double](maxIterations)(0)while (iter < maxIterations) {val start = System.nanoTime()state.next()val elapsedSeconds = (System.nanoTime() - start) / 1e9iterationTimes(iter) = elapsedSecondsiter += 1}state.getLDAModel(iterationTimes)}這段代碼首先調(diào)用initialize方法初始化狀態(tài)信息，然后循環(huán)迭代調(diào)用next方法直到滿(mǎn)足最大的迭代次數(shù)。在我們沒(méi)有指定的情況下，迭代次數(shù)默認(rèn)為20。需要注意的是， ldaOptimizer有兩個(gè)具體的實(shí)現(xiàn)類(lèi)EMLDAOptimizer和OnlineLDAOptimizer，它們分別表示使用EM算法和在線(xiàn)學(xué)習(xí)算法實(shí)現(xiàn)參數(shù)估計(jì)。在未指定的情況下，默認(rèn)使用EMLDAOptimizer。
4.2 變分EM算法的實(shí)現(xiàn)
在spark中，使用GraphX來(lái)實(shí)現(xiàn)EMLDAOptimizer，這個(gè)圖是有兩種類(lèi)型的頂點(diǎn)的二分圖。這兩類(lèi)頂點(diǎn)分別是文檔頂點(diǎn)（Document vertices）和詞頂點(diǎn)（Term vertices）。
- 文檔頂點(diǎn)使用大于0的唯一的指標(biāo)來(lái)索引，保存長(zhǎng)度為k（主題個(gè)數(shù)）的向量
- 詞頂點(diǎn)使用{-1, -2, ..., -vocabSize}來(lái)索引，保存長(zhǎng)度為k（主題個(gè)數(shù)）的向量
- 邊（edges）對(duì)應(yīng)詞出現(xiàn)在文檔中的情況。邊的方向是document -> term，并且根據(jù)document進(jìn)行分區(qū)
我們可以根據(jù)3.1節(jié)中介紹的算法流程來(lái)解析源代碼。
4.2.1 初始化狀態(tài)
spark在EMLDAOptimizer的initialize方法中實(shí)現(xiàn)初始化功能。包括初始化Dirichlet參數(shù)alpha和eta、初始化邊、初始化頂點(diǎn)以及初始化圖。
//對(duì)應(yīng)超參數(shù)alphaval docConcentration = lda.getDocConcentration//對(duì)應(yīng)超參數(shù)etaval topicConcentration = lda.getTopicConcentrationthis.docConcentration = if (docConcentration == -1) (50.0 / k) + 1.0 else docConcentrationthis.topicConcentration = if (topicConcentration == -1) 1.1 else topicConcentration上面的代碼初始化了超參數(shù)alpha和eta，根據(jù)文獻(xiàn)【4】，當(dāng)alpha未指定時(shí)，初始化其為(50.0 / k) + 1.0，其中k表示主題個(gè)數(shù)。當(dāng)eta未指定時(shí)，初始化其為1.1。
//對(duì)于每個(gè)文檔，為每一個(gè)唯一的Term創(chuàng)建一個(gè)(Document->Term)的邊 val edges: RDD[Edge[TokenCount]] = docs.flatMap { case (docID: Long, termCounts: Vector) =>// Add edges for terms with non-zero counts.termCounts.toBreeze.activeIterator.filter(_._2 != 0.0).map { case (term, cnt) =>//文檔id，termindex，詞頻Edge(docID, term2index(term), cnt)} } //term2index將term轉(zhuǎn)為{-1, -2, ..., -vocabSize}索引private[clustering] def term2index(term: Int): Long = -(1 + term.toLong)上面的這段代碼處理每個(gè)文檔，對(duì)文檔中每個(gè)唯一的Term（詞）創(chuàng)建一個(gè)邊，邊的格式為(文檔id，詞索引，詞頻)。詞索引為{-1, -2, ..., -vocabSize}。
//創(chuàng)建頂點(diǎn)val docTermVertices: RDD[(VertexId, TopicCounts)] = {val verticesTMP: RDD[(VertexId, TopicCounts)] =edges.mapPartitionsWithIndex { case (partIndex, partEdges) =>val random = new Random(partIndex + randomSeed)partEdges.flatMap { edge =>val gamma = normalize(BDV.fill[Double](k)(random.nextDouble()), 1.0)//此處的sum是DenseVector，gamma*N_wjval sum = gamma * edge.attr//srcId表示文獻(xiàn)id，dstId表是詞索引Seq((edge.srcId, sum), (edge.dstId, sum))}}verticesTMP.reduceByKey(_ + _)}上面的代碼創(chuàng)建頂點(diǎn)。我們?yōu)槊總€(gè)主題隨機(jī)初始化一個(gè)值，即gamma是隨機(jī)的。sum為gamma * edge.attr，這里的edge.attr即N_wj,所以sum用gamma * N_wj作為頂點(diǎn)的初始值。
this.graph = Graph(docTermVertices, edges).partitionBy(PartitionStrategy.EdgePartition1D)上面的代碼初始化Graph并通過(guò)文檔分區(qū)。
4.2.2 E-步：更新gamma
val eta = topicConcentrationval W = vocabSizeval alpha = docConcentrationval N_k = globalTopicTotalsval sendMsg: EdgeContext[TopicCounts, TokenCount, (Boolean, TopicCounts)] => Unit =(edgeContext) => {// 計(jì)算 N_{wj} gamma_{wjk}val N_wj = edgeContext.attr// E-STEP: 計(jì)算 gamma_{wjk} 通過(guò)N_{wj}來(lái)計(jì)算//此處的edgeContext.srcAttr為當(dāng)前迭代的N_kj , edgeContext.dstAttr為當(dāng)前迭代的N_wk,//后面通過(guò)M-步，會(huì)更新這兩個(gè)值,作為下一次迭代的當(dāng)前值val scaledTopicDistribution: TopicCounts = computePTopic(edgeContext.srcAttr, edgeContext.dstAttr, N_k, W, eta, alpha) *= N_wjedgeContext.sendToDst((false, scaledTopicDistribution))edgeContext.sendToSrc((false, scaledTopicDistribution))}上述代碼中，W表示詞數(shù)，N_k表示所有文檔中，出現(xiàn)在主題k中的詞的詞頻總數(shù)，后續(xù)的實(shí)現(xiàn)會(huì)使用方法computeGlobalTopicTotals來(lái)更新這個(gè)值。N_wj表示詞w出現(xiàn)在文檔j中的詞頻數(shù)，為已知數(shù)。E-步就是利用公式**(3.1.6)**去更新gamma。 代碼中使用computePTopic方法來(lái)實(shí)現(xiàn)這個(gè)更新。edgeContext通過(guò)方法sendToDst將scaledTopicDistribution發(fā)送到目標(biāo)頂點(diǎn)， 通過(guò)方法sendToSrc發(fā)送到源頂點(diǎn)以便于后續(xù)的M-步更新的N_kj和N_wk。下面我們看看computePTopic方法。
private[clustering] def computePTopic(docTopicCounts: TopicCounts,termTopicCounts: TopicCounts,totalTopicCounts: TopicCounts,vocabSize: Int,eta: Double,alpha: Double): TopicCounts = {val K = docTopicCounts.lengthval N_j = docTopicCounts.dataval N_w = termTopicCounts.dataval N = totalTopicCounts.dataval eta1 = eta - 1.0val alpha1 = alpha - 1.0val Weta1 = vocabSize * eta1var sum = 0.0val gamma_wj = new Array[Double](K)var k = 0while (k < K) {val gamma_wjk = (N_w(k) + eta1) * (N_j(k) + alpha1) / (N(k) + Weta1)gamma_wj(k) = gamma_wjksum += gamma_wjkk += 1}// normalizeBDV(gamma_wj) /= sum}這段代碼比較簡(jiǎn)單，完全按照公式**(3.1.6)**表示的樣子來(lái)實(shí)現(xiàn)。val gamma_wjk = (N_w(k) + eta1) * (N_j(k) + alpha1) / (N(k) + Weta1)就是實(shí)現(xiàn)的更新邏輯。
4.2.3 M-步：更新phi和theta
// M-STEP: 聚合計(jì)算新的 N_{kj}, N_{wk} counts. val docTopicDistributions: VertexRDD[TopicCounts] =graph.aggregateMessages[(Boolean, TopicCounts)](sendMsg, mergeMsg).mapValues(_._2)我們由公式**(3.1.7)**可知，更新隱藏變量phi和theta就是更新相應(yīng)的N_kj和N_wk。聚合更新使用aggregateMessages方法來(lái)實(shí)現(xiàn)。請(qǐng)參考文獻(xiàn)【7】來(lái)了解該方法的作用。
4.3 在線(xiàn)變分算法的代碼實(shí)現(xiàn)
4.3.1 初始化狀態(tài)
在線(xiàn)學(xué)習(xí)算法首先使用方法initialize方法初始化參數(shù)值
override private[clustering] def initialize(docs: RDD[(Long, Vector)],lda: LDA): OnlineLDAOptimizer = {this.k = lda.getKthis.corpusSize = docs.count()this.vocabSize = docs.first()._2.sizethis.alpha = if (lda.getAsymmetricDocConcentration.size == 1) {if (lda.getAsymmetricDocConcentration(0) == -1) Vectors.dense(Array.fill(k)(1.0 / k))else {require(lda.getAsymmetricDocConcentration(0) >= 0,s"all entries in alpha must be >=0, got: $alpha")Vectors.dense(Array.fill(k)(lda.getAsymmetricDocConcentration(0)))}} else {require(lda.getAsymmetricDocConcentration.size == k,s"alpha must have length k, got: $alpha")lda.getAsymmetricDocConcentration.foreachActive { case (_, x) =>require(x >= 0, s"all entries in alpha must be >= 0, got: $alpha")}lda.getAsymmetricDocConcentration}this.eta = if (lda.getTopicConcentration == -1) 1.0 / k else lda.getTopicConcentrationthis.randomGenerator = new Random(lda.getSeed)this.docs = docs// 初始化變分分布 q(beta|lambda)this.lambda = getGammaMatrix(k, vocabSize)this.iteration = 0this}根據(jù)文獻(xiàn)【5】，alpha和eta的值大于等于0，并且默認(rèn)為1.0/k。上文使用getGammaMatrix方法來(lái)初始化變分分布q(beta|lambda)。
private def getGammaMatrix(row: Int, col: Int): BDM[Double] = {val randBasis = new RandBasis(new org.apache.commons.math3.random.MersenneTwister(randomGenerator.nextLong()))//初始化一個(gè)gamma分布val gammaRandomGenerator = new Gamma(gammaShape, 1.0 / gammaShape)(randBasis)val temp = gammaRandomGenerator.sample(row * col).toArraynew BDM[Double](col, row, temp).t}getGammaMatrix方法使用gamma分布初始化一個(gè)隨機(jī)矩陣。
4.3.2 更新參數(shù)
override private[clustering] def next(): OnlineLDAOptimizer = {//返回文檔集中采樣的子集//默認(rèn)情況下，文檔可以被采樣多次，且采樣比例是0.05val batch = docs.sample(withReplacement = sampleWithReplacement, miniBatchFraction,randomGenerator.nextLong())if (batch.isEmpty()) return thissubmitMiniBatch(batch)}以上的next方法首先對(duì)文檔進(jìn)行采樣，然后調(diào)用submitMiniBatch對(duì)采樣的文檔子集進(jìn)行處理。下面我們?cè)敿?xì)分解submitMiniBatch方法。
- 1 計(jì)算log(beta)的期望，并將其作為廣播變量廣播到集群中
上述代碼調(diào)用exp(LDAUtils.dirichletExpectation(lambda))方法實(shí)現(xiàn)參數(shù)為lambda的log beta的期望。實(shí)現(xiàn)原理參見(jiàn)公式**(3.2.6)**。
- 2 計(jì)算phi以及gamma，即算法2中的E-步
上面的代碼調(diào)用OnlineLDAOptimizer.variationalTopicInference實(shí)現(xiàn)算法2中的E-步,迭代計(jì)算phi和gamma。
private[clustering] def variationalTopicInference(termCounts: Vector,expElogbeta: BDM[Double],alpha: breeze.linalg.Vector[Double],gammaShape: Double,k: Int): (BDV[Double], BDM[Double]) = {val (ids: List[Int], cts: Array[Double]) = termCounts match {case v: DenseVector => ((0 until v.size).toList, v.values)case v: SparseVector => (v.indices.toList, v.values)}// 初始化變分分布 q(theta|gamma) val gammad: BDV[Double] = new Gamma(gammaShape, 1.0 / gammaShape).samplesVector(k) // K//根據(jù)公式（3.2.6）計(jì)算 E(log theta)val expElogthetad: BDV[Double] = exp(LDAUtils.dirichletExpectation(gammad)) // Kval expElogbetad = expElogbeta(ids, ::).toDenseMatrix // ids * K//根據(jù)公式（3.2.5）計(jì)算phi，這里加1e-100表示并非嚴(yán)格等于val phiNorm: BDV[Double] = expElogbetad * expElogthetad :+ 1e-100 // idsvar meanGammaChange = 1Dval ctsVector = new BDV[Double](cts) // ids// 迭代直至收斂while (meanGammaChange > 1e-3) {val lastgamma = gammad.copy//依據(jù)公式(3.2.5)計(jì)算gammagammad := (expElogthetad :* (expElogbetad.t * (ctsVector :/ phiNorm))) :+ alpha//根據(jù)更新的gamma，計(jì)算E(log theta)expElogthetad := exp(LDAUtils.dirichletExpectation(gammad))// 更新phiphiNorm := expElogbetad * expElogthetad :+ 1e-100//計(jì)算兩次gamma的差值meanGammaChange = sum(abs(gammad - lastgamma)) / k}val sstatsd = expElogthetad.asDenseMatrix.t * (ctsVector :/ phiNorm).asDenseMatrix(gammad, sstatsd)}3 更新lambda
val statsSum: BDM[Double] = stats.map(_._1).reduce(_ += _)val gammat: BDM[Double] = breeze.linalg.DenseMatrix.vertcat(stats.map(_._2).reduce(_ ++ _).map(_.toDenseMatrix): _*)val batchResult = statsSum :* expElogbeta.t// 更新lambda和alphaupdateLambda(batchResult, (miniBatchFraction * corpusSize).ceil.toInt)updateLambda方法實(shí)現(xiàn)算法2中的M-步,更新lambda。實(shí)現(xiàn)代碼如下：
private def updateLambda(stat: BDM[Double], batchSize: Int): Unit = {// 根據(jù)公式（3.2.8）計(jì)算權(quán)重val weight = rho()// 更新lambda，其中stat * (corpusSize.toDouble / batchSize.toDouble)+eta表示rho_caplambda := (1 - weight) * lambda +weight * (stat * (corpusSize.toDouble / batchSize.toDouble) + eta)} // 根據(jù)公式（3.2.8）計(jì)算rho private def rho(): Double = {math.pow(getTau0 + this.iteration, -getKappa)}4 更新alpha
private def updateAlpha(gammat: BDM[Double]): Unit = {//計(jì)算rhoval weight = rho()val N = gammat.rows.toDoubleval alpha = this.alpha.toBreeze.toDenseVector//計(jì)算log p_hatval logphat: BDM[Double] = sum(LDAUtils.dirichletExpectation(gammat)(::, breeze.linalg.*)) / N//計(jì)算梯度為N（-phi(alpha)+log p_hat）val gradf = N * (-LDAUtils.dirichletExpectation(alpha) + logphat.toDenseVector)//計(jì)算公式（3.3.4）中的c，trigamma表示gamma函數(shù)的二階導(dǎo)數(shù)val c = N * trigamma(sum(alpha))//計(jì)算公式（3.3.4）中的qval q = -N * trigamma(alpha)//根據(jù)公式(3.3.7)計(jì)算bval b = sum(gradf / q) / (1D / c + sum(1D / q))val dalpha = -(gradf - b) / qif (all((weight * dalpha + alpha) :> 0D)) {alpha :+= weight * dalphathis.alpha = Vectors.dense(alpha.toArray)}}5 參考文獻(xiàn)
【1】LDA數(shù)學(xué)八卦
【2】通俗理解LDA主題模型
【3】[Latent Dirichlet Allocation](docs/Latent Dirichlet Allocation.pdf)
【4】[On Smoothing and Inference for Topic Models](docs/On Smoothing and Inference for Topic Models.pdf)
【5】[Online Learning for Latent Dirichlet Allocation](docs/Online Learning for Latent Dirichlet Allocation.pdf)
【6】Spark官方文檔
【7】Spark GraphX介紹
【8】Maximum Likelihood Estimation of Dirichlet Distribution Parameters
【9】Variational Bayes
星環(huán)科技：機(jī)器學(xué)習(xí)算法2020線(xiàn)上訓(xùn)練營(yíng) | 首營(yíng)即將開(kāi)班，限額報(bào)名！?zhuanlan.zhihu.com總結(jié)
以上是生活随笔為你收集整理的狄利克雷分布公式_深入机器学习系列11－隐式狄利克雷分布的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。
- 上一篇： springboot中spock如何使用
- 下一篇： 恋与制作人许墨超能力