之前整理过两篇关于主题模型的博客《文本建模之Unigram Model，PLSA与LDA》和《再看LDA主题模型》，主要是整理了主题模型的由来和推导过程，关于模型参数怎么计算没有过多涉及，因此接下来将分两篇博客，分别整理PLSA模型和EM算法求解，LDA模型和Gibbs Sample求解。

PLSA

首先回顾下PLSA，作为生成模型，其在文本生成过程中，引入主题的概念，即先从$K$K个主题中选定一个主题，然后在该主题下生成对应的单词。这里，我们将文档$d_i$di​下每个主题$z_k$zk​的生成概率表示为$p(z_k|d_i)$p(zk​∣di​)，每个主题$z_k$zk​下单词（这里的单词是只词汇表$V$V所对应的每个不同单词）$w_j$wj​的生成概率表示为$p(w_j|z_k)$p(wj​∣zk​)。因此文档$d_i$di​中单词$w_j$wj​的生成概率可以表示为
$p(d_i,w_j)=p(d_i)p(w_j|d_i)$p(di​,wj​)=p(di​)p(wj​∣di​)

根据概率论公式：边缘分布等于联合分布的和有：
$\begin{aligned} p(d_i,w_j)&=p(d_i)\sum_{k=1}^Kp(w_j,z_k|d_i)\\ &=p(d_i)\sum_{k=1}^Kp(w_j|z_k)p(z_k|d_i) \end{aligned}$p(di​,wj​)​=p(di​)k=1∑K​p(wj​,zk​∣di​)=p(di​)k=1∑K​p(wj​∣zk​)p(zk​∣di​)​

单词之间相互独立，因此文档$d_i$di​中所有单词的生成概率为
$p(d_i,\vec{w})=\prod_{j=1}^Np(d_i,w_j)^{n(d_i,w_j)}$p(di​,w)=j=1∏N​p(di​,wj​)n(di​,wj​)

其中，$\vec{w}$w是一个$N$N维向量$\Big(n(d_i,w_1),n(d_i,w_2),\cdots,n(d_i,w_N)\Big )$(n(di​,w1​),n(di​,w2​),⋯,n(di​,wN​))，每个分量是单词$w_j$wj​在文档$d_i$di​中的出现次数，$N$N是词汇表$V$V的大小。
文档之间同样相互独立，因此语料的生成概率为
$\begin{aligned} p(D)&=\prod_{i=1}^Mp(d_i,\vec{w})\\ &=\prod_{i=1}^M\prod_{j=1}^Np(d_i,w_j)^{n(d_i,w_j)} \end{aligned}$p(D)​=i=1∏M​p(di​,w)=i=1∏M​j=1∏N​p(di​,wj​)n(di​,wj​)​

