一、简介

隐含狄利克雷分布(LatentDirichletAllocation，简称LDA)是由DavidM.Blei、AndrewY.Ng、MichaelI.Jordan在2003年提出的，是一种词袋模型，它认为文档是一组词构成的集合，词与词之间是无序的。一篇文档可以包含多个主题，文档中的每个词都是由某个主题生成的，LDA给出文档属于每个主题的概率分布，同时给出每个主题上词的概率分布。LDA是一种无监督学习，在文本主题识别、文本分类、文本相似度计算和文章相似推荐等方面都有应用。

本文将从贝叶公式、Gamma函数、二项分布、Beta分布、多项式分布、Dirichlet分布、共轭先验分布、马氏链及其平稳分布、MCMC、GibbsSampling、EM算法、UnigramModel、贝叶斯UnigramModel、PLSA、LDA几方面介绍LDA模型，需要读者具备一定的概率论和微积分知识。

二、基础知识

▌1.1贝叶公式

贝叶斯学派的最基本的观点是：任一个未知量θ都可看作一个随机变量，应该用一个概率分布去描述对θ的未知状况，这个概率分布是在抽样前就有关于θ的先验信息的概率陈述，这个概率分布被称为先验分布。

从贝叶斯观点看，样本

的产生要分两步进行，首先设想从先验分布p(θ)产生一个样本θ'，这一步是“老天爷”做的，人们是看不到的，故用“设想”二字。第二步是从总体分布p(X|θ')产生一个样本

，这个样本是具体的，人们能看得到的，此样本X发生的概率是与如下联合密度函数成正比。

这个联合密度函数是综合了总体信息和样本信息，常称为似然函数，记为L(θ')。

由于θ'是设想出来的，它仍然是未知的，它是按先验分布p(θ)产生的，要把先验信息进行综合，不能只考虑θ'，而应对θ的一切可能加以考虑，故要用p(θ)参与进一步综合，所以样本X和参数θ的联合分布(三种可用的信息都综合进去了)：

我们的任务是要对未知数θ作出统计推断，在没有样本信息时，人们只能根据先验分布对θ作出推断。在有样本观察值

之后，我们应该依据p(X,θ)对θ作出推断，为此我们把p(X,θ)作如下分解：

其中p(X)是X的边缘密度函数。

它与θ无关，p(X)中不含θ的任何信息。因此能用来对θ作出推断的仅是条件分布p(θ|X)：

这就是贝叶斯公式的密度函数形式，在样本X给定下，θ的条件分布被称为θ的后验分布。它是集中了总体、样本和先验等三种信息中有关θ的一切信息，而又是排除一切与θ无关的信息之后得到的结果，故基于后验分布p(θ|X)对θ进行统计推断是更合理的。

一般说来，先验分布p(θ)是反映人们在抽样前对θ的认识，后验分布p(θ|X)是反映人们在抽样后对θ的认识，之间的差异是由于样本的出现后人们对θ认识的一种调整，所以后验分布p(θ|X)可以看作是人们用总体信息和样本信息(抽样信息)对先验分布p(θ)作调整的结果。下面我们介绍三种估计方法：

1.最大似然估计(ML)

最大似然估计是找到参数θ使得样本X的联合概率最大，并不会考虑先验知识，频率学派和贝叶斯学派都承认似然函数，频率学派认为参数θ是客观存在的，只是未知。求参数θ使似然函数最大，ML估计问题可以用下面公式表示：

通常可以令导数为0求得θ的值。ML估计不会把先验知识考虑进去，很容易出现过拟合的现象。我们举个例子，抛一枚硬币，假设正面向上的概率为p，抛了N次，正面出现

次，反面出现

次，c=1表示正面，c=0表示反面，我们用ML估计：

如果

,

，则

，似乎比我们认知的0.5高了很多。

2.最大后验估计(MAP)

MAP是为了解决ML缺少先验知识的缺点，刚好公式(5)后验概率集中了样本信息和先验信息，所以MAP估计问题可以用下面公式表示：

MAP不仅希望似然函数最大，还希望自己出现的先验概率也最大，加入先验概率，起到正则化的作用，如果θ服从高斯分布，相当于加一个L2范数正则化，如果θ服从拉普拉斯分布，相当于加一个L1范数正则化。我们继续前面抛硬币的例子，大部分人认为应该等于0.5，那么还有少数人认为p取其他值，我们认为p的取值服从Beta分布。

我们取α=5，β=5,即p以最大的概率取0.5，得到

。

3.贝叶斯估计

前面介绍的ML和MAP属于点估计，贝叶斯估计不再把参数θ看成一个未知的确定值，而是看成未知的随机变量，利用贝叶斯定理结合新的样本信息和参数θ的先验分布，来得到θ的新的概率分布(后验分布)。贝叶斯估计的本质是通过贝叶斯决策得到参数θ的最优估计

，使得贝叶斯期望损失最小。贝叶斯期望损失为：

是损失函数，我们希望

最小。如果

，则：

所以贝叶斯估计值为在样本X条件下θ的期望值，贝叶斯估计的步骤为：

确定参数θ的先验分布P(θ)

利用贝叶斯公式，求θ的后验分布：

求出贝叶斯的估计值：

我们继续前面的抛硬币的例子，后验概率：

其中

，所以可以得：

▌1.2Gamma函数

通过分部积分的方法，可以得到一个递归性质。

函数可以当成是阶乘在实数集上的延拓，

。

▌1.3二项分布

在概率论中，试验E只有两个可能结果：A及

，则称E为伯努利(Bernoulli)试验。设p(A)=p，则

。将E独立重复地进行n次，则称这一串重复的独立试验为n重伯努利试验，这里重复是指在每次试验中p(A)=p保持不变，独立是指各次试验的结果互不影响。以X表示n重伯努利试验中事件A发生的次数，称随机变量X服从参数为n,p的二项分布，记为X~B(n,p)。

▌1.4Beta分布

Beta分布是指一组定义在(0,1)区间的连续概率分布，其概率密度函数是：

1)

