【NLP】LDA主题模型详解

前言

本人在开始学习LDA主题模型时，一开始看上去感觉挺简单、挺容易理解的。但是，要想真正理解它不是那么容易，大概是因为里面的数学知识点比较多，让一些数学基础不好的人望而生畏。不过好在，如果不想深究其内部原理的，只想把握LDA主题模型的主题思想，相对不是很难，大家耐心看完就会明白。

一、什么是主题模型？

主题模型(Topic Model)是用来在一系列文档中发现 抽象主题 的一种统计模型。直观来讲，如果一篇文章有一个中心思想，那么一定存在一些特定词语会出现的比较频繁。

比方说，如果现在一篇文章是在讲苹果公司的，那么“乔布斯”和“IPhone”等词语出现的频率会更高一些；如果现在一篇文章是在描述微软公司的，那么“Windows”和“Microsoft”等词语出现的频率会更高一些；
但真实情况下，一篇文章中通常包含多种主题，而且每个主题所占的比例各不相同，比如一篇文章中10%和苹果公司有关，90%和微软公司有关，那么和微软有关的关键字出现的次数应该是苹果关键字出现次数的9倍。

主题模型 就是一种自动分析每个文档，统计文档中的词语，根据统计的信息判断当前文档包含哪些主题以及各个主题所占比例各为多少。

主题模型 是对文本中 隐含主题 的一种建模方法，每个主题其实是词表上单词的概率分布；主题模型是一种生成模型，一篇文章中每个词都是通过“以一定概率选择某个主题，并从这个主题中以一定概率选择某个词语”这样一个过程得到的。

这个过程好比我们在写文章：一般我们写文章的时候都是先给自己的文章定一个主题，然后我们会根据这个主题，用和这个主题有关的词语进行组合最后形成文章。

LDA主题模型之所以为主题模型，自然，它的主体思想便是上面的思想，下面我们先来介绍LDA 主题模型到底是怎么回事。

二、`LDA` 主题模型初步认识

LDA 其实有两种含义，一种是线性判别分析（Linear Discriminant Analysis），一般用来为后续的分类问题做数据的降维处理；另一种便是我们本篇重点介绍的主题模型，它称作为：隐含狄利克雷分布（Latent Dirichlet Allocation，简称LDA），是一种概率主题模型。

LDA由Blei, David M.、Ng, Andrew Y.、Jordan于2003年提出，是一种主题模型，它可以将文档集中每篇文档的主题以概率分布的形式给出，从而通过分析一些文档抽取出它们的主题（分布）出来后，便可以根据主题（分布）进行主题聚类或文本分类。同时，它是一种典型的词袋模型，即一篇文档是由一组词构成，词与词之间没有先后顺序的关系。此外，一篇文档可以包含多个主题，文档中每一个词都由其中的一个主题生成。

2.1 我们人类是怎么写出文章的?

我们引用三位作者在原始论文中的一个简单的例子。比如假设事先给定了这几个主题：Arts、Budgets、Children、Education，然后通过学习训练，获取每个主题 Topic 对应的词语。如下图所示：

然后以一定的概率选取上述某个主题，再以一定的概率选取那个主题下的某个单词，不断的重复这两步，最终生成如下图所示的一篇文章（其中不同颜色的词语分别对应上图中不同主题下的词）：

我以前上高中的时候，文科不好，每次语文考试，总是听到一些某某某作文满分，我就在想为什么我写不出来呢？现在我终于明白了，主要原因有两个：

人家写作文不会跑题；
人家肚子里的词语多；

这就导致，人家在写作文的时候，根据给定的主题，会有很多的词语来创造句子。而我，虽然也没跑题，但是肚子里能用来表达的词太少了，写个作文总是东拼西凑，自然不比人家的。

所以，我们猜想，那些写作比我们优秀的人是先确定自己将要写的文章的几个主题，然后围绕这几个主题遣词造句，写出优秀作文。

LDA 主题模型要干的事就是：根据给定的一篇文档，反推其主题分布。

通俗来说，可以假定认为人类是根据上述文档生成过程写成了各种各样的文章，现在我们想让计算机利用LDA干一件事：你计算机给我推测分析网络上各篇文章分别都写了些啥主题，且各篇文章中各个主题出现的概率大小（主题分布，后面还会提到）是啥。

如果就是按照上面的说法，我们都会觉得 LDA 主题模型也很普通啊，但是探其内部原理，看到那些数学知识点，很多人又觉得好难啊。

2.2 `LDA` 主题模型的主题思想

在 LDA 主题模型中，一篇文档生成的方式如下：

从狄利克雷分布 α\alphaα 中取样生成文档 iii 的主题分布 θi\theta_iθi；
从主题的多项式分布 θi\theta_iθi 中取样生成文档 iii 第 jjj 个词的主题 zi,jz_{i,j}zi,j；
从狄利克雷分布 β\betaβ 中取样生成主题 zi,jz_{i,j}zi,j 对应的词语分布 ϕzi,j\phi_{z_{i,j}}ϕzi,j；
从词语的多项式分布 ϕzi,j\phi_{z_{i,j}}ϕzi,j 中采样最终生成词语 wi,jw_{i,j}wi,j。

其中，类似Beta分布是二项式分布的共轭先验概率分布，而狄利克雷分布（Dirichlet分布）是多项式分布的共轭先验概率分布（这在后面具体介绍）。

LDA 主题模型的图模型结构如下图所示：

以上6句话便是整个 LDA 的主体思想！但是里面重复出现了二项分布、多项式分布、beta分布、狄利克雷分布、共轭先验概率分布、采样，那么，这些都是什么呢？接下来我们便逐一介绍。

三、相关数学知识点

此部分主要介绍LDA主题模型 相关的数学知识点，二项分布、多项式分布、beta分布、狄利克雷分布、共轭先验概率分布。如果只想了解 LDA主题模型 可以跳过此部分直接看第四部分，如果看完后面的又想了解相关数学知识，可以回头再看。

3.1 概率知识回顾

1）先验概率：

在事情尚未发生前，对该事件发生概率的估计。利用过去历史资料计算出来得到的先验概率叫做客观先验概率；凭主观经验来判断而得到的先验概率叫做主观先验概率。

2）后验概率：

通过调查或其它方式获取新的附加信息，利用贝叶斯公式对先验概率进行修正后，而得到的概率。

3）似然函数：

给定模型参数 θ\thetaθ 的条件下，样本数据服从这一概率模型的相似程度。

其中，θ1,θ2,...\theta_1,\theta_2,...θ1,θ2,...表示模型的未知参数，XXX表示样本

4）先验分布：

反映在进行统计试验之前根据其他有关参数知识得到的分布；也就是说在观测获取样本之前，人们对 θ\thetaθ 已经有一些知识，此时这个 θ\thetaθ 的分布函数为 H(θ)H(\theta)H(θ) ， θ\thetaθ 的密度函数为 h(θ)h(\theta)h(θ)，分别称为先验分布函数和先验密度函数，统称先验分布。

5）后验分布：

根据样本 XXX 的分布以及 θ\thetaθ 的先验分布 π(θ)\pi(\theta)π(θ)，使用概率论中求解条件概率的方式可以计算出来已知 XXX 的条件下， θ\thetaθ 的条件分布 π(θ∣X)\pi(\theta|X)π(θ∣X)。因为该分布是在获取样本 XXX 之后计算出来的，所以称为后验分布。

6）共轭分布：

如果先验分布和后验分布具有相同的形式，那么先验分布和似然函数被称为共轭分布。

3.2 二项分布

二项分布是从伯努利分布推导过来的。伯努利分布，又称两点分布或0-1分布，是一个离散型的随机分布，其中的随机变量只有两类取值，非正即负{+,−}\{+,-\}{+,−}。而二项分布即重复 nnn 次的伯努利试验，记为 X∽B(n,p)X \backsim B(n, p)X∽B(n,p) 。简言之，只做一次实验，是伯努利分布，重复做了 nnn 次，是二项分布。二项分布的概率密度函数为：

p(K=k)=Cnkpk(1−p)n−k，k=0,1,2,...,np(K=k) = \mathcal{C}_n^kp^k(1-p)^{n-k}，k = 0, 1, 2,...,np(K=k)=Cnkpk(1−p)n−k，k=0,1,2,...,n

其中 Cnk=n!k!(n−k)!\mathcal{C}_n^k = \frac{n!}{k!(n-k)!}Cnk=k!(n−k)!n!，是二项式系数。

通俗讲：我们做抛硬币试验，正面朝上的概率为 ppp，则反面朝上的概率为1−p1-p1−p，我们记XXX为正面朝上这一事件发生的次数

抛一次硬币，p(X=x)=px(1−p)1−x,x=0,1p(X=x) = p^x(1-p)^{1-x}, x=0,1p(X=x)=px(1−p)1−x,x=0,1，可以看作是x次正面朝上的概率（抛一次硬币也就只有0次和1次了），这就是伯努利分布

如果我们抛多次硬币，求k次正面朝上的概率，我们假设抛3次硬币，则2次朝上的情形有：110、101、011，也就是C32=3!2!(3−2)!=3∗2∗12∗1∗1=3\mathcal{C}_3^2=\frac{3!}{2!(3-2)!} = \frac{3*2*1}{2*1*1} = 3C32=2!(3−2)!3!=2∗1∗13∗2∗1=3 三种情形，于是抛三次硬币，两次正面的概率为：p(X=3)=3∗p2(1−p)1p(X=3) = 3 * p^2(1-p)^1p(X=3)=3∗p2(1−p)1，这就是二项分布

3.3 多项分布

多项分布式二项分布的推广，是指单次试验中随机变量的取值不再是0/1的，而是有多种离散值可能(1,2,3,...,k)(1,2,3,...,k)(1,2,3,...,k)。比如掷骰子，有6个面，N次试验结果服从K=6的多项分布。其中K个离散值的概率满足：

∑i=1kpi=1\sum_{i=1}^kp_i = 1i=1∑kpi=1

多项分布的概率密度为：
p(x1=m1,x2=m2,...,xk=mk,n,p1,p2,...,pk)p(x_1=m_1,x_2=m_2,...,x_k=m_k,n,p_1,p_2,...,p_k)p(x1=m1,x2=m2,...,xk=mk,n,p1,p2,...,pk)
=n!m1!m2!...mk!p1m1p2m2...pkmk= \frac{n!}{m_1!m_2!...m_k!}p_1^{m_1}p_2^{m_2}...p_k^{m_k}=m1!m2!...mk!n!p1m1p2m2...pkmk

其中，随机变量的取值有k个，X=(x1,x2,...,xk)X=(x_1,x_2,...,x_k)X=(x1,x2,...,xk)，并且x1=m1x_1=m_1x1=m1 表示随机变量取 x1x_1x1 时发生的次数，∑i=1kmi=n\sum_{i=1}^km_i = n∑i=1kmi=n 表示所有随机变量取值可能发生次数的总和为 nnn，即 nnn 次试验

通俗地将，我们掷六面骰子，出现每面的概率分别为：p1,p2,p3,..,p6p_1,p_2,p_3,..,p_6p1,p2,p3,..,p6，掷 nnn 次筛子,求 1 出现 m1m_1m1 次，2 出现 m2m_2m2 次，3 出现 m3m_3m3 次，4 出现 m4m_4m4 次，5 出现m5m_5m5 次，6 出现 m6m_6m6 次的概率，且总次数为 nnn：
p(x1=m1,x2=m2,...,x6=m6,n,p1,p2,...,p6)p(x_1=m_1,x_2=m_2,...,x_6=m_6,n,p_1,p_2,...,p_6)p(x1=m1,x2=m2,...,x6=m6,n,p1,p2,...,p6)
=n!m1!m2!...m6!p1m1p2m2...p6m6= \frac{n!}{m_1!m_2!...m_6!}p_1^{m_1}p_2^{m_2}...p_6^{m_6}=m1!m2!...m6!n!p1m1p2m2...p6m6

3.4 `gamma` 函数

gamma 函数定义如下：

Γ(x)=∫0+∞tx−1e−tdt\Gamma(x) = \int_0^{+\infty}t^{x-1}e^{-t}dtΓ(x)=∫0+∞tx−1e−tdt

gamma 函数可以当成是阶乘在实数集上的延拓，对于正整数n，具有如下性质：

Γ(n)=(n−1)!\Gamma(n) = (n-1)!Γ(n)=(n−1)!

我们这里只用到gamma函数的这个性质，所以我们不做详细介绍，只需要知道看到阶乘，可以联想到阶乘在实数上的推广的 gamma 函数

3.5 `Beta` 分布

相信大家学过统计学的都对正态分布、二项分布、 均匀分布 等等很熟悉了，但是却鲜少有人去介绍Beta分布的。

用一句话来说，Beta分布 可以看作一个概率的概率分布，当你不知道一个东西的具体概率是多少时，它可以给出了所有概率出现的可能性大小。

问题1：

假设随机变量X1,X2,...XnX_1,X_2,...X_nX1,X2,...Xn服从(0,1)(0,1)(0,1)的均匀分布，即X∽U(0,1)X \backsim U(0,1)X∽U(0,1)，把这个 nnn 个随机变量排序后得到顺序统计量 X(1),X(2),...,X(n)X_{(1)},X_{(2)},...,X_{(n)}X(1),X(2),...,X(n) ，请问 X(k)X_{(k)}X(k) 的分布是什么？

为了解决这个问题，可以尝试计算 X(k)X_{(k)}X(k) 落在区间 [x,x+Δx][x,x+\Delta x][x,x+Δx] 的概率，即求：