采用最大似然估计，最大化$\log p(D)$logp(D)，对应的$p(z_k|d_i)$p(zk​∣di​)和$p(w_j|z_k)$p(wj​∣zk​)就是我们要找的最优参数。
$\begin{aligned} \mathcal{L}=\log p(D)&=\log \prod_{i=1}^M\prod_{j=1}^Np(d_i,w_j)^{n(d_i,w_j)}\\ &=\sum_{i=1}^M\sum_{j=1}^Nn(d_i,w_j)\log p(d_i,w_j)\\ &=\sum_{i=1}^M\sum_{j=1}^Nn(d_i,w_j)\log [p(d_i)\sum_{k=1}^Kp(w_j|z_k)p(z_k|d_i)]\\ &=\sum_{i=1}^M\sum_{j=1}^Nn(d_i,w_j)\Big[\log p(d_i)+\log\sum_{k=1}^Kp(w_j|z_k)p(z_k|d_i)\Big ]\\ &=\sum_{i=1}^M\Big(\sum_{j=1}^Nn(d_i,w_j)\log p(d_i)+\sum_{j=1}^Nn(d_i,w_j)\log\sum_{k=1}^Kp(w_j|z_k)p(z_k|d_i)\Big)\\ &=\sum_{i=1}^M\Big(n(d_i)\log p(d_i)+\sum_{j=1}^Nn(d_i,w_j)\log\sum_{k=1}^Kp(w_j|z_k)p(z_k|d_i)\Big)\\ &=\sum_{i=1}^Mn(d_i)\Big[\log p(d_i)+\sum_{j=1}^N\frac{n(d_i,w_j)}{n(d_i)}\log\sum_{k=1}^Kp(w_j|z_k)p(z_k|d_i)\Big] \end{aligned}$L=logp(D)​=logi=1∏M​j=1∏N​p(di​,wj​)n(di​,wj​)=i=1∑M​j=1∑N​n(di​,wj​)logp(di​,wj​)=i=1∑M​j=1∑N​n(di​,wj​)log[p(di​)k=1∑K​p(wj​∣zk​)p(zk​∣di​)]=i=1∑M​j=1∑N​n(di​,wj​)[logp(di​)+logk=1∑K​p(wj​∣zk​)p(zk​∣di​)]=i=1∑M​(j=1∑N​n(di​,wj​)logp(di​)+j=1∑N​n(di​,wj​)logk=1∑K​p(wj​∣zk​)p(zk​∣di​))=i=1∑M​(n(di​)logp(di​)+j=1∑N​n(di​,wj​)logk=1∑K​p(wj​∣zk​)p(zk​∣di​))=i=1∑M​n(di​)[logp(di​)+j=1∑N​n(di​)n(di​,wj​)​logk=1∑K​p(wj​∣zk​)p(zk​∣di​)]​

因为文档长度$n(d_i)$n(di​)和文档概率p(d_i)可以单独计算，将其去掉不会影响似然函数的最优化。所以在后面部分，我们统一将$\mathcal{L}$L写成如下形式
$\mathcal{L}=\sum_{i=1}^M\sum_{j=1}^Nn(d_i,w_j)\log\sum_{k=1}^Kp(w_j|z_k)p(z_k|d_i)$L=i=1∑M​j=1∑N​n(di​,wj​)logk=1∑K​p(wj​∣zk​)p(zk​∣di​)

可以看到上式对数函数中包含参数的和，进行求导时将非常的复杂，即直接找到对数似然函数为零的最大值将不那么容易。庆幸的是，EM（Expectation-Maximization，简称EM）算法能够帮助我们解决这个问题。

EM算法

EM算法是一种启发式迭代算法，每次迭代分为两步：E步，求期望（expectation）；M步，极大化对数似然（maximization），这也是算法名称的由来。回到我们的问题，我们将对数函数包含参数和的情况进行一般化，首先目标函数可以写成下式：
$\theta = \arg \max_{\theta}\sum_{i=1}^m\log p(x^{(i)};\theta)$θ=argθmax​i=1∑m​logp(x(i);θ)

$p(x^{(i)};\theta)$p(x(i);θ)是在给定参数$\theta$θ下第$i$i个样本的产生概率，有些情况下这个概率值受到一些隐含因素的影响，如主题模型中文档生成过程中的主题，将其用$z^{(i)}$z(i)表示，此时$p(x^{(i)};\theta)$p(x(i);θ)需要累加所有的$p(x^{(i)},z^{(i)};\theta)$p(x(i),z(i);θ)，而优化目标也变成：
$\theta=\arg\max_\theta\sum_{i=1}^m\log p(x^{(i)};\theta)=\arg\max_\theta\sum_{i=1}^m\log\sum_{z^{(i)}}p(x^{(i)},z^{(i)};\theta)$θ=argθmax​i=1∑m​logp(x(i);θ)=argθmax​i=1∑m​logz(i)∑​p(x(i),z(i);θ)

上式的难办之处就在于对数函数中出现了和的形式，那有没有一个办法将求和符号提到$\log$log外面吗？答案是有的。

Jensen不等式

设$f$f是定义在实数域上的函数，如果对于任意的实数$x$x，都有
$f''\geq0$f′′≥0

那么$f$f是凸函数，此时对于随机变量$x$x，有
$E[f(x)]\geq f(E(x))$E[f(x)]≥f(E(x))

等号成立的条件是随机变量$x$x是常量。Jensen不等式用语言描述或许方便记忆：对于凸函数，函数的期望大于等于期望的函数。 证明如下：

