一、简介

隐含狄利克雷分布（Latent Dirichlet Allocation，简称LDA）是由 David M. Blei、Andrew Y. Ng、Michael I. Jordan 在2003年提出的，是一种词袋模型，它认为文档是一组词构成的集合，词与词之间是无序的。一篇文档可以包含多个主题，文档中的每个词都是由某个主题生成的，LDA给出文档属于每个主题的概率分布，同时给出每个主题上词的概率分布。LDA是一种无监督学习，在文本主题识别、文本分类、文本相似度计算和文章相似推荐等方面都有应用。本文将从贝叶公式、Gamma函数、二项分布、Beta分布、多项式分布、Dirichlet分布、共轭先验分布、马氏链及其平稳分布、MCMC、Gibbs Sampling、EM算法、Unigram Model、贝叶斯Unigram Model、PLSA、LDA 几方面介绍LDA模型，需要读者具备一定的概率论和微积分知识。

二、基础知识

1.1 贝叶公式

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

从贝叶斯观点看，样本 $X = \left( {{x_1}, \cdot \cdot \cdot ,{x_n}} \right)$ 的产生要分两步进行，首先设想从先验分布 $p\left( \theta \right)$ 产生一个样本 $\theta '$ ，这一步是“老天爷”做的，人们是看不到的，故用“设想”二字。第二步是从总体分布 $p\left( {X|\theta '} \right)$ 产生一个样本 $X = \left( {{x_1}, \cdot \cdot \cdot ,{x_n}} \right)$ ，这个样本是具体的，人们能看得到的，此样本 $X$ 发生的概率是与如下联合密度函数成正比。

$\begin{align*} p\left( {X|\theta '} \right) = \prod\limits_{i = 1}^n {p\left( {{x_i}|\theta '} \right)} \tag{1} \end{align*}$