证明：

令t=x+y，当y=0,t=x;y=∞,t=∞，可得：

令x=μt，可得：

2)期望

证明：

▌1.5多项式分布

多项式分布是二项式分布的推广，二项式分布做n次伯努利试验，规定每次试验的结果只有两个，而多项式分布在N次独立试验中结果有K种，且每种结果都有一个确定的概率p，仍骰子是典型的多项式分布。

其中

▌1.6Dirichlet分布

Dirichlet分布是Beta分布在高维度上的推广，概率密度函数是：

▌1.7共轭先验分布

在贝叶斯中，如果后验分布与先验分布属于同类分布，则先验分布与后验分布被称为共轭分布，而先验分布被称为似然函数的共轭先验。

1．Beta-Binomial共轭

1)先验分布

2)二项式似然函数

3)后验分布

即可以表达为

取一个特殊情况理解

Beta(p|1,1)恰好是均匀分布uniform(0,1)，假设有一个不均匀的硬币抛出正面的概率为p，抛出n次后出现正面和反面的次数分别是n1和n2，开始我们对硬币不均匀性一无所知，所以应该假设p~uniform(0,1)，当有了试验样本，我们加入样本信息对p的分布进行修正,p的分布由均匀分布变为Beta分布。

2．Dirichlet-Multinomial共轭

1)先验分布

2)多项分布似然函数

3)后验分布

即可以表达为

▌1.8马氏链及其平稳分布

马氏链的数学定义很简单，状态转移的概率只依赖于前一个状态。

看一个马氏链的具体例子，马氏链表示股市模型，共有三种状态：牛市(Bullmarket)、熊市(Bearmarket)、横盘(Stagnantmarket)，每一个转态都以一定的概率转化到下一个状态，如图1.1所示。

图1.1