如上图所示，对于$(x,f(x))$(x,f(x))，$(y,f(y))$(y,f(y))线段内的任意一点，其横纵坐标可以分别表示为$tx+(1-t)y$tx+(1−t)y，$tf(x)+(1-t)f(y)$tf(x)+(1−t)f(y)，而函数$f(x)$f(x)在$tx+(1-t)y$tx+(1−t)y对应的纵坐标为$f(tf(x)+(1-t)f(y))$f(tf(x)+(1−t)f(y))，因为函数$f(x)$f(x)是凸函数，有
$tf(x)+(1-t)f(y)\geq f(tx+(1-t)y)$tf(x)+(1−t)f(y)≥f(tx+(1−t)y)

因为$t+(1-t)=1$t+(1−t)=1，所以我们可以将$t$t和$1-t$1−t看作随机变量$x,y$x,y对应的概率，即上式表达的是
$E[f(x)]\geq f(E(x))$E[f(x)]≥f(E(x))

若函数$f$f是凹函数，上述不等式符号相反。

EM算法推导

有了Jensen不等式，我们可以对上面的对数似然函数进行缩放如下：
$\begin{aligned} \sum_{i=1}^m\log\sum_{z^{(i)}}p(x^{(i)},z^{(i)};\theta)&=\sum_{i=1}^m\log\sum_{z^{(i)}}Q(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q(z^{(i)})} \tag{1} \end{aligned}$i=1∑m​logz(i)∑​p(x(i),z(i);θ)​=i=1∑m​logz(i)∑​Q(z(i))Q(z(i))p(x(i),z(i);θ)​​(1)$\quad\quad\qquad\qquad\qquad\qquad\quad\geq\sum_{i=1}^m\sum_{z^{(i)}}Q(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\theta)}{Q(z^{(i)})}\tag2$≥i=1∑m​z(i)∑​Q(z(i))logQ(z(i))p(x(i),z(i);θ)​(2)

其中$\sum_{z^{(i)}}Q(z^{(i)})=1$∑z(i)​Q(z(i))=1，即隐变量$z^{(i)}$z(i)的概率分布。在上式中，将$Q(z^{(i)})$Q(z(i))看作随机变量$\frac{p(x^{(i)},z^{(i)};\theta)}{Q(z^{(i)})}$Q(z(i))p(x(i),z(i);θ)​的概率，函数$f$f对应$\log$log函数，根据Jensen不等式就能得到上面不等式。
OK，到这里我们就成功的将求和符号提到了对数函数外面，因而式$(2)$(2)将容易求导。但是式$(2)$(2)和式$(1)$(1)是不等号关系，式$(2)$(2)的最大值不一定是式$(1)$(1)的最大值，那怎么才能得到式$(1)$(1)的最大值呢？我们知道，当给定参数$\theta$θ时，对数似然函数$\mathcal{L}$L的值就确定了，也就是式$(2)$(2)的上界就确定了，同时$p(x^{(i)},z^{(i)})$p(x(i),z(i))的值就确定了，这时我们修改$Q(z^{(i)})$Q(z(i))的值，使下界等于上界。

如果要满足Jensen不等式的等号，需要：
$\begin{aligned} &\quad \frac{p(x^{(i)},z^{(i)};\theta)}{Q(z^{(i)})}=c,\quad c为常量\\ &\rArr p(x^{(i)},z^{(i)};\theta)=cQ(z^{(i)})\\ &\rArr \sum_{z^{(i)}}p(x^{(i)},z^{(i)};\theta)=c\sum_{z^{(i)}}Q(z^{(i)}) \\ &\rArr \sum_{z^{(i)}}p(x^{(i)},z^{(i)};\theta)=c \end{aligned}$​Q(z(i))p(x(i),z(i);θ)​=c,c为常量⇒p(x(i),z(i);θ)=cQ(z(i))⇒z(i)∑​p(x(i),z(i);θ)=cz(i)∑​Q(z(i))⇒z(i)∑​p(x(i),z(i);θ)=c​

