0. latent dirichlet allocation 前言

最近公司分享了這個(gè)topic兔沃， 我自己鉆研了一下寫一下自己對這個(gè)模型的理解橄碾。
因?yàn)楸救嗽陂喿x網(wǎng)上各種資料時(shí)卵沉，發(fā)現(xiàn)寫的比較不適合新手。所以寫下這篇博文法牲。
一是史汗，對自己的學(xué)習(xí)收獲做一個(gè)總結(jié)。另外拒垃，是希望可以幫助新手在學(xué)習(xí)這個(gè)模型的過程中能夠獲得比較好的理解停撞。
我覺得理解一個(gè)事物的最好方式是，是從最高層級開始悼瓮，知道它的全貌戈毒，忽略它的細(xì)枝末節(jié)。隨后看到一個(gè)高粒度抽象的架構(gòu)横堡，去理解它的工作模式埋市。最后才是深入細(xì)節(jié)底層。去逐個(gè)擊破各個(gè)抽象的背后實(shí)現(xiàn)原理命贴。
那么我就用這種方式來寫這篇文章道宅。
所以這篇文章的寫法會和網(wǎng)上大多數(shù)的不同。他們都是從最底層需要的數(shù)學(xué)知識講起胸蛛，做了大量公式推導(dǎo)污茵。開始看的云里霧里，之后串的時(shí)候就會看了后面忘了前面葬项，思維跳躍跨度很大泞当。下面是這篇文章是怎么展開的，好讓讀者心中有數(shù)民珍。

1. 把這個(gè)模型看成黑盒襟士，它需要什么，它能給你什么（最頂層）
2. 理解這個(gè)模型有哪些組件構(gòu)成嚷量，這些組件怎么合作的 （看到高粒度的組件）
3. 這些組件背后的數(shù)學(xué)和算法是什么敌蜂，這些原理間的聯(lián)系是什么 （高粒度組件下用了哪些數(shù)學(xué)算法工具）
4. 我們?nèi)绾位谶@些原理去實(shí)現(xiàn)自己的LDA （不用知道推導(dǎo)，已經(jīng)可以實(shí)現(xiàn)這個(gè)模型了）
5. 背后的數(shù)學(xué)詳細(xì)推導(dǎo) （每個(gè)數(shù)學(xué)結(jié)論是怎么得來的）

1. 這個(gè)模型可以解決什么問題(頂層）

當(dāng)我們有一組文章（如文章A津肛， 文章B章喉， 文章C）； 我們希望知道每個(gè)文章有哪些主題身坐。（如主題A秸脱， 主題B， 主題C）我們就可以使用這個(gè)模型部蛇。與此同時(shí)摊唇，又來了一篇新的文章，我們也可以丟進(jìn)這個(gè)模型里涯鲁，它會告訴我們新文章里有哪些主題巷查。
解決的問題分為2個(gè)階段

1. 訓(xùn)練階段：
   輸入是一組文章（article = list <word>, 所以輸入是list<article>）和 一個(gè)數(shù)字K（表示這組文章里你覺得大概有多少個(gè)TOPIC）有序。 這個(gè)K是個(gè)需要人們給的參數(shù)，會決定這個(gè)模型的表現(xiàn)岛请，所以需要調(diào)參俠旭寿，去對K做調(diào)整，得到最好的模型效果崇败。
   輸出 輸入的文章盅称，每一篇都會給一個(gè)文章到TOPIC的概率分布。 每一個(gè)TOPIC后室，都會告訴你這個(gè)TOPIC下每一個(gè)單詞的概率分布缩膝。

具體例子：

假設(shè)我們有4篇文章, 并且指定K為3（我認(rèn)為這組文章應(yīng)該用3個(gè)TOPIC就可以囊括了）

文章a: ball ball ball planet galaxy 
文章b: referendum planet planet referendum referendum
文章c: planet planet galaxy planet ball
文章d: planet galaxy referendum planet ball

可能的輸出如下：

TOPIC下每一個(gè)單詞的概率分布

topic 1: {ball:70%, galaxy:25%, referendum :1%, planet : 4%}
topic 2: {ball:9%, galaxy:35%, referendum :1%, planet : 55%}
topic 3: {ball:2%, galaxy:3%, referendum :80%, planet : 15%}

根據(jù)觀察我們可以根據(jù)人類的知識，總結(jié)出topic1可能是體育岸霹。因?yàn)橛星蚣膊悖鳪ALAXY可能是一只球隊(duì)名。 topic2可能是科學(xué)贡避，因?yàn)橹v行星和銀河系云芦。 topic3可能是講社會科學(xué)

每一個(gè)文章到TOPIC的概率分布

文章a: {topic1: 80%, topic2: 10%, topic3: 10%} 文章極大概率是體育類
文章b: {topic1: 5%, topic2: 85%, topic3: 10%} 文章極大概率是社科類
文章c: {topic1: 10%, topic2: 2%, topic3: 88%} 文章極大概率是科學(xué)類
文章d: {topic1: 25%, topic2: 15%, topic3: 60%} 文章大概率是科學(xué)類

2. 預(yù)測階段：
   輸入一篇新文章

galaxy ball ball ball ball

輸出這篇文章的概率分布

{topic1: 90%, topic2: 1%, topic3: 9%} 

2. 理解這個(gè)模型有哪些組件構(gòu)成，這些組件怎么合作的 （高粒度的組件）

這個(gè)模型有4個(gè)組件贸桶，2個(gè)旋鈕舅逸，2組齒輪。
當(dāng)我們決定了K是多少之后皇筛。整個(gè)LDA模型就可以想象成一個(gè)機(jī)器如下圖琉历。它有2個(gè)可以調(diào)控的旋鈕，旋鈕調(diào)整的好壞決定了模型的好壞水醋。調(diào)整這2個(gè)旋鈕的作用就是決定內(nèi)部齒輪的工作方式旗笔。

image.png

我們按下紅色按鈕，這個(gè)機(jī)器就會輸出一篇文章拄踪。這個(gè)文章是一個(gè)裝著很多單詞的袋子蝇恶。（就是說單詞之間的順序不重要）
那么模型輸出的文章 越 接近 我們預(yù)先給它的文集，我們就認(rèn)為這個(gè)模型質(zhì)量越好惶桐。

image.png

首先第一個(gè)旋鈕，會按照概率分布D1肠槽，隨機(jī)出一個(gè)多項(xiàng)式分布值M1（eg. topic1:50%擎淤，topic2：30%, topic3:20%）我們告訴第一個(gè)齒輪奢啥。然后第二個(gè)旋鈕按照概率分布D2，隨機(jī)出K(例子中是3)個(gè)多項(xiàng)式分布值
如：

topic 1: {ball:70%, galaxy:25%, referendum :1%, planet : 4%}
topic 2: {ball:9%, galaxy:35%, referendum :1%, planet : 55%}
topic 3: {ball:2%, galaxy:3%, referendum :80%, planet : 15%}