这个概率转化图可以以矩阵的形式表示，如果我们定义矩阵P某一位置(i,j)的值为P(j|i)，表示从状态i转化到状态j的概率，这样我们可以得到马尔科夫链模型的状态转移矩阵为：

假设初始概率分布为

从第60轮开始

的值保持不变，为[0.6250.31250.0625]。我们更改初始概率，

，从55轮开始

的值保持不变，为[0.6250.31250.0625]。两次给定不同的初始概率分布，最终都收敛到概率分布π=[0.6250.31250.0625]，也就是说收敛的行为和初始概率分布π0无关，这个收敛的行为主要是由概率转移矩阵P决定的，可以计算下

。

当n足够大的时候，

矩阵的每一行都是稳定地收敛到π=[0.6250.31250.0625]这个概率分布。这个收敛现象并不是这个马氏链独有的，而是绝大多数马氏链独有的。关于马氏链的收敛有如下定理：

定理1.1如果一个非周期马氏链具有转移概率矩阵P，且它的任何两个状态是连通的，那么

存在且与i无关，我们有：

关于上述定理，给出几点解释：

1)马氏链的状态数可以是有限的，也可以是无限的，因此可以用于连续概率分布和离散概率分布。

2)非周期马氏链：马氏链的状态转化不是循环的，如果是循环的则永远不会收敛，我们遇到的一般都是非周期马氏链。对于任意某一状态i,d为集合

的最大公约数，如果d=1，则该状态为非周期。

3)任何两个状态是连通的：从任意一个状态可以通过有限步到达其他的任意状态，不会出现条件概率一直为0导致不可达的情况。

4)π称为马氏链的平稳分布。

如果从一个具体的初始状态x0开始，沿着马氏链按照概率转移矩阵做跳转，那么可以得到一个转移序列

，由于马氏链的收敛行为，

都将是平稳分布π(x)的样本。

▌1.9MCMC

1.接受-拒绝采样

对于不常见的概率分布π(x)样本，使用接受-拒绝采样对可采样的分布q(x)进行采样得到，如图1.2所示，采样得到Mq(x)的一个样本x0，从均匀分布(0,Mq(x0))中采样得到一个值μ0，如果μ0落在图中灰色区域则拒绝这次采样，否则接受样本x0，重复上面过程得到n个接受的样本，则这些样本服从π(x)分布，具体过程见算法1.1。

图1.2

下面我们来证明下接受-拒绝方法采样得到的样本服从π(x)分布。

证明：acceptx服从π(x)分布，即p(x|accept)=π(x)。

2.MCMC

给定概率分布p(x)，希望能够生成它对应的样本，由于马氏链能收敛到平稳分布，有一个很好的想法：如果我们能构造一个转移矩阵为P的马氏链，使得该马氏链的平稳分布恰好是p(x)，那么我们从任何一个初始状态出发沿着马氏链转移，得到一个转移序列

，如果马氏链在第n步已经收敛了，于是我们可以得到p(x)的样本

，所以关键问题是如何构造转移矩阵，我们是基于下面的定理。

定理1.2(细致平稳条件)如果非周期马氏链的转移矩阵P和分布π(x)满足：

则π(x)是马氏链的平稳分布。

证明很简单，有公式(34)得：

πP=π，满足马氏链的收敛性质。这样我们就有了新的思路寻找转移矩阵P，即只要我们找到矩阵P使得概率分布π(x)满足细致平稳条件即可。

假设有一个转移矩阵为Q的马氏链(Q(i,j)表示从状态i转移到状态j的概率)，通常情况下很难满足细致平稳条件的，即：

我们对公式(36)进行改造，使细致平稳条件成立，引入α(i,j)。

α(i,j)如何取值才能使公式(37)成立？最简单的我们可以取：

Q'(i,j)=Q(i,j)α(i,j),Q'(j,i)=Q(j,i)α(j,i),所以我们有:

转移矩阵Q'满足细致平稳条件，因此马氏链Q'的平稳分布就是π(x)！