这个联合密度函数是综合了总体信息和样本信息，常称为似然函数，记为 $L\left( {\theta '} \right)$ 。

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

$\begin{align*} p\left( {X,\theta } \right) = p\left( {X|\theta } \right)p\left( \theta \right) \tag{2} \end{align*}$

我们的任务是要对未知数 $\theta$ 作出统计推断，在没有样本信息时，人们只能根据先验分布对 $\theta$ 作出推断。在有样本观察值 $X = \left( {{x_1}, \cdot \cdot \cdot ,{x_n}} \right)$ 之后，我们应该依据 $p\left( {X,\theta } \right)$ 对 $\theta$ 作出推断，为此我们把 $p\left( {X,\theta } \right)$ 作如下分解：

$\begin{align*} p\left( {X,\theta } \right){\rm{ = }}p\left( {\theta |X} \right)p\left( X \right) \tag{3} \end{align*}$

其中 $p\left( X \right)$ 是 $X$ 的边缘密度函数。

$\begin{align*} p\left( X \right) = \int_\Theta {p\left( {X,\theta } \right)} d\theta = \int_\Theta {p\left( {X|\theta } \right)p\left( \theta \right)} d\theta \tag{4} \end{align*}$

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

$\begin{align*} p\left( {\theta |X} \right) &= \frac{{p\left( {X,\theta } \right)}}{{p\left( X \right)}} = \frac{{p\left( {X|\theta } \right)p\left( \theta \right)}}{{\int_\Theta {p\left( {X|\theta } \right)p\left( \theta \right)} d\theta }} \\ posterior &= \frac{{{\rm{likelihood}} \cdot {\rm{prior}}}}{{{\rm{evidence}}}} \tag{5} \end{align*}$

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

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

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

$\begin{align*} \hat \theta {\rm{ = }}\mathop {\arg \max }\limits_\theta \sum\limits_{i = 1}^n {\log p\left( {{x_i}|\theta } \right)} \tag{6} \end{align*}$

通常可以令导数为 0 求得 $\theta$ 的值。ML估计不会把先验知识考虑进去，很容易出现过拟合的现象。我们举个例子，抛一枚硬币，假设正面向上的概率为 $p$ ，抛了 $N$ 次，正面出现 ${n^{\left( 1 \right)}}$ 次，反面出现 ${n^{\left( 0 \right)}}$ 次， $c = 1$ 表示正面， $c = 0$ 表示反面，我们用ML估计：

$\begin{align*} L\left( p \right) &= \sum\limits_{i = 1}^N {\log } p\left( {C = {c_i}|p} \right)\\ &= {n^{\left( 1 \right)}}\log p\left( {C = 1|p} \right) + {n^{\left( 0 \right)}}\log p\left( {C = 0|p} \right)\\ &= {n^{\left( 1 \right)}}\log p + {n^{\left( 0 \right)}}\log \left( {1 - p} \right)\\ \frac{{\partial L}}{{\partial p}} &= 0 \Rightarrow {{\hat p}_{ML}} = \frac{{{n^{\left( 1 \right)}}}}{{{n^{\left( 1 \right)}} + {n^{\left( 0 \right)}}}} \tag{7} \end{align*}$

如果 ${n^{\left( 1 \right)}} = 14$ , ${n^{\left( 0 \right)}} = 6$ ，则 ${\hat p_{ML}} = 0.7$ ，似乎比我们认知的 0.5 高了很多。

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

$\begin{align*} \hat \theta &= \mathop {\arg \max }\limits_\theta p\left( {\theta |X} \right)\\ &= \mathop {\arg \max }\limits_\theta \frac{{p\left( {X|\theta } \right)p\left( \theta \right)}}{{p\left( X \right)}}\\ &= \mathop {\arg \max }\limits_\theta p\left( {X|\theta } \right)p\left( \theta \right)\\ &= \mathop {\arg \max }\limits_\theta \left\{ {\sum\limits_{i = 1}^n {\log p\left( {{x_i}|\theta } \right)} + \log p\left( \theta \right)} \right\} \tag{8} \end{align*}$

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

$\begin{align*} p\left( {p|\alpha ,\beta } \right) &= \frac{1}{{B\left( {\alpha ,\beta } \right)}}{p^{\alpha - 1}}{\left( {1 - p} \right)^{\beta - 1}}\\ L\left( p \right) &= {n^{\left( 1 \right)}}\log p + {n^{\left( 0 \right)}}\log \left( {1 - p} \right) + \left( {\alpha - 1} \right)\log p + \left( {\beta - 1} \right)\log \left( {1 - p} \right) - \log B\left( {\alpha ,\beta } \right)\\ \frac{{\partial L}}{{\partial p}} &= \frac{{{n^{\left( 1 \right)}}}}{p} - \frac{{{n^{\left( 0 \right)}}}}{{1 - p}} + \frac{{\alpha - 1}}{p} - \frac{{\beta - 1}}{{1 - p}} \Rightarrow {{\hat p}_{MAP}} = \frac{{{n^{\left( 1 \right)}} + \alpha - 1}}{{{n^{\left( 1 \right)}} + {n^{\left( 0 \right)}} + \alpha {\rm{ + }}\beta - 2}} \tag{9} \end{align*}$

我们取 $\alpha {\rm{ = }}5,\beta = 5$，即 $p$ 以最大的概率取0.5，得到 ${\hat p_{MAP}} = 0.64$。

3. 贝叶斯估计

前面介绍的 ML 和 MAP 属于点估计，贝叶斯估计不再把参数 $\theta$ 看成一个未知的确定值，而是看成未知的随机变量，利用贝叶斯定理结合新的样本信息和参数 $\theta$ 的先验分布，来得到 $\theta$的新的概率分布（后验分布）。贝叶斯估计的本质是通过贝叶斯决策得到参数 $\theta$ 的最优估计 $\hat \theta$ ，使得贝叶斯期望损失最小。贝叶斯期望损失为：

$\begin{align*} R\left( {\hat \theta |X} \right) = \int_\Theta {\lambda \left( {\hat \theta ,\theta } \right)} p\left( {\theta |X} \right)d\theta \tag{10} \end{align*}$

$\lambda \left( {\hat \theta ,\theta } \right)$ 是损失函数，我们希望 $R\left( {\hat \theta |X} \right)$ 最小。如果 $\lambda \left( {\hat \theta ,\theta } \right) = \left\| {\hat \theta - \theta } \right\|_2^2$ ，则：

$\begin{align*} \frac{{\partial R\left( {\hat \theta |X} \right)}}{{\partial \hat \theta }} &= \int_\Theta {2\left( {\hat \theta - \theta } \right)} p\left( {\theta |X} \right)d\theta = 0\\ \Rightarrow \hat \theta &= \int_\Theta \theta p\left( {\theta |X} \right)d\theta \tag{11} \end{align*}$

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

• 确定参数 $\theta$ 的先验分布 $p\left( \theta \right)$
• 利用贝叶斯公式，求$\theta$的后验分布：

$\begin{align*} p\left( {\theta |X} \right) = \frac{{p\left( {X|\theta } \right)p\left( \theta \right)}}{{\int_\Theta {p\left( {X|\theta } \right)p\left( \theta \right)d\theta } }} \tag{12} \end{align*}$

• 求出贝叶斯的估计值：

$\begin{align*} \hat \theta = \int_\Theta \theta p\left( {\theta |X} \right)d\theta \tag{13} \end{align*}$

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

$\begin{align*} p\left( {p|C,\alpha ,\beta } \right) &= \frac{{\prod\nolimits_{i = 1}^N {p\left( {C = {c_i}|p} \right)p\left( {p|\alpha ,\beta } \right)} }}{{\int_0^1 {\prod\nolimits_{i = 1}^N {p\left( {C = {c_i}|p} \right)p\left( {p|\alpha ,\beta } \right)} dp} }}\\ &= \frac{{{p^{n\left( 1 \right)}}{{\left( {1 - p} \right)}^{n\left( 0 \right)}}\frac{1}{{B\left( {\alpha ,\beta } \right)}}{p^{\alpha - 1}}{{\left( {1 - p} \right)}^{\beta - 1}}}}{{\int_0^1 {{p^{n\left( 1 \right)}}{{\left( {1 - p} \right)}^{n\left( 0 \right)}}\frac{1}{{B\left( {\alpha ,\beta } \right)}}{p^{\alpha - 1}}{{\left( {1 - p} \right)}^{\beta - 1}}dp} }}\\ &= \frac{{{p^{n\left( 1 \right) + \alpha - 1}}{{\left( {1 - p} \right)}^{n\left( 0 \right) + \beta - 1}}}}{{\int_0^1 {{p^{n\left( 1 \right) + \alpha - 1}}{{\left( {1 - p} \right)}^{n\left( 0 \right) + \beta - 1}}dp} }}\\ &= Beta\left( {p|n\left( 1 \right) + \alpha ,n\left( 0 \right) + \beta } \right) \tag{14} \end{align*}$

其中 $B\left( {\alpha ,\beta } \right) = \int_0^1 {{p^{\alpha - 1}}{{\left( {1 - p} \right)}^{\beta - 1}}dp}$，所以可以得：

$\begin{align*} \hat p = \frac{{{n^{\left( 1 \right)}} + \alpha }}{{{n^{\left( 1 \right)}} + {n^{\left( 0 \right)}} + \alpha + \beta }} \tag{15} \end{align*}$

1.2 Gamma函数

$\begin{align*} \Gamma \left( x \right) = \int_0^\infty {{t^{x - 1}}{e^{ - t}}dt} \tag{16} \end{align*}$

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

$\begin{align*} \Gamma \left( {x + 1} \right) &= \int_0^\infty {{t^x}{e^{ - t}}dt} \\ &= - \int_0^\infty {{t^x}d{e^{ - t}}} \\ &= - \left[ {{t^x}{e^{ - t}}} \right]_0^\infty + \int_0^\infty {{e^{ - t}}d{t^x}} \\ &= x\int_0^\infty {{t^{x - 1}}{e^{ - t}}d} t = x\Gamma \left( x \right) \tag{17} \end{align*}$

$\Gamma \left( x \right)$ 函数可以当成是阶乘在实数集上的延拓，$\Gamma \left( n \right) = \left( {n - 1} \right)!$。

1.3 二项分布

在概率论中，试验 $E$ 只有两个可能结果：$A$ 及 $\bar A$ ，则称 $E$ 为伯努利(Bernoulli)试验。设 $p\left( A \right) = p$，则 $p\left( {\bar A} \right) = 1 - p$。将 $E$ 独立重复地进行 $n$ 次，则称这一串重复的独立试验为 $n$ 重伯努利试验，这里重复是指在每次试验中 $p\left( A \right) = p$ 保持不变，独立是指各次试验的结果互不影响。以 $X$ 表示 $n$ 重伯努利试验中事件 $A$ 发生的次数，称随机变量 $X$ 服从参数为 $n,p$ 的二项分布，记为 $X \sim B\left( {n,p} \right)$。

$\begin{align*} p\left( {X = k} \right) = \left( {\begin{array}{*{20}{c}} n\\ k \end{array}} \right){p^k}{\left( {1 - p} \right)^{n - k}} \tag{18} \end{align*}$

1.4 Beta分布

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

$\begin{align*} Beta\left( {p|\alpha ,\beta } \right) &= \frac{{{p^{\alpha - 1}}{{\left( {1 - p} \right)}^{\beta - 1}}}}{{\int_0^1 {{p^{\alpha - 1}}{{\left( {1 - p} \right)}^{\beta - 1}}dp} }}\\ &= \frac{{\Gamma \left( {\alpha + \beta } \right)}}{{\Gamma \left( \alpha \right)\Gamma \left( \beta \right)}}{p^{\alpha - 1}}{\left( {1 - p} \right)^{\beta - 1}} \tag{19} \end{align*}$

1）$B\left( {\alpha ,\beta } \right) = \frac{{\Gamma \left( \alpha \right)\Gamma \left( \beta \right)}}{{\Gamma \left( {\alpha + \beta } \right)}} = \int_0^1 {{p^{\alpha - 1}}{{\left( {1 - p} \right)}^{\beta - 1}}dp}$。$\Gamma \left( \alpha \right) = \int_0^\infty {{x^{\alpha - 1}}{e^{ - x}}dx{\rm{ }}} ,\Gamma \left( \beta \right) = \int_0^\infty {{y^{\beta - 1}}{e^{ - y}}dy}$
证明：

$\begin{align*} \Gamma \left( \alpha \right)\Gamma \left( \beta \right) &= \int_0^\infty {{x^{\alpha - 1}}{e^{ - x}}dx{\rm{ }}} \int_0^\infty {{y^{\beta - 1}}{e^{ - y}}dy} \\ &= \int_0^\infty {\int_0^\infty {{x^{\alpha - 1}}{y^{\beta - 1}}{e^{ - \left( {x + y} \right)}}dydx} } \tag{20} \end{align*}$

令$t = x + y$，当$y = 0,t = x$；$y = \infty ,t = \infty$，可得：

$\begin{align*} \Gamma \left( \alpha \right)\Gamma \left( \beta \right) &= \int_0^\infty {\int_x^\infty {{x^{\alpha - 1}}{{\left( {t - x} \right)}^{\beta - 1}}{e^{ - t}}dtdx} } \\ &= \int_0^\infty {\int_0^t {{x^{\alpha - 1}}{{\left( {t - x} \right)}^{\beta - 1}}{e^{ - t}}dxdt} } \tag{21} \end{align*}$

令$x = \mu t,\mu \in \left[ {0,1} \right]$，可得：

$\begin{align*} \Gamma \left( \alpha \right)\Gamma \left( \beta \right) &= \int_0^\infty {\int_0^1 {{{\left( {\mu t} \right)}^{\alpha - 1}}{{\left( {t - \mu t} \right)}^{\beta - 1}}{e^{ - t}}d\mu tdt} } \\ & = \int_0^\infty {{t^{\alpha + \beta - 1}}{e^{ - t}}dt} \int_0^1 {{\mu ^{\alpha - 1}}{{\left( {1 - \mu } \right)}^{\beta - 1}}d\mu } \\ & = \Gamma \left( {\alpha + \beta } \right)\int_0^1 {{\mu ^{\alpha - 1}}{{\left( {1 - \mu } \right)}^{\beta - 1}}d\mu } \\ \frac{{\Gamma \left( \alpha \right)\Gamma \left( \beta \right)}}{{\Gamma \left( {\alpha + \beta } \right)}} &= \int_0^1 {{\mu ^{\alpha - 1}}{{\left( {1 - \mu } \right)}^{\beta - 1}}d\mu } \tag{22} \end{align*}$

2）期望$E\left( p \right) = \frac{\alpha }{{\alpha + \beta }}$
证明：

$\begin{align*} E\left( p \right) &= \int_0^1 {p\frac{{\Gamma \left( {\alpha + \beta } \right)}}{{\Gamma \left( \alpha \right)\Gamma \left( \beta \right)}}{p^{\alpha - 1}}{{\left( {1 - p} \right)}^{\beta - 1}}dp} \\ & = \frac{{\Gamma \left( {\alpha + \beta } \right)}}{{\Gamma \left( \alpha \right)\Gamma \left( \beta \right)}}\int_0^1 {{p^\alpha }{{\left( {1 - p} \right)}^{\beta - 1}}dp} \\ & = \frac{{\Gamma \left( {\alpha + \beta } \right)}}{{\Gamma \left( \alpha \right)\Gamma \left( \beta \right)}}\frac{{\Gamma \left( {\alpha {\rm{ + }}1} \right)\Gamma \left( \beta \right)}}{{\Gamma \left( {\alpha + \beta {\rm{ + }}1} \right)}}\int_0^1 {\frac{{\Gamma \left( {\alpha + \beta {\rm{ + }}1} \right)}}{{\Gamma \left( {\alpha {\rm{ + }}1} \right)\Gamma \left( \beta \right)}}{p^\alpha }{{\left( {1 - p} \right)}^{\beta - 1}}dp} \\ & = \frac{{\Gamma \left( {\alpha + \beta } \right)}}{{\Gamma \left( \alpha \right)\Gamma \left( \beta \right)}}\frac{{\Gamma \left( {\alpha {\rm{ + }}1} \right)\Gamma \left( \beta \right)}}{{\Gamma \left( {\alpha + \beta {\rm{ + }}1} \right)}}\int_0^1 {Beta\left( {p|\alpha + 1,\beta } \right)dp} \\ &= \frac{{\Gamma \left( {\alpha + \beta } \right)}}{{\Gamma \left( \alpha \right)\Gamma \left( \beta \right)}}\frac{{\alpha \Gamma \left( \alpha \right)\Gamma \left( \beta \right)}}{{\left( {\alpha + \beta } \right)\Gamma \left( {\alpha + \beta } \right)}}\\ &= \frac{\alpha }{{\alpha + \beta }} \tag{23} \end{align*}$

1.5 多项式分布

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

$\begin{align*} Mult\left( {\vec n|\vec p,N} \right) = \left( {\begin{array}{*{20}{c}} N\\ {\vec n} \end{array}} \right)\prod\limits_{k = 1}^K {p_k^{{n_k}}} \tag{24} \end{align*}$

其中$\sum\limits_{k = 1}^K {{n_k}} = N,\sum\limits_{k = 1}^K {{p_k}} = 1\left( {\begin{array}{*{20}{c}} N\\ {\vec n} \end{array}} \right) = \frac{{N!}}{{\prod\nolimits_k {{n_k}!} }}$。

1.6 Dirichlet分布

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

$\begin{align*} Dir\left( {\vec p|\vec \alpha } \right) &= \frac{{\Gamma \left( {\sum\nolimits_{k = 1}^K {{\alpha _k}} } \right)}}{{\prod\nolimits_{k = 1}^K {\Gamma \left( {{\alpha _k}} \right)} }}\prod\limits_{k = 1}^K {p_k^{{\alpha _k} - 1}} \\ &= \frac{1}{{\Delta \left( {\vec \alpha } \right)}}\prod\limits_{k = 1}^K {p_k^{{\alpha _k} - 1}} \tag{25} \end{align*}$

1）$\Delta \left( {\vec \alpha } \right){\rm{ = }}\frac{{\prod\nolimits_{k = 1}^K {\Gamma \left( {{\alpha _k}} \right)} }}{{\Gamma \left( {\sum\nolimits_{k = 1}^K {{\alpha _k}} } \right)}}{\rm{ = }}\int_0^1 {\prod\limits_{k = 1}^K {p_k^{{\alpha _k} - 1}} d\vec p}$

2）期望 $E\left( {\vec p} \right) = \left( {\frac{{{\alpha _1}}}{{\sum\nolimits_{k = 1}^K {{\alpha _k}} }},\frac{{{\alpha _2}}}{{\sum\nolimits_{k = 1}^K {{\alpha _k}} }}, \cdot \cdot \cdot ,\frac{{{\alpha _K}}}{{\sum\nolimits_{k = 1}^K {{\alpha _k}} }}} \right)$。

1.7 共轭先验分布

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

1．Beta-Binomial共轭
1）先验分布

$\begin{align*} Beta\left( {p|\alpha ,\beta } \right) = \frac{1}{{B\left( {\alpha ,\beta } \right)}}{p^{\alpha - 1}}{\left( {1 - p} \right)^{\beta - 1}} \tag{26} \end{align*}$

2）二项式似然函数

$\begin{align*} B\left( {{n_1},{n_2}|p} \right) = \left( {\begin{array}{*{20}{c}} n\\ {{n_1}} \end{array}} \right){p^{{n_1}}}{\left( {1 - p} \right)^{{n_2}}} \tag{27} \end{align*}$

3）后验分布

$\begin{align*} \frac{{B\left( {{n_1},{n_2}|p} \right)Beta\left( {p|\alpha ,\beta } \right)}}{{\int_0^1 {B\left( {{n_1},{n_2}|p} \right)Beta\left( {p|\alpha ,\beta } \right)dp} }}&= \frac{{{p^{\alpha + {n_1} - 1}}{{\left( {1 - p} \right)}^{\beta + {n_2} - 1}}}}{{\int_0^1 {{p^{\alpha + {n_1} - 1}}{{\left( {1 - p} \right)}^{\beta + {n_2} - 1}}dp} }}\\ &= \frac{{{p^{\alpha + {n_1} - 1}}{{\left( {1 - p} \right)}^{\beta + {n_2} - 1}}}}{{B\left( {\alpha {\rm{ + }}{n_1},\beta + {n_2}} \right)}} \sim Beta\left( {\alpha {\rm{ + }}{n_1},\beta + {n_2}} \right) \tag{28} \end{align*}$

即可以表达为 $Beta\left( {p|\alpha ,\beta } \right){\rm{ + }}B\left( {{n_1},{n_2}|p} \right){\rm{ = }}Beta\left( {p|\alpha {\rm{ + }}{n_1},\beta + {n_2}} \right)$，取一个特殊情况理解
$Beta\left( {p|1,1} \right){\rm{ + }}B\left( {\alpha - 1,\beta - 1|p} \right){\rm{ = }}Beta\left( {p|\alpha ,\beta } \right)$，$Beta\left( {p|1,1} \right)$ 恰好是均匀分布 $uniform\left( {0,1} \right)$，假设有一个不均匀的硬币抛出正面的概率为 $p$，抛出 $n$ 次后出现正面和反面的次数分别是 ${n_1}$ 和 ${n_2}$ ，开始我们对硬币不均匀性一无所知，所以应该假设 $p \sim uniform\left( {0,1} \right)$ ，当有了试验样本，我们加入样本信息对 $p$ 的分布进行修正, $p$ 的分布由均匀分布变为 $Beta$ 分布。

2．Dirichlet-Multinomial共轭

1）先验分布

$\begin{align*} Dir\left( {\vec p|\vec \alpha } \right) = \frac{1}{{\Delta \left( {\vec \alpha } \right)}}\prod\limits_{k = 1}^K {p_k^{{\alpha _k} - 1}} \tag{29} \end{align*}$

2）多项分布似然函数

$\begin{align*} Mult\left( {\vec n|\vec p,N} \right) = \left( {\begin{array}{*{20}{c}} N\\ {\vec n} \end{array}} \right)\prod\limits_{k = 1}^K {p_k^{{n_k}}} \tag{30} \end{align*}$

3）后验分布

$\begin{align*} \frac{{Dir\left( {\vec p|\vec \alpha } \right)Mult\left( {\vec n|\vec p,N} \right)}}{{\int_0^1 {Dir\left( {\vec p|\vec \alpha } \right)Mult\left( {\vec n|\vec p,N} \right)d\vec p} }}&= \frac{{\prod\limits_{k = 1}^K {p_k^{{\alpha _k} + {n_k} - 1}} }}{{\int_0^1 {\prod\limits_{k = 1}^K {p_k^{{\alpha _k} + {n_k} - 1}} d\vec p} }}\\ &= \frac{{\prod\limits_{k = 1}^K {p_k^{{\alpha _k} + {n_k} - 1}} }}{{\Delta \left( {\vec \alpha + \vec n} \right)}} \sim Dir\left( {\vec p|\vec \alpha + \vec n} \right) \tag{31} \end{align*}$

即可以表达为 $Dir\left( {\vec p|\vec \alpha } \right){\rm{ + }}Mult\left( {\vec n|\vec p,N} \right){\rm{ = }}Dir\left( {\vec p|\vec \alpha {\rm{ + }}\vec n} \right)$

1.8 马氏链及其平稳分布

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

$\begin{align*} P\left( {{X_{t + 1}} = x|{X_t},{X_{t - 1}}, \cdot \cdot \cdot } \right) = P\left( {{X_{t + 1}} = x|{X_t}} \right) \tag{32} \end{align*}$

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

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

$\begin{align*} P = \left( {\begin{array}{*{20}{c}} {0.9}&{0.075}&{0.025}\\ {0.15}&{0.8}&{0.05}\\ {0.25}&{0.25}&{0.5} \end{array}} \right) \end{align*}$

假设初始概率分布为 ${\pi _0}{\rm{ = }}\left[ {\begin{array}{*{20}{c}} {0.3}&{0.4}&{0.3} \end{array}} \right]$，${\pi _1}{\rm{ = }}{\pi _0}P,{\pi _2}{\rm{ = }}{\pi _1}P = {\pi _0}{P^2}, \cdot \cdot \cdot ,{\pi _n}{\rm{ = }}{\pi _{n - 1}}P = {\pi _0}{P^n}$。从第60轮开始 ${\pi _{60}}, \cdot \cdot \cdot ,{\pi _n}$ 的值保持不变，为 $\left[ {\begin{array}{*{20}{c}} {0.625}&{0.3125}&{0.0625} \end{array}} \right]$ 。我们更改初始概率，${\pi _0}{\rm{ = }}\left[ {\begin{array}{*{20}{c}} {0.7}&{0.2}&{0.1} \end{array}} \right]$，从55轮开始 ${\pi _{55}}, \cdot \cdot \cdot ,{\pi _n}$ 的值保持不变，为 $\left[ {\begin{array}{*{20}{c}} {0.625}&{0.3125}&{0.0625} \end{array}} \right]$ 。两次给定不同的初始概率分布，最终都收敛到概率分布 $\pi {\rm{ = }}\left[ {\begin{array}{*{20}{c}} {0.625}&{0.3125}&{0.0625} \end{array}} \right]$ ，也就是说收敛的行为和初始概率分布 ${\pi _0}$ 无关，这个收敛的行为主要是由概率转移矩阵 $P$ 决定的，可以计算下 ${P^n}$。

$\begin{align*} {P^{63}} = {P^{64}} = \cdot \cdot \cdot = \left[ {\begin{array}{*{20}{c}} {0.625}&{0.3125}&{0.0625}\\ {0.625}&{0.3125}&{0.0625}\\ {0.625}&{0.3125}&{0.0625} \end{array}} \right] \end{align*}$

当 $n$ 足够大的时候，${P^n}$ 矩阵的每一行都是稳定地收敛到 $\pi {\rm{ = }}\left[ {\begin{array}{*{20}{c}} {0.625}&{0.3125}&{0.0625} \end{array}} \right]$ 这个概率分布。这个收敛现象并不是这个马氏链独有的，而是绝大多数马氏链独有的。关于马氏链的收敛有如下定理：

定理1.1 如果一个非周期马氏链具有转移概率矩阵 $P$，且它的任何两个状态是连通的，那么 $\mathop {\lim }\limits_{n \to \infty } P_{ij}^n$ 存在且与 $i$ 无关，我们有：
1）$\mathop {\lim }\limits_{n \to \infty } P_{ij}^n{\rm{ = }}\pi \left( j \right)$

2）$\mathop {\lim }\limits_{n \to \infty } {P^n} = \left[ {\begin{array}{*{20}{c}} {\pi \left( 1 \right)}&{\pi \left( 2 \right)}&{ \cdot \cdot \cdot }&{\pi \left( j \right)}&{ \cdot \cdot \cdot }\\ {\pi \left( 1 \right)}&{\pi \left( 2 \right)}&{ \cdot \cdot \cdot }&{\pi \left( j \right)}&{ \cdot \cdot \cdot }\\ { \cdot \cdot \cdot }&{ \cdot \cdot \cdot }&{ \cdot \cdot \cdot }&{ \cdot \cdot \cdot }&{ \cdot \cdot \cdot }\\ {\pi \left( 1 \right)}&{\pi \left( 2 \right)}&{ \cdot \cdot \cdot }&{\pi \left( j \right)}&{ \cdot \cdot \cdot }\\ { \cdot \cdot \cdot }&{ \cdot \cdot \cdot }&{ \cdot \cdot \cdot }&{ \cdot \cdot \cdot }&{ \cdot \cdot \cdot } \end{array}} \right]$

3）$\pi \left( j \right) = \sum\limits_{i = 0}^\infty {\pi \left( i \right){P_{ij}}}$

4）$\pi$ 是方程 $\pi P = \pi$ 的唯一非负解，其中$\pi {\rm{ = }}\left[ {\begin{array}{*{20}{c}} {\pi \left( 1 \right)}&{\pi \left( 2 \right)}&{ \cdot \cdot \cdot }&{\pi \left( j \right)}&{ \cdot \cdot \cdot } \end{array}} \right],\sum\limits_{j = 1}^\infty {\pi \left( j \right)} = 1$。
关于上述定理，给出几点解释：
1） 马氏链的状态数可以是有限的，也可以是无限的，因此可以用于连续概率分布和离散概率分布。
2） 非周期马氏链：马氏链的状态转化不是循环的，如果是循环的则永远不会收敛，我们遇到的一般都是非周期马氏链。对于任意某一状态 $i$，$d$ 为集合$\left\{ {n|n \ge 1,P_{ii}^n > 0} \right\}$ 的最大公约数，如果$d = 1$，则该状态为非周期。
3） 任何两个状态是连通的：从任意一个状态可以通过有限步到达其他的任意状态，不会出现条件概率一直为0导致不可达的情况。
4） $\pi$ 称为马氏链的平稳分布。
如果从一个具体的初始状态 ${x_0}$ 开始，沿着马氏链按照概率转移矩阵做跳转，那么可以得到一个转移序列 ${x_0},{x_1}, \cdot \cdot \cdot ,{x_n},{x_{n + 1}}, \cdot \cdot \cdot$，由于马氏链的收敛行为，${x_n},{x_{n + 1}}, \cdot \cdot \cdot$ 都将是平稳分布 $\pi \left( x \right)$ 的样本。

1.9 MCMC

1. 接受-拒绝采样

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

Algorithm 1.1 接受-拒绝采样算法

1. 目标分布 $\pi \left( x \right)$ ，分布 $q\left( x \right)$ 和常数 $M$，通过对 $q\left( x \right)$ 的采样实现对 $\pi \left( x \right)$ 采样，满足：

   • 对 $q\left( x \right)$ 采样比较容易；
   • $q\left( x \right)$ 的形状接近 $\pi \left( x \right)$，且 $\forall x$，$\pi \left( x \right) \le Mq\left( x \right)$。

2. 采样过程

   • 产生样本 $x \sim q\left( x \right)$，和 $u \sim Uniform\left[ {0,1} \right]$
   • 若$u \le \frac{{\pi \left( x \right)}}{{Mq\left( x \right)}}$，则接受样本$x$
   • 则接受的样本服从$\pi \left( x \right)$分布

下面我们来证明下接受-拒绝方法采样得到的样本服从 $\pi \left( x \right)$ 分布。
证明：accept $x$ 服从 $\pi \left( x \right)$ 分布，即 $p\left( {x|accept} \right) = \pi \left( x \right)$ 。

$\begin{align*} p\left( {x|accept} \right) &= \frac{{p\left( {accept|x} \right)p\left( x \right)}}{{p\left( {accept} \right)}}\\ p\left( {accept|x} \right) &= p\left( {u \le \frac{{\pi \left( x \right)}}{{Mq\left( x \right)}}} \right) = \frac{{\pi \left( x \right)}}{{Mq\left( x \right)}}\\ p\left( x \right) &= q\left( x \right)\\ p\left( {accept} \right) &= \int_x {p\left( {accept|x} \right)p\left( x \right)dx} \\ &= \int_x {\frac{{\pi \left( x \right)}}{{Mq\left( x \right)}}q\left( x \right)dx} \\ &= \frac{1}{M}\\ p\left( {x|accept} \right) &= \pi \left( x \right) \tag{33} \end{align*}$

2. MCMC

给定概率分布 $p\left( x \right)$，希望能够生成它对应的样本，由于马氏链能收敛到平稳分布，有一个很好的想法：如果我们能构造一个转移矩阵为 $P$ 的马氏链，使得该马氏链的平稳分布恰好是 $p\left( x \right)$ ，那么我们从任何一个初始状态出发沿着马氏链转移，得到一个转移序列 ${x_0},{x_1}, \cdot \cdot \cdot ,{x_n},{x_{n + 1}}, \cdot \cdot \cdot$ ，如果马氏链在第 $n$ 步已经收敛了，于是我们可以得到 $p\left( x \right)$ 的样本 ${x_n},{x_{n + 1}}, \cdot \cdot \cdot$ ，所以关键问题是如何构造转移矩阵 $P$，我们是基于下面的定理。

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

$\begin{align*} \pi \left( i \right){P_{ij}} = \pi \left( j \right){P_{ji}} \tag{34} \end{align*}$

则 $\pi \left( x \right)$ 是马氏链的平稳分布。
证明很简单，有公式(34)得：

$\begin{align*} \sum\limits_{i = 1}^\infty {\pi \left( i \right){P_{ij}}} = \sum\limits_{i = 1}^\infty {\pi \left( j \right){P_{ji}}} = \pi \left( j \right)\sum\limits_{i = 1}^\infty {{P_{ji}}} = \pi \left( j \right) \tag{35} \end{align*}$

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

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

$\begin{align*} \pi \left( i \right)Q\left( {i,j} \right) \ne \pi \left( j \right)Q\left( {j,i} \right) \tag{36} \end{align*}$

我们对公式(36)进行改造，使细致平稳条件成立，引入 $\alpha \left( {i,j} \right)$ 。

$\begin{align*} \pi \left( i \right)Q\left( {i,j} \right)\alpha \left( {i,j} \right) = \pi \left( j \right)Q\left( {j,i} \right)\alpha \left( {j,i} \right) \tag{37} \end{align*}$

$\alpha \left( {i,j} \right)$ 如何取值才能使公式(37)成立？最简单的我们可以取：

$\begin{align*} \alpha \left( {i,j} \right) &= \pi \left( j \right)Q\left( {j,i} \right)\\ \alpha \left( {j,i} \right) &= \pi \left( i \right)Q\left( {i,j} \right) \tag{38} \end{align*}$

$Q'\left( {i,j} \right) = Q\left( {i,j} \right)\alpha \left( {i,j} \right)$，$Q'\left( {j,i} \right) = Q\left( {j,i} \right)\alpha \left( {j,i} \right)$ 所以我们有:

$\begin{align*} \pi \left( i \right)Q'\left( {i,j} \right) = \pi \left( j \right)Q'\left( {j,i} \right) \tag{39} \end{align*}$

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

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

总结下MCMC的采样过程。

Algorithm 1.2 MCMC采样算法

1. 初始化马氏链初始状态 ${X_0} = {x_0}$
2. 对 $t = 0,1,2, \cdot \cdot \cdot$ 循环以下过程进行采样
   
   • 第 $t$ 时刻马氏链状态为 ${X_t} = {x_t}$ ，从条件概率 $Q\left( {x|{x_t}} \right)$ 采样得到样本 $x'$
   • 从均匀分布采样 $u \sim uniform\left( {0,1} \right)$
   • 如果 $u < \alpha \left( {{x_t},x'} \right) = \pi \left( {x'} \right)Q\left( {{x_t}|x'} \right)$则接受转移 ${x_t} \to x'$ ，即 ${X_{t + 1}} = x'$
   • 否则不接受转移，即 ${X_{t + 1}} = {x_t}$

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

3. M-H采样

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

$\begin{align*} \pi \left( i \right)Q\left( {i,j} \right) \times 0.1 = \pi \left( j \right)Q\left( {j,i} \right) \times 0.2 \tag{40} \end{align*}$

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

$\begin{align*} \pi \left( i \right)Q\left( {i,j} \right) \times 0.5 = \pi \left( j \right)Q\left( {j,i} \right) \times 1 \tag{41} \end{align*}$

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

$\begin{align*} \alpha \left( {i,j} \right) = \min \left\{ {\frac{{\pi \left( j \right)Q\left( {j,i} \right)}}{{\pi \left( i \right)Q\left( {i,j} \right)}},1} \right\} \tag{42} \end{align*}$

Algorithm 1.3 M-H采样算法

1. 初始化马氏链初始状态 ${X_0} = {x_0}$
2. 对 $t = 0,1,2, \cdot \cdot \cdot$ 循环以下过程进行采样
   
   • 第 [t] 时刻马氏链状态为 ${X_t} = {x_t}$ ，从条件概率 $Q\left( {x|{x_t}} \right)$ 采样得到样本 $x'$
   • 从均匀分布采样 $u \sim uniform\left( {0,1} \right)$
   • 如果 $u < \alpha \left( {{x_t},x'} \right) = \min \left\{ {\frac{{\pi \left( {x'} \right)Q\left( {{x_t}|x'} \right)}}{{\pi \left( {{x_t}} \right)Q\left( {x'|{x_t}} \right)}},1} \right\}$ 则接受转移 ${x_t} \to x'$，即 ${X_{t + 1}} = x'$
   • 否则不接受转移，即 ${X_{t + 1}} = {x_t}$

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