隨后第一個(gè)齒輪拿到M1（eg. topic1:50%嘴拢，topic2：30%, topic3:20%）桩盲， 根據(jù)這個(gè)概率隨機(jī)出一個(gè)topic 告訴第二個(gè)齒輪（eg. topic1）。 第二個(gè)齒輪根據(jù)這個(gè)TOPIC 1和第二個(gè)旋鈕給他的這組分布炊汤，用topic 1: {ball:70%, galaxy:25%, referendum :1%, planet : 4%} 隨機(jī)出一個(gè)詞=》比如ball
循環(huán)上述這個(gè)過程，又可以得到一個(gè)詞弊攘。這樣一直輸出到文章要求的詞數(shù)N抢腐，那么就產(chǎn)生了一篇文章。
這2個(gè)旋鈕控制的是一個(gè)叫狄里克雷分布的函數(shù)襟交。而2個(gè)齒輪則是多項(xiàng)式分布迈倍。這里你只要知道這2個(gè)名詞。也不妨礙理解第二層的結(jié)構(gòu)捣域。我會在第3層講這些組件背后的數(shù)學(xué)和算法是什么的時(shí)候啼染，展開介紹你需要理解第3層要知道的這2個(gè)分布的意義。那么我們就有了下圖焕梅。

image.png

image.png

下面我們直觀的來感受下這個(gè)公式，用可視化的方式知道這個(gè)模型的運(yùn)作機(jī)制牛哺。
在這里我們把一個(gè)2維的狄里克雷分布（又叫beta分布）表示為一條線陋气。一個(gè)樣本點(diǎn)可以落在這條線的任何一個(gè)地方。這條線的范圍區(qū)間是【0,1】引润。
小伙伴這邊可以打開這個(gè)網(wǎng)站和我一起動(dòng)手玩一些巩趁，一旦你理解了beta的分布和二項(xiàng)分布的聯(lián)系。那么狄里克雷只不過是把beta分布的2維變到K維淳附。二項(xiàng)分布的2項(xiàng)變成多項(xiàng)分布的K項(xiàng)议慰。

2.1 理解beta分布和二項(xiàng)分布

我們從上面的模型介紹可以知道，狄里克雷分布可以按一定的概率分布來隨機(jī)出一個(gè)多項(xiàng)分布奴曙。那么對應(yīng)到2維的情況别凹。等價(jià)于一個(gè)beta分布 按一定的概率分布可以隨機(jī)出一個(gè)二項(xiàng)分布。

我們控制beta分布的旋鈕有2個(gè)變量（a, b）代表2個(gè)維度洽糟。比如丟硬幣的話番川，a代表正面向上，b代表背面向上脊框。

beta分布颁督，試玩https://keisan.casio.com/exec/system/1180573226； 我們先把A=1, B=1; 觀察一下beta分布會長成什么樣浇雹？(下面藍(lán)色的區(qū)域是對X坐標(biāo)軸縮放用的沉御，我們暫時(shí)用不到。我們只要關(guān)于調(diào)整a, b就可以了昭灵。)

image.png

我們發(fā)現(xiàn)圖像變成了一個(gè)直線吠裆，并且圍成的面積是1. 上面圖中的面積是1，所以要求一個(gè)二項(xiàng)分布產(chǎn)生的概率烂完，我們就要對一個(gè)X的區(qū)間做積分试疙。比如a=1,b=1的beta分布產(chǎn)生 二項(xiàng)分布在0.5~0.6的概率 為下圖這個(gè)面積的大小。

image.png

如果我們想更精確到0.5的概率抠蚣。我們就讓delta X比較小即可祝旷。比如0.5~0.500001
二項(xiàng)分布在0.5~0.6的概率 的意思是 我們從這個(gè)beta分布隨機(jī)出一個(gè)2項(xiàng)分布。
這個(gè)二項(xiàng)分布，（正面向上的概率為(0.5~0.6)中的一個(gè)數(shù)怀跛。背面向上的概率為1-正面向上的概率）的概率

比如我們這次用這個(gè)beta分布隨機(jī)出0.55這個(gè)值距贷。 就等價(jià)我們隨機(jī)出2個(gè)topic下，每個(gè)topic的概率是多少的多項(xiàng)分布吻谋。為（topic1:55%, topic2:45%）

接下來會引入忠蝗，共軛的概念，beta分布是二項(xiàng)分布的共軛分布漓拾。

大白話的意思是阁最，我們有一個(gè)beta分布的圖像。然后我們觀察到一批二項(xiàng)分布的樣本骇两。我們把樣本的數(shù)據(jù)丟進(jìn)這個(gè) beta分布速种，他會根據(jù)這組數(shù)據(jù)重新生成一個(gè)同時(shí)貼近數(shù)據(jù)和原來beta分布的新beta分布。
比如這個(gè)硬幣脯颜，我們丟了50次（在二項(xiàng)分布上采樣了50個(gè)樣本）哟旗；那么我們就得到50個(gè)數(shù)據(jù)贩据。其中40次正面向上栋操，10次反面向上。
按照beta分布和二項(xiàng)分布的公式（式2饱亮， 數(shù)學(xué)推導(dǎo)在最后矾芙，現(xiàn)在只要記住結(jié)論理解過程即可）

image.png

image.png

image.png

值得注意的是岖是，因?yàn)閲饋淼目偯娣e始終要等于1（因?yàn)楦怕屎蜑?）.因?yàn)閄的可能范圍不斷縮小帮毁。所以Y會增大。所以我們可以看縱坐標(biāo)軸的大小是在變大的豺撑。

那么總結(jié)一下：
1. 把一個(gè)二項(xiàng)分布產(chǎn)生出的數(shù)據(jù)結(jié)果 丟給一個(gè)beta分布烈疚，我們會得到一個(gè)新的beta分布。在專業(yè)術(shù)語里聪轿，就是把一個(gè)似然函數(shù) 丟給一個(gè) 先驗(yàn)分布爷肝，會得到一個(gè)后驗(yàn)分布，并且后驗(yàn)分布和先驗(yàn)分布 是 同一種分布。

2. 數(shù)據(jù)量越大阶剑，模型就越可以聚攏跃巡。當(dāng)數(shù)量趨近于無窮。beta分布退化至二項(xiàng)分布牧愁。比如上圖素邪，就會成為在0.8上的一條線，Y在【0.8,0.8+無窮小】為無窮大猪半。

我們理解了2維的情況兔朦。我們可以以此類推到3維。就引出了狄里克雷分布和多項(xiàng)式分布磨确。其實(shí)本質(zhì)是一樣的沽甥。為了更好的用圖像表示，我們把一個(gè)Y值交高的區(qū)域乏奥，用更深的顏色去表示摆舟。就是概率密度高的地方顏色深，然后低的地方淺邓了。

image.png

這樣我們就可以用畫出一個(gè)3維的狄里克雷分布的分布情況恨诱。不過要強(qiáng)調(diào)的是，如果顏色的深度用一個(gè)值可以表示的話骗炉，那么區(qū)域里顏色的總和是個(gè)恒定值照宝。就和上面的1一樣。也就是顏色總量不變句葵。我們可以讓他均勻分配厕鹃。也可以有些地方深（出現(xiàn)的概率高），有些地方淺（出現(xiàn)的概率低）
因?yàn)槭?維了乍丈，所以我們由原來的beta分布里的a, b 2個(gè)參數(shù)剂碴。變?yōu)閍1, a2, a3 3個(gè)參數(shù)。 那么我們就可以從下面的等邊三角形里轻专，隨機(jī)出一個(gè)點(diǎn)忆矛。這個(gè)點(diǎn)到3個(gè)頂點(diǎn)的距離的遠(yuǎn)近分別代表對應(yīng)3種可能的概率。越近代表概率越高铭若。
如下圖洪碳，

image.png

當(dāng)然我們的分布也可以是叼屠，

image.png