所以，我们可以得到，等号成立需要满足的条件是
$Q(z^{(i)})=\frac{p(x^{(i)},z^{(i)};\theta)}{\sum_{z^{(i)}}p(x^{(i)},z^{(i)};\theta)}=p(z^{(i)}|x^{(i)};\theta)$Q(z(i))=∑z(i)​p(x(i),z(i);θ)p(x(i),z(i);θ)​=p(z(i)∣x(i);θ)

即隐变量的后验概率。知道了$Q(z^{(i)})$Q(z(i))，将其带入式$(2)$(2)，然后找出使其最大的$\theta$θ，这时我们就确定了在该$\theta$θ值处的一个下界，同样我们调整$Q(z^{(i)})$Q(z(i))使这个下界最大化，同时也在尝试提升我们的对数似然函数，即我们需要最大化下式：
$\begin{aligned} &\quad\enspace\arg\max_{\theta}\sum_{i=1}^m\sum_{z^{(i)}}Q(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\theta)}{Q(z^{(i)})}\\ &\rArr\arg\max_{\theta}\sum_{i=1}^m\sum_{z^{(i)}}Q(z^{(i)})\log p(x^{(i)},z^{(i)};\theta)\\ &=\arg\max_{\theta}\sum_{i=1}^m\sum_{z^{(i)}}p(z^{(i)}|x^{(i)};\theta)\log p(x^{(i)},z^{(i)};\theta) \end{aligned}$​argθmax​i=1∑m​z(i)∑​Q(z(i))logQ(z(i))p(x(i),z(i);θ)​⇒argθmax​i=1∑m​z(i)∑​Q(z(i))logp(x(i),z(i);θ)=argθmax​i=1∑m​z(i)∑​p(z(i)∣x(i);θ)logp(x(i),z(i);θ)​

上式可以理解为对数似然函数$\log p(x^{(i)},z^{(i)};\theta)$logp(x(i),z(i);θ)基于条件概率$p(z^{(i)}|x^{(i)};\theta)$p(z(i)∣x(i);θ)的期望，正好对应EM算法的E步，而最大化过程，正好对应M步。重复上述过程，直到收敛到对数似然函数$\mathcal{L}$L的最大处$\theta^*$θ∗，至于为什么会收敛？证明可以参考李航《统计学习方法》等。注意，与所有迭代算法一样，EM算法不能保证得到全局最优解。

EM算法的形象过程如下图：

总结EM算法的流程如下：
输入：观测数据$x=(x^{(1)},x^{(2)},\cdots,x^{(m)})$x=(x(1),x(2),⋯,x(m))，联合分布$p(x^{(i)},z^{(i)};\theta)$p(x(i),z(i);θ)，条件概率$p(z^{(i)}|x^{(i)};\theta)$p(z(i)∣x(i);θ)，最大迭代次数$J$J。

1. 随机初始化模型参数$\theta^0$θ0;
2. E步：计算联合分布的条件概率期望：
   $Q_j(z^{(i)})=p(z^{(i)}|x^{(i)};\theta^j)$Qj​(z(i))=p(z(i)∣x(i);θj)$E=\sum_{i=1}^m\sum_{z^{(i)}}Q(z^{(i)})\log p(x^{(i)},z^{(i)};\theta)$E=i=1∑m​z(i)∑​Q(z(i))logp(x(i),z(i);θ)
3. M步：最大化期望，得到$\theta^{j+1}$θj+1
   $\theta^{j+1}=\arg\max_{\theta}\sum_{i=1}^m\sum_{z^{(i)}}Q_j(z^{(i)})\log p(x^{(i)},z^{(i)};\theta)$θj+1=argθmax​i=1∑m​z(i)∑​Qj​(z(i))logp(x(i),z(i);θ)
4. 重复第（2）步和第（3）步，直至收敛或达到最大迭代次数$J$J。

PLSA模型的EM算法步骤