• 当 $i \ne j$，$Q'\left( {i,j} \right) = Q\left( {i,j} \right)\alpha \left( {i,j} \right)$；
• 当 $i = j$，考虑拒绝转移概率，$Q'\left( {i,i} \right) = Q\left( {i,i} \right)\alpha \left( {i,i} \right) + \sum\limits_j {Q\left( {i,j} \right)\left( {1 - \alpha \left( {i,j} \right)} \right)}$。
  $\begin{align*} \sum\limits_j {Q'\left( {i,j} \right)} &= Q'\left( {i,i} \right) + \sum\limits_{i \ne j} {Q'\left( {i,j} \right)} \\ &= Q\left( {i,i} \right)\alpha \left( {i,i} \right) + \sum\limits_j {Q\left( {i,j} \right)\left( {1 - \alpha \left( {i,j} \right)} \right)} + \sum\limits_{i \ne j} {Q\left( {i,j} \right)\alpha \left( {i,j} \right)} \\ &= \sum\limits_j {Q\left( {i,j} \right)} = 1 \end{align*}$

1.10 Gibbs Sampling

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

$\begin{align*} p\left( {{x_1},{y_1}} \right)p\left( {{y_2}|{x_1}} \right) = p\left( {{x_1}} \right)p\left( {{y_1}|{x_1}} \right)p\left( {{y_2}|{x_1}} \right)\\ p\left( {{x_1},{y_2}} \right)p\left( {{y_1}|{x_1}} \right) = p\left( {{x_1}} \right)p\left( {{y_2}|{x_1}} \right)p\left( {{y_1}|{x_1}} \right) \tag{43} \end{align*}$

所以可得：

$\begin{align*} p\left( {{x_1},{y_1}} \right)p\left( {{y_2}|{x_1}} \right){\rm{ = }}p\left( {{x_1},{y_2}} \right)p\left( {{y_1}|{x_1}} \right) \tag{44} \end{align*}$

也就是：

$\begin{align*} p\left( A \right)p\left( {{y_2}|{x_1}} \right){\rm{ = }}p\left( B \right)p\left( {{y_1}|{x_1}} \right) \tag{45} \end{align*}$

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

$\begin{align*} p\left( A \right)p\left( {{x_2}|{y_1}} \right){\rm{ = }}p\left( C \right)p\left( {{x_1}|{y_1}} \right) \tag{46} \end{align*}$

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

$\begin{align*} Q\left( {A \to B} \right) &= p\left( {{y_B}|{x_1}} \right) &if\left( {{x_A} = {x_B} = {x_1}} \right) \\ Q\left( {A \to C} \right) &= p\left( {{x_C}|{y_1}} \right) &if\left( {{y_A} = {y_C} = {y_1}} \right)\\ Q\left( {A \to D} \right) &= 0 &others \tag{47} \end{align*}$

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

$\begin{align*} P\left( X \right)Q\left( {X \to Y} \right) = P\left( Y \right)Q\left( {Y \to X} \right) \tag{48} \end{align*}$

所以这个二维空间上的马氏链将收敛到平稳分布 $p\left( {x,y} \right)$，称为Gibbs Sampling算法。

Algorithm 1.4 Gibbs Sampling算法

1. 随机初始化 ${X_0} = {x_0},{Y_0} = {y_0}$
2. 对 $t = 0,1,2, \cdot \cdot \cdot$ 循环以下过程进行采样
   
   • 从条件概率分布 $p\left( {y|{x_t}} \right)$ 中采样得到 ${y_{t + 1}} \sim p\left( {y|{x_t}} \right)$
   • 从条件概率分布 $p\left( {x|{y_{t + 1}}} \right)$ 中采样得到 ${x_{t + 1}} \sim p\left( {x|{y_{t + 1}}} \right)$

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

二维可以很容易推广到高维的情况，在 $n$ 维空间中对于概率分布 $p\left( {{x_1},{x_2}, \cdot \cdot \cdot ,{x_n}} \right)$。

Algorithm 1.5 n维Gibbs Sampling算法

1. 随机初始化 $\left\{ {{x_i}:i = 1, \cdot \cdot \cdot ,n} \right\}$
2. 对 $t = 0,1,2, \cdot \cdot \cdot$ 循环以下过程进行采样
   
   • $x_1^{\left( {t + 1} \right)} \sim p\left( {{x_1}|x_2^{\left( t \right)},x_3^{\left( t \right)}, \cdot \cdot \cdot ,x_n^{\left( t \right)}} \right)$
   • $x_2^{\left( {t + 1} \right)} \sim p\left( {{x_2}|x_1^{\left( {t{\rm{ + }}1} \right)},x_3^{\left( t \right)}, \cdot \cdot \cdot ,x_n^{\left( t \right)}} \right)$
   • $\cdot \cdot \cdot$
   • $x_j^{\left( {t + 1} \right)} \sim p\left( {{x_j}|x_1^{\left( {t{\rm{ + }}1} \right)}, \cdot \cdot \cdot ,x_{j - 1}^{\left( {t + 1} \right)},x_j^{\left( t \right)} \cdot \cdot \cdot ,x_n^{\left( t \right)}} \right)$
   • $\cdot \cdot \cdot$
   • $x_n^{\left( {t + 1} \right)} \sim p\left( {{x_n}|x_1^{\left( {t + 1} \right)}, \cdot \cdot \cdot ,x_{n - 1}^{\left( {t + 1} \right)}} \right)$

1.11 EM算法

我们先介绍凸函数的概念，$f$ 的定义域是实数集，若 $x \in R$ 且 $f''\left( x \right) \ge 0$ ，则 $f$ 是凸函数，若 $f''\left( x \right) > 0$ ，则 $f$ 是严格凸函数；若${\rm{x}}$ 是向量且hessian矩阵 $H$ 是半正定矩阵，则 $f$ 是凸函数，若 $H$ 是正定矩阵，则 $f$ 是严格凸函数。
定理1.3（Jensen不等式）$f$ 的定义域是实数集，且是凸函数，则有：

$\begin{align*} E\left[ {f\left( X \right)} \right] \ge f\left( {E\left[ X \right]} \right) \tag{49} \end{align*}$

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

假设训练集 $\left\{ {{x^{\left( 1 \right)}},{x^{\left( 2 \right)}}, \cdot \cdot \cdot ,{x^{\left( m \right)}}} \right\}$，每个样本相互独立，我们需要估计模型 $p\left( {x,z} \right)$ 的参数 $\theta$，由于含有隐变量 $z$，所以很难直接用最大似然求解，如果 $z$ 已知，那么就可以用最大似然求解。

$\begin{align*} l\left( \theta \right) &= \sum\limits_{i = 1}^m {\log p\left( {x;\theta } \right)} \\ &= \sum\limits_{i = 1}^m {\log \sum\limits_z {p\left( {x,z;\theta } \right)} } \tag{50} \end{align*}$

其实我们的目标是找到 $z$ 和 $\theta$ 使 $l\left( \theta \right)$ 最大，也就是分别对 $z$ 和 $\theta$ 求偏导，然后令其为0，理想是美好的，现实是残酷的，公式(49)求偏导后变的很复杂，求导前要是能把求和符号从对数函数中提出来就好了。EM算法可以有效地解决这个问题，引入 ${Q_i}\left( {{z^{\left( i \right)}}} \right)$ 表示 ${z^{\left( i \right)}}$ 的概率分布（$\sum\nolimits_z {{Q_i}\left( {{z^{\left( i \right)}}} \right)} = 1,{Q_i}\left( {{z^{\left( i \right)}}} \right) \ge 0$）。由公式(50)可得：

$\begin{align*} \sum\limits_i {\log p\left( {{x^{\left( i \right)}};\theta } \right)} &= \sum\limits_i {\log \sum\limits_{{z^{\left( i \right)}}} {p\left( {{x^{\left( i \right)}},{z^{\left( i \right)}};\theta } \right)} } \\ &= \sum\limits_i {\log \sum\limits_{{z^{\left( i \right)}}} {{Q_i}\left( {{z^{\left( i \right)}}} \right)\frac{{p\left( {{x^{\left( i \right)}},{z^{\left( i \right)}};\theta } \right)}}{{{Q_i}\left( {{z^{\left( i \right)}}} \right)}}} } \\ & \ge \sum\limits_i {\sum\limits_{{z^{\left( i \right)}}} {{Q_i}\left( {{z^{\left( i \right)}}} \right)\log \frac{{p\left( {{x^{\left( i \right)}},{z^{\left( i \right)}};\theta } \right)}}{{{Q_i}\left( {{z^{\left( i \right)}}} \right)}}} } \tag{51} \end{align*}$

最后一步是利用Jensen不等式，$f = \log \left( x \right)$，$f'' = - \frac{1}{{{x^2}}} < 0$，所以 $f$ 是凹函数，$\sum\limits_{{z^{\left( i \right)}}} {{Q_i}\left( {{z^{\left( i \right)}}} \right)\frac{{p\left( {{x^{\left( i \right)}},{z^{\left( i \right)}};\theta } \right)}}{{{Q_i}\left( {{z^{\left( i \right)}}} \right)}}}$是 $\frac{{p\left( {{x^{\left( i \right)}},{z^{\left( i \right)}};\theta } \right)}}{{{Q_i}\left( {{z^{\left( i \right)}}} \right)}}$ 的期望，所以有：

$\begin{align*} f\left( {E\left[ {\frac{{p\left( {{x^{\left( i \right)}},{z^{\left( i \right)}};\theta } \right)}}{{{Q_i}\left( {{z^{\left( i \right)}}} \right)}}} \right]} \right) \ge E\left[ {f\left( {\frac{{p\left( {{x^{\left( i \right)}},{z^{\left( i \right)}};\theta } \right)}}{{{Q_i}\left( {{z^{\left( i \right)}}} \right)}}} \right)} \right] \tag{52} \end{align*}$

由公式(51)可知，我们可以不断地最大化下界，以提高 $l\left( \theta \right)$，最终达到最大值。如果固定 $\theta$，那么 $l\left( \theta \right)$ 的下界就取决于 ${Q_i}\left( {{z^{\left( i \right)}}} \right)$，可以通过调整这个概率，使得下界不断地上升逼近 $l\left( \theta \right)$，最终相等，然后固定 ${Q_i}\left( {{z^{\left( i \right)}}} \right)$，调整 $\theta$，使下界达到最大值，此时 $\theta$ 为新的值，再固定 $\theta$，调整 ${Q_i}\left( {{z^{\left( i \right)}}} \right)$，反复直到收敛到 $l\left( \theta \right)$ 的最大值。现在我们有两个问题需要证明，1. 下界何时等于 $l\left( \theta \right)$；2. 为什么可以收敛到最大值。

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

$\begin{align*} \frac{{p\left( {{x^{\left( i \right)}},{z^{\left( i \right)}};\theta } \right)}}{{{Q_i}\left( {{z^{\left( i \right)}}} \right)}} = c \tag{53} \end{align*}$

再由 $\sum\nolimits_z {{Q_i}\left( {{z^{\left( i \right)}}} \right)} = 1$ 得：

$\begin{align*} \sum\nolimits_z {\frac{{p\left( {{x^{\left( i \right)}},{z^{\left( i \right)}};\theta } \right)}}{c}} &= 1 \Rightarrow c = \sum\nolimits_z {p\left( {{x^{\left( i \right)}},{z^{\left( i \right)}};\theta } \right)} \\ {Q_i}\left( {{z^{\left( i \right)}}} \right) &= \frac{{p\left( {{x^{\left( i \right)}},{z^{\left( i \right)}};\theta } \right)}}{{\sum\nolimits_z {p\left( {{x^{\left( i \right)}},{z^{\left( i \right)}};\theta } \right)} }}{\rm{ = }}\frac{{p\left( {{x^{\left( i \right)}},{z^{\left( i \right)}};\theta } \right)}}{{p\left( {{x^{\left( i \right)}};\theta } \right)}}{\rm{ = }}p\left( {{z^{\left( i \right)}}|{x^{\left( i \right)}},\theta } \right) \tag{54} \end{align*}$

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

Algorithm 1.6 EM算法

1. 初始化$\theta$
2. Repeat until convergence
   
   • (E-step) For each i,set
     ${Q_i}\left( {{z^{\left( i \right)}}} \right){\rm{ = }}p\left( {{z^{\left( i \right)}}|{x^{\left( i \right)}},\theta } \right)$
   • (M-step) Set
     $\theta : = \mathop {\arg \max }\limits_\theta \sum\limits_i {\sum\limits_{{z^{\left( i \right)}}} {{Q_i}\left( {{z^{\left( i \right)}}} \right)\log \frac{{p\left( {{x^{\left( i \right)}},{z^{\left( i \right)}};\theta } \right)}}{{{Q_i}\left( {{z^{\left( i \right)}}} \right)}}} }$

现在我们开始讨论第二个问题，${\theta ^{\left( t \right)}}$ 和 ${\theta ^{\left( {t{\rm{ + }}1} \right)}}$ 是EM迭代过程的参数估计，我们需要证明$l\left( {{\theta ^{\left( t \right)}}} \right) \le l\left( {{\theta ^{\left( {t + 1} \right)}}} \right)$ ，也就是EM算法是单调地提高 $l\left( \theta \right)$ ，$Q_i^{\left( t \right)}\left( {{z^{\left( i \right)}}} \right){\rm{ = }}p\left( {{z^{\left( i \right)}}|{x^{\left( i \right)}},{\theta ^{\left( t \right)}}} \right)$。

$\begin{align*} l\left( {{\theta ^{\left( t \right)}}} \right) = \sum\limits_i {\sum\limits_{{z^{\left( i \right)}}} {Q_i^{\left( t \right)}\left( {{z^{\left( i \right)}}} \right)\log \frac{{p\left( {{x^{\left( i \right)}},{z^{\left( i \right)}};{\theta ^{\left( t \right)}}} \right)}}{{Q_i^{\left( t \right)}\left( {{z^{\left( i \right)}}} \right)}}} } \tag{55} \end{align*}$

$\begin{align*} l\left( {{\theta ^{\left( {t{\rm{ + }}1} \right)}}} \right) &\ge \sum\limits_i {\sum\limits_{{z^{\left( i \right)}}} {Q_i^{\left( t \right)}\left( {{z^{\left( i \right)}}} \right)\log \frac{{p\left( {{x^{\left( i \right)}},{z^{\left( i \right)}};{\theta ^{\left( {t{\rm{ + }}1} \right)}}} \right)}}{{Q_i^{\left( t \right)}\left( {{z^{\left( i \right)}}} \right)}}} } \\ &\ge \sum\limits_i {\sum\limits_{{z^{\left( i \right)}}} {Q_i^{\left( t \right)}\left( {{z^{\left( i \right)}}} \right)\log \frac{{p\left( {{x^{\left( i \right)}},{z^{\left( i \right)}};{\theta ^{\left( t \right)}}} \right)}}{{Q_i^{\left( t \right)}\left( {{z^{\left( i \right)}}} \right)}}} } \\ &= l\left( {{\theta ^{\left( t \right)}}} \right) \tag{56} \end{align*}$

第一个不等式是因为：

$\begin{align*} l\left( \theta \right) \ge \sum\limits_i {\sum\limits_{{z^{\left( i \right)}}} {{Q_i}\left( {{z^{\left( i \right)}}} \right)\log \frac{{p\left( {{x^{\left( i \right)}},{z^{\left( i \right)}};\theta } \right)}}{{{Q_i}\left( {{z^{\left( i \right)}}} \right)}}} } \tag{57} \end{align*}$

公式(57)中，$\theta = {\theta ^{\left( {t + 1} \right)}}$，${Q_i} = Q_i^{\left( t \right)}$。

第二个不等式是因为${\theta ^{\left( {t + 1} \right)}}$是为了

$\begin{align*} \mathop {\arg \max }\limits_\theta \sum\limits_i {\sum\limits_{{z^{\left( i \right)}}} {{Q_i}\left( {{z^{\left( i \right)}}} \right)\log \frac{{p\left( {{x^{\left( i \right)}},{z^{\left( i \right)}};\theta } \right)}}{{{Q_i}\left( {{z^{\left( i \right)}}} \right)}}} } \tag{58} \end{align*}$

三、LDA

2.1 Unigram Model

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

Game 2.1 Unigram Model

1. 上帝只有一个骰子，这个骰子有 $V$ 个面，每个面对应一个词，各个面概率不同；
2. 每抛一次骰子，抛出的面就对应的产生一个词，如果一篇文档中有 $n$ 个词，那么就独立地抛 $n$ 次骰子产生这 $n$ 个词

骰子各个面的概率记为 $\vec p = \left( {{p_1},{p_2}, \cdot \cdot \cdot ,{p_V}} \right)$，对于一篇文档 $\vec w = \left( {{w_1},{w_2}, \cdot \cdot \cdot ,{w_n}} \right)$，生成该文档的概率为：

$\begin{align*} p\left( {{w_1},{w_2}, \cdot \cdot \cdot ,{w_n}} \right) = \prod\limits_{i = 1}^n {p\left( {{w_i}} \right)} \tag{59} \end{align*}$

假设我们预料是由 $m$ 篇文档组成即 $W = \left( {{{\vec w}_1},{{\vec w}_2}, \cdot \cdot \cdot ,{{\vec w}_m}} \right)$，每篇文档是相互独立的，则该预料的概率为：

$\begin{align*} p\left( {W|\vec p} \right) = \prod\limits_{i = 1}^n {p\left( {{{\vec w}_i}|\vec p} \right)} \tag{60} \end{align*}$

假设预料中总共有 $N$ 个词，每个词 ${w_i}$ 的词频为 ${n_i}$ ，那么$\vec n = \left( {{n_1},{n_2}, \cdot \cdot \cdot ,{n_V}} \right)$ 服从多项式分布，可参考1.5节的多项式分布概念。

$\begin{align*} p\left( {\vec n|\vec p} \right) = Mult\left( {\vec n|\vec p,N} \right) = \left( {\begin{array}{*{20}{c}} N\\ {\vec n} \end{array}} \right)\prod\limits_{k = 1}^V {p_k^{{n_k}}} \tag{61} \end{align*}$

此时公式（60）为：

$\begin{align*} p\left( {W|\vec p} \right) = \prod\limits_{i = 1}^n {p\left( {{{\vec w}_i}|\vec p} \right)} = \prod\limits_{k = 1}^V {p_k^{{n_k}}} \tag{62} \end{align*}$

我们需要估计模型中的参数 $\vec p$，可以用最大似然估计：

$\begin{align*} {{\hat p}_{ML}} = \mathop {\arg \max }\limits_{\vec p} \sum\limits_{k = 1}^V {{n_k}\ln \left( {{p_k}} \right)} \\ s.t.{\rm{ }}\sum\limits_{k = 1}^V {{p_k}} = 1,\sum\limits_{k = 1}^V {{n_k}} = N \tag{63} \end{align*}$

于是参数 ${p_k}$ 的估计值就是：

$\begin{align*} L\left( {\vec p} \right) &= \sum\limits_{k = 1}^V {{n_k}\ln \left( {{p_k}} \right)} + \lambda \left( {\sum\limits_{k = 1}^V {{p_k}} - 1} \right)\\ \frac{{\partial L\left( {\vec p} \right)}}{{\partial {p_k}}} &= \frac{{{n_k}}}{{{p_k}}} + \lambda = 0 \Rightarrow {p_k} = - \frac{{{n_k}}}{\lambda }\\ \sum\limits_{k = 1}^V {{n_k}} &= N \Rightarrow \lambda = - N\\ {p_k} &= \frac{{{n_k}}}{N} \tag{64} \end{align*}$

2.2 贝叶斯Unigram Model

对于以上模型，统计学家中贝叶斯学派就不同意了，为什么上帝只有一个固定的筛子呢，在贝叶斯学派看来，一切参数都是随机变量，模型中 $\vec p$ 不是唯一固定的，而是服从一个分布，所以贝叶斯Unigram Model游戏规则变为：

Game 2.2 贝叶斯Unigram Model

1. 上帝有一个装有无穷多个骰子的坛子，每个骰子有 $V$ 个面，每个面对应一个词，各个面概率不同；
2. 上帝从坛子中抽一个骰子出来，然后用这个骰子不断地抛，产生预料中的所有词。

上帝这个坛子里面有些骰子数量多，有些骰子数量少，所以从概率分布的角度看，坛子里面的骰子 $\vec p$ 服从一个概率分布 $p\left( {\vec p} \right)$，这个分布称为参数 $\vec p$ 的先验分布。先验分布 $p\left( {\vec p} \right)$ 可以有多种选择，注意到 $\vec n$ 是服从多项式分布的，$p\left( {\vec n|\vec p} \right) = Mult\left( {\vec n|\vec p,N} \right)$，回顾1.7节可知，$p\left( {\vec p} \right)$ 最好的选择是Dirichlet分布：

$\begin{align*} p\left( {\vec p|\vec \alpha } \right) = Dir\left( {\vec p|\vec \alpha } \right) = \frac{1}{{\Delta \left( {\vec \alpha } \right)}}\prod\limits_{k = 1}^V {p_k^{{\alpha _k} - 1}} ,\vec \alpha = \left( {{\alpha _1},{\alpha _2}, \cdot \cdot \cdot ,{\alpha _V}} \right) \tag{65} \end{align*}$

于是，在给定了参数 $\vec p$ 的先验分布 $Dir\left( {\vec p|\vec \alpha } \right)$ 时候，语料中各个词出现的次数服从多项式分布 $\vec n \sim Mult\left( {\vec n|\vec p,N} \right)$，所以后验分布为：

$\begin{align*} p\left( {\vec p|\vec n,\vec \alpha } \right) = \frac{{p\left( {\vec n|\vec p} \right)p\left( {\vec p|\vec \alpha } \right)}}{{\int {p\left( {\vec n|\vec p} \right)p\left( {\vec p|\vec \alpha } \right)d\vec p} }} = \frac{1}{{\Delta \left( {\vec \alpha + \vec n} \right)}}\prod\limits_{k = 1}^V {p_k^{{\alpha _k} + {n_k} - 1}} \tag{66} \end{align*}$

对参数 $\vec p$ 采用贝叶斯估计，假设参数 $\vec p$ 服从 $Dir\left( {\vec p|\vec \alpha } \right)$ 分布，我们利用样本信息对 $\vec p$ 的先验分布进行修正，得到 $\vec p$ 的后验分布也是服从 $Dir\left( {\vec p|\vec \alpha + \vec n} \right)$ 分布。可以用 $\vec p$ 的期望值作为参数 $\vec p$ 的估计值。由1.6节可知，$\vec p$ 的期望值为：

$\begin{align*} E\left( {\vec p} \right) = \left( {\frac{{{n_1} + {\alpha _1}}}{{\sum\nolimits_{k = 1}^V {\left( {{n_k} + {\alpha _k}} \right)} }},\frac{{{n_2} + {\alpha _2}}}{{\sum\nolimits_{k = 1}^V {\left( {{n_k} + {\alpha _k}} \right)} }}, \cdot \cdot \cdot ,\frac{{{n_V} + {\alpha _V}}}{{\sum\nolimits_{k = 1}^V {\left( {{n_k} + {\alpha _k}} \right)} }}} \right) \tag{67} \end{align*}$

接下来我们计算语料产生的概率，开始并不知道上帝到底用哪个骰子，所以每个骰子都有可能被使用，使用的概率由 $p\left( {\vec p|\vec \alpha } \right)$ 决定的，对于每个具体的骰子，由该骰子产生预料的概率为 $p\left( {W|\vec p} \right)$，所以语料产生的概率为：

$\begin{align*} p\left( {W|\vec \alpha } \right) &= \int {p\left( {W|\vec p} \right)p\left( {\vec p|\vec \alpha } \right)d\vec p} \\ &= \int {\prod\limits_{k = 1}^V {p_k^{{n_k}}Dir\left( {\vec p|\vec \alpha } \right)} d\vec p} \\ &= \int {\prod\limits_{k = 1}^V {p_k^{{n_k}}\frac{1}{{\Delta \left( {\vec \alpha } \right)}}\prod\limits_{k = 1}^V {p_k^{{\alpha _k} - 1}} } d\vec p} \\ &= \frac{1}{{\Delta \left( {\vec \alpha } \right)}}\int {\prod\limits_{k = 1}^V {p_k^{{n_k} + {\alpha _k} - 1}} d\vec p} \\ &= \frac{{\Delta \left( {\vec \alpha + \vec n} \right)}}{{\Delta \left( {\vec \alpha } \right)}} \tag{68} \end{align*}$

2.3 PLSA

1. PLSA Model

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

Game 2.3 PLSA Topic Model

1. 上帝有两种类型的骰子，一类是doc-topic，每个骰子有 $K$ 个面，每个面是一个topic编号；一类是topic-word，每个骰子有 $V$ 个面，每个面对应一个词。
2. 上帝一共有 $K$ 个topic-word骰子，每个骰子有一个编号，从1到 $K$。
3. 生成每篇文档之前，上帝都先为这篇文档制造一个特定的doc-topic骰子，然后重复下面过程生成文档的词：
   
   • 投掷这个doc-topic骰子，得到一个编号为 $z$ 的topic
   • 选择 $K$ 个topic-word骰子中编号为 $z$ 的骰子，投掷这个骰子，得到一个词。

2. EM算法推导PLSA

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

$D = \left\{ {{d_1},{d_2}, \cdot \cdot \cdot ,{d_N}} \right\}$：表示语料中 $N$ 篇文档；$W = \left\{ {{w_1},{w_2}, \cdot \cdot \cdot ,{w_M}} \right\}$：表示语料中 $M$ 个词组；$n\left( {{d_i},{w_j}} \right)$：表示词 ${w_j}$ 在文档 ${d_i}$ 中出现的频次，${\rm{N}} = {\left( {n\left( {{d_i},{w_j}} \right)} \right)_{ij}} \in {R^{N \times M}}$；$Z = \left\{ {{z_1},{z_2}, \cdot \cdot \cdot ,{z_K}} \right\}$：表示 $K$ 个主题，每篇文档可以有多个主题；$p\left( {{w_j}|d{}_i} \right)$：表示词 ${w_j}$ 在给定文档 ${d_i}$ 中出现的概率；$p\left( {{z_k}|d{}_i} \right)$：表示主题 ${z_k}$ 在给定文档 ${d_i}$ 下出现的概率；$p\left( {{w_j}|{z_k}} \right)$：表示词 ${w_j}$ 在给定主题 ${z_k}$ 下出现的概率。

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

• select a document ${d_i}$ with probability $p\left( {{d_i}} \right)$
• pick a latent class ${z_k}$ with probability $p\left( {{z_k}|{d_i}} \right)$
• generate a word ${w_j}$ with probability $p\left( {{w_j}|{z_k}} \right)$

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

$\begin{align*} L &= \prod\limits_{i = 1}^N {\prod\limits_{j = 1}^M {p{{\left( {{d_i},{w_j}} \right)}^{n\left( {{d_i},{w_j}} \right)}}} } \\ \log L &= \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {n\left( {{d_i},{w_j}} \right)\log p\left( {{d_i},{w_j}} \right)} } \\ &= \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {n\left( {{d_i},{w_j}} \right)\log \left[ {\sum\limits_{k = 1}^K {p\left( {{d_i}} \right)p\left( {{z_k}|{d_i}} \right)p\left( {{w_j}|{z_k}} \right)} } \right]} } \\ &= \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {n\left( {{d_i},{w_j}} \right)\left[ {\log p\left( {{d_i}} \right) + \log \sum\limits_{k = 1}^K {p\left( {{z_k}|{d_i}} \right)p\left( {{w_j}|{z_k}} \right)} } \right]} } \\ &= \sum\limits_{i = 1}^N {n\left( {{d_i}} \right)\log p\left( {{d_i}} \right)} + \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {n\left( {{d_i},{w_j}} \right)\log \sum\limits_{k = 1}^K {p\left( {{z_k}|{d_i}} \right)p\left( {{w_j}|{z_k}} \right)} } } \tag{69} \end{align*}$