你們就會發(fā)現(xiàn)瞳腌，3個(gè)角的顏色很深，就代表點(diǎn)出現(xiàn)在角落的概率就很高镜雨。這個(gè)就好比一個(gè)派對嫂侍。每個(gè)角落都有一個(gè)正向的刺激，那么人們就會傾向于聚集到各個(gè)角落。

image.png

我們可以看到三角形中黃色的點(diǎn)碎浇，就是隨機(jī)出來的一個(gè)3項(xiàng)分布临谱。這里比如某個(gè)點(diǎn)可能是（{topic1: 90%, topic2: 1%, topic3: 9%} ）， 在圖中具體的位置奴璃，就是離角落1很近悉默。離角落2和角落3都很遠(yuǎn)。

image.png

如果我們有了大量樣本數(shù)據(jù)：如這組文集里苟穆。有1000個(gè)TOPIC1抄课， 10個(gè)TOPIC2， 10個(gè)TOPIC3. 我們可以讓原本顏色均勻的三角形雳旅，會根據(jù)新到來的數(shù)據(jù)調(diào)整自己的顏色分布跟磨。

這個(gè)調(diào)整的目的就是為了這個(gè)三角形下次隨機(jī)出來的結(jié)果，能夠更加接近于樣本集岭辣。

這里就和上面串起來了吱晒，訓(xùn)練模型甸饱，其實(shí)就是把一組文章當(dāng)成輸入沦童，而獲得一個(gè)更接近樣本的三角形（狄里克雷分布）

根據(jù)我們之前的例子，我們設(shè)定了3個(gè)主題（體育叹话，社科偷遗，科學(xué)）， 我們有4個(gè)詞（planet, galaxy, ball, referendum）那么從主題到詞的那個(gè)狄里克雷分布驼壶，我們就可以用下圖來表示氏豌。這張圖恒定只有3個(gè)點(diǎn)，因?yàn)槲覀冇?個(gè)主題热凹。

image.png

最后我們把那個(gè)公式可視化就變成了下圖泵喘。后面2個(gè)方盒就是多項(xiàng)分布。3個(gè)顏色的球各1個(gè)再箱子里般妙，代表{1/3纪铺， 1/3, 1/3}的多項(xiàng)分布

image.png

所以整個(gè)模型的運(yùn)作，就是現(xiàn)在第一個(gè)三角形李根據(jù)顏色分布的概率隨機(jī)出一個(gè)球的位置碟渺，得到topic一個(gè)分布數(shù)值
第二個(gè)正四面體鲜锚，根據(jù)顏色密度，隨機(jī)出3個(gè)球，得到TOPIC到詞的概率芜繁。

image.png

上面這2個(gè)概率值旺隙，就能構(gòu)建出對應(yīng)的BOX。然后針對每次要隨機(jī)的詞骏令，我們先選TOPIC蔬捷，然后在到對應(yīng)的TOPIC->詞的BOX里選詞。
最后構(gòu)建出了整個(gè)文章榔袋。

上面就是LDA模型的全過程抠刺。

下面我們就是要了解我們?nèi)绾斡?xùn)練模型，找到那個(gè)best setting摘昌？

image.png

在深入介紹訓(xùn)練之前我們先

2.2 暫退傺回顧一下

到目前為止，我們介紹了LDA的模型聪黎。
希望你可以

1. 自己畫出我們上面貼過的那張概率圖模型的樣子罕容。也就是從alpha, beta為源頭 最后生成到w的過程。并且可以自己寫一下每個(gè)詞的概率的計(jì)算公式稿饰。
2. 在思考一下锦秒， beta分布和二項(xiàng)分布的關(guān)系。beta分布如何隨機(jī)出一個(gè)二項(xiàng)分布喉镰。二項(xiàng)分布的數(shù)據(jù)是怎么反向更新beta分布的旅择。那么擴(kuò)展到K項(xiàng)，本質(zhì)是一樣的侣姆。

3. 結(jié)合我們的第一節(jié)的頂層設(shè)計(jì)生真，再理解一下，我們的輸入其實(shí)扮演的是多項(xiàng)分布的數(shù)據(jù)作用捺宗。但是他不是直接能反向更新2個(gè)狄里克雷分布的柱蟀。因?yàn)樗谋壤磻?yīng)的是文章到詞的多項(xiàng)分布蚜厉。我們需要的是文章到主題的多項(xiàng)分布长已，和 主題到詞的多項(xiàng)分布。這中間有一個(gè)主題在這里昼牛。那么我們應(yīng)該怎么根據(jù)我們先有的（文章到詞的樣本） 轉(zhuǎn)換 為我們需要的 （文章到主題的樣本术瓮，主題到詞的樣本）才是我們真正的目標(biāo)。
   而根據(jù)狄里克雷分布和多項(xiàng)分布的公式贰健。有了樣本轴术，我們只要加回去胃惜，就能得到新的參數(shù)了。這不是難點(diǎn)汹忠。

   
   
   image.png

下一章我就會介紹LDA訓(xùn)練里最核心的吉布斯采樣，他是怎么根據(jù)
（文章到詞的樣本），來達(dá)成對（文章到主題的樣本，主題到詞的樣本）。 而有了這些樣本甘苍，我們又是如何得到我們頂層設(shè)計(jì)里提到的輸出。
輸出：（每一篇都會給一個(gè)文章到TOPIC的概率分布烘豌。 每一個(gè)TOPIC载庭，都會告訴你這個(gè)TOPIC下每一個(gè)單詞的概率分布。）

3. 這些組件背后的數(shù)學(xué)和算法是什么廊佩，這些原理間的聯(lián)系是什么 （高粒度組件下用了哪些數(shù)學(xué)算法工具）

我們先來看一下我們手上有哪些組件了囚聚。

1. 輸入（一組文章，一個(gè)K）
2. 輸出 (每篇文章到主題的分布标锄， 每個(gè)主題到單詞的分布）
3. 文章到主題的多項(xiàng)分布生成器-狄里克雷分布alpha（旋鈕1）
4. 主題到詞匯的多項(xiàng)分布生成器-狄里克雷分布beta (旋鈕2）
5. m個(gè)文章到主題的多項(xiàng)式分布（齒輪1）
6. k個(gè)主題到詞的多項(xiàng)式分布（齒輪2）
7. 吉布斯采樣算法（內(nèi)置算法顽铸，用于獲得文章到主題的樣本，主題到詞的樣本）

3料皇，4號組件背后的數(shù)學(xué)： 狄里克雷分布 可以 根據(jù)概率密度 生成 多項(xiàng)式分布

5,6 號組件背后的數(shù)學(xué)： 如果我們有文章下主題的樣本數(shù)據(jù)谓松，我們可以優(yōu)化3號組件。如果我們有主題下詞的樣本數(shù)據(jù)践剂，我們可以優(yōu)化4號組件鬼譬。

如果3，4號組件優(yōu)化的越好逊脯，那么根據(jù)概率公式优质，我們產(chǎn)生的文檔所構(gòu)成的詞袋（順序無關(guān)的詞的集合）和原始輸入本章集合的詞袋越接近。

因?yàn)閠opic的定義是由原始輸入集定義的军洼。所以來了一篇新的文章巩螃，我們的模型參數(shù)更加符合原始的詞袋，也就意味著預(yù)測它基于文章里的詞歉眷，屬于的topic有哪些就更準(zhǔn)確牺六。

3.1下面就是如何訓(xùn)練的問題

如果沒有吉布斯采樣颤枪，我們會怎么解決這個(gè)問題汗捡。

我們現(xiàn)在手里有的是，文章-》單詞的樣本數(shù)據(jù)畏纲。我們希望有的是文章-》主題扇住，和 主題-》單詞的樣本數(shù)據(jù)。因?yàn)槲覀兿Ｍ聪蚋碌男o的定義是文章-》主題盗胀，主題-》單詞艘蹋。