P(x≤X(k)≤x+Δx)=?P(x \leq X_{(k)} \leq x+ \Delta x) = ?P(x≤X(k)≤x+Δx)=?

首先，把 [0,1][0,1][0,1] 区间分成三段 [0,x)，[x,x+Δx]，(x+Δx,1][0,x)，[x,x+Δx]，(x+Δx,1][0,x)，[x,x+Δx]，(x+Δx,1]，然后考虑下简单的情形：即假设 nnn 个数中只有1个落在了区间 [x,x+Δx][x,x+Δx][x,x+Δx]内，由于这个区间内的数 X(k)X_{(k)}X(k) 是第k大的，所以[0,x)[0,x)[0,x)中应该有 k−1k−1k−1 个数，(x+Δx,1](x+Δx,1](x+Δx,1] 这个区间中应该有 n−kn−kn−k 个数。如下图所示：

从而问题转换为下述事件EEE：

E={X1∈[x,x+Δx];Xi∈[0,x),(i=2,3,...,k);Xj∈(x+Δx,1],(j=k+1,...,n)}E = \{X_1 \in [x,x+\Delta x];X_i \in [0,x), \ (i=2,3,...,k);X_j \in (x+Δx,1], \ (j = k+1,...,n)\}E={X1∈[x,x+Δx];Xi∈[0,x), (i=2,3,...,k);Xj∈(x+Δx,1], (j=k+1,...,n)}

对于上述事件EEE，有：

P(E)=∏i=1nP(Xi)=xk−1(1−x−Δx)n−kΔxP(E) = \prod_{i=1}^n P(X_i) = x^{k-1}(1-x-\Delta x)^{n-k}\Delta xP(E)=i=1∏nP(Xi)=xk−1(1−x−Δx)n−kΔx
=xk−1(1−x)n−kΔx+o(Δx)= x^{k-1}(1-x)^{n-k}\Delta x + o(\Delta x)=xk−1(1−x)n−kΔx+o(Δx)

其中，o(Δx)o(\Delta x)o(Δx) 表示 ΔxΔxΔx 的高阶无穷小。显然，由于不同的排列组合，即 nnn 个数中有一个落在 [x,x+Δx][x,x+Δx][x,x+Δx] 区间的有 nnn 种取法，余下n−1n−1n−1个数中有 k−1k−1k−1 个落在[0,x)[0,x)[0,x) 的有 Cn−1k−1\mathcal{C}_{n-1}^{k-1}Cn−1k−1 种组合，所以和事件EEE 等价的事件一共 nCn−1k−1n\mathcal{C}_{n-1}^{k-1}nCn−1k−1有个。

如果有2个数落在区间 [x,x+Δx][x,x+Δx][x,x+Δx] 呢？如下图所示：

类似于事件EEE，对于2个数落在区间 [x,x+Δx][x,x+Δx][x,x+Δx] 的事件 E′E'E′：

E′={X1,X2∈[x,x+Δx];Xi∈[0,x),(i=3,...,k);Xj∈(x+Δx,1],(j=k+1,...,n)}E' = \{X_1,X_2 \in [x,x+\Delta x];X_i \in [0,x), \ (i=3,...,k);X_j \in (x+Δx,1], \ (j = k+1,...,n)\}E′={X1,X2∈[x,x+Δx];Xi∈[0,x), (i=3,...,k);Xj∈(x+Δx,1], (j=k+1,...,n)}

则有：