我们可以得到一个非常好的结论，转移矩阵Q'可以通过任意一个马氏链转移矩阵Q乘以α(i,j)得到，α(i,j)一般称为接受率，其取值范围为[0,1]，可以理解为一个概率值，在原来的马氏链上，从状态i以Q(i,j)的概率跳转到状态j的时候，我们以一定的概率α(i,j)接受这个转移，很像前面介绍的接受-拒绝采样，那里以一个常见的分布通过一定的接受-拒绝概率得到一个不常见的分布，这里以一个常见的马氏链状态转移矩阵Q通过一定的接受-拒绝概率得到新的马氏链状态转移矩阵Q'。

图1.3

总结下MCMC的采样过程。

MCMC采样算法有一个问题，如果接受率α(xt,x')比较小，马氏链容易原地踏步，拒绝大量的跳转，收敛到平稳分布π(x)的速度很慢，有没有办法可以使α(xt,x')变大？

3.M-H采样

M-H采样可以解决MCMC采样接受概率过低问题，回到公式(37)，若α(i,j)=0.1,α(j,i)=0.2，即：

公式(40)两边同时扩大5倍，仍然满足细致平稳条件，即：

所以我们可以把公式(37)中的α(i,j)和α(j,i)同比例放大，使得其中最大的放大到1，这样提高了采样中的接受率，细致平稳条件也没有打破，所以可以取：

提出一个问题：按照MCMC中介绍的方法把Q→Q'，是否可以保证Q'每行加和为1？

▌1.10GibbsSampling

对于高维的情形，由于接受率α≤1，M-H算法效率不够高，我们能否找到一个转移矩阵Q使得接受率α=1呢？从二维分布开始，假设p(x,y)是一个二维联合概率分布，考察某个特征维度相同的两个点A(x1,y1)和B(x1,y2)，容易发现下面等式成立：

所以可得：

也就是：

观察细致平稳条件公式，我们发现在x=x1这条直线上，如果用条件分布p(y|x1)作为任何两点之间的转移概率，那么任何两点之间的转移都满足细致平稳条件。同样的，在y=y1这条直线上任取两点A(x1,y1)和C(x2,y1)，我们可以得到：

图1.4

基于上面的发现，我们可以构造分布p(x,y)的马氏链的状态转移矩阵Q。

有了上面的转移矩阵Q，很容易验证对于平面任意两点X,Y，都满足细致平稳条件。

所以这个二维空间上的马氏链将收敛到平稳分布p(x,y)，称为GibbsSampling算法。

整个采样过程中，我们通过轮换坐标轴，得到样本(x0,y0),(x0,y1),(x1,y1),...，马氏链收敛后，最终得到的样本就是p(x,y)的样本。当然坐标轴轮换不是必须的，我们也可以每次随机选择一个坐标轴进行采样，在t时刻，可以在x轴和y轴之间随机的选择一个坐标轴，然后按照条件概率做转移。

图1.5

二维可以很容易推广到高维的情况，在n维空间中对于概率分布p(x1,x2,...xn)。

▌1.11EM算法

我们先介绍凸函数的概念，f的定义域是实数集，若x∈R且f''(x)≥0，则f是凸函数，若f''(x)>0，则f是严格凸函数；若x是向量且hessian矩阵H是半正定矩阵，则f是凸函数，若H是正定矩阵，则f是严格凸函数。

定理1.3(Jensen不等式)f的定义域是实数集，且是凸函数，则有：

如果f是严格凸函数，只有当X是常量，公式(49)等式成立即E[f(X)]=f(E[X])。

图1.6

假设训练集

，每个样本相互独立，我们需要估计模型p(x,z)的参数θ，由于含有隐变量z，所以很难直接用最大似然求解，如果z已知，那么就可以用最大似然求解。

其实我们的目标是找到z和θ使l(θ)最大，也就是分别对Z和θ求偏导，然后令其为0，理想是美好的，现实是残酷的，公式(49)求偏导后变的很复杂，求导前要是能把求和符号从对数函数中提出来就好了。EM算法可以有效地解决这个问题，引入