那么我們的目標(biāo)就是可以估算出
p(z | w) 的概率。 這個(gè)式子的定義是票灰，已知詞是什么的情況下女阀，它是z topic 的條件概率是多少宅荤。
比如我們有一個(gè)文章 （ball, galaxy）。 我們希望知道 第一個(gè)ball 它是（topic1: x%, topic2: y%, topic3: d%）

p(z | w) 如果word 是這篇文章的第一個(gè)單詞,z=2浸策，他的值就是y%
隨后我們在W和Z上加一個(gè)剪頭冯键，表示W(wǎng)是向量，那么Z也是一個(gè)向量了庸汗。也就是已知多個(gè)單詞惫确，那么每個(gè)單詞分別是Z向量里對應(yīng)的z的概率是多少。

我們希望得到這個(gè)是因?yàn)?有了這個(gè)概率分布蚯舱，我們就可以用這個(gè)分布 對每個(gè)單詞去采樣改化。那么每篇文章的每個(gè)單詞 就會得到 一個(gè)他是屬于哪個(gè)TOPIC的樣本數(shù)據(jù)。比如上面的例子文章 （ball, galaxy）枉昏， 我們可以運(yùn)用p(z|w) 采樣得到 (ball:topic1, galaxy: topic 3)
那么我們有M篇文章陈肛，假設(shè)每篇文章N個(gè)詞。我們就可以統(tǒng)計(jì)出M*N個(gè)詞兄裂，哪些屬于TOPIC1燥爷， 哪些屬于TOPIC2， 哪些屬于TOPIC3懦窘。然后看屬于TOPIC1的詞前翎，分別BALL占幾個(gè)，GALAXY占幾個(gè)畅涂。 這樣我們就有了TOPIC 到 詞的多項(xiàng)分布采樣數(shù)據(jù)了港华。

同時(shí)我們可以分別統(tǒng)計(jì)每篇文章下，這個(gè)詞是TOPIC1午衰，那么就TOPIC1 的計(jì)數(shù)+1. 下一個(gè)詞是TOPIC2立宜，那么TOPIC2的計(jì)數(shù)+1. 最后我們可以得到這篇文章如果100個(gè)單詞，那么會統(tǒng)計(jì)得到（topic1: 90, topic2: 5, topic3: 5）的文章到主題的單詞分布臊岸。

而這2個(gè)采樣數(shù)據(jù)正是我們需要的橙数。
所以核心是要知道p(z | w) 。

如果我們直接去求 我們會怎么做帅戒？
首先我們把前面的公式

image.png

用向量來表示得到如下公式

image.png

p(w|z) 表示了我們前一個(gè)公式的第二個(gè)旋鈕和第二個(gè)齒輪做的事情。就是根據(jù)beta的狄里克雷分布瞎访，隨機(jī)出phi的多項(xiàng)分布腻贰，并且在知道Z是哪個(gè)的條件下，隨機(jī)出word是哪個(gè)的過程扒秸。

那么現(xiàn)在我們的式子里有了p(w|z) 我們想要的p(z | w). 那么我就想到了貝葉斯公式播演，可以把已知的和要求的做一個(gè)交換冀瓦。

image.png

image.png

上面應(yīng)該還是比較好理解的肄程，文章中某一個(gè)單詞的概率，因?yàn)镵個(gè)TOPIC 都有一定概率會生成這個(gè)單詞选浑。所以需要分別求每個(gè)TOPIC下生成這個(gè)的概率蓝厌，然后求和。
我們假設(shè)第一個(gè)單詞的概率（p1k1 + p1k2 + p1k3）古徒。 那么如果一共這篇文章有2個(gè)單詞拓提，就會有
（p1k1 + p1k2 + p1k3）* (p2k1 + p2k2 + p2k3)
我們展開就會有9個(gè)項(xiàng)。3^2
如果有N個(gè)單詞的話隧膘，K個(gè)TOPIC的話代态。就可以有K^N 的項(xiàng)。這個(gè)離散狀態(tài)空間太大疹吃，無法列舉蹦疑。

我們試了正統(tǒng)方法去算p（z|w）遇到困難。
所以我們需要用一個(gè)別的方法來采樣到p(z|w)
下面就是MCMC方法的介紹萨驶，以及MCMC算法家族中吉布斯采樣具體解決這個(gè)問題的思路歉摧。

3.2 MCMC 的工作原理

這塊部分其他介紹LDA的文章，都有比較好的說明腔呜。我就簡要講解最核心的部分叁温，不詳細(xì)展開。如果需要深入掌握這塊的核畴，可以參考別的文章膝但。

整個(gè)邏輯思路是這樣的。首先在計(jì)算機(jī)中我們用線性同余法可以很容易生成均勻分布的隨機(jī)數(shù)谤草「可是有一些分布，我們沒有方法可以很好的采樣（從均勻分布變換過去比較困難）咖刃。這時(shí)就需要用復(fù)雜的隨機(jī)模擬的方法來生成樣本泳炉。

其核心1就是馬爾科夫鏈，這種鏈的特性就是基于當(dāng)前狀態(tài)嚎杨，就可以推出下一狀態(tài)。并且大多數(shù)馬氏鏈在滿足平穩(wěn)條件下氧腰，最終是會收斂的枫浙。

核心2就是蒙特卡羅方法刨肃，我們通過自己給出一個(gè)建議概率轉(zhuǎn)移矩陣，使得這個(gè)建議概率產(chǎn)生的新狀態(tài)在根據(jù)接受率有些采納箩帚，有些放棄真友。來獲得符合原來分布的樣本集。
這就是MCMC的核心2個(gè)概念紧帕。

3.2.1 理解馬氏鏈會收斂的性質(zhì)（理解可跳過）

image.png

上面的矩陣就是概率轉(zhuǎn)移矩陣盔然。這個(gè)矩陣的作用就是我們當(dāng)前在某一個(gè)狀態(tài)。比如我現(xiàn)在在I狀態(tài)是嗜，我們可以找到矩陣中對應(yīng)的那一行愈案，就能知道我下一個(gè)狀態(tài)是I,J,K,L的概率分別是多少。根據(jù)這個(gè)概率鹅搪，我就可以去到下一個(gè)狀態(tài)站绪。

有這么一種情況，我的位置向量在若干次轉(zhuǎn)移后會達(dá)到一個(gè)穩(wěn)定狀態(tài)丽柿。
舉個(gè)例子恢准，我最開始在一個(gè)（25%是i, 25% 是j, 25%是k, 25%是l）經(jīng)過前面的矩陣轉(zhuǎn)換。我下一個(gè)狀態(tài)可能是（30%甫题， 15%馁筐，25%， 30%）坠非。在經(jīng)過前面的矩陣眯漩，會繼續(xù)變化自己的狀態(tài)÷槎ィ可是經(jīng)過了足夠多次之后赦抖，我的狀態(tài)會收斂到一個(gè)固定的值上。

下面給一段python代碼辅肾，小伙伴可以自己做下實(shí)驗(yàn)队萤，如果轉(zhuǎn)移矩陣不變。你們隨便初始化初始狀態(tài)矫钓，最后的結(jié)果都會是個(gè)定值要尔。