对于给定训练预料，希望公式 (69) 最大化。$p\left( {{z_k}|{d_i}} \right)$ 和 $p\left( {{w_j}|{z_k}} \right)$ 是 PLSA 模型需要求解的参数，按照通常的做法是令偏导数 为0，但是参数是以求和的形式出现在对数函数里面，求导后会变得很复杂。$n\left( {{d_i}} \right)$ 表示第 $i$ 篇文档的词数，所以当预料固定，公式(69)中第一项可以看作常量，所以只要最大化(69)中的第二项即可，如公式(70)所示。

$\begin{align*} \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {n\left( {{d_i},{w_j}} \right)\log \sum\limits_{k = 1}^K {p\left( {{z_k}|{d_i}} \right)p\left( {{w_j}|{z_k}} \right)} } } \tag{70} \end{align*}$

引入 ${Q_k}\left( {{z_k}} \right)$ 表示 ${z_k}$ 的概率分布（$\sum\limits_{k = 1}^K {{Q_k}\left( {{z_k}} \right)} = 1,{Q_k}\left( {{z_k}} \right) \ge 0$），根据Jensen不等式得：

$\begin{align*} \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {n\left( {{d_i},{w_j}} \right)\log \sum\limits_{k = 1}^K {{Q_k}\left( {{z_k}} \right)\frac{{p\left( {{z_k}|{d_i}} \right)p\left( {{w_j}|{z_k}} \right)}}{{{Q_k}\left( {{z_k}} \right)}}} } } \\ {\rm{ }} \ge \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {n\left( {{d_i},{w_j}} \right)\sum\limits_{k = 1}^K {{Q_k}\left( {{z_k}} \right)\log \frac{{p\left( {{z_k}|{d_i}} \right)p\left( {{w_j}|{z_k}} \right)}}{{{Q_k}\left( {{z_k}} \right)}}} } } \tag{71} \end{align*}$