表示

的概率分布

。由公式(50)可得：

最后一步是利用Jensen不等式

所以f是凹函数，

是

的期望，所以有：

由公式(51)可知，我们可以不断地最大化下界，以提高l(θ)，最终达到最大值。如果固定θ，那么l(θ)的下界就取决于

，可以通过调整这个概率，使得下界不断地上升逼近l(θ)，最终相等，然后固定

，调整θ，使下界达到最大值，此时θ为新的值，再固定θ，调整

，反复直到收敛到l(θ)的最大值。现在我们有两个问题需要证明，1.下界何时等于l(θ)；2.为什么可以收敛到最大值。

第一个问题，由Jensen不等式定理中等式成立条件可知，X为常量，即：

再由

得：

下面我们先给出EM算法，然后再讨论第二个问题，E步：固定θ，根据公式(53)选择Qi使得下界等于l(θ)，M步：最大化下界，得到新的θ值。EM算法如下：

在我们开始讨论第二个问题，

是EM迭代过程的参数估计，我们需要证明

，也就是EM算法是单调地提高

。

第一个不等式是因为：

公式(57)中，

。

第二个不等式是因为

是为了

三、LDA

▌2.1UnigramModel

假设我们的词典中一共有V个词，UnigramModel就是认为上帝按照下面游戏规则产生文本的。

Game2.1UnigramModel

骰子各个面的概率记为

，对于一篇文档

，生成该文档的概率为：

假设我们预料是由m篇文档组成即

，每篇文档是相互独立的，则该预料的概率为：

假设预料中总共有N个词，每个词wi的词频为ni，那么

服从多项式分布，可参考1.5节的多项式分布概念。

此时公式(60)为：

我们需要估计模型中的参数

，可以用最大似然估计：

于是参数pk的估计值就是：

▌2.2贝叶斯UnigramModel

对于以上模型，统计学家中贝叶斯学派就不同意了，为什么上帝只有一个固定的筛子呢，在贝叶斯学派看来，一切参数都是随机变量，模型中

不是唯一固定的，而是服从一个分布，所以贝叶斯UnigramModel游戏规则变为：

Game2.2贝叶斯UnigramModel

上帝这个坛子里面有些骰子数量多，有些骰子数量少，所以从概率分布的角度看，坛子里面的骰子

服从一个概率分布

，这个分布称为参数

的先验分布。先验分布

是服从多项式分布的，

，回顾1.7节可知，

于是，在给定了参数

时候，语料中各个词出现的次数服从多项式分布

，所以后验分布为：

对参数

分布。可以用的期望值作为参数

接下来我们计算语料产生的概率，开始并不知道上帝到底用哪个骰子，所以每个骰子都有可能被使用，使用的概率由

决定的，对于每个具体的骰子，由该骰子产生预料的概率为

，所以语料产生的概率为：

▌2.3PLSA

1.PLSAModel

概率隐语义分析，是主题模型的一种。上面介绍的UnigramModel相对简单，没有考虑文档有多个主题的情况，一般一篇文档可以由多个主题(Topic)组成，文档中的每个词都是由一个固定的Topic生成的，所以PLSA的游戏规则为：

2.EM算法推导PLSA

PLSA模型中doc-topic和topic-word的每个面的概率值是固定的，所以属于点估计，但是PLSA模型既含有观测变量di,wj，又含有隐变量zk，就不能简单地直接使用极大似然估计法估计模型参数，我们可以采用EM算法估计参数。我们先介绍推导过程用到的符号含义：

：表示语料中N篇文档；

：表示语料中M个词组；

：表示词wj在文档di中出现的频次，

；

：表示K个主题，每篇文档可以有多个主题；

wj在给定文档di中出现的概率；

：表示主题zk在给定文档di下出现的概率；

：表示词wj在给定主题zk下出现的概率。

一般给定语料di,wj是可以观测的，zk是隐变量，不可以直观地观测到。我们定义“doc-word”的生成模型，如图1.8所示。

图2.3