import numpy as np

input = np.array([[0.05, 0.70, 0.25]])
matrix = np.matrix([[0.9, 0.075, 0.025],
                    [0.15, 0.8, 0.05],
                    [0.25, 0.25, 0.5]])
for i in range(100):
    input = np.dot(input, matrix)
    print(input)

....
[[ 0.625 0.3125 0.0625]]
[[ 0.625 0.3125 0.0625]]
[[ 0.625 0.3125 0.0625]]
[[ 0.625 0.3125 0.0625]]
[[ 0.625 0.3125 0.0625]]
[[ 0.625 0.3125 0.0625]]
[[ 0.625 0.3125 0.0625]]

那么insight就是，轉(zhuǎn)移矩陣決定最后的概率分布和初始值無關(guān)新娜。那么如果一個(gè)概率分布不好求赵辕，我們可以通過構(gòu)造一個(gè)馬爾科夫鏈的轉(zhuǎn)移矩陣，使得最后收斂的這個(gè)概率分布正好等于我們想求的這個(gè)概率分布概龄，即可还惠。
重點(diǎn)，整體概率分布不好求私杜，但是具體給定某個(gè)狀態(tài)蚕键，我們還是要能夠求這個(gè)狀態(tài)的產(chǎn)生的概率

其實(shí)到這里救欧，我們依然會問，原始概率分布不好求锣光，那收斂到這個(gè)概率分布的轉(zhuǎn)移矩陣就好求了嗎笆怠？

3.2.2 理解蒙特卡羅思想（理解可跳過）

這里需要用到1個(gè)技巧，這就是第二個(gè)MC誊爹，也就是蒙特卡羅算法蹬刷。我們通過自己給出一個(gè)建議轉(zhuǎn)移矩陣。再根據(jù)原來的當(dāng)前狀態(tài)的概率值和自己的建議的轉(zhuǎn)移概率值频丘，可以算出一個(gè)接受率办成。就可以知道我們隨機(jī)出來的新狀態(tài)，是保留還是放棄椎镣。

而蒙特卡羅的思想我們可以舉個(gè)例子诈火，我們想求圓周率pi的值是多少。直接求困難的話状答，我們可以用均勻分布隨機(jī)出很多在范圍【0,1】【0,1】的（x,y）冷守。只對滿足（x^2 + y^2 <= 1）的點(diǎn)留下來。最后比如我們丟了1W個(gè)點(diǎn)惊科。計(jì)數(shù)的只有X個(gè)
那么4 * X / 10000就是PI的值拍摇。因?yàn)橛?jì)數(shù)的是半徑為1的1/4圓的面積 ，總量是面積變長為1的正方形面積

python代碼如下：

import random

cnt = 0
samples = 1000000
for i in range(samples):
    x = random.random()
    y = random.random()
    if x*x + y*y <= 1:
        cnt += 1
print(4 * cnt / samples)

3.2.3 理解MH算法（理解可跳過）

那么我們要構(gòu)造的就是一個(gè)采樣算法馆截。那么構(gòu)造這樣的算法的洞察在于只有滿足一個(gè)叫detailed balance 的馬爾科夫鏈才會收斂到一個(gè)平穩(wěn)分布充活。而我們需要這個(gè)平穩(wěn)分布。所以我們需要去構(gòu)造蜡娶。MH是構(gòu)造平穩(wěn)分布的其中一種算法

image.png

[[ 0.625 0.3125 0.0625]]
所以pi1 = 0.625, pi2 = 0.3125, pi3 = 0.0625
轉(zhuǎn)移矩陣是
[[0.9, 0.075, 0.025],
[0.15, 0.8, 0.05],
[0.25, 0.25, 0.5]]

那么假設(shè)i = 1, j = 2（1開始計(jì)數(shù)）

image.png

對于任意的I,J。如果detailed balance 都滿足窖张。我們就可以得到一個(gè)式子

image.png

那么我們先從一個(gè)一般的算法講起叫MH算法（吉布斯采樣是它的特例）

MH 算法的目的：是根據(jù)一個(gè)需求的(desired distribution)概率分布 P(x)生成一系列樣本狀態(tài)點(diǎn)（因此睦霎，這個(gè)算法可以生成任意的概率分布）梢卸。為了達(dá)到這個(gè)目的，該算法使用馬爾可夫過程去到達(dá)一個(gè)平穩(wěn)分布pi(x), 以使pi(x)=P(x)副女。
根據(jù)下面的做法蛤高，可以滿足平穩(wěn)分布（式3， 數(shù)學(xué)推導(dǎo)在最后，現(xiàn)在只要記住結(jié)論理解過程即可）

image.png

還有一個(gè)容易混淆的點(diǎn)，上文提及的丟棄懊昨。不是說完全放棄這次迭代窄潭。而是如果沒有在接受率范圍里。就使用舊的樣本作為當(dāng)前這次采樣的樣本酵颁。所以迭代1000次嫉你，接受了600次，還是會有1000個(gè)樣本

3.2.4 理解吉布斯采樣

而吉布斯算法躏惋，就是這個(gè)接受率變?yōu)榱擞肋h(yuǎn)的1. 也就是不再會放棄一些新樣本幽污。當(dāng)然這就要求這個(gè)算法的建議概率為原概率分布的完全條件概率。

image.png

其實(shí)用大白話說很容易簿姨，如果一個(gè)狀態(tài)是N維的距误，這個(gè)完全條件概率，就是只看其中1維扁位，在別的維度都已知的條件下准潭，這個(gè)我看中的這維的概率分布會是怎么樣。
然后再這個(gè)基礎(chǔ)上對這維的數(shù)據(jù)進(jìn)行采樣域仇，并且更新刑然。
把所有維度都求一遍，才得到一個(gè)新的樣本
FOR里面所有行即采樣一個(gè)樣本要做的事情殉簸。N為采樣N個(gè)樣本
算法如下

image.png

對我們的影響就是：

1. 不再需要自己給一個(gè)建議概率g
2. 不再需要計(jì)算原有的概率分布下的概率值P（也是平穩(wěn)分布pi）
3. 需要我們可以計(jì)算完全條件概率分布值般卑，并采樣武鲁。

3.2.5 吉布斯采樣完全條件概率如何計(jì)算

你理解完這段，你就把前面所有知識點(diǎn)串起來了 這里面涉及很多數(shù)學(xué)推導(dǎo)蝠检，我都只給結(jié)論沐鼠。我們先拋開細(xì)節(jié)，關(guān)于結(jié)論和結(jié)論之間的邏輯關(guān)系。最后在回到每個(gè)結(jié)論的推導(dǎo)證明上饲梭。

步驟1. 找到LDA下的要求的完全條件概率乘盖。

我們的初心是得到 p(z | w) 也就是在已知單詞的情況下知道每個(gè)單詞的主題分布是什么。

那么我們繼續(xù)思考憔涉，LDA下的維度是什么订框。因?yàn)橥耆珬l件概率，就是在維度上推演的兜叨。
我們可以發(fā)現(xiàn)M篇文章穿扳，每篇文章的N個(gè)單詞。都是可以變換TOPIC的国旷。所以這里的維度是M * N維矛物。

那么基于上述2點(diǎn)， 我們

image.png

我們知道p(w, z) 我們根據(jù)LDA模型的定義是可以算出來的跪但。那么這個(gè)概率應(yīng)該怎么算呢履羞？

步驟二， 想辦法計(jì)算這個(gè)概率