当${Q_k}\left( {{z_k}} \right) = \frac{{p\left( {{z_k}|{d_i}} \right)p\left( {{w_j}|{z_k}} \right)}}{{\sum\limits_{l = 1}^K {p\left( {{z_l}|{d_i}} \right)p\left( {{w_j}|{z_l}} \right)} }} = \frac{{p\left( {{w_j},{z_k}|{d_i}} \right)}}{{\sum\limits_{l = 1}^K {p\left( {{w_j},{z_l}|{d_i}} \right)} }} = p\left( {{z_k}|{w_j},{d_i}} \right)$时，公式(71)不等式中等号成立，所以只需要最大化：

$\begin{align*} & \mathop {\min }\limits_{p\left( {{z_k}|{d_i}} \right),p\left( {{w_j}|{z_k}} \right)} {\rm{ }}\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {n\left( {{d_i},{w_j}} \right)\sum\limits_{k = 1}^K {{Q_k}\left( {{z_k}} \right)\log \frac{{p\left( {{z_k}|{d_i}} \right)p\left( {{w_j}|{z_k}} \right)}}{{{Q_k}\left( {{z_k}} \right)}}} } } \\ & s.t.{\rm{ }}\sum\limits_{j = 1}^M {p\left( {{w_j}|{z_k}} \right)} = 1,{\rm{ }}\sum\limits_{k = 1}^K {p\left( {{z_k}|{d_i}} \right)} = 1 \tag{72} \end{align*}$