下面进入正题，用EM算法进行模型参数估计，似然函数为：

对于给定训练预料，希望公式(69)最大化。

引入

表示zk的概率分布

，根据Jensen不等式得：

当

时，

公式(71)不等式中等号成立，所以只需要最大化：

根据拉格朗日乘子法

所以可得：

总结EM算法为：

1.E-step随机初始化变量，

计算隐变量后验概率。

2.M-step最大化似然函数，更新变量

，

3.重复1、2两步，直到收敛。

▌2.4LDA

对于PLSA模型，贝叶斯学派表示不同意，为什么上帝只有一个doc-topic骰子，为什么上帝只有固定K个topic-word骰子？

是模型的参数，一切参数都是随机变量，模型中

不是唯一固定的，类似2.2节贝叶斯UnigramModel和2.1节UnigramModel的关系。所以LDA游戏规则为：

假设我们训练语料有M篇doc，词典中有V个word，K个topic。对于第m篇文档有Nm个词。

，第m篇文档的主题分布概率，

;

，主题为k的词的概率分布，

；

：第m篇文档中属于topick的词的个数，

；

：topick产生词t的个数，

；

：

先验分布超参数；

：

先验分布超参数；

：第m篇文档中第n个词的主题；

：第m篇文档中第n个词。

LDA的概率图模型表示如图2.4所示。

图2.4

1.联合概率分布

1)

：第一步对

分布进行采样得到样本

(也就是从第一个坛子中抽取doc-topic骰子)；第二步doc-topic骰子有K个面，每个面表示一个主题，那么在一次投掷骰子过程中，每个主题的概率为

，第m篇文档有Nm个词，所以需要投掷Nm次骰子，为该篇文档中的每个词生成一个主题，第n个词对应的主题为

，整篇文档的主题表示为

。在Nm次投掷过程中，每个主题出现的次数为

，那么

服从多项式分布

(只生成每个词的主题，并未由主题产生具体的词)。可以采用贝叶斯估计对参数

进行估计。

的先验分布为

后验分布为(推导过程可以参考1.7节)

的贝叶斯估计值为

下面我们计算第m篇文档的主题概率分布：

整个语料中的M篇文档是相互独立的，所以可以得到语料中主题的概率为：

2)

：第一步对

分布进行K采样得到样本

(从第二个坛子中独立地抽取了K个topic-word骰子

)；第二步根据之前得到的主题

，为每个

生成对应的词

，

，第k个主题有

个词，所以需要投掷

，那么

服从多项式分布

，可以采用贝叶斯估计对参数

进行估计。

的先验分布为

后验分布为(推导过程可以参考1.7节)

的贝叶斯估计值为

下面我们计算第k个主题的词概率分布：

整个语料中的K个主题是相互独立的，所以可以得到语料中词的概率为：

由公式(74)、(78)、(82)可得联合概率分布为：

2.GibbsSampling

上面我们已经推导出参数的贝叶斯估计公式，但是仍然存在一个问题，公式中的

无法根据语料直接得到，如果我们知道语料中的每个词的主题，即得到

，那么就可以推断出

，进一步就可以得出贝叶斯的参数估计。我们需要利用GibbsSampling对

进行采样来得到

。先介绍一些符号定义。

下标索引；

：表示去除下标为i的词；

：第m篇文档中第n个词为t；

：第m篇文档中第n个词的主题为k；

：除去下标为i这个词，剩下的所有词中，词t属于主题k的统计次数，

(这里假设

)；

：除去下标为i的这个词，第m篇文档中主题m产生词的个数，

(这里假设

)；

：语料的主题；

：语料的单词。

1)

的计算过程类似

，仅仅在计算的时候不考虑下标为i的这个词，我们假设

；当已知语料时，

可以从语料中统计出来，所以可以认为是常量。

2)我们是推断i=(m,n)词t的主题为k的条件概率

我们再利用另外一种方法推导条件概率：

已经推导出条件概率，可以用GibbsSampling公式进行采样了。