（式5屡久， 數(shù)學(xué)推導(dǎo)在最后忆首， 只給結(jié)論之間的關(guān)系，具體證明在最后）
給出狄里克雷的概率密度函數(shù)

image.png

基于上述概率分布涂身。我們可以得到p(w,z)形式如下：（式6雄卷， 數(shù)學(xué)推導(dǎo)在最后）

image.png

有了p(w, z)的公式，我們就可以 擴(kuò)展到 p(w, z（除去某一維度））
維度是某一篇文章下的某一個(gè)單詞蛤售。

最后可以根據(jù)p(w,z)推得

image.png

那么這個(gè)概率已經(jīng)變成我們可以計(jì)算的形式了丁鹉。
因?yàn)橄闰?yàn)狄里克雷分布的alpha和beta不影響最后模型的結(jié)果，所以我們可以把它們?nèi)咳橥粋€(gè)值悴能，這樣我們就可以把上面的式子轉(zhuǎn)換為下面的式子
拿beta分布來舉例揣钦，其實(shí)就是alpha 和 beta取等值。在k維狄里克雷分布下漠酿，就是alpha 1 ~k 都是一個(gè)值冯凹。這種所有超參數(shù)都為同一個(gè)值稱為對稱超參數(shù)

image.png

步驟三， 有了這個(gè)概率值炒嘲，我們就可以采樣了宇姚。

因?yàn)轳R爾科夫鏈和初始狀態(tài)是什么無關(guān)，所以我們可以用均勻分布隨機(jī)給所有維度一個(gè)初始的狀態(tài)（topic 號）

那么現(xiàn)在采樣的步驟就是夫凸，比如我們在處理文章1的第一個(gè)單詞浑劳，假設(shè)是ball。

我們先要計(jì)算他在topic1下的概率值是多少夭拌。
那就是看除了它之外的其他所有的單詞的4個(gè)統(tǒng)計(jì)量魔熏。

分別是乘號左邊的2個(gè)衷咽。
就是在全文本集中，topic1有多少個(gè)ball詞是分子蒜绽。 和 topic1下一共有多少個(gè)詞是分母镶骗。

乘號右邊的2個(gè)。
在第一篇文章中躲雅，topic1的詞有多少個(gè)鼎姊。 第一篇文章一共有多少個(gè)詞。

那么整個(gè)采樣過程如下(java代碼）：

private int multiSamplingAcc(int m, int n) {
        // do multinomial sampling via cumulative method:
        double[] p = new double[K];
        for (int k = 0; k < K; k++) {
            // nw[w][k] 第k個(gè)topic下有多少個(gè)w號單詞吏夯， nwsum[k] 第k個(gè)topic下一共有多少個(gè)單詞
            // nd[m][k] 第m篇文章里有多少個(gè)屬于topic k的單詞此蜈， ndsum[m] 第m篇文章里有多少個(gè)單詞
            p[k] = (nw[docs[m][n]][k] + beta) / (nwsum[k] + V * beta)
                    * (nd[m][k] + alpha) / (ndsum[m] + K * alpha);
        }
        for (int k = 1; k < p.length; k++) {
            p[k] += p[k - 1];
        }
        double u = Math.random() * p[K - 1];
        int k_new;
        for (k_new = 0; k_new < p.length; k_new++) {
            if (u <= p[k_new])
                break;
        }
        z[m][n] = k_new;
        return k_new;
    }

步驟四即横，那么我們就可以一直采樣噪生，直到收斂。得到我們想要的輸出

還記得頂層我們程序需要的輸出嗎东囚？

每一篇都會給一個(gè)文章到TOPIC的概率分布跺嗽。 每一個(gè)TOPIC，都會告訴你這個(gè)TOPIC下每一個(gè)單詞的概率分布页藻。

那么吉布斯采樣在收斂完以后桨嫁，我們手里有的是。一套數(shù)據(jù)樣本份帐。這個(gè)樣本是每個(gè)詞分別屬于哪個(gè)topic璃吧。

我們在來看一下，一開始我們有的是文章-》單詞的樣本數(shù)據(jù)废境。

我們希望有的是文章-》主題畜挨，和 主題-》單詞的樣本數(shù)據(jù)。因?yàn)槲覀兿Ｍ聪蚋碌男o的定義是文章-》主題噩凹，主題-》單詞巴元。

那么現(xiàn)在每個(gè)詞都有一個(gè)主題。所以我們可以統(tǒng)計(jì)這個(gè)文章下驮宴，每個(gè)主題所占的比例（根據(jù)文章內(nèi)的詞所占的比例）逮刨； 已經(jīng)在所有文集里，這個(gè)主題下的所有單詞（每個(gè)單詞所占的比例）那么我們就得到了主題到單詞的多項(xiàng)分布堵泽。

這樣我們就得到了我們的輸出修己。沒錯(cuò)此時(shí)我們已經(jīng)通過這2個(gè)輸出，作為似然函數(shù)的數(shù)據(jù)值迎罗，把我們最開始的LDA模型的2個(gè)狄里克雷分布的旋鈕給訓(xùn)練好了睬愤。我們不需要更新alpha和beta的數(shù)值，因?yàn)橛邢闰?yàn)和數(shù)據(jù)也是等價(jià)的佳谦，并且我們不需要模型真的去隨機(jī)產(chǎn)生文章戴涝。而是要對新的文章去做預(yù)測。

上述是一個(gè)直觀的理解。其實(shí)這背后也有數(shù)學(xué)上的推導(dǎo)和證明啥刻。我們把這個(gè)稱為（式7奸鸯，在最后統(tǒng)一處理）

那么我們來看一下這塊對應(yīng)的代碼

// nw[w][k] 第k個(gè)topic下有多少個(gè)w號單詞， nwsum[k] 第k個(gè)topic下一共有多少個(gè)單詞
// nd[m][k] 第m篇文章里有多少個(gè)屬于topic k的單詞可帽， ndsum[m] 第m篇文章里有多少個(gè)單詞
// V為distinct的文集單詞集合數(shù)量
// K為指定的主題數(shù)量
// M為文章數(shù)量

// 得到每個(gè)主題下娄涩，每個(gè)詞的概率
    public double[][] getPhi() {
        double[][] phi = new double[K][V];
        for (int k = 0; k < K; k++) {
            for (int t = 0; t < V; t++) {
                phi[k][t] = (nw[t][k] + beta) / (nwsum[k] + V * beta);
            }
        }
        return phi;
    }

// 得到每個(gè)文章下，每個(gè)主題的概率分布
    public double[][] getTheta() {
        double[][] theta = new double[M][K];
        for (int m = 0; m < M; m++) {
            for (int k = 0; k < K; k++) {
                theta[m][k] = (nd[m][k] + alpha) / (ndsum[m] + K * alpha);
            }
        }
        return theta;
    }

3.2.6 回顧一下全貌

我們現(xiàn)在已經(jīng)掌握可以編寫LDA代碼的全部邏輯鏈條和組件映跟。我們來回顧一下蓄拣。