根据拉格朗日乘子法

$\begin{align*} &l = \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {n\left( {{d_i},{w_j}} \right)\sum\limits_{k = 1}^K {{Q_k}\left( {{z_k}} \right)\log \frac{{p\left( {{z_k}|{d_i}} \right)p\left( {{w_j}|{z_k}} \right)}}{{{Q_k}\left( {{z_k}} \right)}}} } } \\&+ {\lambda _1}\left( {\sum\limits_{j = 1}^M {p\left( {{w_j}|{z_k}} \right)} - 1} \right) + {\lambda _2}\left( {\sum\limits_{k = 1}^K {p\left( {{z_j}|{d_i}} \right)} - 1} \right)\\ &\frac{{\partial l}}{{\partial p\left( {{w_j}|{z_k}} \right)}} = \sum\limits_{i = 1}^N {\frac{{n\left( {{d_i},{w_j}} \right){Q_k}\left( {{z_k}} \right)}}{{p\left( {{w_j}|{z_k}} \right)}}} + {\lambda _1},\\&{\rm{ }}\frac{{\partial l}}{{\partial p\left( {{z_k}|{d_i}} \right)}} = \sum\limits_{j = 1}^M {\frac{{n\left( {{d_i},{w_j}} \right){Q_k}\left( {{z_k}} \right)}}{{p\left( {{z_k}|{d_i}} \right)}}} + {\lambda _2} \end{align*}$

所以可得：

$\begin{align*} p\left( {{w_j}|{z_k}} \right) &= - \frac{{\sum\limits_{i = 1}^N {n\left( {{d_i},{w_j}} \right){Q_k}\left( {{z_k}} \right)} }}{{{\lambda _1}}},{\rm{ }}\sum\limits_{j = 1}^M {p\left( {{w_j}|{z_k}} \right)} = 1\\ & \Rightarrow {\lambda _1} = - \sum\limits_{j = 1}^M {\sum\limits_{i = 1}^N {n\left( {{d_i},{w_j}} \right){Q_k}\left( {{z_k}} \right)} } \\ & \Rightarrow p\left( {{w_j}|{z_k}} \right) = \frac{{\sum\limits_{i = 1}^N {n\left( {{d_i},{w_j}} \right){Q_k}\left( {{z_k}} \right)} }}{{\sum\limits_{j = 1}^M {\sum\limits_{i = 1}^N {n\left( {{d_i},{w_j}} \right){Q_k}\left( {{z_k}} \right)} } }}\\ p\left( {{z_k}|{d_i}} \right) &= - \frac{{\sum\limits_{j = 1}^M {n\left( {{d_i},{w_j}} \right){Q_k}\left( {{z_k}} \right)} }}{{{\lambda _2}}},{\rm{ }}\sum\limits_{k = 1}^K {p\left( {{z_k}|{d_i}} \right)} = 1\\ & \Rightarrow {\lambda _2} = - \sum\limits_{k = 1}^K {\sum\limits_{j = 1}^M {n\left( {{d_i},{w_j}} \right){Q_k}\left( {{z_k}} \right)} = } n\left( {{d_i}} \right)\\ & \Rightarrow p\left( {{z_k}|{d_i}} \right) = \frac{{\sum\limits_{j = 1}^M {n\left( {{d_i},{w_j}} \right){Q_k}\left( {{z_k}} \right)} }}{{n\left( {{d_i}} \right)}} \tag{73} \end{align*}$

总结EM算法为：

1. E-step 随机初始化变量 $p\left( {{z_k}|{d_i}} \right)$，$p\left( {{w_j}|{z_k}} \right)$，计算隐变量后验概率。

${Q_k}\left( {{z_k}} \right) = \frac{{p\left( {{z_k}|{d_i}} \right)p\left( {{w_j}|{z_k}} \right)}}{{\sum\limits_{l = 1}^K {p\left( {{z_l}|{d_i}} \right)p\left( {{w_j}|{z_l}} \right)} }} = \frac{{p\left( {{w_j},{z_k}|{d_i}} \right)}}{{\sum\limits_{l = 1}^K {p\left( {{w_j},{z_l}|{d_i}} \right)} }} = p\left( {{z_k}|{w_j},{d_i}} \right)$

1. M-step 最大化似然函数，更新变量$p\left( {{z_k}|{d_i}} \right)$，$p\left( {{w_j}|{z_k}} \right)$

$p\left( {{w_j}|{z_k}} \right) = \frac{{\sum\limits_{i = 1}^N {n\left( {{d_i},{w_j}} \right){Q_k}\left( {{z_k}} \right)} }}{{\sum\limits_{j = 1}^M {\sum\limits_{i = 1}^N {n\left( {{d_i},{w_j}} \right){Q_k}\left( {{z_k}} \right)} } }},{\rm{ }}p\left( {{z_k}|{d_i}} \right) = \frac{{\sum\limits_{j = 1}^M {n\left( {{d_i},{w_j}} \right){Q_k}\left( {{z_k}} \right)} }}{{n\left( {{d_i}} \right)}}$

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

2.4 LDA

对于 PLSA 模型，贝叶斯学派表示不同意，为什么上帝只有一个 doc-topic 骰子，为什么上帝只有固定 $K$ 个topic-word骰子？$p\left( {{z_k}|{d_i}} \right)$ 和 $p\left( {{w_j}|{z_k}} \right)$ 是模型的参数，一切参数都是随机变量，模型中 $p\left( {{z_k}|{d_i}} \right)$ 和 $p\left( {{w_j}|{z_k}} \right)$ 不是唯一固定的，类似 2.2 节贝叶斯 Unigram Model 和 2.1 节 Unigram Model 的关系。所以 LDA 游戏规则为：

Game 2.4 LDA Topic Model

1. 上帝有两坛的骰子，第一个坛子装的是 doc-topic 骰子，第二个坛子装的是 topic-word 骰子。
2. 上帝随机地从第二个坛子中独立的抽取 $K$ 个 topic-word 骰子，编号为1到 $K$。
3. 生成每篇文档之前，上帝先从第一个坛子中随机抽取一个 doc-topic 骰子，然后重复地投掷这个骰子，为文档中每个词生成一个 topic 编号为 $z$ （这里只是生成每个词的主题，词并未生成），重复生成文档中每个词对应的 topic。
4. 对语料中每篇文档的每个 topic 编号为 $z$ 的主题，选择 $K$ 个 topic-word 骰子中编号为 $z$ 骰子，投掷这个 topic-word 骰子，生成对应的 word。

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

• ${\vec \vartheta _m}$：$p\left( {z|{d_m}} \right)$，第 $m$ 篇文档的主题分布概率，$\vec \Theta = \left\{ {{{\vec \vartheta }_m}} \right\}_{m = 1}^M \in {R^{M \times K}}$;
• ${\vec \varphi _k}$：$p\left( {w|{z_k}} \right)$，主题为 $k$ 的词的概率分布，$\vec \Phi = \left\{ {{{\vec \varphi }_k}} \right\}_{k = 1}^K \in {R^{K \times V}}$；
• $n_m^{\left( k \right)}$：第 $m$ 篇文档中属于 topic $k$ 的词的个数，${\vec n_m} = \left( {n_m^{\left( 1 \right)},n_m^{\left( 2 \right)}, \cdot \cdot \cdot ,n_m^{\left( K \right)}} \right)$；
• $n_k^{\left( t \right)}$：topic $k$ 产生词 $t$ 的个数，${\vec n_k} = \left( {n_k^{\left( 1 \right)},n_k^{\left( 2 \right)}, \cdot \cdot \cdot ,n_k^{\left( V \right)}} \right)$；
• $\vec \alpha \in {R^{K \times 1}}$：${\vec \vartheta _m}$先验分布超参数；
• $\vec \beta \in {R^{V \times 1}}$：${\vec \varphi _k}$ 先验分布超参数；
• ${z_{m,n}}$：第 $m$ 篇文档中第 $n$ 个词的主题；
• ${w_{m,n}}$：第 $m$ 篇文档中第 $n$ 个词。

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

