[墙裂推荐]自然语言处理(序章):我爱自然语言处理(II)
本文紧接上一篇自然语言处理(序章):我爱自然语言处理(I),由于文章篇幅过长导致编辑器响应速度变慢,所以将其拆分为两篇,本文即为第二部分。(本博文引用内容版权属我爱自然语言博客作者及其引用文章作者,特此再次声明)。
小编推荐:本部分是我爱自然语言博客里写的最好的几篇文章了,墙裂推荐阅读;通过阅读本部分的文章,我对Metropolis Hastings算法以及Gibbs Sampling有了更为深入的了解。BTW,数学史真的很好玩啊！
MCMC 和 Gibbs Sampling
随机模拟
随机模拟(或者统计模拟)方法有一个很酷的别名是蒙特卡罗方法(Monte Carlo Simulation)。这个方法的发展始于20世纪40年代，和原子弹制造的曼哈顿计划密切相关，当时的几个大牛，包括乌拉姆、冯.诺依曼、费米、费曼、Nicholas Metropolis，在美国洛斯阿拉莫斯国家实验室研究裂变物质的中子连锁反应的时候，开始使用统计模拟的方法,并在最早的计算机上进行编程实现。
现代的统计模拟方法最早由数学家乌拉姆提出，被Metropolis命名为蒙特卡罗方法，蒙特卡罗是著名的赌场，赌博总是和统计密切关联的，所以这个命名风趣而贴切，很快被大家广泛接受。被不过据说费米之前就已经在实验中使用了，但是没有发表。说起蒙特卡罗方法的源头，可以追溯到18世纪，布丰当年用于计算π的著名的投针实验就是蒙特卡罗模拟实验。统计采样的方法其实数学家们很早就知道，但是在计算机出现以前，随机数生成的成本很高，所以该方法也没有实用价值。随着计算机技术在二十世纪后半叶的迅猛发展，随机模拟技术很快进入实用阶段。对那些用确定算法不可行或不可能解决的问题，蒙特卡罗方法常常为人们带来希望。
统计模拟中有一个重要的问题就是给定一个概率分布$p(x)$，我们如何在计算机中生成它的样本。一般而言均匀分布 $Uniform(0,1)$的样本是相对容易生成的。通过线性同余发生器可以生成伪随机数，我们用确定性算法生成$[0,1]$之间的伪随机数序列后，这些序列的各种统计指标和均匀分布$Uniform(0,1)$的理论计算结果非常接近。这样的伪随机序列就有比较好的统计性质，可以被当成真实的随机数使用。
而我们常见的概率分布，无论是连续的还是离散的分布，都可以基于$Uniform(0,1)$的样本生成。例如正态分布可以通过著名的Box-Muller变换得到
[Box-Muller变换]
如果随机变量$U_1$,$U_2$独立且$U_1,U_2 \sim\ Uniform[0,1]$,
\begin{equation} \begin{split} Z_0 & = \sqrt{-2\ln U_1} cos(2\pi U_2) \\ Z_1 & = \sqrt{-2\ln U_1} sin(2\pi U_2) \end{split} \end{equation}
则$Z_0,Z_1$独立且服从标准正态分布。
其它几个著名的连续分布，包括指数分布、Gamma分布、t分布、F分布、Beta分布、Dirichlet分布等等,也都可以通过类似的数学变换得到；离散的分布通过均匀分布更加容易生成。更多的统计分布如何通过均匀分布的变换生成出来，大家可以参考统计计算的书，其中 Sheldon M. Ross 的《统计模拟》是写得非常通俗易懂的一本。
不过我们并不是总是这么幸运的，当$p(x)$的形式很复杂，或者$p(x)$是个高维的分布的时候，样本的生成就可能很困难了。 譬如有如下的情况:
- $p(x) = \frac{\tilde{p}(x)}{\int \tilde{p}(x) dx}$,而$\tilde{p}(x)$我们是可以计算的，但是底下的积分式无法显式计算。
- $p(x,y)$是一个二维的分布函数，这个函数本身计算很困难，但是条件分布$p(x|y),p(y|x)$的计算相对简单;如果$p(x)$是高维的，这种情形就更加明显。
此时就需要使用一些更加复杂的随机模拟的方法来生成样本。而本节中将要重点介绍的 MCMC(Markov Chain Monte Carlo) 和 Gibbs Sampling算法就是最常用的一种，这两个方法在现代贝叶斯分析中被广泛使用。要了解这两个算法，我们首先要对马氏链的平稳分布的性质有基本的认识。
马氏链及其平稳分布
马氏链的数学定义很简单:
\begin{equation} P(X_{t+1}=x|X_t, X_{t-1}, \cdots) =P(X_{t+1}=x|X_t) \end{equation}
也就是状态转移的概率只依赖于前一个状态。
我们先来看马氏链的一个具体的例子。社会学家经常把人按其经济状况分成3类：下层(lower-class)、中层(middle-class)、上层(upper-class)，我们用1,2,3分别代表这三个阶层。社会学家们发现决定一个人的收入阶层的最重要的因素就是其父母的收入阶层。如果一个人的收入属于下层类别，那么他的孩子属于下层收入的概率是 0.65, 属于中层收入的概率是 0.28, 属于上层收入的概率是 0.07。事实上，从父代到子代，收入阶层的变化的转移概率如下:
使用矩阵的表示方式，转移概率矩阵记为:
\begin{equation} P=\left[ \begin{array}{cc} 0.65 & 0.28 & 0.07 \\ 0.15 & 0.67 & 0.18 \\ 0.12 & 0.36 & 0.52 \\ \end{array} \right] \end{equation}
假设当前这一代人处在下层、中层、上层的人的比例是概率分布向量 $\pi_0=[\pi_0(1),\pi_0(2),\pi_0(3)]$，那么他们的子女的分布比例将是$\pi_1=\pi_0P$, 他们的孙子代的分布比例将是 $\pi_2=\pi_1P=\pi_0P^2$, ……, 第n代子孙的收入分布比例将是$\pi_n=\pi_{n−1}P=\pi_0P^n$。
假设初始概率分布为$\pi_0=[0.21,0.68,0.11]$，则我们可以计算前$n$代人的分布状况如下
我们发现从第7代人开始，这个分布就稳定不变了，这个是偶然的吗？我们换一个初始概率分布$\pi_0=[0.75,0.15,0.1]$.试试看，继续计算前$n$代人的分布状况如下
我们发现，到第9代人的时候, 分布又收敛了。最为奇特的是，两次给定不同的初始概率分布，最终都收敛到概率分布 $\pi=[0.286,0.489,0.225]$，也就是说收敛的行为和初始概率分布$\pi_0$无关。这说明这个收敛行为主要是由概率转移矩阵$P$决定的。我们计算一下$P^n$.
\begin{equation} P^{20} = P^{21} = \cdots = P^{100} = \cdots = \begin{bmatrix} 0.286 & 0.489 & 0.225 \\ 0.286 & 0.489 & 0.225 \\ 0.286 & 0.489 & 0.225 \\ \end{bmatrix} \end{equation}
我们发现，当$n$足够大的时候，这个$P^n$矩阵的每一行都是稳定地收敛到$\pi=[0.286,0.489,0.225]$这个概率分布。自然的，这个收敛现象并非是我们这个马氏链独有的，而是绝大多数马氏链的共同行为，关于马氏链的收敛我们有如下漂亮的定理：
马氏链定理
如果一个非周期马氏链具有转移概率矩阵$P$,且它的任何两个状态是连通的，那么$\lim_{n\rightarrow\infty}P_{ij}^n$存在且与$i$无关，记$\lim_{n\rightarrow\infty}P_{ij}^n = \pi(j)$, 我们有:
\begin{equation} \begin{split} \lim_{n \rightarrow \infty} P^n =\begin{bmatrix} \pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \\ \pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ \pi(1) & \pi(2) & \cdots & \pi(j) & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ \end{bmatrix} \\ \pi(j) = \sum_{i=0}^{\infty}\pi(i)P_{ij} \\ \end{split} \end{equation}
$\pi$是方程$\pi P=\pi$的唯一非负解。其中,
\begin{equation} \pi = [\pi(1), \pi(2), \cdots, \pi(j),\cdots ], \quad \sum_{i=0}^{\infty} \pi_i = 1 \end{equation}
$\pi$称为马氏链的平稳分布。
这个马氏链的收敛定理非常重要，所有的 MCMC(Markov Chain Monte Carlo) 方法都是以这个定理作为理论基础的。 定理的证明相对复杂，一般的随机过程课本中也不给证明，所以我们就不用纠结它的证明了，直接用这个定理的结论就好了。我们对这个定理的内容做一些解释说明：
- 该定理中马氏链的状态不要求有限，可以是有无穷多个的；
- 定理中的“非周期“这个概念我们不打算解释了，因为我们遇到的绝大多数马氏链都是非周期的；
- 两个状态$i,j$是连通并非指$i$可以直接一步转移到$j$($P_{ij}>0)$,而是指$i$可以通过有限的$n$步转移到达$j$($P^n_{ij}>0$)。马氏链的任何两个状态是连通的含义是指存在一个$n$,使得矩阵$P^n$中的任何一个元素的数值都大于零。
- 我们用$X_i$表示在马氏链上跳转第$i$步后所处的状态，如果$\lim_{n\rightarrow\infty}P_{ij}^n=\pi(j)$存在，很容易证明以上定理的第二个结论。由于
\begin{equation} \begin{split} P(X_{n+1}=j) & = \sum_{i=0}^\infty P(X_n=i) P(X_{n+1}=j|X_n=i) \\ & = \sum_{i=0}^\infty P(X_n=i) P_{ij} \end{split} \end{equation}
上式两边取极限就得到$\pi(j) = \sum_{i=0}^{\infty}\pi(i)P_{ij}$.
从初始概率分布$\pi_0$出发，我们在马氏链上做状态转移，记$X_i$的概率分布为$\pi_i$, 则有:
\begin{equation} \begin{split} X_0 & \sim \pi_0(x) \ X_i & \sim \pi_i(x), \quad\quad \pi_i(x) = \pi_{i-1}(x)P = \pi_0(x)P^n \end{split} \end{equation}
由马氏链收敛的定理, 概率分布$\pi_i(x)$将收敛到平稳分布$\pi_(x)$。假设到第$n$步的时候马氏链收敛，则有
\begin{equation} \begin{split} X_0 & \sim \pi_0(x) \\ X_1 & \sim \pi_1(x) \\ & \cdots \\ X_n & \sim \pi_n(x)=\pi(x) \\ X_{n+1} & \sim \pi(x) \\ X_{n+2}& \sim \pi(x) \\ & \cdots \end{split} \end{equation}
所以$X_n,X_{n+1},X_{n+2},\cdots \sim \pi(x)$都是同分布的随机变量，当然他们并不独立。如果我们从一个具体的初始状态$x_0$开始,沿着马氏链按照概率转移矩阵做跳转，那么我们得到一个转移序列$x_0,x_1,x_2,\cdots,x_n,x_{n+1},\cdots$, 由于马氏链的收敛行为，$x_n,x_{n+1},...$ 都将是平稳分布$\pi(x)$的样本。
Markov Chain Monte Carlo
对于给定的概率分布$p(x)$,我们希望能有便捷的方式生成它对应的样本。由于马氏链能收敛到平稳分布，于是一个很漂亮的想法是：如果我们能构造一个转移矩阵为$P$的马氏链，使得该马氏链的平稳分布恰好是$p(x)$,那么我们从任何一个初始状态$x_0$出发沿着马氏链转移, 得到一个转移序列$x_0,x_1,x_2,\cdots,x_n,x_{n+1}\cdots$，如果马氏链在第$n$步已经收敛了，于是我们就得到了$p(x)$的样本$x_n,x_{n+1},\cdots$。
这个绝妙的想法在1953年被Metropolis想到了，为了研究粒子系统的平稳性质，Metropolis考虑了物理学中常见的波尔兹曼分布的采样问题，首次提出了基于马氏链的蒙特卡罗方法，即Metropolis算法，并在最早的计算机上编程实现。Metropolis算法是首个普适的采样方法，并启发了一系列MCMC方法，所以人们把它视为随机模拟技术腾飞的起点。Metropolis的这篇论文被收录在《统计学中的重大突破》中，Metropolis算法也被遴选为二十世纪的十个最重要的算法之一。
我们接下来介绍的MCMC算法是Metropolis算法的一个改进变种，即常用的Metropolis-Hastings算法。由上一节的例子和定理我们看到了，马氏链的收敛性质主要由转移矩阵$P$决定,所以基于马氏链做采样的关键问题是如何构造转移矩阵$P$,使得平稳分布恰好是我们要的分布$p(x)$。如何能做到这一点呢？我们主要使用如下的定理。
定理：[细致平稳条件]如果非周期马氏链的转移矩阵$P$和分布$\pi(x)$满足:对于任意$i$和$j$,有:
\begin{equation} \pi(i)P_{ij} = \pi(j)P_{ji} \end{equation}
则$\pi(x)$是马氏链的平稳分布，上式被称为细致平稳条件(Detailed balance condition)。
其实这个定理是显而易见的，因为细致平稳条件的物理含义就是对于任何两个状态$i$,$j$,从$i$转移出去到$j$而丢失的概率质量，恰好会被从$j$转移回$i$的概率质量补充回来，所以状态$i$上的概率质量$\pi(i)$是稳定的，从而$\pi(x)$是马氏链的平稳分布。数学上的证明也很简单，由细致平稳条件可得:
\begin{equation} \begin{split} & \sum_{i=1}^\infty \pi(i)P_{ij} = \sum_{i=1}^\infty \pi(j)P_{ji} = \sum_{i=1}^\infty \pi(j)P_{ji} = \pi(j) \\ & \rightarrow \pi P = \pi \end{split} \end{equation}
由于$\pi$是方程$\pi P =\pi$的解，所以$\pi$是平稳分布。
假设我们已经有一个转移矩阵为$Q$的马氏链($q(i,j)$表示从状态$i$转移到状态$j$的概率，也可以写为$q(j|i)$或者$q(i \to j)$), 显然，通常情况下$p(i) q(i,j)\neq p(j)q(j,i)$也就是细致平稳条件不成立，所以$p(x)$不太可能是这个马氏链的平稳分布。我们可否对马氏链做一个改造，使得细致平稳条件成立呢？譬如，我们引入一个$\alpha(i,j)$, 我们希望:
\begin{equation} p(i)q(i,j)\alpha(i,j)=p(j)q(j,i)\alpha(j,i) \quad (∗) \end{equation}
取什么样的$\alpha(i,j)$以上等式能成立呢？最简单的，按照对称性，我们可以取:
\begin{equation} \alpha(i,j)=p(j)q(j,i)，\alpha(j,i)=p(i)q(i,j) \end{equation}
于是(*)式就成立了。所以有:
于是我们把原来具有转移矩阵$Q$的一个很普通的马氏链，改造为了具有转移矩阵$Q\prime$的马氏链，而 $Q\prime$恰好满足细致平稳条件，由此马氏链$Q\prime$的平稳分布就是$p(x)$了!
在改造$Q$的过程中引入的$\alpha(i,j)$称为接受率，物理意义可以理解为在原来的马氏链上，从状态$i$以$q(i,j)$的概率转跳转到状态$j$的时候，我们以$\alpha(i,j)$的概率接受这个转移，于是得到新的马氏链$Q\prime$的转移概率为$q(i,j)\alpha(i,j)$。
假设我们已经有一个转移矩阵$Q$(对应元素为$q(i,j)$), 把以上的过程整理一下，我们就得到了如下的用于采样概率分布$p(x)$的算法。
上述过程中$p(x),q(x|y)$说的都是离散的情形，事实上即便这两个分布是连续的，以上算法仍然是有效，于是就得到更一般的连续概率分布$p(x)$的采样算法，而$q(x|y)$就是任意一个连续二元概率分布对应的条件分布。
以上的MCMC采样算法已经能很漂亮的工作了，不过它有一个小的问题：马氏链$Q$在转移的过程中的接受率 $\alpha(i,j)$可能偏小，这样采样过程中马氏链容易原地踏步，拒绝大量的跳转，这使得马氏链遍历所有的状态空间要花费太长的时间，收敛到平稳分布$p(x)的$速度太慢。有没有办法提升一些接受率呢?
假设$\alpha(i,j)=0.1,\alpha(j,i)=0.2$, 此时满足细致平稳条件，于是
\begin{equation} p(i)q(i,j)×0.1=p(j)q(j,i)×0.2 \end{equation}
上式两边扩大5倍，我们改写为
\begin{equation} p(i)q(i,j)×0.5=p(j)q(j,i)×1 \end{equation}
看，我们提高了接受率，而细致平稳条件并没有打破！这启发我们可以把细致平稳条件(**) 式中的$\alpha(i,j),\alpha(j,i)$同比例放大，使得两数中最大的一个放大到1，这样我们就提高了采样中的跳转接受率。所以我们可以取:
于是，经过对上述MCMC采样算法中接受率的微小改造，我们就得到了如下教科书中最常见的Metropolis-Hastings算法。
对于分布$p(x)$,我们构造转移矩阵$Q\prime$使其满足细致平稳条件:
\begin{equation} p(x)Q\prime (x→y)=p(y)Q\prime(y→x) \end{equation}
此处$x$并不要求是一维的，对于高维空间的$p(\mathbf{x})$，如果满足细致平稳条件:
\begin{equation} p(\mathbf{x}) Q’(\mathbf{x}\rightarrow \mathbf{y}) = p(\mathbf{y}) Q’(\mathbf{y}\rightarrow \mathbf{x}) \end{equation}
那么以上的Metropolis-Hastings算法一样有效。
Gibbs Sampling
对于高维的情形，由于接受率$\alpha$的存在(通常$\alpha$<1), 以上Metropolis-Hastings算法的效率不够高。能否找到一个转移矩阵$Q$使得接受率$\alpha=1$呢？我们先看看二维的情形，假设有一个概率分布 $p(x,y)$, 考察$x$坐标相同的两个点$A(x_1,y_1),B(x_1,y_2)$,我们发现:
\begin{equation} \begin{split} p(x_1,y_1)p(y_2|x_1) &= p(x_1)p(y_1|x_1)p(y_2|x_1) \\ p(x_1,y_2)p(y_1|x_1) &= p(x_1)p(y_2|x_1)p(y_1|x_1) \end{split} \end{equation}
所以得到:
\begin{equation} p(x_1,y_1)p(y_2|x_1) = p(x_1,y_2)p(y_1|x_1) \quad (***) \end{equation}
即:
\begin{equation} p(A)p(y_2|x_1) = p(B)p(y_1|x_1) \end{equation}
基于以上等式，我们发现，在$x=x_1$这条平行于$y$轴的直线上，如果使用条件分布$p(y|x_1)$做为任何两个点之间的转移概率，那么任何两个点之间的转移满足细致平稳条件。同样的，如果我们在$y=y_1$这条直线上任意取两个点$A(x_1,y_1),C(x_2,y_1)$,也有如下等式:
\begin{equation} p(A)p(x_2|y_1)=p(C)p(x_1|y_1) \end{equation}
于是我们可以如下构造平面上任意两点之间的转移概率矩阵$Q$:
\begin{equation} \begin{split} Q(A\rightarrow B) & = p(y_B|x_1) & \text{如果} \quad x_A=x_B=x_1 & \\ Q(A\rightarrow C) & = p(x_C|y_1) & \text{如果} \quad y_A=y_C=y_1 & \\ Q(A\rightarrow D) & = 0 & \text{其它} & \end{split} \end{equation}
有了如上的转移矩阵$Q$, 我们很容易验证对平面上任意两点$X,Y$, 满足细致平稳条件:
\begin{equation} p(X)Q(X→Y)=p(Y)Q(Y→X) \end{equation}
于是这个二维空间上的马氏链将收敛到平稳分布$p(x,y)$。而这个算法就称为Gibbs Sampling算法,是 Stuart Geman 和Donald Geman 这两兄弟于1984年提出来的，之所以叫做Gibbs Sampling 是因为他们研究了Gibbs random field, 这个算法在现代贝叶斯分析中占据重要位置。
以上采样过程中，如上图所示，马氏链的转移只是轮换的沿着坐标轴$x$轴和$y$轴做转移，于是得到样本 $(x_0,y_0),(x_0,y_1),(x_1,y_1),(x_1,y_2),(x_2,y_2),\cdots$,马氏链收敛后，最终得到的样本就是$p(x,y)$的样本，而收敛之前的阶段称为 burn-in period。额外说明一下，我们看到教科书上的 Gibbs Sampling 算法大都是坐标轴轮换采样的，但是这其实是不强制要求的。最一般的情形可以是，在$t$时刻，可以在$x$轴和$y$轴之间随机的选一个坐标轴，然后按条件概率做转移，马氏链也是一样收敛的。轮换两个坐标轴只是一种方便的形式。
以上的过程我们很容易推广到高维的情形，对于(***)式，如果$x_1$变为多维情形$x_1$,可以看出推导过程不变，所以细致平稳条件同样是成立的.
\begin{equation} p(\mathbf{x_1},y_1)p(y_2|\mathbf{x_1}) = p(\mathbf{x_1},y_2)p(y_1|\mathbf{x_1}) \end{equation}
此时转移矩阵$Q$由条件分布$p(y|x_1)$定义。上式只是说明了一根坐标轴的情形，和二维情形类似，很容易验证对所有坐标轴都有类似的结论。所以$n$维空间中对于概率分布$p(x_1,x_2,\cdots,x_n)$可以如下定义转移矩阵:
- 如果当前状态为$(x_1,x_2,⋯,x_n)$，马氏链转移的过程中，只能沿着坐标轴做转移。沿着$x_i$这根坐标轴做转移的时候，转移概率由条件概率$p(x_i|x_1,⋯,x_{i−1},x_{i+1},⋯,x_n)$定义;
- 其它无法沿着单根坐标轴进行的跳转，转移概率都设置为 0。
于是我们可以把Gibbs Smapling算法从采样二维的$p(\mathbf{x},\mathbf{y})$推广到采样$n$维的 $p(x_1,x_2,⋯,x_n)$.
以上算法收敛后，得到的就是概率分布$p(x_1,x_2,⋯,x_n)$的样本，当然这些样本并不独立，但是我们此处要求的是采样得到的样本符合给定的概率分布，并不要求独立。同样的，在以上算法中，坐标轴轮换采样不是必须的，可以在坐标轴轮换中引入随机性，这时候转移矩阵$Q$中任何两个点的转移概率中就会包含坐标轴选择的概率，而在通常的 Gibbs Sampling 算法中，坐标轴轮换是一个确定性的过程，也就是在给定刻$t$，在一根固定的坐标轴上转移的概率是1。
文本建模
文本建模
我们日常生活中总是产生大量的文本，如果每一个文本存储为一篇文档，那每篇文档从人的观察来说就是有序的词的序列$d=(w_1,w_2,\cdots,w_n)$。
包含$M$篇文档的语料库统计文本建模的目的就是追问这些观察到语料库中的的词序列是如何生成的。统计学被人们描述为猜测上帝的游戏，人类产生的所有的语料文本我们都可以看成是一个伟大的上帝在天堂中抛掷骰子生成的，我们观察到的只是上帝玩这个游戏的结果 —— 词序列构成的语料，而上帝玩这个游戏的过程对我们是个黑盒子。所以在统计文本建模中，我们希望猜测出上帝是如何玩这个游戏的，具体一点，最核心的两个问题是:
- 上帝都有什么样的骰子；
- 上帝是如何抛掷这些骰子的；
第一个问题就是表示模型中都有哪些参数，骰子的每一个面的概率都对应于模型中的参数；第二个问题就表示游戏规则是什么，上帝可能有各种不同类型的骰子，上帝可以按照一定的规则抛掷这些骰子从而产生词序列。
Unigram Model
假设我们的词典中一共有$V$个词$v_1,v_2,⋯v_V$，那么最简单的Unigram Model就是认为上帝是按照如下的游戏规则产生文本的。
上帝的这个唯一的骰子各个面的概率记为$\vec{p} = (p_1, p_2, \cdots, p_V)$, 所以每次投掷骰子类似于一个抛钢镚时候的贝努利实验， 记为$w\sim Mult(w|\vec{p})$.
上帝投掷$V$个面的骰子对于一篇文档$d=\vec w=(w_1, w_2, \cdots, w_n)$, 该文档被生成的概率就是:
\begin{equation} p(\vec w) = p(w_1, w_2, \cdots, w_n) = p(w_1)p(w_2) \cdots p(w_n) \end{equation}
而文档和文档之间我们认为是独立的， 所以如果语料中有多篇文档$\mathcal{W}=(\vec{w_1}, \vec{w_2},…,\vec{w_m})$,则该语料的概率是:
\begin{equation} p(\mathcal{W})= p(\vec{w_1})p(\vec{w_2}) \cdots p(\vec{w_m}) \end{equation}
在Unigram Model中， 我们假设了文档之间是独立可交换的，而文档中的词也是独立可交换的，所以一篇文档相当于一个袋子，里面装了一些词，而词的顺序信息就无关紧要了，这样的模型也称为词袋模型(Bag-of-words)。
假设语料中总的词数是$N$, 在所有的$N$个词中,如果我们关注每个词$v_i$的发生次数$n_i$，那么$ n=(n_1, n_2,\cdots, n_V)$正好是一个多项分布:
\begin{equation} p( n) = Mult( n|\vec{p}, N) = \binom{N}{ n} \prod_{k=1}^V p_k^{n_k} \end{equation}
此时， 语料的概率是:
\begin{equation} p(\mathcal{W})= p(\vec{w_1})p(\vec{w_2}) \cdots p(\vec{w_m}) = \prod_{k=1}^V p_k^{n_k} \end{equation}
当然，我们很重要的一个任务就是估计模型中的参数$\vec{p}$，也就是问上帝拥有的这个骰子的各个面的概率是多大，按照统计学家中频率派的观点，使用最大似然估计最大化$P(\mathcal{W})$，于是参数$p_i$的估计值就是$\hat{p_i}=\frac{n_i}{N}$.
对于以上模型，贝叶斯统计学派的统计学家会有不同意见，他们会很挑剔的批评只假设上帝拥有唯一一个固定的骰子是不合理的。在贝叶斯学派看来，一切参数都是随机变量，以上模型中的骰子$\vec{p}$不是唯一固定的，它也是一个随机变量。所以按照贝叶斯学派的观点，上帝是按照以下的过程在玩游戏的:
上帝的这个坛子里面，骰子可以是无穷多个，有些类型的骰子数量多，有些类型的骰子少，所以从概率分布的角度看，坛子里面的骰子$\vec{p}$服从一个概率分布$p(\vec{p})$，这个分布称为参数$\vec{p}$的先验分布。
以上贝叶斯学派的游戏规则的假设之下，语料$\mathcal{W}$产生的概率如何计算呢？由于我们并不知道上帝到底用了哪个骰子$\vec{p}$,所以每个骰子都是可能被使用的，只是使用的概率由先验分布$p(\vec{p})$来决定。对每一个具体的骰子$\vec{p}$,由该骰子产生数据的概率是$p(\mathcal{W}|\vec{p})$, 所以最终数据产生的概率就是对每一个骰子$\vec{p}$上产生的数据概率进行积分累加求和:
\begin{equation} p(\mathcal{W}) = \int p(\mathcal{W}|\vec{p}) p(\vec{p})d\vec{p} \end{equation}
在贝叶斯分析的框架下，此处先验分布$p(\vec{p})$就可以有很多种选择了，注意到:
\begin{equation} p( n) = Mult( n|\vec{p}, N) \end{equation}
实际上是在计算一个多项分布的概率，所以对先验分布的一个比较好的选择就是多项分布对应的共轭分布,即 Dirichlet 分布:
\begin{equation} Dir(\vec{p}|\vec\alpha)= \frac{1}{\Delta(\vec\alpha)} \prod_{k=1}^V p_k^{\alpha_k -1}， \quad \vec\alpha=(\alpha_1, \cdots, \alpha_V) \end{equation}
此处，$\Delta(\vec\alpha)$就是归一化因子$Dir(\vec\alpha)$，即:
\begin{equation} \Delta(\vec\alpha) = \int \prod_{k=1}^V p_k^{\alpha_k -1} d\vec{p} \end{equation}
回顾前一个小节介绍的Drichlet分布的一些知识，其中很重要的一点就是:Dirichlet 先验 + 多项分布的数据 → 后验分布为 Dirichlet 分布.
\begin{equation} Dir(\vec{p}|\vec\alpha) + MultCount( n)= Dir(\vec{p}|\vec\alpha+ n) \end{equation}
于是，在给定了参数$\vec{p}$的先验分布$Dir(\vec{p}|\vec\alpha)$的时候，各个词出现频次的数据$ n \sim Mult( n|\vec{p},N)$为多项分布, 所以无需计算，我们就可以推出后验分布是:
\begin{equation} p(\vec{p}|\mathcal{W},\vec\alpha) = Dir(\vec{p}| n+ \vec\alpha) = \frac{1}{\Delta( n+\vec\alpha)} \prod_{k=1}^V p_k^{n_k + \alpha_k -1} d\vec{p} \end{equation}
在贝叶斯的框架下，参数$\vec{p}$如何估计呢？由于我们已经有了参数的后验分布，所以合理的方式是使用后验分布的极大值点，或者是参数在后验分布下的平均值。在该文档中，我们取平均值作为参数的估计值。使用上个小节中的结论，由于$\vec{p}$的后验分布为$Dir(\vec{p}| n + \vec\alpha)$，于是:
\begin{equation} E(\vec{p}) = \Bigl(\frac{n_1 + \alpha_1}{\sum_{i=1}^V(n_i + \alpha_i)}, \frac{n_2 + \alpha_2}{\sum_{i=1}^V(n_i + \alpha_i)}, \cdots, \frac{n_V + \alpha_V}{\sum_{i=1}^V(n_i + \alpha_i)} \Bigr) \end{equation}
也就是说对每一个$p_i$, 我们用下式做参数估计:
\begin{equation} \hat{p_i} = \frac{n_i + \alpha_i}{\sum_{i=1}^V(n_i + \alpha_i)} \end{equation}
考虑到$\alpha_i$在Dirichlet 分布中的物理意义是事件的先验的伪计数，这个估计式子的含义是很直观的：每个参数的估计值是其对应事件的先验的伪计数和数据中的计数的和在整体计数中的比例。
进一步，我们可以计算出文本语料的产生概率为:
\begin{equation} \begin{split} p(\mathcal{W}|\vec\alpha) & = \int p(\mathcal{W}|\vec{p}) p(\vec{p}|\vec\alpha)d\vec{p} \notag \\ & = \int \prod_{k=1}^V p_k^{n_k} Dir(\vec{p}|\vec\alpha) d\vec{p} \notag \\ & = \int \prod_{k=1}^V p_k^{n_k} \frac{1}{\Delta(\vec\alpha)} \prod_{k=1}^V p_k^{\alpha_k -1} d\vec{p} \notag \\ & = \frac{1}{\Delta(\vec\alpha)} \int \prod_{k=1}^V p_k^{n_k + \alpha_k -1} d\vec{p} \notag \\ & = \frac{\Delta( n+\vec\alpha)}{\Delta(\vec\alpha)} \end{split} \end{equation}
Topic Model 和 PLSA
以上 Unigram Model 是一个很简单的模型，模型中的假设看起来过于简单，和人类写文章产生每一个词的过程差距比较大，有没有更好的模型呢？
我们可以看看日常生活中人是如何构思文章的。如果我们要写一篇文章，往往是先确定要写哪几个主题。譬如构思一篇自然语言处理相关的文章，可能$40\%$会谈论语言学,$30\%$谈论概率统计,$20\%$谈论计算机、还有$10\%$谈论其它的主题：
- 说到语言学，我们容易想到的词包括：语法、句子、乔姆斯基、句法分析、主语…；
- 谈论概率统计，我们容易想到以下一些词: 概率、模型、均值、方差、证明、独立、马尔科夫链、…；
- 谈论计算机，我们容易想到的词是： 内存、硬盘、编程、二进制、对象、算法、复杂度…；
我们之所以能马上想到这些词，是因为这些词在对应的主题下出现的概率很高。我们可以很自然的看到，一篇文章通常是由多个主题构成的、而每一个主题大概可以用与该主题相关的频率最高的一些词来描述。
以上这种直观的想法由Hoffman 于 1999 年给出的PLSA(Probabilistic Latent Semantic Analysis) 模型中首先进行了明确的数学化。Hoffman 认为一篇文档(Document) 可以由多个主题(Topic) 混合而成， 而每个Topic 都是词汇上的概率分布，文章中的每个词都是由一个固定的 Topic 生成的。下图是英语中几个Topic 的例子。
所有人类思考和写文章的行为都可以认为是上帝的行为，我们继续回到上帝的假设中，那么在 PLSA 模型中，Hoffman 认为上帝是按照如下的游戏规则来生成文本的。
以上PLSA 模型的文档生成的过程可以图形化的表示为:
PLSA模型的文档生成过程我们可以发现在以上的游戏规则下，文档和文档之间是独立可交换的，同一个文档内的词也是独立可交换的，还是一个bag-of-words模型。游戏中的$K$个topic-word骰子，我们可以记为$\vec\varphi_1, \cdots, \vec\varphi_K$, 对于包含$M$篇文档的语料$C=(d_1, d_2, \cdots, d_M)$中的每篇文档$d_m$，都会有一个特定的doc-topic骰子$\vec{\theta}_m$，所有对应的骰子记为$\vec{\theta}_1, \cdots, \vec{\theta}_M$。为了方便，我们假设每个词$w$ 都是一个编号，对应到topic-word骰子的面。于是在 PLSA 这个模型中，第$m$篇文档$d_m$中的每个词的生成概率为:
\begin{equation} p(w|d_m) = \sum_{z=1}^K p(w|z)p(z|d_m) = \sum_{z=1}^K \varphi_{zw} \theta_{mz} \end{equation}
所以整篇文档的生成概率为:
\begin{equation} p(\vec w|d_m) = \prod_{i=1}^n \sum_{z=1}^K p(w_i|z)p(z|d_m) = \prod_{i=1}^n \sum_{z=1}^K \varphi_{zw_i} \theta_{mz} \end{equation}
由于文档之间相互独立，我们也容易写出整个语料的生成概率。求解PLSA这个Topic Model的过程汇总，模型参数并容易求解，可以使用著名的EM算法进行求得局部最优解，由于该模型的求解并不是本文的介绍要点，有兴趣的同学参考Hoffman的原始论文，此处略去不讲。
LDA文本建模
游戏规则
对于上述的 PLSA 模型，贝叶斯学派显然是有意见的，doc-topic 骰子$\vec{\theta}_m$和topic-word骰子$\vec\varphi_k$都是模型中的参数，参数都是随机变量，怎么能没有先验分布呢？于是，类似于对Unigram Model的贝叶斯改造， 我们也可以如下在两个骰子参数前加上先验分布从而把PLSA对应的游戏过程改造为一个贝叶斯的游戏过程。由于$\vec\varphi_k$和$\vec{\theta}_m$都对应到多项分布，所以先验分布的一个好的选择就是Drichlet分布，于是我们就得到了LDA(Latent Dirichlet Allocation)模型。
在 LDA 模型中, 上帝是按照如下的规则玩文档生成的游戏的.
假设语料库中有$M$ 篇文档，所有的的word和对应的 topic 如下表示:
\begin{equation} \begin{split} \vec{\mathbf{w}} & = (\vec w_1, \cdots, \vec w_M) \\ \vec{\mathbf{z}} & = (\vec z_1, \cdots, \vec z_M) \end{split} \end{equation}
其中，$\vec w_m$表示第$m$篇文档中的词，$\vec z_m$表示这些词对应的topic编号。
物理过程分解
使用概率图模型表示， LDA模型的游戏过程如图所示。
这个概率图可以分解为两个主要的物理过程：
理解LDA最重要的就是理解这两个物理过程。LDA模型在基于$K$个topic生成语料中的$M$篇文档的过程中,由于是bag-of-words模型，有一些物理过程是相互独立可交换的。由此,LDA生成模型中,$M$篇文档会对应于$M$个独立的Dirichlet-Multinomial共轭结构；$K$个topic会对应于$K$个独立的Dirichlet-Multinomial 共轭结构.所以理解LDA所需要的所有数学就是理解Dirichlet-Multiomail共轭，其它就是理解物理过程。现在我们进入细节，来看看LDA模型是如何被分解为$M+K$个Dirichlet-Multinomial共轭结构的。
由第一个物理过程，我们知道$\vec\alpha \to \vec\theta_m \to \vec z_m$表示生成第$m$篇文档中的所有词对应的topics，显然$\vec\alpha \to \vec\theta_m$对应于Dirichlet分布，$\vec\theta_m \to \vec z_{m}$对应于Multinomial分布,所以整体是一个Dirichlet-Multinomial共轭结构；
前文介绍Bayesian Unigram Model的小节中我们对Dirichlet-Multinomial共轭结构做了一些计算。借助于该小节中的结论，我们可以得到:
\begin{equation} p(\vec z_m |\vec\alpha) = \frac{\Delta(\vec n_m+\vec\alpha)}{\Delta(\vec\alpha)} \end{equation}
其中$\vec n_m = (n_{m}^{(1)}, \cdots, n_{m}^{(K)})$表示第$m$篇文档中第$k$个topic产生的词的个数。进一步，利用Dirichlet-Multiomial共轭结构，我们得到参数$\vec{\theta}_m$的后验分布恰好是$Dir(\vec{\theta}_m| n_m + \vec\alpha)$.
由于语料中$M$篇文档的 topics生成过程相互独立，所以我们得到$M$个相互独立的Dirichlet-Multinomial共轭结构，从而我们可以得到整个语料中topics生成概率:
\begin{equation} \begin{split} p(\vec{\mathbf{z}} |\vec\alpha) & = \prod_{m=1}^M p(\vec z_m |\vec\alpha) \notag \\ &= \prod_{m=1}^M \frac{\Delta( n_m+\vec\alpha)}{\Delta(\vec\alpha)} \quad\quad (*) \end{split} \end{equation}
目前为止，我们由$M$篇文档得到了$M$个Dirichlet-Multinomial共轭结构，还有额外$K$个Dirichlet-Multinomial共轭结构在哪儿呢？在上帝按照之前的规则玩LDA游戏的时候，上帝是先完全处理完成一篇文档，再处理下一篇文档。文档中每个词的生成都要抛两次骰子，第一次抛一个doc-topic骰子得到topic, 第二次抛一个topic-word骰子得到word，每次生成每篇文档中的一个词的时候这两次抛骰子的动作是紧邻轮换进行的。如果语料中一共有$N$个词，则上帝一共要抛$2N$次骰子，轮换的抛doc-topic骰子和topic-word骰子。但实际上有一些抛骰子的顺序是可以交换的，我们可以等价的调整$2N$次抛骰子的次序：前$N$次只抛doc-topic骰子得到语料中所有词的topics,然后基于得到的每个词的topic编号，后$N$次只抛topic-word骰子生成$N$个word。于是上帝在玩 LDA 游戏的时候，可以等价的按照如下过程进行：
以上游戏是先生成了语料中所有词的topic, 然后对每个词在给定topic的条件下生成 word.在语料中所有词的 topic已经生成的条件下，任何两个word的生成动作都是可交换的。于是我们把语料中的词进行交换，把具有相同topic的词放在一起:
\begin{equation} \begin{split} \vec{\mathbf{w}}’ &= (\vec w_{(1)}, \cdots, \vec w_{(K)}) \\ \vec{\mathbf{z}}’ &= (\vec z_{(1)}, \cdots, \vec z_{(K)}) \end{split} \end{equation}
其中，$\vec w_{(k)}$表示这些词都是由第$k$个topic生成的，$\vec z_{(k)}$对应于这些词的topic编号，所以$\vec z_{(k)}$中的分量都是$k$。
对应于概率图中的第二个物理过程$\vec\beta \rightarrow \vec\varphi_k \rightarrow w_{m,n} | k=z_{m,n}$，在$k=z_{m,n}$的限制下，语料中任何两个由 topic $k$生成的词都是可交换的，即便他们不再同一个文档中，所以我们此处不再考虑文档的概念，转而考虑由同一个topic生成的词。考虑如下过程 $\vec\beta \rightarrow \vec\varphi_k \rightarrow \vec w_{(k)}$，容易看出， 此时$\vec\beta \rightarrow \vec\varphi_k$对应于 Dirichlet分布， $\vec\varphi_k \rightarrow \vec w_{(k)}$对应于 Multinomial 分布， 所以整体也还是一个Dirichlet-Multinomial共轭结构；
同样的，我们可以得到:
\begin{equation} p(\vec w_{(k)} |\vec\beta) = \frac{\Delta( n_k+\vec\beta)}{\Delta(\vec\beta)} \end{equation}
其中$ n_k = (n_{k}^{(1)}, \cdots, n_{k}^{(V)})$， $n_{k}^{(t)}$表示第$k$个topic产生的词中 word $t$的个数。进一步，利用Dirichlet-Multiomial共轭结构，我们得到参数$\vec\varphi_k$的后验分布恰好是$Dir( \vec\varphi_k| n_k + \vec\beta)$.
而语料中$K$个topics生成words的过程相互独立，所以我们得到$K$个相互独立的Dirichlet-Multinomial共轭结构，从而我们可以得到整个语料中词生成概率:
\begin{equation} \begin{split} p(\vec{\mathbf{w}} |\vec{\mathbf{z}},\vec\beta) &= p(\vec{\mathbf{w}}’ |\vec{\mathbf{z}}’,\vec\beta) \notag \\ &= \prod_{k=1}^K p(\vec w_{(k)} | \vec z_{(k)}, \vec\beta) \notag \\ &= \prod_{k=1}^K \frac{\Delta( n_k+\vec\beta)}{\Delta(\vec\beta)} \quad\quad (**) \end{split} \end{equation}
结合(*)和(**)于是我们得到:
\begin{equation} \begin{split} p(\vec{\mathbf{w}},\vec{\mathbf{z}} |\vec\alpha, \vec\beta) &= p(\vec{\mathbf{w}} |\vec{\mathbf{z}}, \vec\beta) p(\vec{\mathbf{z}} |\vec\alpha) \notag \\ &= \prod_{k=1}^K \frac{\Delta( n_k+\vec\beta)}{\Delta(\vec\beta)} \prod_{m=1}^M \frac{\Delta( n_m+\vec\alpha)}{\Delta(\vec\alpha)} \quad\quad (***) \end{split} \end{equation}
此处的符号表示稍微不够严谨, 向量$ n_k$, $ n_m$都用$n$表示， 主要通过下标进行区分， $k$下标为topic编号, $m$下标为文档编号。
Gibbs Sampling
有了联合分布$p(\vec{\mathbf{w}},\vec{\mathbf{z}})$, 万能的MCMC算法就可以发挥作用了！于是我们可以考虑使用Gibbs Sampling算法对这个分布进行采样。当然由于$\vec{\mathbf{w}}$是观测到的已知数据，只有$\vec{\mathbf{z}}$是隐含的变量，所以我们真正需要采样的是分布$p(\vec{\mathbf{z}}|\vec{\mathbf{w}})$。在Gregor Heinrich 那篇很有名的LDA 模型科普文章Parameter estimation for text analysis中，是基于(***) 式推导Gibbs Sampling 公式的。此小节中我们使用不同的方式，主要是基于Dirichlet-Multinomial共轭来推导 Gibbs Sampling 公式，这样对于理解采样中的概率物理过程有帮助。
语料库$\vec{\mathbf{z}}$中的第$i$个词我们记为$z_i$, 其中$i=(m,n)$是一个二维下标，对应于第$m$篇文档的第$n$个词，我们用$\neg i$表示去除下标为$i$的词。那么按照 Gibbs Sampling 算法的要求，我们要求得任一个坐标轴$i$对应的条件分布$p(z_i = k|\vec{\mathbf{z}}_{\neg i}, \vec{\mathbf{w}})$。假设已经观测到的词$w_i=t$, 则由贝叶斯法则，我们容易得到:
由于$z_i=k,w_i=t$只涉及到第$m$篇文档和第$k$个topic，所以上式的条件概率计算中, 实际上也只会涉及到如下两个Dirichlet-Multinomial 共轭结构:
- $\vec\alpha \rightarrow \vec\theta_m \rightarrow \vec z_{m}$;
- $\vec\beta \rightarrow \vec\varphi_k \rightarrow \vec w_{(k)}$;
其它的$M+K−2$个Dirichlet-Multinomial共轭结构和$z_i=k,w_i=t$是独立的。由于在语料去掉第$i$个词对应的 $(z_i,w_i)$，并不改变我们之前讨论的$M+K$个Dirichlet-Multinomial共轭结构，只是某些地方的计数会减少。所以$\vec{\theta}_m, \vec\varphi_k$的后验分布都是Dirichlet:
使用上面两个式子，把以上想法综合一下，我们就得到了如下的Gibbs Sampling公式的推导:
以上推导估计是整篇文章中最复杂的数学了，表面上看上去复杂，但是推导过程中的概率物理意义是简单明了的：$z_i=k,w_i=t$的概率只和两个Dirichlet-Multinomial共轭结构关联。而最终得到的$\hat{\theta}{mk}, \hat{\varphi}{kt}$就是对应的两个Dirichlet后验分布在贝叶斯框架下的参数估计。借助于前面介绍的Dirichlet 参数估计的公式 ，我们有:
\begin{equation} \begin{split} \hat\theta_{mk} &= \frac{n_{m,\neg i}^{(k)} + \alpha_k}{\sum_{k=1}^K (n_{m,\neg i}^{(k)} + \alpha_k)} \\ \hat\varphi_{kt} &= \frac{n_{k,\neg i}^{(t)} + \beta_t}{\sum_{t=1}^V (n_{k,\neg i}^{(t)} + \beta_t)} \end{split} \end{equation}
于是，我们最终得到了LDA模型的Gibbs Sampling公式:
这个公式是很漂亮的， 右边其实就是$p(topic|doc)⋅p(word|topic)$,这个概率其实是doc→topic→word的路径概率，由于topic 有$K$个，所以Gibbs Sampling 公式的物理意义其实就是在这$K$条路径中进行采样。
Training and Inference
有了LDA模型，当然我们的目标有两个:
- 估计模型中的参数$\vec\varphi_1, \cdots, \vec\varphi_K$和$\vec{\theta}_1, \cdots, \vec{\theta}_M$；
- 对于新来的一篇文档$doc_{new}$，我们能够计算这篇文档的topic分布$\vec{\theta}_{new}$。
有了Gibbs Sampling公式， 我们就可以基于语料训练LDA模型，并应用训练得到的模型对新的文档进行topic 语义分析。训练的过程就是获取语料中的$(z,w)$的样本，而模型中的所有的参数都可以基于最终采样得到的样本进行估计。训练的流程很简单:
对于Gibbs Sampling算法实现的细节，请参考Gregor Heinrich的 Parameter estimation for text analysis 中对算法的描述，以及PLDA(http://code.google.com/p/plda)的代码实现，此处不再赘述。
由这个topic-word频率矩阵我们可以计算每一个$p(word|topic)$概率，从而算出模型参数$\vec\varphi_1, \cdots, \vec\varphi_K$, 这就是上帝用的$K$个topic-word骰子。当然，语料中的文档对应的骰子参数$\vec{\theta}_1, \cdots, \vec{\theta}_M$在以上训练过程中也是可以计算出来的，只要在Gibbs Sampling收敛之后，统计每篇文档中的topic的频率分布，我们就可以计算每一个$p(topic|doc)$概率，于是就可以计算出每一个$\vec{\theta}_m$。由于参数$\vec{\theta}_m$是和训练语料中的每篇文档相关的，对于我们理解新的文档并无用处，所以工程上最终存储LDA模型时候一般没有必要保留。通常，在LDA模型训练的过程中，我们是取Gibbs Sampling收敛之后的$n$个迭代的结果进行平均来做参数估计，这样模型质量更高。
有了LDA的模型，对于新来的文档 $doc_{new}$, 我们如何做该文档的topic语义分布的计算呢？基本上inference的过程和training的过程完全类似。对于新的文档， 我们只要认为Gibbs Sampling公式中的$\hat\varphi_{kt}$部分是稳定不变的，是由训练语料得到的模型提供的，所以采样过程中我们只要估计该文档的topic分布$\vec{\theta}_{new}$就好了。
后记LDA
对于专业做机器学习的兄弟而言，只能算是一个简单的Topic Model。但是对于互联网中做数据挖掘、语义分析的工程师，LDA 的门槛并不低。 LDA 典型的属于这样一种机器学习模型：要想理解它，需要比较多的数学背景，要在工程上进行实现，却相对简单。 Gregor Heinrich 的LDA 模型科普文章 Parameter estimation for text analysis 写得非常的出色，这是学习 LDA 的必看文章。不过即便是这篇文章，对于工程师也是有门槛的。我写的这个科普最好对照 Gregor Heinrich 的这篇文章来看， 我用的数学符号也是尽可能和这篇文章保持一致。这份LDA 科普是基于给组内兄弟做报告的 ppt 整理而成的，说是科普其实也不简单，涉及到的数学还是太多。在工业界也混了几年，经常感觉到工程师对于学术界的玩的模型有很强的学习和尝试的欲望，只是学习成本往往太高。所以我写 LDA 的初衷就是写给工业界的工程师们看的，希望把学术界玩的一些模型用相对通俗的方式介绍给工程师；如果这个科普对于读研究生的一些兄弟姐妹也有所启发，只能说那是一个 side effect :-)。我个人很喜欢LDA ，它是在文本建模中一个非常优雅的模型，相比于很多其它的贝叶斯模型， LDA 在数学推导上简洁优美。学术界自 2003 年以来也输出了很多基于LDA 的 Topic Model 的变体，要想理解这些更加高级的 Topic Model, 首先需要很好的理解标准的 LDA 模型。在工业界， Topic Model 在 Google、Baidu 等大公司的产品的语义分析中都有着重要的应用；所以Topic Model 对于工程师而言，这是一个很有应用价值、值得学习的模型。我接触 Topic Model 的时间不长，主要是由于2年前和 PLDA 的作者 Wangyi 一起合作的过程中，从他身上学到了很多 Topic Model 方面的知识。关于 LDA 的相关知识，其实可以写的还有很多：如何提高 LDA Gibbs Sampling 的速度、如何优化超参数、如何做大规模并行化、LDA 的应用、LDA 的各种变体…… 不过我的主要目标还是科普如何理解标准的LDA模型。学习一个模型的时候我喜欢追根溯源，常常希望把模型中的每一个数学推导的细节搞明白，把公式的物理意义想清楚，不过数学推导本身并不是我想要的，把数学推导还原为物理过程才是我乐意做的事。最后引用一下物理学家费曼的名言结束 LDA 的数学科普：
What I cannot create, I do not understand. — Richard Feynman