1. 最頂層上，我們明確了輸入是list<list<words>>, K; 輸出是double[K][V] phi (代表主題到詞的概率分布）
   double[M][K] theta(代表文章到主題的分布 ）
2. 模型的好壞 由2個(gè)旋鈕決定努隙，分別是文章到主題的alpha的狄里克雷分布球恤，和主題到詞的beta狄里克雷分布；隨后會由2組齒輪去產(chǎn)生單詞荸镊，分別是文章到主題的多項(xiàng)分布咽斧，主題到詞的多項(xiàng)分布）
3. 根據(jù)上述4個(gè)組件，我們可以算出在已知p(w, z|a,b) 的概率躬存；我們需要的是p(z |w,a,b) 的概率张惹， 因?yàn)槲覀円烂總€(gè)詞具體的主題才可以采樣得到概率期望值作為頂層模型的輸出
4. 在直接計(jì)算的時(shí)候，遇到了p(w) 項(xiàng)數(shù)指數(shù)爆炸不好計(jì)算岭洲。所以我們曲線救國宛逗，引入了MCMC的gibbs 采樣算法
5. gibbs采樣算法，需要用到完全條件概率盾剩。這個(gè)可以由p(w,z) 去推導(dǎo)出來雷激。
6. 有了完全條件概率的值的計(jì)算方法，我們就可以知道當(dāng)前每一個(gè)維度彪腔，轉(zhuǎn)換到下一個(gè)狀態(tài)侥锦，這個(gè)維度的新的值是基于什么概率轉(zhuǎn)移方程去做的。那么我們就能采樣到一個(gè)新的值對這個(gè)維度來說德挣。
7. 這樣依次對每一個(gè)維度恭垦，基于其他維度的信息去采樣。做完之后格嗅，我們就得到一個(gè)全新的采樣樣本番挺。這樣做直到gibbs收斂。
8. gibbs收斂后屯掖，那么現(xiàn)在的每個(gè)詞上的TOPIC分布玄柏，就滿足了和文章輸入最密切相關(guān)的聯(lián)合概率分布 p(w, z). 那么我們可以根據(jù)狄里克雷概率的期望公式拿到頂層設(shè)計(jì)的2個(gè)輸出。

4. 我們?nèi)绾位谶@些原理去實(shí)現(xiàn)自己的LDA

4.1 根據(jù)上面的算法贴铜， 定義出我們需要的數(shù)據(jù)結(jié)構(gòu)粪摘。

image.png

public class GibbsSampler {
    int K, V, M;
    double beta = 1, alpha = 1;
    // nw[:][k] that of term observation counts for topic k
    // topic–term count
    int[][] nw;
    // topic–term sum
    int[] nwsum;

    // document–topic count
    int[][] nd;
    // document–topic sum
    int[] ndsum;

    // Z(m,n) => m doc, n word 's topic idx
    int[][] z;
    // docs[m][n]
    int[][] docs;

    Random random = new Random();

4.2 寫構(gòu)造函數(shù)處理輸入

這個(gè)時(shí)候我們要傳進(jìn)來輸入了瀑晒。這里需要對每個(gè)單詞做一個(gè)唯一編號。也就是給每個(gè)word一個(gè)遞增id徘意，來加速運(yùn)算苔悦。（否則要用map<string>）， 在外層椎咧，我們需要得到這個(gè)編號ID的最大值+1之后玖详，就是distinct 單詞集合的數(shù)量，作為V傳進(jìn)來

public GibbsSampler(int[][] docs, int k, int v) {
this.docs = docs;
K = k;
M = docs.length;
V = v;
}

4.3 gibbs 采樣框架書寫

init 是隨機(jī)初始化勤讽，因?yàn)槲覀兊淖詈蟮慕Y(jié)果和初始的概率分布無關(guān)蟋座。所以隨機(jī)初始化，把幾個(gè)值統(tǒng)計(jì)一下就好了脚牍。

在gibbs采樣的每一個(gè)維度下向臀，都要先排除這個(gè)維度的數(shù)據(jù)統(tǒng)計(jì)值，這個(gè)時(shí)候需要調(diào)用
decrementCntAndSum

multiSamplingAcc 就是上文提到了完全條件概率具體值的計(jì)算和運(yùn)用莫矗。返回出一個(gè)這個(gè)維度下的新topic

在得到一個(gè)這個(gè)維度下的新的topic的時(shí)候飒硅，要把統(tǒng)計(jì)值更新回來。

偽代碼為

init()
for 一個(gè)足夠大使得gibbs可以收斂的迭代次數(shù)
  for 每一個(gè)維度
      排除當(dāng)前維度的統(tǒng)計(jì)信息
      用完全條件概率公式的計(jì)算方式去采樣得到一個(gè)新的topic
      更新當(dāng)前維度的topic和統(tǒng)計(jì)信息

把模型寫到磁盤里

public void gibbs()  {
        init();
        for (int i = 0; i < ITERATIONS; i++) {
            for (int m = 0; m < M; m++) {
                int N = docs[m].length;
                for (int n = 0; n < N; n++) {
                    int k = z[m][n];
                    decrementCntAndSum(m, k, n);
                    int k_new = multiSamplingAcc(m, n);
                    incrementCntAndSum(m, k_new, n);
                }
            }
        }
    }

4.4 實(shí)現(xiàn)輔助函數(shù)

private void init() {
        nw = new int[V][K];
        nwsum = new int[K];
        nd = new int[M][K];
        ndsum = new int[M];

        z = new int[M][];
        for (int m = 0; m < M; m++) {
            int N = docs[m].length;
            z[m] = new int[N];
            for (int n = 0; n < N; n++) {
                int k = random.nextInt(K);
                z[m][n] = k;
                incrementCntAndSum(m, k, n);
            }
        }
    }
    private void incrementCntAndSum(int m, int k, int n) {
        nd[m][k]++;
        ndsum[m]++;
        int t = docs[m][n];
        nw[t][k]++;
        nwsum[k]++;
    }
    private void decrementCntAndSum(int m, int k, int n) {
        nd[m][k]--;
        ndsum[m]--;
        int t = docs[m][n];
        nw[t][k]--;
        nwsum[k]--;
    }

4.5 實(shí)現(xiàn)核心完全條件概率函數(shù)采樣算法

private int multiSamplingAcc(int m, int n) {
        // do multinomial sampling via cumulative method:
        double[] p = new double[K];
        for (int k = 0; k < K; k++) {
            // nw[w][k] 第k個(gè)topic下有多少個(gè)w號單詞作谚， nwsum[k] 第k個(gè)topic下一共有多少個(gè)單詞
            // nd[m][k] 第m篇文章里有多少個(gè)屬于topic k的單詞， ndsum[m] 第m篇文章里有多少個(gè)單詞
            p[k] = (nw[docs[m][n]][k] + beta) / (nwsum[k] + V * beta)
                    * (nd[m][k] + alpha) / (ndsum[m] + K * alpha);
        }
        for (int k = 1; k < p.length; k++) {
            p[k] += p[k - 1];
        }
        double u = Math.random() * p[K - 1];
        int k_new;
        for (k_new = 0; k_new < p.length; k_new++) {
            if (u <= p[k_new])
                break;
        }
        z[m][n] = k_new;
        return k_new;
    }

4.6 計(jì)算輸出的2個(gè)數(shù)組

// 得到每個(gè)主題下庵芭，每個(gè)詞的概率
    public double[][] getPhi() {
        double[][] phi = new double[K][V];
        for (int k = 0; k < K; k++) {
            for (int t = 0; t < V; t++) {
                phi[k][t] = (nw[t][k] + beta) / (nwsum[k] + V * beta);
            }
        }
        return phi;
    }

    // 得到每個(gè)文章下妹懒，每個(gè)主題的概率分布
    public double[][] getTheta() {
        double[][] theta = new double[M][K];
        for (int m = 0; m < M; m++) {
            for (int k = 0; k < K; k++) {
                theta[m][k] = (nd[m][k] + alpha) / (ndsum[m] + K * alpha);
            }
        }
        return theta;
    }