1. 联合概率分布

$\begin{align*}p\left( {\vec w,\vec z|\vec \alpha ,\vec \beta } \right) = p\left( {\vec w|\vec z,\vec \beta } \right)p\left( {\vec z|\vec \alpha } \right)\tag{74}\end{align*}$

1）$\vec \alpha \underbrace {\xrightarrow{{\begin{array}{*{20}{c}} {}&{} \end{array}}}}_{Dir}{\vec \vartheta _m}\underbrace {\xrightarrow{{\begin{array}{*{20}{c}} {}&{} \end{array}}}}_{Mult}{\vec z_m}$：第一步对 $Dir$ 分布进行采样得到样本 ${\vec \vartheta _m}$（也就是从第一个坛子中抽取 doc-topic 骰子 ${\vec \vartheta _m}$）；第二步 doc-topic 骰子有 $K$ 个面，每个面表示一个主题，那么在一次投掷骰子过程中，每个主题的概率为 ${\vec \vartheta _m}{\rm{ = }}\left( {\vartheta _m^{\left( 1 \right)},\vartheta _m^{\left( 2 \right)}, \cdot \cdot \cdot ,\vartheta _m^{\left( K \right)}} \right)$，第 $m$ 篇文档有 ${N_m}$ 个词，所以需要投掷 ${N_m}$ 次骰子，为该篇文档中的每个词生成一个主题， 第 $n$ 个词对应的主题为 ${z_{m,n}}$，整篇文档的主题表示为 ${\vec z_m}$。在 ${N_m}$ 次投掷过程中，每个主题出现的次数为${\vec n_m} = \left( {n_m^{\left( 1 \right)},n_m^{\left( 2 \right)}, \cdot \cdot \cdot ,n_m^{\left( K \right)}} \right)$，那么 ${\vec n_m}$ 服从多项式分布$Mult\left( {{{\vec n}_m}|{{\vec \vartheta }_m},{N_m}} \right)$（只生成每个词的主题，并未由主题产生具体的词）。可以采用贝叶斯估计对参数 ${\vec \vartheta _m}$ 进行估计。

• ${\vec \vartheta _m}$ 的先验分布为 $Dir\left( {{{\vec \vartheta }_m}|\vec \alpha } \right)$

• 后验分布为（推导过程可以参考1.7节）
  $\begin{align*}p\left( {{{\vec \vartheta }_m}|{{\vec n}_m},\vec \alpha } \right) &= \frac{{p\left( {{{\vec n}_m}|{{\vec \vartheta }_m}} \right)p\left( {{{\vec \vartheta }_m}|\vec \alpha } \right)}}{{\int {p\left( {{{\vec n}_m}|{{\vec \vartheta }_m}} \right)p\left( {{{\vec \vartheta }_m}|\vec \alpha } \right)d{{\vec \vartheta }_m}} }}\\ &= \frac{{mult\left( {{{\vec n}_m}|{{\vec \vartheta }_m},{N_m}} \right)Dir\left( {{{\vec \vartheta }_m}|\vec \alpha } \right)}}{{\int {mult\left( {{{\vec n}_m}|{{\vec \vartheta }_m},{N_m}} \right)Dir\left( {{{\vec \vartheta }_m}|\vec \alpha } \right)d{{\vec \vartheta }_m}} }}\\ &= Dir\left( {{{\vec \vartheta }_m}|\vec \alpha + {{\vec n}_m}} \right)\tag{75}\end{align*}$

• ${\vec \vartheta _m}$ 的贝叶斯估计值为

$\begin{align*}E\left( {{{\vec \vartheta }_m}} \right) = \left( {\frac{{{\alpha _1} + n_m^{\left( 1 \right)}}}{{\sum\nolimits_{k = 1}^K {\left( {{\alpha _k} + n_m^{\left( k \right)}} \right)} }}, \cdot \cdot \cdot ,\frac{{{\alpha _K} + n_m^{\left( K \right)}}}{{\sum\nolimits_{k = 1}^K {\left( {{\alpha _k} + n_m^{\left( k \right)}} \right)} }}} \right) \tag{76} \end{align*}$

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

$\begin{align*}p\left( {{{\vec z}_m}|\vec \alpha } \right) &= \int {p\left( {{{\vec z}_m}|{{\vec \vartheta }_m}} \right)p\left( {{{\vec \vartheta }_m}|\vec \alpha } \right)d{{\vec \vartheta }_m}} \\ &= \int {\prod\limits_{k = 1}^K {{{\left( {\vartheta _m^{\left( k \right)}} \right)}^{n_m^{\left( k \right)}}}\frac{1}{{\Delta \left( {\vec \alpha } \right)}}\prod\limits_{k = 1}^K {{{\left( {\vartheta _m^{\left( k \right)}} \right)}^{{\alpha _k} - 1}}} } d{{\vec \vartheta }_m}} \\ &= \frac{1}{{\Delta \left( {\vec \alpha } \right)}}\int {\prod\limits_{k = 1}^K {{{\left( {\vartheta _m^{\left( k \right)}} \right)}^{n_m^{\left( k \right)} + {\alpha _k} - 1}}} d{{\vec \vartheta }_m}} \\ &= \frac{{\Delta \left( {\vec \alpha + {{\vec n}_m}} \right)}}{{\Delta \left( {\vec \alpha } \right)}} \tag{77}\end{align*}$

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

$\begin{align*}p\left( {\vec z|\vec \alpha } \right) = \prod\limits_{m = 1}^M {p\left( {{{\vec z}_m}|\vec \alpha } \right)} = \prod\limits_{m = 1}^M {\frac{{\Delta \left( {\vec \alpha + {{\vec n}_m}} \right)}}{{\Delta \left( {\vec \alpha } \right)}}}\tag{78}\end{align*}$

2）$\vec \beta \underbrace {\xrightarrow{{\begin{array}{*{20}{c}} {}&{} \end{array}}}}_{Dir}{\vec \varphi _k}\underbrace {\xrightarrow{{\begin{array}{*{20}{c}} {}&{} \end{array}}}}_{Mult}{\vec w_k}$：第一步对 $Dir$ 分布进行 $K$ 次采样得到样本$\left\{ {{{\vec \varphi }_k}} \right\}_{k = 1}^K$（从第二个坛子中独立地抽取了 $K$ 个topic-word骰子 $\left\{ {{{\vec \varphi }_k}} \right\}_{k = 1}^K$）；第二步根据之前得到的主题 $\vec z$，为每个 ${z_{m,n}}$ 生成对应的词 ${w_{m,n}}$，$\vec z$ 中的值有 $K$ 种不同的取值（因为我们假设语料有 $K$ 个主题），所以可以将 $\vec z$ 中的元素分为 $K$ 类。我们现在为第 $k$ 个主题生成对应的词，那么需要选择编号为 $k$ 的topic-word骰子，该骰子有 $V$ 个面，每个面表示一个词，那么在一次投掷骰子过程中，每个词的概率为 ${\vec \varphi _k} = \left( {\varphi _k^{\left( 1 \right)},\varphi _k^{\left( 2 \right)}, \cdot \cdot \cdot ,\varphi _k^{\left( V \right)}} \right)$ ，第 $k$ 个主题有 ${N_k}$ 个词，所以需要投掷 ${N_k}$ 次骰子，为该主题生成 ${N_k}$ 个词。在 ${N_k}$ 次投掷过程中，每个词出现的次数为 ${\vec n_k} = \left( {n_k^{\left( 1 \right)},n_k^{\left( 2 \right)}, \cdot \cdot \cdot ,n_k^{\left( V \right)}} \right)$，那么 ${\vec n_k}$ 服从多项式分布 $Mult\left( {{{\vec n}_k}|{{\vec \varphi }_k},{N_k}} \right)$，可以采用贝叶斯估计对参数 ${\vec \varphi _k}$ 进行估计。

• ${\vec \varphi _k}$ 的先验分布为 $Dir\left( {{{\vec \varphi }_k}|\vec \beta } \right)$

• 后验分布为（推导过程可以参考1.7节）

$\begin{align*} p\left( {{{\vec \varphi }_k}|{{\vec n}_k},\vec \beta } \right) &= \frac{{p\left( {{{\vec n}_k}|{{\vec \varphi }_k}} \right)p\left( {{{\vec \varphi }_k}|\vec \beta } \right)}}{{\int {p\left( {{{\vec n}_k}|{{\vec \varphi }_k}} \right)p\left( {{{\vec \varphi }_k}|\vec \beta } \right)d{{\vec \varphi }_k}} }}\\ &= \frac{{mult\left( {{{\vec n}_k}|{{\vec \varphi }_k},{N_k}} \right)Dir\left( {{{\vec \varphi }_k}|\vec \beta } \right)}}{{\int {mult\left( {{{\vec n}_k}|{{\vec \varphi }_k},{N_k}} \right)Dir\left( {{{\vec \varphi }_k}|\vec \beta } \right)d{{\vec \varphi }_k}} }}\\ &= Dir\left( {{{\vec \varphi }_k}|\vec \beta + {{\vec n}_k}} \right) \tag{79}\end{align*}$

• ${\vec \varphi _k}$ 的贝叶斯估计值为

$\begin{align*}E\left( {{{\vec \varphi }_k}} \right) = \left( {\frac{{{\beta _1} + n_k^{\left( 1 \right)}}}{{\sum\nolimits_{t = 1}^V {\left( {{\beta _t} + n_k^{\left( t \right)}} \right)} }}, \cdot \cdot \cdot ,\frac{{{\beta _V} + n_k^{\left( V \right)}}}{{\sum\nolimits_{t = 1}^V {\left( {{\beta _t} + n_k^{\left( t \right)}} \right)} }}} \right) \tag{80}\end{align*}$

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

$\begin{align*}p\left( {{{\vec w}_k}|\vec \beta } \right) &= \int {p\left( {{{\vec w}_k}|{{\vec \varphi }_k}} \right)p\left( {{{\vec \varphi }_k}|\vec \beta } \right)d{{\vec \varphi }_k}} \\ &= \int {\prod\limits_{t = 1}^V {{{\left( {\varphi _k^{\left( t \right)}} \right)}^{n_k^{\left( t \right)}}}\frac{1}{{\Delta \left( {\vec \beta } \right)}}\prod\limits_{t = 1}^V {{{\left( {\varphi _k^{\left( t \right)}} \right)}^{{\beta _t} - 1}}} } d{{\vec \varphi }_k}} \\ &= \frac{1}{{\Delta \left( {\vec \beta } \right)}}\int {\prod\limits_{t = 1}^V {{{\left( {\varphi _k^{\left( t \right)}} \right)}^{n_k^{\left( t \right)} + {\beta _t} - 1}}} d{{\vec \varphi }_k}} \\ &= \frac{{\Delta \left( {\vec \beta + {{\vec n}_k}} \right)}}{{\Delta \left( {\vec \beta } \right)}} \tag{81} \end{align*}$

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

$\begin{align*}p\left( {\vec w|\vec z,\vec \beta } \right) = \prod\limits_{k = 1}^K {p\left( {{{\vec w}_k}|\vec \beta } \right)} = \prod\limits_{k = 1}^K {\frac{{\Delta \left( {\vec \beta + {{\vec n}_k}} \right)}}{{\Delta \left( {\vec \beta } \right)}}} \tag{82} \end{align*}$

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

$\begin{align*}p\left( {\vec w,\vec z|\vec \alpha ,\vec \beta } \right) &= p\left( {\vec w|\vec z,\vec \beta } \right)p\left( {\vec z|\vec \alpha } \right)\\ &= \prod\limits_{k = 1}^K {\frac{{\Delta \left( {\vec \beta + {{\vec n}_k}} \right)}}{{\Delta \left( {\vec \beta } \right)}}\prod\limits_{m = 1}^M {\frac{{\Delta \left( {\vec \alpha + {{\vec n}_m}} \right)}}{{\Delta \left( {\vec \alpha } \right)}}} } \tag{83} \end{align*}$

2. Gibbs Sampling

上面我们已经推导出参数的贝叶斯估计公式，但是仍然存在一个问题，公式中的 ${\vec n_k}$ 和 ${\vec n_m}$ 无法根据语料直接得到，如果我们知道语料中的每个词的主题，即得到 $\vec z$，那么就可以推断出 ${\vec n_k}$ 和 ${\vec n_m}$，进一步就可以得出贝叶斯的参数估计。
我们需要利用 Gibbs Sampling 对 $p\left( {\vec z|\vec w} \right)$ 进行采样来得到 $\vec z$。根据1.10节 Gibbs Sampling 的原理可知，我们首先需要推导条件概率 $p\left( {{z_i} = k|{{\vec z}_{ - i}},\vec w} \right)$。 先介绍一些符号定义。