P(E′)=xk−2(1−x−Δx)n−k(Δx)2=o(Δx)P(E') = x^{k-2}(1-x-\Delta x)^{n - k}(\Delta x)^2 = o(\Delta x)P(E′)=xk−2(1−x−Δx)n−k(Δx)2=o(Δx)

从上述的事件 EEE、事件 E′E'E′ 中，可以看出，只要落在 [x,x+Δx][x,x+Δx][x,x+Δx] 内的数字超过一个，则对应的事件的概率就是 o(Δx)o(Δx)o(Δx)。于是乎有：

P(x≤X(k)≤x+Δx)=nCn−1k−1P(E)+o(Δx)P(x \leq X_{(k)} \leq x+ \Delta x) = n\mathcal{C}_{n-1}^{k-1}P(E) + o(\Delta x)P(x≤X(k)≤x+Δx)=nCn−1k−1P(E)+o(Δx)
=nCn−1k−1xk−1(1−x)n−kΔx+o(Δx)= n\mathcal{C}_{n-1}^{k-1}x^{k-1}(1-x)^{n-k}\Delta x + o(\Delta x)=nCn−1k−1xk−1(1−x)n−kΔx+o(Δx)

从而得到 X(k)X_{(k)}X(k) 的概率密度函数 f(x)f(x)f(x) 为：

f(x)=lim⁡Δx→0P(x≤X(k)≤x+Δx)Δxf(x) = \lim_{\Delta x \rightarrow 0}\frac{P(x \leq X_{(k)} \leq x+ \Delta x) }{\Delta x}f(x)=Δx→0limΔxP(x≤X(k)≤x+Δx)
=nCn−1k−1xk−1(1−x)n−k=n\mathcal{C}_{n-1}^{k-1}x^{k-1}(1-x)^{n-k}=nCn−1k−1xk−1(1−x)n−k
=n!(k−1)!(n−k)!xk−1(1−x)n−k,x∈[0,1]=\frac{n!}{(k-1)!(n-k)!}x^{k-1}(1-x)^{n-k}, x \in [0,1]=(k−1)!(n−k)!n!xk−1(1−x)n−k,x∈[0,1]

此处我们考虑gamma函数的性质，Γ(n)=(n−1)!\Gamma(n) = (n-1)!Γ(n)=(n−1)!

得到：

f(x)=Γ(n+1)Γ(k)Γ(n−k+1)xk−1(1−x)n−kf(x) = \frac{\Gamma(n+1)}{\Gamma(k)\Gamma(n-k+1)}x^{k-1}(1-x)^{n-k}f(x)=Γ(k)Γ(n−k+1)Γ(n+1)xk−1(1−x)n−k

取 α=k,β=n−k+1\alpha = k,\beta= n-k+1α=k,β=n−k+1，最终得到：

f(x)=Γ(α+β)Γ(α)Γ(β)xα−1(1−x)β−1f(x) = \frac{\Gamma (\alpha+\beta)}{\Gamma (\alpha)\Gamma (\beta)}x^{\alpha-1}(1-x)^{\beta-1}f(x)=Γ(α)Γ(β)Γ(α+β)xα−1(1−x)β−1

这便是 Beta分布，其中：
Γ(α+β)Γ(α)Γ(β)=1B(α,β)\frac{\Gamma (\alpha+\beta)}{\Gamma (\alpha)\Gamma (\beta)}=\frac{1}{B(\alpha,\beta)} Γ(α)Γ(β)Γ(α+β)=B(α,β)1
B(α,β)B(\alpha,\beta)B(α,β) 是 Beta函数

最终 Beta分布 可表示为：

f(x)=1B(α,β)xα−1(1−x)β−1f(x) = \frac{1}{B(\alpha,\beta)} x^{\alpha-1}(1-x)^{\beta-1}f(x)=B(α,β)1xα−1(1−x)β−1

那么Beta分布 有什么用处呢？

我们举个例子，熟悉棒球运动的都知道有一个指标就是棒球击球率(batting average)，就是用一个运动员击中的球数除以击球的总数，我们一般认为0.266是正常水平的击球率，而如果击球率高达0.3就被认为是非常优秀的。

现在有一个棒球运动员，我们希望能够预测他在这一赛季中的棒球击球率是多少。你可能就会直接计算棒球击球率，用击中的数除以击球数，但是如果这个棒球运动员只打了一次，而且还命中了，那么他就击球率就是100%了，这显然是不合理的，因为根据棒球的历史信息，我们知道这个击球率应该是0.215到0.36之间才对啊。

对于这个问题，我们可以用一个二项分布表示（一系列成功或失败），一个最好的方法来表示这些经验（在统计中称为先验信息）就是用 Beta分布 ，这表示在我们没有看到这个运动员打球之前，我们就有了一个大概的范围。

Beta分布是指一组定义在(0,1)区间的连续概率分布，具有两个参数：α,β>0\alpha,\beta>0α,β>0

接下来我们将这些先验信息转换为 Beta分布 的参数，我们知道一个运动员的击球率应该是平均0.27左右，而他的范围是0.21到0.35，那么根据这个信息，我们可以取α=81,β=219\alpha=81,\beta=219α=81,β=219：

之所以取这两个参数是因为：

Beta分布 的均值是 0.27;
从图中可以看到这个分布主要落在了(0.2,0.35)间，这是从经验中得出的合理的范围。

其中，x轴就表示各个击球率的取值，x对应的y值就是这个击球率所对应的概率。也就是说beta分布可以看作一个概率的概率分布。

那么有了先验信息后，现在我们考虑一个运动员只打一次球，那么他现在的数据就是”1中;1击”。这时候我们就可以更新我们的分布了，让这个曲线做一些移动去适应我们的新信息。Beta分布 在数学上就给我们提供了这一性质，他与二项分布是共轭先验的（Conjugate_prior）。所谓共轭先验就是先验分布是 Beta分布 ，而后验分布同样是 Beta分布。结果很简单：

其中和是一开始的参数，在这里是81和219。所以在这一例子里，α\alphaα 增加了1(击中了一次)。β\betaβ 没有增加(没有漏球)。这就是我们的新的beta分布，我们跟原来的比较一下：

可以看到这个分布其实没多大变化，这是因为只打了1次球并不能说明什么问题。但是如果我们得到了更多的数据，假设一共打了300次，其中击中了100次，200次没击中，那么这一新分布就是：

注意到这个曲线变得更加尖，并且平移到了一个右边的位置，表示比平均水平要高。

一个有趣的事情是，根据这个新的 Beta分布，我们可以得出他的数学期望为：0.3，这一结果要比直接的估计要小。你可能已经意识到，我们事实上就是在这个运动员在击球之前可以理解为他已经成功了81次，失败了219次这样一个先验信息。

因此，对于一个我们不知道概率是什么，而又有一些合理的猜测时，Beta分布 能很好的作为一个表示概率的概率分布

3.6 `Beta分布`和`二项分布`共轭

回顾前面的问题1：假设随机变量X1,X2,...XnX_1,X_2,...X_nX1,X2,...Xn服从(0,1)(0,1)(0,1)的均匀分布，即X∽U(0,1)X \backsim U(0,1)X∽U(0,1)，把这个 nnn 个随机变量排序后得到顺序统计量 X(1),X(2),...,X(n)X_{(1)},X_{(2)},...,X_{(n)}X(1),X(2),...,X(n) ，请问 X(k)X_{(k)}X(k) 的分布是什么？

如果，现在我们要在问题1的基础上增加一些观测数据，变成问题2。

问题2：

随机变量X∽U(0,1)X \backsim U(0,1)X∽U(0,1)，对于顺序统计量为 X(1),X(2),...,X(n)X_{(1)},X_{(2)},...,X_{(n)}X(1),X(2),...,X(n)，需要猜测 p=X(k)p = X_{(k)}p=X(k);
观测数据 Y∽U(0,1)Y \backsim U(0,1)Y∽U(0,1)，YYY 中有 m1m_1m1 个比 ppp 小，有 m2m_2m2 个比 ppp 大；
那么请问：P(p∣Y)P(p|Y)P(p∣Y)的分布是什么？

根据 YYY 中有 m1m_1m1 个比 ppp 小，有 m2m_2m2 个比 ppp 大，可知，YYY 中有 m1m_1m1 个比 X(k)X_{(k)}X(k) 小，有 m2m_2m2 个比 X(k)X_{(k)}X(k) 大，所以 X(k)X_{(k)}X(k) 是 X1,X2,...,Xn,Y1,Y2,...,Ym∽U(0,1)X_1,X_2,...,X_n,Y_1,Y_2,...,Y_m \backsim U(0,1)X1,X2,...,Xn,Y1,Y2,...,Ym∽U(0,1)中的第 k+m1k+m_1k+m1 个数。

根据前面的结论：只要落在 [x,x+Δx][x,x+Δx][x,x+Δx] 内的数字超过一个，则对应的事件的概率就是 o(Δx)o(Δx)o(Δx)，继而推出事件服从 Beta分布 ，从而可知 p=X(k)p = X_{(k)}p=X(k) 的概率密度函数为：

Beta(p∣k+m1,n−k+1+m2)Beta(p|k+m_1,n-k+1+m_2)Beta(p∣k+m1,n−k+1+m2)

上面过程描述如下：

为了猜测 p=X(k)p = X_{(k)}p=X(k)，在获得一定的观测数据前，我们对 ppp 的认知是：f(p)=Beta(p∣k,n−k+1)f(p)=Beta(p|k,n-k+1)f(p)=Beta(p∣k,n−k+1)，此称为 ppp 的先验分布；
然后为了获得这个结果，YYY 中有 m1m_1m1 个比 ppp 小，有 m2m_2m2 个比 ppp 大，针对 YYY 做了 mmm 次伯努利试验，所以 mmm 服从二项分布 B(m,p)B(m,p)B(m,p)；
在给定了来自数据提供的 (m1,m2)(m_1,m_2)(m1,m2) 的知识后，ppp 的后验分布变为： f(p∣m1,m2)=Beta(p∣k+m1,n−k+1+m2)f(p|m_1,m_2)=Beta(p|k+m_1,n-k+1+m_2)f(p∣m1,m2)=Beta(p∣k+m1,n−k+1+m2)

熟悉贝叶斯方法的人会发现这一过程不正是贝叶斯方法的思考过程吗？

先验分布 π(θ)\pi(\theta)π(θ) + 样本信息 xxx ⇒\Rightarrow⇒ 后验分布 π(θ∣x)\pi(\theta|x)π(θ∣x)

根据上述，这种观测到的数据符合二项分布，参数的先验分布和后验分布都是Beta分布的情况，就是Beta分布-二项分布共轭，换言之， Beta分布是二项分布的共轭先验概率分布

二项分布和Beta分布是共轭分布意味着，如果我们为二项分布的参数p选取的先验分布是Beta分布，那么以p为参数的二项分布用贝叶斯估计得到的后验分布仍然服从Beta分布。

3.7 共轭先验分布

在贝叶斯概率理论中，如果后验概率 p(θ∣x)p(θ|x)p(θ∣x) 和先验概率 p(θ)p(θ)p(θ) 满足同样的分布律，那么，先验分布和后验分布被叫做共轭分布，同时，先验分布叫做似然函数的共轭先验分布。

投掷一个非均匀硬币，可以使用参数为 θ\thetaθ 的伯努利模型， θ\thetaθ 为硬币为正面的概率，那么结果 xxx 的分布（似然函数）形式为（这个在前面的二项分布中介绍过）：

p(x∣θ)=θx×(1−θ)1−xp(x|\theta) = \theta^x \times (1-\theta)^{1-x}p(x∣θ)=θx×(1−θ)1−x

其共轭先验为 Beta分布，具有两个参数 α,β\alpha,\betaα,β，称为超参数。且这两个参数决定了 θ\thetaθ 参数的分布，为：

p(θ∣α,β)=θα−1(1−θ)β−1∫01θα−1(1−θ)β−1dθp(\theta|\alpha,\beta)=\frac{\theta^{\alpha-1}(1-\theta)^{\beta-1}}{\int_0^1 \theta^{\alpha-1}(1-\theta)^{\beta-1}d\theta}p(θ∣α,β)=∫01θα−1(1−θ)β−1dθθα−1(1−θ)β−1

然后根据贝叶斯公式计算后验分布的：

p(θ∣x)∝p(x∣θ)×p(θ)=θx×(1−θ)1−xθα−1(1−θ)β−1=θx+α−1(1−θ)1−x+β−1p(\theta|x) \propto p(x|\theta)\times p(\theta) = \theta^x \times (1-\theta)^{1-x}\theta^{\alpha-1}(1-\theta)^{\beta-1} = \theta^{x+\alpha-1}(1-\theta)^{1-x+\beta-1}p(θ∣x)∝p(x∣θ)×p(θ)=θx×(1−θ)1−xθα−1(1−θ)β−1=θx+α−1(1−θ)1−x+β−1

我们发现得到的结果是另一个Beta分布，因此先验分布和后验分布都是Beta分布，两者共轭。从而说明了先验分布 Beta分布 确实是似然函数 二项分布 的共轭先验分布。

3.8 狄利克雷分布（`Dirichlet` 分布）

前面我们讲到 Beta分布 是 二项分布 的共轭先验分布，Dirichlet分布 则是 多项分布 的共轭先验分布。Dirichlet（狄利克雷）同时可以看做是将 Beta分布 推广到多变量的情形。概率密度函数定义如下:

f(x1,x2,...xK−1;α1,α2,...,αK)=1B(α)∏i=1Kxiαi−1f(x_1,x_2,...x_{K-1};\alpha_1,\alpha_2,...,\alpha_K) = \frac{1}{B(\alpha)}\prod_{i=1}^Kx_i^{\alpha_i-1}f(x1,x2,...xK−1;α1,α2,...,αK)=B(α)1i=1∏Kxiαi−1

其中，B(α)B(\alpha)B(α) 相当于是多项 Beta函数：

B(α)=∏i=1KΓ(αi)Γ(∑i=1Kαi)B(\alpha) = \frac{\prod_{i=1}^K\Gamma(\alpha_i)}{\Gamma(\sum_{i=1}^K\alpha_i)}B(α)=Γ(∑i=1Kαi)∏i=1KΓ(αi)

且 α=(α1,α2,...,αK)\alpha = (\alpha_1,\alpha_2,...,\alpha_K)α=(α1,α2,...,αK)

问题3：

我们由问题2引出问题3：

随机变量X∽U(0,1)X \backsim U(0,1)X∽U(0,1)
排序后的顺序统计量为 X(1),X(2),...,X(n)X_{(1)},X_{(2)},...,X_{(n)}X(1),X(2),...,X(n)
问 (X(k1),X(k1+k2))(X_{(k_1)},X_{(k_1+k_2)})(X(k1),X(k1+k2)) 的联合分布是什么？

为了简化计算，取 x3x_3x3满足 x1+x2+x3=1x_1+x_2+x_3=1x1+x2+x3=1,但只有x1,x2x_1,x_2x1,x2是变量，如下图所示：

从而有：

P(X(k1)∈(x1,x1+Δx),X(k1+k2)∈(x2,x2+Δx))P(X_{(k_1)} \in (x_1,x_1+\Delta x), X_{(k_1+k_2)} \in (x_2, x_2+\Delta x))P(X(k1)∈(x1,x1+Δx),X(k1+k2)∈(x2,x2+Δx))
=n(n−1)Ck1−1,k2−1n−2x1k1−1x2k2−1x3n−k1−k2(Δx)2=n(n-1)\mathcal{C}_{k_1-1,k_2-1}^{n-2}x_1^{k_1-1}x_2^{k_2-1 }x_3^{n-k_1-k_2}(\Delta x)^2=n(n−1)Ck1−1,k2−1n−2x1k1−1x2k2−1x3n−k1−k2(Δx)2
=n!(k1−1)!(k2−1)!(n−k1−k2)!x1k1−1x2k2−1x3n−k1−k2(Δx)2=\frac{n!}{(k_1-1)!(k_2-1)!(n-k_1-k_2)!}x_1^{k_1-1}x_2^{k_2-1 }x_3^{n-k_1-k_2}(\Delta x)^2=(k1−1)!(k2−1)!(n−k1−k2)!n!x1k1−1x2k2−1x3n−k1−k2(Δx)2

继而得到(X(k1),X(k1+k2))(X_{(k_1)},X_{(k_1+k_2)})(X(k1),X(k1+k2)) 的联合分布为：

f(x1,x2,x3)=n!(k1−1)!(k2−1)!(n−k1−k2)!x1k1−1x2k2−1x3n−k1−k2f(x_1,x_2,x_3) = \frac{n!}{(k_1-1)!(k_2-1)!(n-k_1-k_2)!}x_1^{k_1-1}x_2^{k_2-1 }x_3^{n-k_1-k_2}f(x1,x2,x3)=(k1−1)!(k2−1)!(n−k1−k2)!n!x1k1−1x2k2−1x3n−k1−k2
=Γ(n+1)Γ(k1)Γ(k2)Γ(n−k1−k2+1)x1k1−1x2k2−1x3n−k1−k2= \frac{\Gamma(n+1)}{\Gamma(k_1)\Gamma(k_2)\Gamma(n-k_1-k_2+1)}x_1^{k_1-1}x_2^{k_2-1 }x_3^{n-k_1-k_2}=Γ(k1)Γ(k2)Γ(n−k1−k2+1)Γ(n+1)x1k1−1x2k2−1x3n−k1−k2

观察上述式子的最终结果，可以看出上面这个分布其实就是3维形式的 Dirichlet 分布:

Dir(x1,x2,x3∣k1,k2,n−k1−k2+1)Dir(x_1,x_2,x_3|k_1,k_2,n-k_1-k_2+1)Dir(x1,x2,x3∣k1,k2,n−k1−k2+1)

令α=k1,α2=k2,α3=n−k1−k2+1\alpha=k_1,\alpha_2=k_2,\alpha_3=n-k_1-k_2+1α=k1,α2=k2,α3=n−k1−k2+1，得到：

f(x1,x2,x3)=Γ(α1+α2+α3)Γ(α1)Γ(α2)Γ(α3)x1k1−1x2k2−1x3n−k1−k2f(x_1,x_2,x_3) = \frac{\Gamma(\alpha_1+\alpha_2+\alpha_3)}{\Gamma(\alpha_1)\Gamma(\alpha_2)\Gamma(\alpha_3)}x_1^{k_1-1}x_2^{k_2-1 }x_3^{n-k_1-k_2}f(x1,x2,x3)=Γ(α1)Γ(α2)Γ(α3)Γ(α1+α2+α3)x1k1−1x2k2−1x3n−k1−k2

这个就是一般形式的3维 Dirichlet 分布。

3.9 `Dirichlet分布` 和 `多项分布` 共轭

为了论证Dirichlet分布是多项式分布的共轭先验概率分布，下面咱们继续在上述问题3的基础上再进一步，提出问题4。

问题4：

随机变量X∽U(0,1)X \backsim U(0,1)X∽U(0,1)，排序后的顺序统计量为 X(1),X(2),...,X(n)X_{(1)},X_{(2)},...,X_{(n)}X(1),X(2),...,X(n);
令p1=X(k1),p2=X(k1+k2),p3=1−p1−p2p_1=X_{(k_1)},p_2=X_{(k_1+k_2)},p_3=1-p_1-p_2p1=X(k1),p2=X(k1+k2),p3=1−p1−p2，(此处的p3p_3p3非变量，只是为了表达方便)，现在要猜测 p⃗=(p1,p2,p3)\vec{p} = (p_1,p_2,p_3)p=(p1,p2,p3);
观测数据 Y∽U(0,1)Y \backsim U(0,1)Y∽U(0,1)，YYY中落到[0,p1),[p1,p2),[p2,1][0,p_1),[p_1,p_2),[p_2,1][0,p1),[p1,p2),[p2,1]三个区间的个数分别为 m1,m2,m3,m=m1+m2+m3m_1,m_2,m_3,m=m_1+m_2+m_3m1,m2,m3,m=m1+m2+m3;
问后验分别 P(p⃗∣Y)P(\vec{p}|Y)P(p∣Y)的分布是什么？

为了方便讨论，我们记 m⃗=(m1,m2,m3),k⃗=(k1,k2,n−k1−k2+1)\vec{m} = (m_1,m_2,m_3), \vec{k} = (k_1,k_2,n-k_1-k_2+1)m=(m1,m2,m3),k=(k1,k2,n−k1−k2+1)，根据已知条件，YYY中落到[0,p1),[p1,p2),[p2,1][0,p_1),[p_1,p_2),[p_2,1][0,p1),[p1,p2),[p2,1]三个区间的个数分别为m1,m2,m3m_1,m_2,m_3m1,m2,m3；可知 p1,p2p_1,p_2p1,p2分别是这 m+nm+nm+n 个数中第 k1+m1k_1+m_1k1+m1 大，第 k2+m2k_2+m_2k2+m2 大的数，于是，后验分布：

P(p⃗∣Y)=Dir(k⃗∣k1+m1,k2+m2,n−k1−k2+1+m3)P(\vec{p}|Y) = Dir(\vec{k}|k_1+m_1,k_2+m_2,n-k_1-k_2+1+m_3)P(p∣Y)=Dir(k∣k1+m1,k2+m2,n−k1−k2+1+m3)

即：Dir(p⃗∣k⃗+m⃗)Dir(\vec{p}|\vec{k}+\vec{m})Dir(p∣k+m)

同样，按照贝叶斯推理的逻辑，可将上述过程整理如下：

我们要猜测参数p⃗=(p1,p2,p3)\vec{p} = (p_1,p_2,p_3)p=(p1,p2,p3)，其先验分布为 Dir(p⃗∣k⃗)Dir(\vec{p}|\vec{k})Dir(p∣k)
数据 YYY中落到[0,p1),[p1,p2),[p2,1][0,p_1),[p_1,p_2),[p_2,1][0,p1),[p1,p2),[p2,1]三个区间的个数分别为 m1,m2,m3m_1,m_2,m_3m1,m2,m3，所以 m⃗=(m1,m2,m3)\vec{m} = (m_1,m_2,m_3)m=(m1,m2,m3)服从多项分布 Mult(m⃗∣p⃗)Mult(\vec{m}|\vec{p})Mult(m∣p)
在给定了来自数据提供的知识 m⃗\vec{m}m 后，p⃗\vec{p}p 的后验分布变为：Dir(p⃗∣k⃗+m⃗)Dir(\vec{p}|\vec{k}+\vec{m})Dir(p∣k+m)

针对于这种观测到的数据符合多项分布，参数的先验分布和后验分布都是Dirichlet 分布的情况，就是Dirichlet分布-多项分布共轭。换言之，至此已经证明了Dirichlet分布的确就是多项式分布的共轭先验概率分布。

至此，我们已经介绍了gamma函数、二项分布、多项分布、Beta分布、Dirichlet分布，接下来我们就要正式介绍LDA模型

四、`LDA` 主题模型

首先我们对前面的数学知识点进行总结，便于跳过前面数学部分的人可以直接从这里看。

4.1 知识点回顾

知识点1：

Beta分布 是 二项分布 的共轭先验概率分布；

对于α,β>0\alpha, \beta > 0α,β>0，有如下关系：
Beta(p∣α,β)+Count(m1,m2)=Beta(p∣α+m1,β+m2)Beta(p|\alpha,\beta) + Count(m_1,m_2) =Beta(p|\alpha+m_1,\beta+m_2) Beta(p∣α,β)+Count(m1,m2)=Beta(p∣α+m1,β+m2)
其中：(m1,m2)(m_1,m_2)(m1,m2) 对应的是二项分布 B(m1+m2,p)B(m_1+m_2,p)B(m1+m2,p) 的计数。针对于这种观测到的数据符合二项分布，参数的先验分布和后验分布都是Beta分布的情况，就是Beta-Binomial 共轭。

知识点2：

Dirichlet分布 是 多项式分布 的共轭先验概率分布;

将α⃗\vec{\alpha}α 从整数集合延拓到实数集合，得到更一般的表达式：
Dir(p⃗∣α⃗)+MultCount(m⃗)=Dir(p⃗∣α+m⃗⃗)Dir(\vec{p}|\vec{\alpha}) + MultCount(\vec{m}) = Dir(\vec{p}|\vec{\alpha+\vec{m}})Dir(p∣α)+MultCount(m)=Dir(p∣α+m)
针对于这种观测到的数据符合 多项分布 ，参数的先验分布和后验分布都是 Dirichlet 分布 的情况，就是 Dirichlet-Multinomial 共轭。

知识点3：

贝叶斯派思考问题的固定模式

先验分布 π(θ)\pi(\theta)π(θ) + 样本信息 xxx ⇒\Rightarrow⇒ 后验分布 π(θ∣x)\pi(\theta|x)π(θ∣x)

上述思考模式意味着，新观察到的样本信息将修正人们以前对事物的认知。换言之，在得到新的样本信息之前，人们对 θ\thetaθ 的认知是先验分布 π(θ)\pi(\theta)π(θ)，在得到新的样本信息 X\mathcal{X}X后，人们对的认知为 π(θ∣x)\pi(\theta|x)π(θ∣x)。

下面我们将通过由基础模型到LDA模型的渐进方式进行介绍。

为了便于下面描述，首先定义一些变量：

www 表示词， VVV表示所有单词的个数（固定值）
zzz 表示主题，KKK是主题的个数（预先给定，固定值）
D=(W1,W2,...,WM)D = (W_1,W_2,...,W_M)D=(W1,W2,...,WM)表示语料库，其中的MMM是语料库中的文档数（固定值）
W=(w1,w2,..,wN)W = (w_1,w_2,..,w_N)W=(w1,w2,..,wN)表示文档，其中 NNN 表示一个文档中的词数（随机变量）

4.2 基础模型1：`Unigram model`

对于文档W=(w1,w2,..,wN)W = (w_1,w_2,..,w_N)W=(w1,w2,..,wN)，用p(wn)p(w_n)p(wn)表示词wnw_nwn的先验概率，生成文档 WWW 的概率为：

p(W)=∏n=1Np(wn)p(W) = \prod_{n=1}^Np(w_n)p(W)=n=1∏Np(wn)

其图模型为（图中被涂色的w表示可观测变量，N表示一篇文档中总共N个单词，M表示M篇文档）：

或者表示为：

unigram model 假设文本中的词服从Multinomial分布，而我们已经知道Multinomial分布的先验分布为Dirichlet分布。

上图中的表示在文本中观察到的第n个词，n∈[1,N]n∈[1,N]n∈[1,N]表示该文本中一共有N个单词。加上方框表示重复，即一共有N个这样的随机变量。其中，ppp 和 ααα 是隐含未知变量：

ppp是词服从的Multinomial分布的参数
ααα 是Dirichlet分布（即Multinomial分布的先验分布）的参数。

一般 ααα 由经验事先给定，ppp 由观察到的文本中出现的词学习得到，表示文本中出现每个词的概率。

4.3 基础模型2：`Mixture of unigrams model`

该模型的生成过程是：给某个文档先选择一个主题，再根据该主题生成文档，该文档中的所有词都来自一个主题。假设主题有 z1,z2,...,zkz_1,z_2,...,z_kz1,z2,...,zk，生成文档 WWW 的概率为：

p(W)=p(z1)∏n=1Np(wn∣z1)+...+p(zk)∏n=1Np(wn∣zk)=∑zp(z)∏n=1Np(wn∣z)p(W) = p(z_1)\prod_{n=1}^Np(w_n|z_1)+...+p(z_k)\prod_{n=1}^Np(w_n|z_k) = \sum_zp(z)\prod_{n=1}^Np(w_n|z)p(W)=p(z1)n=1∏Np(wn∣z1)+...+p(zk)n=1∏Np(wn∣zk)=z∑p(z)n=1∏Np(wn∣z)

其图模型为（图中被涂色的w表示可观测变量，未被涂色的z表示未知的隐变量，N表示一篇文档中总共N个单词，M表示M篇文档）：

4.4 基础模型4：`PLSA模型`

PLSA模型 是最接近 LDA模型，所以理解PLSA模型有助于我们理解 LDA模型。

（1）pLSA模型下生成文档过程

在上面的Mixture of unigrams model中，我们假定一篇文档只有一个主题生成，可实际中，一篇文章往往有多个主题，只是这多个主题各自在文档中出现的概率大小不一样。比如介绍一个国家的文档中，往往会分别从教育、经济、交通等多个主题进行介绍。那么在pLSA中，文档是怎样被生成的呢？

假设你要写M篇文档，由于一篇文档由各个不同的词组成，所以你需要确定每篇文档里每个位置上的词。
再假定你一共有K个可选的主题，有V个可选的词，咱们来玩一个扔骰子的游戏。

第一步：

假设你每写一篇文档会制作一颗K面的“文档-主题”骰子（扔此骰子能得到K个主题中的任意一个），和K个V面的“主题-词项” 骰子（每个骰子对应一个主题，K个骰子对应之前的K个主题，且骰子的每一面对应要选择的词项，V个面对应着V个可选的词）。

比如可令K=3，即制作1个含有3个主题的“文档-主题”骰子，这3个主题可以是：教育、经济、交通。然后令V = 3，制作3个有着3面的“主题-词项”骰子，其中，教育主题骰子的3个面上的词可以是：大学、老师、课程，经济主题骰子的3个面上的词可以是：市场、企业、金融，交通主题骰子的3个面上的词可以是：高铁、汽车、飞机。

第二步：

每写一个词，先扔该“文档-主题”骰子选择主题，得到主题的结果后，使用和主题结果对应的那颗“主题-词项”骰子，扔该骰子选择要写的词。

先扔“文档-主题”的骰子，假设（以一定的概率）得到的主题是 教育，所以下一步便是扔 教育 主题筛子，（以一定的概率）得到教育主题筛子对应的某个词：大学。
上面这个投骰子产生词的过程简化下便是：“先以一定的概率选取主题，再以一定的概率选取词”。事实上，一开始可供选择的主题有3个：教育、经济、交通，那为何偏偏选取教育这个主题呢？其实是随机选取的，只是这个随机遵循一定的概率分布。比如可能选取教育主题的概率是0.5，选取经济主题的概率是0.3，选取交通主题的概率是0.2，那么这3个主题的概率分布便是 {教育：0.5，经济：0.3，交通：0.2}，我们把各个主题z在文档d中出现的概率分布称之为 主题分布，且是一个多项分布。
同样的，从主题分布中随机抽取出教育主题后，依然面对着3个词：大学、老师、课程，这3个词都可能被选中，但它们被选中的概率也是不一样的。比如大学这个词被选中的概率是0.5，老师这个词被选中的概率是0.3，课程被选中的概率是0.2，那么这3个词的概率分布便是 {大学：0.5，老师：0.3，课程：0.2}，我们把各个词语w在主题z下出现的概率分布称之为 词分布，这个词分布也是一个多项分布。
所以，选主题和选词都是两个随机的过程，先从主题分布{教育：0.5，经济：0.3，交通：0.2}中抽取出主题：教育，然后从该教育主题对应的词分布{大学：0.5，老师：0.3，课程：0.2}中抽取出词：大学。

第三步：

最后，你不停的重复扔“文档-主题”骰子和”主题-词项“骰子，重复N次（产生N个词），完成一篇文档，重复这产生一篇文档的方法M次，则完成M篇文档。

上述过程抽象出来即是PLSA的文档生成模型。在这个过程中，我们并未关注词和词之间的出现顺序，所以pLSA是一种词袋方法。具体说来，该模型假设一组共现(co-occurrence)词项关联着一个隐含的主题类别 zk∈{z1,z2,...,zK}z_k \in \{z_1,z_2,...,z_K\}zk∈{z1,z2,...,zK}。同时定义：
1、P(di)P(d_i)P(di)表示海量文档中某篇文档被选中的概率。

2、P(wj∣di)P(w_j|d_i)P(wj∣di) 表示词 wjw_jwj 在给定文档 did_idi中出现的概率。

怎么计算得到呢？针对海量文档，对所有文档进行分词后，得到一个词汇列表，这样每篇文档就是一个词语的集合。对于每个词语，用它在文档中出现的次数除以文档中词语总的数目便是它在文档中出现的概率。

3、P(zk∣di)P(z_k|d_i)P(zk∣di)表示具体某个主题 zkz_kzk 在给定文档 did_idi下出现的概率。

4、P(wj∣zk)P(w_j|z_k)P(wj∣zk)表示具体某个词 wjw_jwj 在给定主题 zkz_kzk 下出现的概率，与主题关系越密切的词，其条件概率 p(wj∣zk)p(w_j|z_k)p(wj∣zk)越大。

利用上述的第1、3、4个概率，我们便可以按照如下的步骤得到“文档-词项”的生成模型：

按照概率 P(di)P(d_i)P(di) 选择一篇文档 did_idi;
选定文档 did_idi 后，从主题分布中按照概率 p(zk∣di)p(z_k|d_i)p(zk∣di)选择一个隐含的主题类别 zkz_kzk;
选定主题 zkz_kzk 后，从词分布中按照概率 p(wj∣zk)p(w_j|z_k)p(wj∣zk) 选择一个词 wjw_jwj

所以pLSA中生成文档的整个过程便是选定文档生成主题，确定主题生成词。

（2）根据文档反推其主题分布

反过来，既然文档已经产生，那么如何根据已经产生好的文档反推其主题呢？这个利用看到的文档推断其隐藏的主题（分布）的过程（其实也就是产生文档的逆过程），便是主题建模的目的：自动地发现文档集中的主题（分布）。

换言之，人类根据文档生成模型写成了各类文章，然后丢给了计算机，相当于计算机看到的是一篇篇已经写好的文章。现在计算机需要根据一篇篇文章中看到的一系列词归纳出当篇文章的主题，进而得出各个主题各自不同的出现概率：主题分布。即文档 ddd 和单词 www 是可被观察到的，但主题 zzz 却是隐藏的。

如下图所示（图中被涂色的d、w表示可观测变量，未被涂色的z表示未知的隐变量，N表示一篇文档中总共N个单词，M表示M篇文档）：

上图中，文档d和词w是我们得到的样本（样本随机，参数虽未知但固定，所以pLSA属于频率派思想。区别于下文要介绍的LDA中：样本固定，参数未知但不固定，是个随机变量，服从一定的分布，所以 LDA 属于贝叶斯派思想），可观测得到，所以对于任意一篇文档，其 P(wj∣di)P(w_j|d_i)P(wj∣di)是已知的。

频率派与贝叶斯派各自不同的思考方式：

频率派把需要推断的参数θ看做是固定的未知常数，即概率虽然是未知的，但最起码是确定的一个值，同时，样本X 是随机的，所以频率派重点研究样本空间，大部分的概率计算都是针对样本X 的分布；

而贝叶斯派的观点则截然相反，他们认为待估计的参数是随机变量，服从一定的分布，而样本X 是固定的，由于样本是固定的，所以他们重点研究的是参数的分布。

从而可以根据大量已知的文档-词项信息 P(wj∣di)P(w_j|d_i)P(wj∣di) ，训练出文档-主题 P(zk∣di)P(z_k|d_i)P(zk∣di) 和主题-词项 P(wj∣zk)P(w_j|z_k)P(wj∣zk)，如下公式所示：

P(wj∣zk)=∑k=1KP(wj∣zk)P(zk∣di)P(w_j|z_k) = \sum_{k=1}^KP(w_j|z_k)P(z_k|d_i)P(wj∣zk)=k=1∑KP(wj∣zk)P(zk∣di)

故得到文档中每个词的生成概率为：

P(di,wj)=P(di)P(wj∣di)=P(di)∑k=1KP(wj∣zk)P(zk∣di)P(d_i,w_j)=P(d_i)P(w_j|d_i)=P(d_i)\sum_{k=1}^KP(w_j|z_k)P(z_k|d_i)P(di,wj)=P(di)P(wj∣di)=P(di)k=1∑KP(wj∣zk)P(zk∣di)

由于 P(di)P(d_i)P(di) 可事先计算求出，而 P(wj∣zk),P(zk∣di)P(w_j|z_k),P(z_k|d_i)P(wj∣zk),P(zk∣di)未知，所以 θ=(P(wj∣zk),P(zk∣di))\theta = (P(w_j|z_k),P(z_k|d_i))θ=(P(wj∣zk),P(zk∣di))就是我们要估计的参数值，通俗地将，就是要最大化这个 θ\thetaθ。

用什么方法进行估计呢，常用的参数估计方法有极大似然估计MLE、最大后验证估计MAP、贝叶斯估计等等。因为该待估计的参数中含有隐变量z，所以我们可以考虑EM算法。

关于EM算法可以参考：EM算法

（3）EM算法估计pLSA的两个未知参数

首先尝试从矩阵的角度来描述待估计的两个未知变量P(wj∣zk)P(w_j|z_k)P(wj∣zk)和 P(zk∣di)P(z_k|d_i)P(zk∣di)。

假定用 ϕk\phi_kϕk 表示词表 VVV 在主题 zkz_kzk 上的一个多项分布，则 ϕk\phi_kϕk 可以表示成一个向量，每个元素 ϕk,j\phi_{k,j}ϕk,j表示词项 wjw_jwj 出现在主题 zkz_kzk 中的概率，即：

P(wj∣zk)=ϕk,j,∑wj∈Vϕk,j=1P(w_j|z_k) = \phi_{k,j} , \sum_{w_j \in V}\phi_{k,j} = 1 P(wj∣zk)=ϕk,j,wj∈V∑ϕk,j=1

用θi\theta_iθi 表示所有主题 ZZZ 在文档 did_idi 上的一个多项分布，则 θi\theta_iθi 可以表示成一个向量，每个元素 θi,k\theta_{i,k}θi,k表示主题 zkz_kzk 出现在文档 did_idi 中的概率，即：

P(zk∣di)=θi,k,∑zk∈Zθi,k=1P(z_k|d_i) = \theta_{i,k}, \sum_{z_k \in Z} \theta_{i,k} = 1P(zk∣di)=θi,k,zk∈Z∑θi,k=1

这样，巧妙的把 P(wj∣zk)P(w_j|z_k)P(wj∣zk) 和 P(zk∣di)P(z_k|d_i)P(zk∣di) 转换成了两个矩阵。换言之，最终我们要求解的参数是这两个矩阵：

Φ=[ϕ1,ϕ2,...,ϕK],zk∈Z\Phi = [\phi_1,\phi_2,...,\phi_K], z_k \in ZΦ=[ϕ1,ϕ2,...,ϕK],zk∈Z
Θ=[θ1,θ2,...,θM],di∈D\Theta = [\theta_1,\theta_2,...,\theta_M], d_i \in DΘ=[θ1,θ2,...,θM],di∈D

由于词与词之间相互独立，所以整篇文档 N 个词的分布为：

P(W∣di)=∏j=1NP(di,wj)n(di,wj)P(W|d_i) = \prod_{j=1}^N P(d_i,w_j)^{n(d_i,w_j)}P(W∣di)=j=1∏NP(di,wj)n(di,wj)

再由于文档和文档之间也是相互独立的，所以整个语料库中词的分布为（整个语料库M篇文档，每篇文档N个词）：

P(W∣D)=∏i=1M∏j=1NP(di,wj)n(di,wj)P(W|D)=\prod_{i=1}^M\prod_{j=1}^NP(d_i,w_j)^{n(d_i,w_j)}P(W∣D)=i=1∏Mj=1∏NP(di,wj)n(di,wj)

其中，n(di,wj)n(d_i,w_j)n(di,wj) 表示词项 wjw_jwj 在文档 did_idi 中的词频，n(di)n(d_i)n(di) 表示文档 did_idi 中词的总数，显然，n(di)=∑wj∈Vn(di,wj)n(d_i) = \sum_{w_j \in V}n(d_i, w_j)n(di)=∑wj∈Vn(di,wj)

从而得到整个语料库的词分布的对数似然函数：

ℓ(Φ,Θ)=∑i=1M∑j=1Nn(di,wj)logP(di,wj)\ell(\Phi, \Theta) = \sum_{i=1}^M\sum_{j=1}^Nn(d_i,w_j)logP(d_i,w_j)ℓ(Φ,Θ)=i=1∑Mj=1∑Nn(di,wj)logP(di,wj)
=∑i=1Mn(di)(logP(di)+∑j=1Nn(di,wj)n(di)log∑k=1KP(wj∣zk)P(zk∣di))= \sum_{i=1}^Mn(d_i)\Big(logP(d_i)+\sum_{j=1}^N\frac{n(d_i,w_j)}{n(d_i)}log\sum_{k=1}^KP(w_j|z_k)P(z_k|d_i)\Big)=i=1∑Mn(di)(logP(di)+j=1∑Nn(di)n(di,wj)logk=1∑KP(wj∣zk)P(zk∣di))
=∑i=1Mn(di)(logP(di)+∑j=1Nn(di,wj)n(di)log∑k=1Kϕk,jθi,k)=\sum_{i=1}^Mn(d_i)\Big(logP(d_i)+\sum_{j=1}^N\frac{n(d_i,w_j)}{n(d_i)}log\sum_{k=1}^K\phi_{k,j}\theta_{i,k}\Big)=i=1∑Mn(di)(logP(di)+j=1∑Nn(di)n(di,wj)logk=1∑Kϕk,jθi,k)

现在，我们需要最大化上述这个对数似然函数来求解参数 ϕk,j,θi,k\phi_{k,j},\theta_{i,k}ϕk,j,θi,k，对于这种含有隐变量的最大似然估计，可以使用EM算法，EM算法分为两个步骤：先E-step，后M-step

E-step：假定参数已知，计算此时隐变量的后验概率

利用贝叶斯法则，可以得到：

P(zk∣di,wj)=P(zk,di,wj)∑l=1KP(zl,di,wj)P(z_k|d_i,w_j) = \frac{P(z_k,d_i,w_j)}{\sum_{l=1}^KP(z_l,d_i,w_j)}P(zk∣di,wj)=∑l=1KP(zl,di,wj)P(zk,di,wj)
=P(wj∣di,zk)P(zk∣di)P(di)∑l=1KP(wj∣di,zl)P(zl∣di)P(di)=\frac{P(w_j|d_i,z_k)P(z_k|d_i)P(d_i)}{\sum_{l=1}^KP(w_j|d_i,z_l)P(z_l|d_i)P(d_i)}=∑l=1KP(wj∣di,zl)P(zl∣di)P(di)P(wj∣di,zk)P(zk∣di)P(di)
=P(wj∣zk)P(zk∣di)∑l=1KP(wj∣zk)P(zl∣di)=\frac{P(w_j|z_k)P(z_k|d_i)}{\sum_{l=1}^KP(w_j|z_k)P(z_l|d_i)}=∑l=1KP(wj∣zk)P(zl∣di)P(wj∣zk)P(zk∣di)
=ϕk,jθi,k∑l=1Kϕl,jθi,l=\frac{\phi_{k,j}\theta_{i,k}}{\sum_{l=1}^K\phi_{l,j}\theta_{i,l}}=∑l=1Kϕl,jθi,lϕk,jθi,k

M-step：带入隐变量的后验概率，最大化样本分布的对数似然函数，求解相应的参数。

观察之前得到的对数似然函数 ℓ(Φ,Θ)\ell(\Phi, \Theta)ℓ(Φ,Θ) 的结果，由于文档长度 P(di)∝n(di)P(d_i) \propto n(d_i)P(di)∝n(di)可以单独计算，所以去掉它不影响最大化似然函数。此外，根据 E-step 的计算结果，把 ϕk,jθi,k=P(zk∣di,wj)∑l=1Kϕl,jθi,l\phi_{k,j}\theta_{i,k} = P(z_k|d_i,w_j)\sum_{l=1}^K\phi_{l,j}\theta_{i,l}ϕk,jθi,k=P(zk∣di,wj)∑l=1Kϕl,jθi,l 代入 ℓ(Φ,Θ)\ell(\Phi, \Theta)ℓ(Φ,Θ) ，于是我们只要最大化下面这个函数 ℓ\ellℓ 即：
ℓ=∑i=1M∑j=1Nn(di,wj)∑k=1KP(zk∣di,wj)log[ϕk,jθi,k]\ell = \sum_{i=1}^M\sum_{j=1}^Nn(d_i,w_j)\sum_{k=1}^KP(z_k|d_i,w_j)log[\phi_{k,j}\theta_{i,k}]ℓ=i=1∑Mj=1∑Nn(di,wj)k=1∑KP(zk∣di,wj)log[ϕk,jθi,k]

这是一个多元函数求极值问题，并且已知有如下约束条件：

∑j=1Nϕk,j=1\sum_{j=1}^N\phi_{k,j} = 1j=1∑Nϕk,j=1
∑k=1Kθi,k=1\sum_{k=1}^K\theta_{i,k} = 1k=1∑Kθi,k=1

熟悉凸优化，一般处理这种带有约束条件的极值问题，常用的方法便是 拉格朗日乘数法，即通过引入拉格朗日乘子将约束条件和多元（目标）函数融合到一起，转化为无约束条件的极值问题。

这里我们引入两个拉格朗日乘子 γ\gammaγ 和 ρ\rhoρ，从而写出拉格朗日函数：

H=ℓ+∑k=1Kγk(1−∑j=1Nϕk,j)+∑i=1Mρi(1−∑k=1Kθi,k)\mathcal{H} = \ell + \sum_{k=1}^K \gamma_k\Big(1-\sum_{j=1}^N\phi_{k,j} \Big) + \sum_{i=1}^M\rho_i\Big(1-\sum_{k=1}^K\theta_{i,k}\Big)H=ℓ+k=1∑Kγk(1−j=1∑Nϕk,j)+i=1∑Mρi(1−k=1∑Kθi,k)

因为我们要求解的参数是 ϕk,j\phi_{k,j}ϕk,j 和 θi,k\theta_{i,k}θi,k，所以对 ϕk,j\phi_{k,j}ϕk,j 和 θi,k\theta_{i,k}θi,k 求偏导，然后令偏导结果为 0，得到：

∂H∂ϕk,j=∑i=1Mn(di,wj)P(zk∣di,wj)−γkϕk,j=0,1≤j≤<N,1≤k≤K\frac{\partial \mathcal{H}}{\partial \phi_{k,j}} = \sum_{i=1}^Mn(d_i,w_j)P(z_k|d_i,w_j) -\gamma_k\phi_{k,j} = 0, 1 \leq j \leq < N,1\leq k\leq K∂ϕk,j∂H=i=1∑Mn(di,wj)P(zk∣di,wj)−γkϕk,j=0,1≤j≤<N,1≤k≤K

∂H∂θi,k=∑j=1Nn(di,wj)P(zk∣di,wj)−ρiθi,k=0,1≤i≤<M,1≤k≤K\frac{\partial \mathcal{H}}{\partial \theta_{i,k}} = \sum_{j=1}^Nn(d_i,w_j)P(z_k|d_i,w_j) -\rho_i\theta_{i,k} = 0, 1 \leq i \leq < M,1\leq k\leq K∂θi,k∂H=j=1∑Nn(di,wj)P(zk∣di,wj)−ρiθi,k=0,1≤i≤<M,1≤k≤K

消去拉格朗日乘子，最终可估计出参数 ϕk,j\phi_{k,j}ϕk,j 和 θi,k\theta_{i,k}θi,k：

ϕk,j=∑i=1Mn(di,wj)P(zk∣di,wj)∑l=1N∑i=1Mn(di,wl)P(zk∣di,wl)\phi_{k,j} =\frac{ \sum_{i=1}^Mn(d_i,w_j)P(z_k|d_i,w_j)}{\sum_{l=1}^N \sum_{i=1}^Mn(d_i,w_l)P(z_k|d_i,w_l)}ϕk,j=∑l=1N∑i=1Mn(di,wl)P(zk∣di,wl)∑i=1Mn(di,wj)P(zk∣di,wj)

θi,k=∑j=1Nn(di,wj)P(zk∣di,wj)n(di)\theta_{i,k}=\frac{\sum_{j=1}^Nn(d_i,w_j)P(z_k|d_i,w_j)}{n(d_i)}θi,k=n(di)∑j=1Nn(di,wj)P(zk∣di,wj)

综上，在pLSA 中：

由于P(wj∣zk)P(w_j|z_k)P(wj∣zk)和 P(zk∣di)P(z_k|d_i)P(zk∣di)未知，所以我们用EM算法去估计 θ=(P(wj∣zk),P(zk∣di))\theta = (P(w_j|z_k),P(z_k|d_i))θ=(P(wj∣zk),P(zk∣di)) 这个参数的值。
而后，用 ϕk,j\phi_{k,j}ϕk,j表示词项 wjw_jwj 出现在主题 zkz_kzk 中的概率，即 P(wj∣zk)=ϕk,jP(w_j|z_k) = \phi_{k,j}P(wj∣zk)=ϕk,j，用 θi,k\theta_{i,k}θi,k表示主题 zkz_kzk 出现在文档 did_idi 中的概率，即P(zk∣di)=θi,kP(z_k|d_i) = \theta_{i,k}P(zk∣di)=θi,k，从而把 P(wj∣zk)P(w_j|z_k)P(wj∣zk)转换成了“主题-词项”矩阵 Φ\PhiΦ（主题生成词），把 P(zk∣di)P(z_k|d_i)P(zk∣di) 转换成了“文档-主题”矩阵Θ\ThetaΘ（文档生成主题）。
最终求解出 ϕk,j\phi_{k,j}ϕk,j 、θi,k\theta_{i,k}θi,k。

4.5 `LDA` 模型

事实上，理解了pLSA模型，也就差不多快理解了LDA模型，因为LDA就是在pLSA的基础上加层贝叶斯框架，即 LDA 就是 pLSA 的贝叶斯版本（正因为LDA被贝叶斯化了，所以才需要考虑历史先验知识，才加的两个先验参数）。

（1）pLSA跟LDA的对比：生成文档与参数估计

在pLSA模型中，我们按照如下的步骤得到“文档-词项”的生成模型：

按照概率 P(di)P(d_i)P(di) 选择一篇文档 did_idi
选定文档 did_idi 后，确定文章的主题分布
从主题分布中按照概率 P(zk∣di)P(z_k|d_i)P(zk∣di) 选择一个隐含的主题类别 zkz_kzk
选定zkz_kzk 后，确定主题下的词分布
从词分布中按照概率 P(wj∣zk)P(w_j|z_k)P(wj∣zk)选择一个词 wjw_jwj

下面，咱们对比下本文开头所述的LDA模型中一篇文档生成的方式是怎样的：

按照先验概率 P(di)P(d_i)P(di) 选择一篇文档 did_idi
从狄利克雷分布（即Dirichlet分布）α\alphaα 中取样生成文档 did_idi 的主题分布 θi\theta_iθi，换言之，主题分布 θi\theta_iθi由超参数为 α\alphaα 的 Dirichlet分布 生成
从主题的多项式分布 θi\theta_iθi 中取样生成文档 did_idi 第 jjj个词的主题 zi,jz_{i,j}zi,j
从狄利克雷分布（即Dirichlet分布） β\betaβ 中取样生成主题 zi,jz_{i,j}zi,j 对应的词语分布 ϕzi,j\phi_{z_{i,j}}ϕzi,j，换言之，词语分布ϕzi,j\phi_{z_{i,j}}ϕzi,j 由参数为 β\betaβ 的 Dirichlet分布 生成
从词语的多项式分布 ϕzi,j\phi_{z_{i,j}}ϕzi,j 中采样最终生成词语 wi,jw_{i,j}wi,j

从上面两个过程可以看出，LDA在PLSA的基础上，为主题分布和词分布分别加了两个 Dirichlet先验。

继续拿之前讲解PLSA的例子进行具体说明

如前所述，在PLSA中，选主题和选词都是两个随机的过程，先从主题分布 {教育：0.5，经济：0.3，交通：0.2} 中抽取出主题：教育，然后从该主题对应的词分布 {大学：0.5，老师：0.3，课程：0.2} 中抽取出词：大学。

而在LDA中，选主题和选词依然都是两个随机的过程，依然可能是先从主题分布 {教育：0.5，经济：0.3，交通：0.2} 中抽取出主题：教育，然后再从该主题对应的词分布 {大学：0.5，老师：0.3，课程：0.2} 中抽取出词：大学。

（2）那PLSA跟LDA的区别在于什么地方呢？

PLSA中，主题分布和词分布是唯一确定的，能明确的指出主题分布可能就是{教育：0.5，经济：0.3，交通：0.2}，词分布可能就是{大学：0.5，老师：0.3，课程：0.2}。
但在LDA中，主题分布和词分布不再唯一确定不变，即无法确切给出。例如主题分布可能是{教育：0.5，经济：0.3，交通：0.2}，也可能是{教育：0.6，经济：0.2，交通：0.2}，到底是哪个我们不再确定（即不知道），因为它是随机的可变化的。但再怎么变化，也依然服从一定的分布，即主题分布跟词分布由Dirichlet先验随机确定。

前面我们提到的是，有很多个主题或词，各个主题或词被抽中的概率不一样，所以抽取主题或词是随机抽取了。

但是现在主题分布和词分布也不确定了，也就是说，在我们抽取主题的时候，之前是从一个固定的主题分布中，现在是这个主题分布也是不确定的，可能是{教育：0.5，经济：0.3，交通：0.2}，也可能是{教育：0.6，经济：0.2，交通：0.2}。正因为LDA是PLSA的贝叶斯版本，所以主题分布跟词分布本身由先验知识随机给定。

pLSA中：

主题分布和词分布确定后，以一定的概率 (P(zk,di),P(wj∣zk))(P(z_k,d_i),P(w_j|z_k))(P(zk,di),P(wj∣zk)) 分别选取具体的主题和词项，生成好文档。而后根据生成好的文档反推其主题分布、词分布时，最终用EM算法（极大似然估计思想）求解出了两个未知但固定的参数的值：ϕk,j\phi_{k,j}ϕk,j（由P(wj∣zk)P(w_j|z_k)P(wj∣zk)转换而来）和 θi,k\theta_{i,k}θi,k（由 (P(zk,di)(P(z_k,d_i)(P(zk,di)转换而来）。

文档 ddd 产生主题 zzz 的概率，主题 zzz 产生单词 www 的概率都是两个固定的值。

举个文档 ddd 产生主题 zzz 的例子:

给定一篇文档 ddd，主题分布是一定的，比如{P(zi∣d),i=1,2,3}\{P(z_i|d), i = 1,2,3\}{P(zi∣d),i=1,2,3}可能就是{0.4,0.5,0.1}\{0.4,0.5,0.1\}{0.4,0.5,0.1}，表示z1、z2、z3z_1、z_2、z_3z1、z2、z3，这3个主题被文档 ddd 选中的概率都是个固定的值：P(z1∣d)=0.4、P(z2∣d)=0.5、P(z3∣d)=0.1P(z_1|d) = 0.4、P(z_2|d) = 0.5、P(z_3|d) = 0.1P(z1∣d)=0.4、P(z2∣d)=0.5、P(z3∣d)=0.1，如下图所示：

LDA中：

我们不再认为主题分布（各个主题在文档中出现的概率分布）和词分布（各个词语在某个主题下出现的概率分布）是唯一确定的（而是随机变量），而是有很多种可能。但一篇文档总得对应一个主题分布和一个词分布吧，怎么办呢？LDA为它们弄了两个Dirichlet先验参数，这个Dirichlet先验为某篇文档随机抽取出某个主题分布和词分布。

文档 ddd 产生主题 zzz（准确的说，其实是Dirichlet先验为文档 ddd 生成主题分布Θ\ThetaΘ，然后根据主题分布Θ\ThetaΘ产生主题 zzz）的概率，主题 zzz 产生单词 www 的概率都不再是某两个确定的值，而是随机变量。

举下文档 ddd 具体产生主题 zzz 的例子:

给定一篇文档 ddd，现在有多个主题z1、z2、z3z_1、z_2、z_3z1、z2、z3，它们的主题分布 {P(zi∣d),i=1,2,3}\{ P(z_i|d), i = 1,2,3\}{P(zi∣d),i=1,2,3}可能是{0.4,0.5,0.1}\{0.4,0.5,0.1\}{0.4,0.5,0.1}，也可能是{0.2,0.2,0.6}\{0.2,0.2,0.6\}{0.2,0.2,0.6}，即这些主题被 ddd 选中的概率都不再认为是确定的值，可能是P(z1∣d)=0.4、P(z2∣d)=0.5、P(z3∣d)=0.1P(z_1|d) = 0.4、P(z_2|d) = 0.5、P(z_3|d) = 0.1P(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.6P(z_1|d) = 0.2、P(z_2|d) = 0.2、P(z_3|d) = 0.6P(z1∣d)=0.2、P(z2∣d)=0.2、P(z3∣d)=0.6 等等，而主题分布到底是哪个取值集合我们不确定（为什么？这就是贝叶斯派的核心思想，把未知参数当作是随机变量，不再认为是某一个确定的值），但其先验分布是dirichlet 分布，所以可以从无穷多个主题分布中按照dirichlet 先验随机抽取出某个主题分布出来。如下图所示：

换言之，LDA在pLSA的基础上给这两参数 (P(zk∣di),P(wj∣zk))(P(z_k|d_i),P(w_j|z_k))(P(zk∣di),P(wj∣zk)) 加了两个先验分布的参数（贝叶斯化）：一个主题分布的先验分布Dirichlet分布 α\alphaα ，和一个词语分布的先验分布Dirichlet分布 β\betaβ。

综上，LDA真的只是pLSA的贝叶斯版本，文档生成后，两者都要根据文档去推断其主题分布和词语分布（即两者本质都是为了估计给定文档生成主题，给定主题生成词语的概率），只是用的参数推断方法不同，在pLSA中用极大似然估计的思想去推断两未知的固定参数，而LDA则把这两参数弄成随机变量，且加入dirichlet先验。

所以，pLSA跟LDA的本质区别就在于它们去估计未知参数所采用的思想不同，前者用的是频率派思想，后者用的是贝叶斯派思想。

好比，我去一朋友家：

按照频率派的思想，我估计他在家的概率是1/2，不在家的概率也是1/2，是个定值。
而按照贝叶斯派的思想，他在家不在家的概率不再认为是个定值1/2，而是随机变量。比如按照我们的经验（比如当天周末），猜测他在家的概率是0.6，但这个0.6不是说就是完全确定的，也有可能是0.7。如此，贝叶斯派没法确切给出参数的确定值（0.3,0.4,0.6,0.7,0.8,0.9都有可能），但至少明白在哪个范围或哪些取值（0.6,0.7,0.8,0.9）更有可能，哪个范围或哪些取值（0.3,0.4）不太可能。进一步，贝叶斯估计中，参数的多个估计值服从一定的先验分布，而后根据实践获得的数据（例如周末不断跑他家），不断修正之前的参数估计，从先验分布慢慢过渡到后验分布。

（3）LDA生成文档过程的进一步理解

上面说，LDA中，主题分布 —— 比如{P(zi),i=1,2,3}\{P(z_i), i =1,2,3\}{P(zi),i=1,2,3} 等于{0.4,0.5,0.1}\{0.4,0.5,0.1\}{0.4,0.5,0.1}或{0.2,0.2,0.6}\{0.2,0.2,0.6\}{0.2,0.2,0.6} —— 是由dirichlet先验给定的，不是根据文档产生的。所以，LDA生成文档的过程中，先从dirichlet先验中“随机”抽取出主题分布，然后从主题分布中“随机”抽取出主题，最后从确定后的主题对应的词分布中“随机”抽取出词。

那么，dirichlet先验 到底是如何“随机”抽取主题分布的呢？

事实上，从dirichlet分布 中随机抽取主题分布，这个过程不是完全随机的。为了说清楚这个问题，咱们得回顾下dirichlet分布。事实上，如果我们取3个事件的话，可以建立一个三维坐标系，类似xyz三维坐标系，这里，我们把3个坐标轴弄为p1、p2、p3p_1、p_2、p_3p1、p2、p3，如下图所示：

在这个三维坐标轴所划分的空间里，每一个坐标点(p1,p2,p3)(p_1,p_2,p_3)(p1,p2,p3)就对应着一个主题分布，且某一个点(p1,p2,p3)(p_1,p_2,p_3)(p1,p2,p3) 的大小表示3个主题z1、z2、z3z_1、z_2、z_3z1、z2、z3出现的概率大小（因为各个主题出现的概率和为1，所以p1+p2+p3=1p_1+p_2+p_3 = 1p1+p2+p3=1，且p1、p2、p3p_1、p_2、p_3p1、p2、p3 这3个维度最大取值为1）。比如(p1,p2,p3)=(0.4,0.5,0.1)(p_1,p_2,p_3) = (0.4,0.5,0.1)(p1,p2,p3)=(0.4,0.5,0.1)便对应着主题分布{P(zi),i=1,2,3}={0.4,0.5,0.1}\{ P(z_i), i =1,2,3 \} = \{0.4,0.5,0.1\}{P(zi),i=1,2,3}={0.4,0.5,0.1}

可以想象到，空间里有很多这样的点(p1,p2,p3)(p_1,p_2,p_3)(p1,p2,p3)，意味着有很多的主题分布可供选择，那dirichlet分布如何选择主题分布呢？把上面的斜三角形放倒，映射到底面的平面上，便得到如下所示的一些彩图（3个彩图中，每一个点对应一个主题分布，高度代表某个主题分布被dirichlet分布选中的概率，且选不同的dirichlet 分布会偏向不同的主题分布）：

我们来看上图中左边这个图，高度就是代表dirichlet分布选取某个坐标点(p1,p2,p3)(p_1,p_2,p_3)(p1,p2,p3)（这个点就是一个主题分布）的概率大小。如下图所示，平面投影三角形上的三个顶点上的点：A=(0.9,0.05,0.05)、B=(0.05,0.9,0.05)、C=(0.05,0.05,0.9)A=(0.9,0.05,0.05)、B=(0.05,0.9,0.05)、C=(0.05,0.05,0.9)A=(0.9,0.05,0.05)、B=(0.05,0.9,0.05)、C=(0.05,0.05,0.9)各自对应的主题分布被dirichlet分布选中的概率值很大，而平面三角形内部的两个点：D、E对应的主题分布被dirichlet分布选中的概率值很小。

所以虽然说dirichlet分布是随机选取任意一个主题分布的，但依然存在着P(A)=P(B)=P(C)>>P(D)=P(E)P(A) = P(B) = P(C) >> P(D) = P(E)P(A)=P(B)=P(C)>>P(D)=P(E)，即dirichlet分布还是“偏爱”某些主题分布的。至于dirichlet分布的参数是如何决定dirichlet分布的形状的，可以从dirichlet分布的定义和公式思考。

此外，就算说“随机”选主题也是根据主题分布来“随机”选取，这里的随机不是完全随机的意思，而是根据各个主题出现的概率值大小来抽取。比如当dirichlet先验为文档 ddd 生成的主题分布{P(zi),i=1,2,3}\{ P(z_i), i =1,2,3 \}{P(zi),i=1,2,3}是{0.4,0.5,0.1}\{0.4,0.5,0.1\}{0.4,0.5,0.1}时，那么主题z2z_2z2在文档 ddd 中出现的概率便是0.5。所以，从主题分布中抽取主题，这个过程也不是完全随机的，而是按照各个主题出现的概率值大小进行抽取。

五、LDA参数估计：Gibbs采样

在pLSA 中，我们使用 EM算法 去估计“主题-词项”矩阵 Φ\PhiΦ（由 P(wj∣zk)P(w_j|z_k)P(wj∣zk) 转换得到）和“文档-主题”矩阵Θ\ThetaΘ（由 P(zk,di)P(z_k,d_i)P(zk,di) 转换得到）这两个参数，而且这两参数都是个固定的值，只是未知，使用的思想其实就是极大似然估计MLE。

那么在 LDA 中，又是如何估计的呢？

在 LDA 中，估计 Φ\PhiΦ 和 Θ\ThetaΘ 这两未知参数可以用 变分(Variational inference)-EM算法，也可以用gibbs采样，前者的思想是最大后验估计MAP（MAP与MLE类似，都把未知参数当作固定的值），后者的思想是贝叶斯估计。贝叶斯估计是对MAP的扩展，但它与MAP有着本质的不同，即贝叶斯估计把待估计的参数看作是服从某种先验分布的随机变量。

那么什么是贝叶斯估计呢？其实前面我们有介绍过，这里我们再举个例子：

假设中国的大学只有两种：理工科和文科，这两种学校数量的比例是1:1，其中，理工科男女比例7:1，文科男女比例1:7。某天你被外星人随机扔到一个校园，问你该学校可能的男女比例是多少？然后，你实际到该校园里逛了一圈，看到的5个人全是男的，这时候再次问你这个校园的男女比例是多少？

因为刚开始时，有先验知识，所以该学校的男女比例要么是7:1，要么是1:7，即 P(比例为7:1)=1/2，P(比例为1:7)=1/2P(比例为7:1) = 1/2，P(比例为1:7) = 1/2P(比例为7:1)=1/2，P(比例为1:7)=1/2。
然后看到5个男生后重新估计男女比例，其实就是求P(比例7:1∣5个男生）=？，P(比例1:7∣5个男生)=？P(比例7:1|5个男生）= ？，P(比例1:7|5个男生) = ？P(比例7:1∣5个男生）=？，P(比例1:7∣5个男生)=？
用贝叶斯公式 :
P(A∣B)=P(B∣A)P(A)P(B)P(A|B)=\frac{P(B|A)P(A)}{P(B)}P(A∣B)=P(B)P(B∣A)P(A)
可得：
P(比例7:1∣5个男生)=P(比例7:1)∗P(5个男生∣比例7:1)P(5个男生)P(比例7:1|5个男生) = \frac{P(比例7:1)*P(5个男生|比例7:1)} {P(5个男生)}P(比例7:1∣5个男生)=P(5个男生)P(比例7:1)∗P(5个男生∣比例7:1)

P(5个男生)P(5个男生)P(5个男生) 是5个男生的先验概率，与学校无关，所以是个常数；

类似的：
P(比例1:7∣5个男生)=P((比例1:7)∗P(5个男生∣比例1:7)P(5个男生)P(比例1:7|5个男生) = \frac{P((比例1:7)*P(5个男生|比例1:7)}{P(5个男生)}P(比例1:7∣5个男生)=P(5个男生)P((比例1:7)∗P(5个男生∣比例1:7)
最后将上述两个等式比一下，可得：
P(比例7:1∣5个男生)P(比例1:7∣5个男生)=P((比例7:1)∗P(5个男生∣比例7:1)P(比例1:7)∗P(5个男生∣比例1:7)\frac{P(比例7:1|5个男生)}{P(比例1:7|5个男生)} = \frac{P((比例7:1)*P(5个男生|比例7:1)}{ P(比例1:7)*P(5个男生|比例1:7)}P(比例1:7∣5个男生)P(比例7:1∣5个男生)=P(比例1:7)∗P(5个男生∣比例1:7)P((比例7:1)∗P(5个男生∣比例7:1)

由于LDA 把要估计的主题分布和词分布看作是其先验分布是Dirichlet分布的随机变量，所以，在LDA这个估计主题分布、词分布的过程中，它们的先验分布（即Dirichlet分布）事先由人为给定，那么LDA就是要去求它们的后验分布（LDA中可用 gibbs采样 去求解它们的后验分布，得到期望 θ^mk,φ^kt\hat{\theta}_{mk},\hat{\varphi}_{kt}θ^mk,φ^kt）

在LDA中，主题分布和词分布本身都是多项分布，而前面介绍的可知 “Dirichlet分布是多项式分布的共轭先验概率分布”，因此选择Dirichlet 分布作为它们的共轭先验分布。意味着为多项分布的参数p选取的先验分布是Dirichlet分布，那么以p为参数的多项分布用贝叶斯估计得到的后验分布仍然是Dirichlet分布。

5.1 Gibbs采样流程

类似于 pLSA，LDA 的原始论文中是用的 变分-EM算法估计 未知参数，后来发现另一种估计 LDA 未知参数的方法更好，这种方法就是：Gibbs Sampling，有时叫Gibbs采样或Gibbs抽样，都一个意思。Gibbs抽样 是马尔可夫链蒙特卡尔理论（MCMC）中用来获取一系列近似等于指定多维概率分布（比如2个或者多个随机变量的联合概率分布）观察样本的算法。

对于一个 nnn 维的概率分布 π(x1,x2,...,xn)\pi(x_1,x_2,...,x_n)π(x1,x2,...,xn)，可以通过在 nnn 个坐标上轮换采样，来得到新的样本，对于轮换到任意一个坐标 xix_ixi 上的转移，马尔科夫链的的状态转移概率为：p(xi∣x1,x2,...,xi−1,xi+1,...,xn)p(x_i|x_1,x_2,...,x_{i-1},x_{i+1},...,x_n)p(xi∣x1,x2,...,xi−1,xi+1,...,xn)，即固定 n−1n-1n−1 个坐标轴，在某一个坐标上移动。

Gibbs采样算法在高维空间采样的时候具有比较高的优势，Gibbs采样的过程比较类似这个坐标轴下降法。

Gibbs采样算法流程如下：

输入稳定的分布 π(x1,x2,...,xn)\pi(x_1,x_2,...,x_n)π(x1,x2,...,xn) 或者对应特征的条件概率分布，设定状态转移次数阀值 n1n_1n1，需要的样本数 n2n_2n2；
随机初始化状态值 (x11,x21,...,xn1)(x_1^1,x_2^1,...,x_n^1)(x11,x21,...,xn1)；
进行迭代数据采样（迭代 n1+n+2−1n_1+n+2-1n1+n+2−1 次）；

从条件概率分布中采样得到对应的样本：
xjt+1→p(xj∣x1t+1,x2t+1,...,xj−1t+1,xj+1t,...,xnt)x_j^{t+1} \rightarrow p(x_j|x_1^{t+1},x_2^{t+1},...,x_{j-1}^{t+1},x_{j+1}^t,...,x_n^t)xjt+1→p(xj∣x1t+1,x2t+1,...,xj−1t+1,xj+1t,...,xnt)

最终得到样本集

{(x1n1,x2n1,...,xnn1),...,(x1n1+n2−1,x2n1+n2−1,...,xnn1+n2−1)}\Big\{\big(x_1^{n_1},x_2^{n_1},...,x_n^{n_1}\big),...,\big(x_1^{n_1+n_2-1},x_2^{n_1+n_2-1},...,x_n^{n_1+n_2-1}\big)\Big\}{(x1n1,x2n1,...,xnn1),...,(x1n1+n2−1,x2n1+n2−1,...,xnn1+n2−1)}

5.2 LDA参数学习-Gibbs采样

给定一个文档集合，www 是可以观察到的已知变量，α\alphaα 和 β\betaβ 是根据经验给定的先验参数，其他的变量 zzz，θ\thetaθ 和 ϕ\phiϕ 都是未知的隐含变量，需要根据观察到的变量来学习估计的。

具体来讲，所有文档联合起来形成的词向量 www 是已知数据，但是不知道语料库的主题 zzz 的分布。假设可以先求解出 w、zw、zw、z 的联合分布 p(w,z)p(w,z)p(w,z)，进而就可以求出某个词 wiw_iwi对应主题特征 ziz_izi的条件概率分布p(zi=k∣w,z−i)p(z_i=k|w,z_{-i})p(zi=k∣w,z−i)，其中 ziz_izi表示去掉下标为 iii 后的主题分布，有了条件概率，那么就可以使用 Gibbs采样，最终可以得到第 iii 个词的主题。

如果通过采样得到所有词的主题，那么可以通过统计所有词的主题数，从而得到各个主题的词分布。接着统计各个文档对应词的主题数，从而可以得到各个文档的主题分布。

（1）LDA 的图模型如下:

因为 α\alphaα 产生主题分布 θi\theta_iθi，主题分布θi\theta_iθi确定具体主题 zi,jz_{i,j}zi,j，且 β\betaβ 产生词分布 ϕzi,j\phi_{z_{i,j}}ϕzi,j、词分布 ϕzi,j\phi_{z_{i,j}}ϕzi,j 确定具体词 wi,jw_{i,j}wi,j，所以联合概率分布可表示为：

p(w⃗,z⃗∣α⃗,β⃗)=p(w⃗∣z⃗,β⃗)p(z⃗∣α⃗)p(\vec{w},\vec{z}|\vec{\alpha},\vec{\beta}) = p(\vec{w}|\vec{z},\vec{\beta})p(\vec{z}|\vec{\alpha})p(w,z∣α,β)=p(w∣z,β)p(z∣α)

其中

第一项 p(z⃗∣α⃗)p(\vec{z}|\vec{\alpha})p(z∣α) 表示根据主题分布的先验分布参数 α\alphaα 采样主题的过程

第二项 p(w⃗∣z⃗,β⃗)p(\vec{w}|\vec{z},\vec{\beta})p(w∣z,β) 表示的是根据确定的主题 z⃗\vec{z}z 和词分布的先验分布参数 β\betaβ 采样词的过程；
这两项因子是需要计算的两个未知参数。

由于这两个过程是独立的，所以下面可以分别处理，各个击破。

第一项：p(z⃗∣α⃗)p(\vec{z}|\vec{\alpha})p(z∣α)

可以根据主题分布的先验分布参数 α\alphaα 采样得到的主题分布 θ\thetaθ 产生：

p(z⃗∣θ)=∏i=1Wp(zi∣di)=∏m=1M∏k=1Kp(zi=k∣di=m)=∏m=1M∏k=1Kθm,knm(k)p(\vec{z}|\theta) = \prod_{i=1}^Wp(z_i|d_i)=\prod_{m=1}^M\prod_{k=1}^Kp(z_i=k|d_i=m) =\prod_{m=1}^M\prod_{k=1}^K\theta_{m,k}^{n_m^{(k)}} p(z∣θ)=i=1∏Wp(zi∣di)=m=1∏Mk=1∏Kp(zi=k∣di=m)=m=1∏Mk=1∏Kθm,knm(k)

其中：

did_idi 表示单词 iii 所属的文档；

nm(k)n_m^{(k)}nm(k) 表示主题 kkk 在文章 mmm 中出现的次数

对主题分布 θ\thetaθ 积分：

p(z⃗∣α⃗)=∫p(z⃗∣θ)p(θ∣α⃗)dθp(\vec{z}|\vec{\alpha}) = \int p(\vec{z}|\theta)p(\theta|\vec{\alpha})d\thetap(z∣α)=∫p(z∣θ)p(θ∣α)dθ
=∫∏m=1M1Δ(α⃗)∏k=1Kθm,knm(k)+αk−1dθ⃗m= \int \prod_{m=1}^M\frac{1}{\Delta(\vec{\alpha})}\prod_{k=1}^K\theta_{m,k}^{n_m^{(k)}+\alpha_k-1}d\vec{\theta}_m=∫m=1∏MΔ(α)1k=1∏Kθm,knm(k)+αk−1dθm
=∏m=1MΔ(n⃗m+α⃗)Δ(α⃗)=\prod_{m=1}^M\frac{\Delta(\vec{n}_m + \vec{\alpha})}{\Delta(\vec{\alpha})}=m=1∏MΔ(α)Δ(nm+α)

其中：

n⃗m={nm(1),nm(2),...,nm(k)}\vec{n}_m = \{n_m^{(1)},n_m^{(2)},...,n_m^{(k)}\}nm={nm(1),nm(2),...,nm(k)}

第二项：p(w⃗∣z⃗,β⃗)p(\vec{w}|\vec{z},\vec{\beta})p(w∣z,β)

可以根据确定的主题 z⃗\vec{z}z 和从先验分布参数为 β\betaβ 的词分布 ϕ\phiϕ 产生：

p(w⃗∣z⃗,ϕ)=∏i=1Wp(wi∣zi)=∏i=1Wφzi,wip(\vec{w}|\vec{z},\phi) = \prod_{i=1}^Wp(w_i|z_i) = \prod_{i=1}^W\varphi_{z_i,w_i}p(w∣z,ϕ)=i=1∏Wp(wi∣zi)=i=1∏Wφzi,wi

由于样本中的词服从参数为 β\betaβ，主题为 ziz_izi 的独立多项分布，这意味着可以把上面对词的乘积分解成分别对主题和对词的两层乘积：

p(w⃗∣z⃗,ϕ)=∏k=1K∏t=1Vp(wi=t∣zi=k)=∏k=1K∏t=1Vφk,tnk(t)p(\vec{w}|\vec{z},\phi) = \prod_{k=1}^K\prod_{t=1}^Vp(w_i = t|z_i=k) = \prod_{k=1}^K\prod_{t=1}^V \varphi_{k,t}^{n_k^{(t)}}p(w∣z,ϕ)=k=1∏Kt=1∏Vp(wi=t∣zi=k)=k=1∏Kt=1∏Vφk,tnk(t)

其中：

nk(t)n_k^{(t)}nk(t) 表示词 ttt 在主题 kkk 中出现的次数

目标函数p(w⃗∣z⃗,β⃗)p(\vec{w}|\vec{z},\vec{\beta})p(w∣z,β)需要对词分布 ϕ\phiϕ 积分，并且结合Dirichlet 分布的归一化系数Δ(α⃗)\Delta(\vec{\alpha})Δ(α) 的公式：

Δ(α⃗)=∫∏k=1Vpkαk−1dp⃗\Delta(\vec{\alpha}) = \int \prod_{k=1}^Vp_k^{\alpha_k-1}d\vec{p}Δ(α)=∫k=1∏Vpkαk−1dp

于是：

p(w⃗∣z⃗,β⃗)=∫p(w⃗∣z⃗,ϕ)p(ϕ∣β⃗)dϕp(\vec{w}|\vec{z},\vec{\beta}) = \int p(\vec{w}|\vec{z},\phi)p(\phi|\vec{\beta})d \phip(w∣z,β)=∫p(w∣z,ϕ)p(ϕ∣β)dϕ
=∫∏z=1K1Δ(β⃗)∏t=1Vφz,tnz(t)+βt−1dφ⃗z= \int \prod_{z=1}^K \frac{1}{\Delta(\vec{\beta})}\prod_{t=1}^V\varphi_{z,t}^{n_z^{(t)} +\beta_t-1}d\vec{\varphi}_z=∫z=1∏KΔ(β)1t=1∏Vφz,tnz(t)+βt−1dφz
=∏z=1KΔ(n⃗z+β⃗)Δ(β⃗)=\prod_{z=1}^K\frac{\Delta(\vec{n}_z + \vec{\beta})}{\Delta(\vec{\beta})}=z=1∏KΔ(β)Δ(nz+β)

其中：

n⃗z={nz(1),nz(2),...,nz(V)}\vec{n}_z =\{n_z^{(1)},n_z^{(2)},...,n_z^{(V)}\}nz={nz(1),nz(2),...,nz(V)}

nz(t)n_z^{(t)}nz(t) 表示词 ttt 在主题 zzz 中出现的次数

结合第一项和第二项，最终得到主题和词向量的联合分布为：

p((w⃗,z⃗)∝p(w⃗,z⃗∣α⃗,β⃗)=p(w⃗∣z⃗,β⃗)p(z⃗∣α⃗)=∏z=1KΔ(n⃗z+β⃗)Δ(β⃗)×∏m=1MΔ(n⃗m+α⃗)Δ(α⃗)p((\vec{w},\vec{z}) \propto p(\vec{w},\vec{z}|\vec{\alpha},\vec{\beta}) = p(\vec{w}|\vec{z},\vec{\beta})p(\vec{z}|\vec{\alpha})=\prod_{z=1}^K\frac{\Delta(\vec{n}_z + \vec{\beta})}{\Delta(\vec{\beta})}\times \prod_{m=1}^M\frac{\Delta(\vec{n}_m + \vec{\alpha})}{\Delta(\vec{\alpha})}p((w,z)∝p(w,z∣α,β)=p(w∣z,β)p(z∣α)=z=1∏KΔ(β)Δ(nz+β)×m=1∏MΔ(α)Δ(nm+α)

基于联合分布，就可以使用求解 Gibbs采样 所需要的条件分布p(zi=k∣w,z¬i)p(z_i=k|w,z_{\neg i})p(zi=k∣w,z¬i)；对于下标 iii，由于它对应的词wiw_iwi 是可以观察到的，因此有公式如下：

p(zi=k∣w⃗,z⃗¬i)∝p(zi=k,wi=t∣w⃗¬i,z⃗¬i)p(z_i=k|\vec{w},\vec{z}_{\neg i}) \propto p(z_i=k,w_i=t|\vec{w}_{\neg i},\vec{z}_{\neg i})p(zi=k∣w,z¬i)∝p(zi=k,wi=t∣w¬i,z¬i)

其中：

w⃗¬i\vec{w}_{\neg i}w¬i 表示除去下标为 iii 的词；

z⃗¬i\vec{z}_{\neg i}z¬i 表示除去下标为 iii 的主题

对于 zi=k,wi=tz_i=k,w_i=tzi=k,wi=t ，只涉及到第 ddd 篇文档和第 kkk 个主题两个 Dirichlet共轭，即：

α⃗→θ⃗d→z⃗d\vec{\alpha} \rightarrow \vec{\theta}_d \rightarrow \vec{z}_dα→θd→zd
β⃗→φ⃗k→w⃗(k)\vec{\beta} \rightarrow \vec{\varphi}_k \rightarrow \vec{w}_{(k)}β→φk→w(k)

至于其他的 Dirichlet共轭 和这两个是互相独立的，也就是说从语料库中去掉 ziz_izi 和 wiw_iwi 后，并不会改变共轭结构。所以对应的后验分布为：

p(θ⃗d∣w⃗¬i,z⃗¬i)=Dir(θ⃗d∣n⃗d,¬i+α⃗)p(\vec{\theta}_d| \vec{w}_{\neg i},\vec{z}_{\neg i}) = Dir(\vec{\theta}_d | \vec{n}_{d, \neg i} + \vec{\alpha})p(θd∣w¬i,z¬i)=Dir(θd∣nd,¬i+α)
p(φ⃗k∣w⃗¬i,z⃗¬i)=Dir(φ⃗k∣n⃗k,¬i+β⃗)p(\vec{\varphi}_k| \vec{w}_{\neg i},\vec{z}_{\neg i}) = Dir(\vec{\varphi}_k | \vec{n}_{k, \neg i} + \vec{\beta})p(φk∣w¬i,z¬i)=Dir(φk∣nk,¬i+β)

计算Gibbs采样的条件概率：

p(zi=k∣w⃗,z⃗¬i)∝p(zi=k,wi=t∣w⃗¬i,z⃗¬i)p(z_i=k|\vec{w},\vec{z}_{\neg i}) \propto p(z_i=k,w_i=t|\vec{w}_{\neg i},\vec{z}_{\neg i})p(zi=k∣w,z¬i)∝p(zi=k,wi=t∣w¬i,z¬i)
=∫p(zi=k,wi=t,θ⃗d,φ⃗k∣w⃗¬i,z⃗¬i)dθ⃗ddφ⃗k= \int p(z_i=k,w_i=t,\vec{\theta}_d,\vec{\varphi}_k|\vec{w}_{\neg i},\vec{z}_{\neg i})d\vec{\theta}_dd\vec{\varphi}_k=∫p(zi=k,wi=t,θd,φk∣w¬i,z¬i)dθddφk
=∫p(zi=k,θ⃗d∣w⃗¬i,z⃗¬i)p(wi=t,φ⃗k∣w⃗¬i,z⃗¬i)dθ⃗ddφ⃗k=\int p(z_i=k,\vec{\theta}_d|\vec{w}_{\neg i},\vec{z}_{\neg i})p(w_i=t,\vec{\varphi}_k|\vec{w}_{\neg i},\vec{z}_{\neg i})d\vec{\theta}_dd\vec{\varphi}_k=∫p(zi=k,θd∣w¬i,z¬i)p(wi=t,φk∣w¬i,z¬i)dθddφk
=∫p(zi=k∣θ⃗d)p(θ⃗d∣w⃗¬i,z⃗¬i)p(wi=t∣φ⃗k)p(φ⃗k∣w⃗¬i,z⃗¬i)dθ⃗ddφ⃗k=\int p(z_i=k|\vec{\theta}_d)p(\vec{\theta}_d|\vec{w}_{\neg i},\vec{z}_{\neg i})p(w_i=t|\vec{\varphi}_k)p(\vec{\varphi}_k|\vec{w}_{\neg i},\vec{z}_{\neg i})d\vec{\theta}_dd\vec{\varphi}_k=∫p(zi=k∣θd)p(θd∣w¬i,z¬i)p(wi=t∣φk)p(φk∣w¬i,z¬i)dθddφk
=∫p(zi=k∣θ⃗d)Dir(θ⃗d∣n⃗d,¬i+α⃗)dθ⃗d∗∫p(wi=t∣φ⃗k)Dir(φ⃗k∣n⃗k,¬i+β⃗)dφ⃗k=\int p(z_i=k|\vec{\theta}_d)Dir(\vec{\theta}_d|\vec{n}_{d, \neg i} + \vec{\alpha})d\vec{\theta}_d *\int p(w_i=t|\vec{\varphi}_k)Dir(\vec{\varphi}_k | \vec{n}_{k, \neg i} + \vec{\beta})d\vec{\varphi}_k=∫p(zi=k∣θd)Dir(θd∣nd,¬i+α)dθd∗∫p(wi=t∣φk)Dir(φk∣nk,¬i+β)dφk
=∫θdkDir(θ⃗d∣n⃗d,¬i+α⃗)dθ⃗d∫φktDir(φ⃗k∣n⃗k,¬i+β⃗)dφ⃗k=\int \theta_{dk}Dir(\vec{\theta}_d|\vec{n}_{d, \neg i} + \vec{\alpha})d\vec{\theta}_d \int \varphi_{kt}Dir(\vec{\varphi}_k | \vec{n}_{k, \neg i} + \vec{\beta})d\vec{\varphi}_k=∫θdkDir(θd∣nd,¬i+α)dθd∫φktDir(φk∣nk,¬i+β)dφk
=EDir(θd)(θdk)EDir(φk)(φkt)=E_{Dir(\theta_d)}(\theta_{dk})E_{Dir(\varphi_k)}(\varphi_{kt})=EDir(θd)(θdk)EDir(φk)(φkt)

Dirichlet分布的期望公式如下：

EDir(θd)(θdk)=nd,¬ik+αk∑s=1Knd,¬is+αsE_{Dir(\theta_d)}(\theta_{dk}) = \frac{n_{d, \neg i}^{k}+\alpha_k}{\sum_{s=1}^Kn_{d, \neg i}^{s}+\alpha_s}EDir(θd)(θdk)=∑s=1Knd,¬is+αsnd,¬ik+αk

EDir(φk)(φkt)=nk,¬it+βt∑f=1Vnk,¬if+βfE_{Dir(\varphi_k)}(\varphi_{kt}) = \frac{n_{k, \neg i}^{t}+\beta_t}{\sum_{f=1}^Vn_{k, \neg i}^{f}+\beta_f}EDir(φk)(φkt)=∑f=1Vnk,¬if+βfnk,¬it+βt

带入条件概率中，可以得到最终的条件概率公式：

p(zi=k∣w⃗,z⃗¬i)=nd,¬ik+αk∑s=1Knd,¬is+αs⋅nk,¬it+βt∑f=1Vnk,¬if+βfp(z_i=k|\vec{w},\vec{z}_{\neg i}) = \frac{n_{d, \neg i}^{k}+\alpha_k}{\sum_{s=1}^Kn_{d, \neg i}^{s}+\alpha_s} \cdot \frac{n_{k, \neg i}^{t}+\beta_t}{\sum_{f=1}^Vn_{k, \neg i}^{f}+\beta_f}p(zi=k∣w,z¬i)=∑s=1Knd,¬is+αsnd,¬ik+αk⋅∑f=1Vnk,¬if+βfnk,¬it+βt

仔细观察上述结果，可以发现，式子的右半部分便是 p(topic∣doc)⋅(word∣topic)p(topic|doc)\cdotp(word|topic)p(topic∣doc)⋅(word∣topic)，这个概率的值对应着 doc→topic→worddoc \rightarrow topic \rightarrow worddoc→topic→word 的路径概率。如此，KKK 个 topictopictopic 对应着 KKK 条路径，Gibbs Sampling 便在这 KKK 条路径中进行采样，如下图所示：

就这样，Gibbs Sampling通过求解出主题分布和词分布的后验分布，从而成功解决主题分布和词分布这两参数未知的问题。

5.3 Gibbs采样总结

（1）Gibbs采样训练流程

选择合适的主题数 KKK ，选择合适的超参数 α,β\alpha,\betaα,β
对于语料库中每一篇文档的每一个词，随机的赋予一个主题编号 zzz
重新扫描语料库，对于每一个词，利用 Gibbs采样 公式更新它的 topictopictopic 的编号，并更新语料库中该词的编号
重复第三步中基于坐标轴轮换的Gibbs采样，直到 Gibbs采样 收敛。
统计语料库中各个文档各个词的主题，得到文档主题分布；然后统计语料库中各个主题词的分布，得到主题与词的分布。

（2）Gibbs采样预测流程

对当前文档的每一个词，随机的赋予一个主题编号 zzz
重新扫描当前文档，对于每一个词，利用 Gibbs采样 算法更新它的 topictopictopic 编号
重复第二步的基于坐标轴轮换的Gibbs采样，直到Gibbs采样收敛
统计文档中各个词的主题，得到该文档主题分布。

六、总结

本篇可谓是长篇大作，由于博主也在学习中，所以其中难免有些错误，对于一些细节性的，博主也会在后续的深入学习中进一步改进，本篇主要是想记录学习 LDA 一个过程，按照本篇顺序一一深入，我觉得是学懂 LDA 的一个很好的方式。

本篇主要内容如下：

1. 介绍什么是主题模型；
1. 初步了解LDA模型
1. 相关数学知识点介绍
1. 由基础模型到LDA模型（重点）
1. LDA模型参数估计：Gibbs采样
1. 案例演示（链接）

内容比较多，可能存在表达不清楚的地方，希望大家多多指教，谢谢！！！