所以整個(gè)訓(xùn)練的代碼 也就100多行，代碼量其實(shí)不大双吆。

4.7 寫預(yù)測文章的主題分布代碼

預(yù)測其實(shí)就是還是在同樣的LDA模型下做預(yù)測眨唬。那么我們就可以假定原來訓(xùn)練出來的參數(shù)是最佳的。也就是主題到詞的概率好乐。所以我們需要的輸出是針對這篇文章的文章到主題的概率分布匾竿。
而我們有的輸入，是這篇文章的文章到詞蔚万。
那么我們就單獨(dú)針對這篇文章用gibbs采樣岭妖，來獲得這篇文章的p(z|w). 其中算法的一個(gè)乘法因子（文章到主題的）是根據(jù)文章自己的詞語分布特性重新計(jì)算的。另一個(gè)因子（主題到詞的）反璃，我們就用模型算好的phi數(shù)組即可昵慌。

那么模型收斂以后，我們就可以用樣本數(shù)據(jù)求期望獲得對文章到主題的分布了淮蜈。
也就是在最上面的我們畫出來的狄里克雷三角形中斋攀，再找到一個(gè)最合適的黃色小球在三角形中的位置。

public class LdaPredicator {
    public static final int ITERATIONS = 1000;

    public static double[] predict(double[][] phi, double alpha, double beta, int[] doc) {
        int K = phi.length;
        int V = phi[0].length;

        int[][] nw = new int[V][K];
        int[] nd = new int[K];
        int[] nwsum = new int[K];
        int ndsum = 0;

        // init
        int N = doc.length;
        int[] z = new int[N];   // z_i := 1到K之間的值梧田，表示馬氏鏈的初始狀態(tài)
        for (int n = 0; n < N; n++) {
            int k = (int) (Math.random() * K);
            z[n] = k;
            nw[doc[n]][k]++;
            nd[k]++;
            nwsum[k]++;
        }
        ndsum = N;

        for (int i = 0; i < ITERATIONS; i++) {
            for (int n = 0; n < N; n++) {

                int topic = z[n];
                nw[doc[n]][topic]--;
                nd[topic]--;
                nwsum[topic]--;
                ndsum--;

                // do multinomial sampling via cumulative method: 通過多項(xiàng)式方法采樣多項(xiàng)式分布
                double[] p = new double[K];
                for (int k = 0; k < K; k++) {
                    p[k] = phi[k][doc[n]]
                            * (nd[k] + alpha) / (ndsum + K * alpha);
                }
                // cumulate multinomial parameters  累加多項(xiàng)式分布的參數(shù)
                for (int k = 1; k < p.length; k++) {
                    p[k] += p[k - 1];
                }
                // scaled sample because of unnormalised p[] 正則化
                double u = Math.random() * p[K - 1];
                for (topic = 0; topic < p.length; topic++) {
                    if (u < p[topic])
                        break;
                }

                // add newly estimated z_i to count variables   將重新估計(jì)的該詞語加入計(jì)數(shù)器
                nw[doc[n]][topic]++;
                nd[topic]++;
                nwsum[topic]++;
                ndsum++;
                z[n] = topic;
            }
        }

        double[] theta = new double[K];

        for (int k = 0; k < K; k++) {
            theta[k] = (nd[k] + alpha) / (ndsum + K * alpha);
        }
        return theta;
    }
}

最后寫一個(gè)main函數(shù)淳蔼，大家自己來玩一下吧

import java.util.Arrays;
import java.util.HashMap;
import java.util.Map;

public class LdaMain {
    private int[][] docs;
    private Map<String, Integer> w2id = new HashMap<>();
    private Map<Integer,String> id2w = new HashMap<>();
    private GibbsSampler gibbsSampler;
    public LdaMain(String[] oriDocs, int k) {
        int m = oriDocs.length;
        docs = new int[m][];
        for (int i = 0; i < m; i++) {
            String[] ws = oriDocs[i].split(" ");
            docs[i] = new int[ws.length];
            int j = 0;
            for (String w : ws) {
                w2id.putIfAbsent(w, w2id.size());
                id2w.putIfAbsent(w2id.get(w), w);
                docs[i][j++] = w2id.get(w);
            }
        }
        gibbsSampler = new GibbsSampler(docs, k, w2id.size());
    }

    public static void main(String[] args) {
        String[] docs = {
                "ball ball ball planet galaxy ball ball ball ball ball ball ball ball ball ball",
                "referendum planet planet referendum referendum referendum referendum referendum referendum referendum referendum referendum referendum referendum referendum",
                "planet planet galaxy planet ball planet galaxy planet galaxy planet galaxy planet galaxy planet planet planet planet",
                "planet galaxy referendum planet ball planet galaxy planet galaxy planet galaxy planet galaxy planet galaxy planet planet planet planet"};
        LdaMain main = new LdaMain(docs, 3);
        main.train();
        main.predict("planet galaxy planet galaxy ball planet galaxy planet galaxy planet galaxy planet galaxy planet planet planet planet planet planet");
    }

    private void predict(String doc) {
        double[][] phi = gibbsSampler.getPhi();
        String[] ws = doc.split(" ");
        int[] ddoc = new int[ws.length];
        for (int i = 0; i < ws.length; i++) {
            ddoc[i] = w2id.get(ws[i]);
        }
        printTopic2Word(phi);
        double[] theta = LdaPredicator.predict(phi, gibbsSampler.alpha, gibbsSampler.beta, ddoc);
        System.out.println(Arrays.toString(theta));
    }

    private void printTopic2Word(double[][] phi) {
        for (int i = 0; i < phi.length; i++) {
            System.out.println("topic " + i + ": ");
            for (int j = 0; j < w2id.size(); j++) {
                System.out.print(id2w.get(j) + ":" + phi[i][j] + " ");
            }
            System.out.println();
        }
    }

    private void train() {
        gibbsSampler.gibbs();
    }
}

5. 背后的數(shù)學(xué)詳細(xì)推導(dǎo) （每個(gè)數(shù)學(xué)結(jié)論是怎么得來的）

這篇文章的重點(diǎn)我已經(jīng)講完了侧蘸。我們需要數(shù)學(xué)推導(dǎo)的地方，我都留了（式X）的字樣鹉梨， 應(yīng)該是有7個(gè)闺魏，需要數(shù)學(xué)推導(dǎo)一下。
我們把這塊俯画，當(dāng)做讀后作業(yè)析桥，交給讀者自己去尋找。所有證明都可以在LDA漫游指南里找到艰垂。
當(dāng)然有余力的小伙伴泡仗，還可以看看LDA數(shù)學(xué)八卦。
看懂?dāng)?shù)學(xué)推導(dǎo)需要一定的微積分基礎(chǔ)和概率統(tǒng)計(jì)基礎(chǔ)猜憎。小伙伴們加油娩怎！

感謝

看到最后的小伙伴，非常不容易胰柑。這個(gè)模型的原理還是有不小的門檻截亦。我也不確定哪里有筆誤，或者理解不準(zhǔn)確的地方柬讨。
不過把整個(gè)邏輯鏈條理清楚了還是很暢快的效拭。
你們投入了時(shí)間希望你們也可以收獲到了這個(gè)LDA模型的魅力。如果有任何疑問剃斧，可以在文章下留言碱呼。thank you for your time and suggestion~