$i = \left( {m,n} \right)$：下标索引；$- i$：表示去除下标为 $i$ 的词；${w_i} = t$：第 $m$ 篇文档中第 $n$ 个词为 $t$；${z_i} = k$：第 $m$ 篇文档中第 $n$ 个词的主题为 $k$；$n_{k, - i}^{\left( t \right)}$：除去下标为 $i$ 这个词，剩下的所有词中，词 $t$ 属于主题 $k$ 的统计次数，${\vec n_{k, - i}} = \left( {n_k^{\left( 1 \right)},n_k^{\left( 2 \right)}, \cdot \cdot \cdot ,n_k^{\left( t \right)} - 1, \cdot \cdot \cdot ,n_k^{\left( V \right)}} \right)$（这里假设${w_i} = t,{z_i} = k$）；$n_{m, - i}^{\left( k \right)}$：除去下标为 $i$ 的这个词，第 $m$ 篇文档中主题 $k$ 产生词的个数， ${\vec n_{m, - i}} = \left( {n_m^{\left( 1 \right)},n_m^{\left( 2 \right)}, \cdot \cdot \cdot ,n_m^{\left( k \right)} - 1, \cdot \cdot \cdot ,n_m^{\left( K \right)}} \right)$ （这里假设 ${z_i} = k$ ）；$\vec z = \left\{ {{z_i} = k,{{\vec z}_{ - i}}} \right\}$：语料的主题；$\vec w = \left\{ {{w_i} = t,{{\vec w}_{ - i}}} \right\}$：语料的单词。

$\begin{align*}p\left( {{z_i} = k|{{\vec z}_{ - i}},\vec w} \right) &= \frac{{p\left( {\vec w,\vec z} \right)}}{{p\left( {\vec w,{{\vec z}_{ - i}}} \right)}}\\ &= \frac{{p\left( {\vec w|\vec z} \right)}}{{p\left( {{{\vec w}_{ - i}}|{{\vec z}_{ - i}}} \right)p\left( {{w_i}} \right)}} \cdot \frac{{p\left( {\vec z} \right)}}{{p\left( {{{\vec z}_{ - i}}} \right)}}\\ &\propto \frac{{\Delta \left( {{{\vec n}_k} + \vec \beta } \right)}}{{\Delta \left( {{{\vec n}_{k, - i}} + \vec \beta } \right)}} \cdot \frac{{\Delta \left( {{{\vec n}_m} + \vec \alpha } \right)}}{{\Delta \left( {{{\vec n}_{m, - i}} + \vec \alpha } \right)}}\\ &= \frac{{n_{k, - i}^{\left( t \right)} + {\beta ^{\left( t \right)}}}}{{\sum\nolimits_{v = 1}^V {\left( {n_{k, - i}^{\left( v \right)} + {\beta ^{\left( v \right)}}} \right)} }} \cdot \frac{{n_{m, - i}^{\left( k \right)} + {\alpha ^{\left( k \right)}}}}{{\sum\nolimits_{j = 1}^K {\left( {n_{m, - i}^{\left( j \right)} + {\alpha ^{\left( j \right)}}} \right)} }} \tag{84} \end{align*}$

1）$p\left( {{{\vec w}_{ - i}}|{{\vec z}_{ - i}}} \right)$ 和 $p\left( {{{\vec z}_{ - i}}} \right)$ 的计算过程类似 $p\left( {\vec w|\vec z} \right)$ 和 $p\left( {\vec z} \right)$，仅仅在计算的时候不考虑下标为 $i$ 的这个词，我们假设 ${w_i} = t,{z_i} = k$ ；当已知语料时，$p\left( {{w_i}} \right)$可以从语料中统计出来，所以可以认为是常量。

$\begin{align*}p\left( {\vec w|\vec z} \right) &= \Delta \left( {\vec \beta + {{\vec n}_k}} \right)\prod\limits_{z = 1,z \ne k}^K {\frac{{\Delta \left( {\vec \beta + {{\vec n}_z}} \right)}}{{\Delta \left( {\vec \beta } \right)}}} \\ p\left( {{{\vec w}_{ - i}}|{{\vec z}_{ - i}}} \right) &= \Delta \left( {\vec \beta + {{\vec n}_{k, - i}}} \right)\prod\limits_{z = 1,z \ne k}^K {\frac{{\Delta \left( {\vec \beta + {{\vec n}_z}} \right)}}{{\Delta \left( {\vec \beta } \right)}}} \\ p\left( {\vec z} \right) &= \Delta \left( {\vec \alpha + {{\vec n}_m}} \right)\prod\limits_{i = 1,i \ne m}^M {\frac{{\Delta \left( {\vec \alpha + {{\vec n}_i}} \right)}}{{\Delta \left( {\vec \alpha } \right)}}} \\ p\left( {{{\vec z}_{ - i}}} \right) &= \Delta \left( {\vec \alpha + {{\vec n}_{m, - i}}} \right)\prod\limits_{i = 1,i \ne m}^M {\frac{{\Delta \left( {\vec \alpha + {{\vec n}_i}} \right)}}{{\Delta \left( {\vec \alpha } \right)}}} \tag{85} \end{align*}$

2）我们是推断 $i = \left( {m,n} \right)$ 词 $t$ 的主题为 $k$ 的条件概率

$\begin{align*}\Delta \left( {{{\vec n}_k} + \vec \beta } \right) &= \frac{{\prod\nolimits_{v = 1}^V {\Gamma \left( {n_k^{\left( v \right)} + {\beta ^{\left( v \right)}}} \right)} }}{{\Gamma \left( {\sum\nolimits_{v = 1}^V {\left( {n_k^{\left( v \right)} + {\beta ^{\left( v \right)}}} \right)} } \right)}}\\ \Delta \left( {{{\vec n}_{k, - i}} + \vec \beta } \right) &= \frac{{\Gamma \left( {n_k^{\left( t \right)} + {\beta ^{\left( t \right)}} - 1} \right)\prod\nolimits_{v = 1,v \ne t}^V {\Gamma \left( {n_k^{\left( v \right)} + {\beta ^{\left( v \right)}}} \right)} }}{{\Gamma \left( {\sum\nolimits_{v = 1,v \ne t}^V {\left( {n_k^{\left( v \right)} + {\beta ^{\left( v \right)}}} \right)} + n_k^{\left( t \right)} + {\beta ^{\left( t \right)}} - 1} \right)}}\\ \Delta \left( {{{\vec n}_m} + \vec \alpha } \right) &= \frac{{\prod\nolimits_{j = 1}^K {\Gamma \left( {n_m^{\left( j \right)} + {\alpha ^{\left( j \right)}}} \right)} }}{{\Gamma \left( {\sum\nolimits_{j = 1}^K {\left( {n_m^{\left( j \right)} + {\alpha ^{\left( j \right)}}} \right)} } \right)}}\\ \Delta \left( {{{\vec n}_{m, - i}} + \vec \alpha } \right) &= \frac{{\Gamma \left( {n_m^{\left( k \right)} + {\alpha ^{\left( k \right)}} - 1} \right)\prod\nolimits_{j = 1,j \ne k}^K {\Gamma \left( {n_m^{\left( j \right)} + {\alpha ^{\left( j \right)}}} \right)} }}{{\Gamma \left( {\sum\nolimits_{j = 1,j \ne k}^K {\left( {n_m^{\left( j \right)} + {\alpha ^{\left( j \right)}}} \right)} + n_m^{\left( k \right)} + {\alpha ^{\left( k \right)}} - 1} \right)}} \end{align*}$

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

$\begin{align*}p\left( {{z_i} = k|{{\vec z}_{ - i}},\vec w} \right) &\propto p\left( {{z_i} = k,{w_i} = t|{{\vec z}_{ - i}},{{\vec w}_{ - i}}} \right)\\ &= \int {p\left( {{z_i} = k,{w_i} = t,{{\vec \vartheta }_m},{{\vec \varphi }_k}|{{\vec z}_{ - i}},{{\vec w}_{ - i}}} \right)} d{{\vec \vartheta }_m}d{{\vec \varphi }_k}\\ &= \int {p\left( {{z_i} = k,{{\vec \vartheta }_m}|{{\vec z}_{ - i}},{{\vec w}_{ - i}}} \right) \cdot } p\left( {{w_i} = t,{{\vec \varphi }_k}|{{\vec z}_{ - i}},{{\vec w}_{ - i}}} \right)d{{\vec \vartheta }_m}d{{\vec \varphi }_k}\\ &= \int {p\left( {{z_i} = k|{{\vec \vartheta }_m}} \right)p\left( {{{\vec \vartheta }_m}|{{\vec z}_{ - i}},{{\vec w}_{ - i}}} \right) \cdot } p\left( {{w_i} = t|{{\vec \varphi }_k}} \right)p\left( {{{\vec \varphi }_k}|{{\vec z}_{ - i}},{{\vec w}_{ - i}}} \right)d{{\vec \vartheta }_m}d{{\vec \varphi }_k}\\ &= \int {p\left( {{z_i} = k|{{\vec \vartheta }_m}} \right)Dir\left( {{{\vec \vartheta }_m}|{{\vec n}_{m, - i}} + \vec \alpha } \right) \cdot } p\left( {{w_i} = t|{{\vec \varphi }_k}} \right)Dir\left( {{{\vec \varphi }_k}|{{\vec n}_{k, - i}} + \vec \beta } \right)d{{\vec \vartheta }_m}d{{\vec \varphi }_k}\\ &= \int {\vartheta _m^{\left( k \right)}Dir\left( {{{\vec \vartheta }_m}|{{\vec n}_{m, - i}} + \vec \alpha } \right)d{{\vec \vartheta }_m}} \cdot \int {\varphi _k^{\left( t \right)}Dir\left( {{{\vec \varphi }_k}|{{\vec n}_{k, - i}} + \vec \beta } \right)d{{\vec \varphi }_k}} \\ &= E\left( {\vartheta _m^{\left( k \right)}} \right) \cdot E\left( {\varphi _k^{\left( t \right)}} \right)\\ &= \frac{{n_{m, - i}^{\left( k \right)} + {\alpha ^{\left( k \right)}}}}{{\sum\nolimits_{j = 1}^K {\left( {n_{m, - i}^{\left( j \right)} + {\alpha ^{\left( j \right)}}} \right)} }} \cdot \frac{{n_{k, - i}^{\left( t \right)} + {\beta ^{\left( t \right)}}}}{{\sum\nolimits_{v = 1}^V {\left( {n_{k, - i}^{\left( v \right)} + {\beta ^{\left( v \right)}}} \right)} }}\end{align*}$
已经推导出条件概率，可以用Gibbs Sampling公式进行采样了。

Algorithm 2.1 LDA Gibbs sampling

1. Initialisation
   $n_m^{\left( k \right)} = 0,{n_m} = 0,n_k^{\left( t \right)} = 0,{n_k} = 0$
   For all documents $m \in \left[ {1,M} \right]$ do:
   For all words $n \in \left[ {1,{N_m}} \right]$ in document m do:

   • Sample ${z_{m,n}} = k \sim Mult\left( {\frac{1}{K}} \right)$,${w_{m,n}} = t$

   • $\begin{align*}n_m^{\left( k \right)}{\rm{ = }}n_m^{\left( k \right)}{\rm{ + }}1,{n_m} = {n_m} + 1\\ n_k^{\left( t \right)} = n_k^{\left( t \right)} + 1,{n_k} = {n_k} + 1\end{align*}$

2. Gibbs Sampling
   While not finished do:
   For all documents $m \in \left[ {1,M} \right]$ do:
   For all words $n \in \left[ {1,{N_m}} \right]$ in document m do:

   • ${z_{m,n}} = k,{w_{m,n}} = t \Rightarrow \begin{array}{*{20}{c}} {n_m^{\left( k \right)}{\rm{ = }}n_m^{\left( k \right)} - 1,{n_m} = {n_m} - 1}\\ {n_k^{\left( t \right)} = n_k^{\left( t \right)} - 1,{n_k} = {n_k} - 1} \end{array}$
   • Mult. sampling acc. to Eq. 85 sample ${z_{m,n}} = \tilde k \sim p\left( {{z_i}|{{\vec z}_{ - i}},\vec w} \right)$
     $\begin{array}{l} n_m^{\left( {\tilde k} \right)}{\rm{ = }}n_m^{\left( {\tilde k} \right)}{\rm{ + }}1,{n_m} = {n_m} + 1\\ n_{\tilde k}^{\left( t \right)} = n_{\tilde k}^{\left( t \right)} + 1,{n_{\tilde k}} = {n_{\tilde k}} + 1 \end{array}$

3. Parameter estimation
   according to Eq. 76 and Eq. 80 estimate ${\vec \vartheta _m},{\vec \varphi _k}$.

参考文献
[1] Parameter estimation for text analysis
[2] Probabilistic Latent Semantic Analysis
[3] Latent Dirichlet Allocation
[4] The EM algorithm