$\begin{aligned} \mathcal{L}&=\log p(D)=\sum_{i=1}^M\sum_{j=1}^Nn(d_i,w_j)\log\sum_{k=1}^Kp(w_j|z_k)p(z_k|d_i)\\ &\geq\sum_{i=1}^M\sum_{j=1}^Nn(d_i,w_j)\sum_{k=1}^KQ(z_k)\log\frac{p(w_j|z_k)p(z_k|d_i)}{Q(z_k)} \end{aligned}$L​=logp(D)=i=1∑M​j=1∑N​n(di​,wj​)logk=1∑K​p(wj​∣zk​)p(zk​∣di​)≥i=1∑M​j=1∑N​n(di​,wj​)k=1∑K​Q(zk​)logQ(zk​)p(wj​∣zk​)p(zk​∣di​)​​

根据上面的推导得
$Q(z_k)=p(z_k|d_i,w_j)$Q(zk​)=p(zk​∣di​,wj​)

然后开始EM算法迭代过程。
E步：计算联合分布的条件概率期望
（1）计算隐变量的后验概率
$\begin{aligned} p(z_k|d_i,w_j)&=\frac{p(z_k,d_i,w_j)}{\sum_{k=1}^Kp(z_k,d_i,w_j)}\\ &=\frac{p(w_j|z_k)p(z_k|d_i)}{\sum_{k=1}^Kp(w_j|z_k)p(z_k|d_i)} \end{aligned}$p(zk​∣di​,wj​)​=∑k=1K​p(zk​,di​,wj​)p(zk​,di​,wj​)​=∑k=1K​p(wj​∣zk​)p(zk​∣di​)p(wj​∣zk​)p(zk​∣di​)​​

（2）将隐变量的后验概率带入，得到期望函数
$\begin{aligned} \mathcal{L}^c&=\sum_{i=1}^M\sum_{j=1}^Nn(d_i,w_j)\sum_{k=1}^Kp(z_k|d_i,w_j)\log p(z_k,d_i,w_j)\\ &=\sum_{i=1}^M\sum_{j=1}^Nn(d_i,w_j)\sum_{k=1}^Kp(z_k|d_i,w_j)\log p(w_j|z_k)p(z_k|d_i) \end{aligned}$Lc​=i=1∑M​j=1∑N​n(di​,wj​)k=1∑K​p(zk​∣di​,wj​)logp(zk​,di​,wj​)=i=1∑M​j=1∑N​n(di​,wj​)k=1∑K​p(zk​∣di​,wj​)logp(wj​∣zk​)p(zk​∣di​)​

M步：最大化期望函数，这是一个带约束的最优化问题：
$\sum_{j=1}^Np(w_j|z_k)=1\\ \sum_{k=1}^Kp(z_k|d_i)=1$j=1∑N​p(wj​∣zk​)=1k=1∑K​p(zk​∣di​)=1

可以采用拉格朗日乘子法求解，拉格朗日函数如下：
$\mathcal{H}=\mathcal{L}^c+\sum_{k=1}^K\tau_k\Big(1-\sum_{j=1}^Np(w_j|z_k)\Big)+\sum_{i=1}^M\rho\Big(1-\sum_{k=1}^Kp(z_k|d_i)\Big)$H=Lc+k=1∑K​τk​(1−j=1∑N​p(wj​∣zk​))+i=1∑M​ρ(1−k=1∑K​p(zk​∣di​))

分别求偏导，令其结果等于0，最终可以估计出参数$p(w_j|z_k)$p(wj​∣zk​)和$p(z_k|d_i)$p(zk​∣di​)
$p(w_j|z_k)=\frac{\sum_{i=1}^Mn(d_i,w_j)p(z_k|d_i,w_j)}{\sum_{j=1}^N\sum_{i=1}^Mn(d_i,w_j)p(z_k|d_i,w_j)}$p(wj​∣zk​)=∑j=1N​∑i=1M​n(di​,wj​)p(zk​∣di​,wj​)∑i=1M​n(di​,wj​)p(zk​∣di​,wj​)​$p(z_k|d_i)=\frac{\sum_{j=1}^Nn(d_i,w_j)p(z_k|d_i,w_j)}{n(d_i)}$p(zk​∣di​)=n(di​)∑j=1N​n(di​,wj​)p(zk​∣di​,wj​)​

参考文献

EM算法原理总结
从最大似然到EM算法浅解