机器学习算法笔记：EM期望最大

直观案例

案例 a
假设有两硬币 A 和 B ，现在想要知道两枚硬币各自为正面的概率 ${\theta }_A$ 和 ${\theta }_B$，

• 总共做了 5 次试验，每次随机取出 A 或 B 抛10次，并记录 A 和 B 的正反情况，现在 A 和 B 抛的总次数和正面的次数都知道，很容易估计得到：$$\begin{gathered} {\hat{\theta }}_A=\frac{\text{A正面次数}}{A正面+反面总次数}=\frac{24}{24+6}=0.80 \\ {\hat{\theta }}_B=\frac{\text{B正面次数}}{B正面+反面总次数}=\frac{9}{9+11}=0.45 \end{gathered}$$
  

案例 b
与案例 a相比：现在只知道正反面的记录情况，但不知道每次试验取的是 A 还是 B （A 和 B 取出概率相等），那么现在该如何估计 ${\theta }_A$ 和 ${\theta }_B$ 的值呢？
EM算法：

• 1、 ${\hat{\theta }}_A^{(0)}=0.60$ 和 ${\hat{\theta }}_B^{(0)}=0.50$

• 2、 第一次试验 5 正 5 反 ，那么 ： $$\begin{array}{rl}& {}_{\text{ A出现5H5T的概率}}{P}_A=({\hat{\theta }}_A^{(0)}{)}^5(1-{\hat{\theta }}_A^{(0)}{)}^5=(0.60{)}^5(0.40{)}^5\\ & {}_{\text{ B出现5H5T的概率}}{P}_B=({\theta }_B^{(0)}{)}^5(1-{\theta }_B^{(0)}{)}^5=(0.50{)}^5(0.50{)}^5\\ & {}_{\text{ 第一次试验取出的是A的概率}}P(A)=\frac{P_A}{P_A+P_B}=0.45\\ & {}_{\text{ 第一次试验取出的是B的概率}}P(B)=\frac{P_B}{P_A+P_B}=0.55\end{array}$$所以第一次试验可以看作是 A 抛了 $10$次$\times P(A)=10\times 0.45=4.5次$，其中正面为 2.25 次，反面为 2.25 次。同理，B 抛了 $10$次$\times P(A)=10\times 0.55=5.5次$，其中正面为 2.75 次，反面为 2.75 次。

• 3、 同步骤 2 中一样，分别计算第二至第五次得到下表（注意：下表只保留了一位小数）：
  

• 4、 得到上表后，与案例 a 一样：相当于 A 和 B 抛的总次数和正面的次数都已知，就可以对 ${\theta }_A$ 和 ${\theta }_B$进行估计：$$\begin{gathered} {\hat{\theta }}_A^{(1)}=\frac{\text{A正面次数}}{A正面+反面总次数}=\frac{21.3}{21.3+8.6}=0.71 \\ {\hat{\theta }}_B^{(1)}=\frac{\text{B正面次数}}{B正面+反面总次数}=\frac{11.7}{11.7+8.4}=0.58 \end{gathered}$$

• 5、 将 ${\hat{\theta }}_A^{(1)}$ 和 ${\hat{\theta }}_B^{(1)}$ 作为初始值重复 $1\to 4$ ，直到迭代 $t$ 次之后 ${\hat{\theta }}_A^{(t)}$ 和 ${\hat{\theta }}_B^{(t)}$ 趋于稳定,，即可将其作为最终估计值。本例迭代十次后 ${\hat{\theta }}_A^{(10)}=0.80$ 和 ${\hat{\theta }}_B^{(10)}=0.52$，与案例 a中的估计值已经较为接近了

E-step：步骤 $1\to 3$ 根据模型参数 ${\theta }_A$ 和 ${\theta }_B$ 计算期望（隐变量分布），即 A 和 B 分布情况
M-step：步骤 $4$ 根据 A 和 B 分布情况重新估算参数 ${\hat{\theta }}_A$ 和 ${\hat{\theta }}_B$ ，可作为下一次迭代输入

补充

$M-step$ 的 $M$ 由 $MLE$ 体现，令 $i$ 表示第 $i$ 次抛硬币，$x_i=1$ 表示第 $i$ 次为正面 $H$，否则为反面 $F$，则每次抛硬币的概率可表示为：
$$p(x_i,\theta )={\theta }^{x_i}(1-\theta {)}^{1-x_i}=\{\begin{array}{rl}& \theta x_i=1\\ & 1-\theta x_i=0\end{array}$$

每次试验 $N$ 次，用极大似然估计出参数 $\theta$：
$$\begin{array}{rl}\hat{\theta }& =\underset{\theta }{\mathrm{argmax}}\mathcal{L}(\theta )=\underset{\theta }{\mathrm{argmax}}\prod _{i=1}^{N}p(x_i,\theta )\text{MLE}\\ & \simeq \underset{\theta }{\mathrm{argmin}}-\sum _{i=1}^N\log p(x_i,\theta )=\underset{\theta }{\mathrm{argmin}}-\sum _{i=1}^N{\theta }^{x_i}(1-\theta {)}^{1-x_i}\text{NLL}\\ & \frac{\partial \mathcal{L}(\theta )}{\partial \theta }=-\sum _{i=1}^N\left(\frac{x_i}{\theta }+(1-x_i)\frac{-1}{1-\theta }\right)=0\text{令偏导为0}\\ & ⟹\hat{\theta }=\frac{\sum _{i=1}^N{x}_i}{N}\end{array}$$

由此可知 案例 b 中第 $4$ 步已经蕴含了极大似然的思想

EM算法

期望最大算法的目的是解决具有隐变量的混合模型的参数估计，隐变量可以是每个样本类别$z_i$（如案例 b 中的 A 和 B），该变量属于位置状态。

MLE 对 $p(x|\theta )$ 参数的估计记为：${\theta }_{MLE}=\underset{\theta }{\mathrm{argmax}}\log p(x|\theta )$。但多了一个未知隐变量不好求解，于是就要用到 EM 算法对其迭代求解： $${\theta }^{t+1}=\underset{\theta }{\mathrm{argmax}}{\int }_z\log [p(x,z|\theta )]\underset{类似{\hat{\theta }}_A^{(0)}{\hat{\theta }}_A^{(0)}}{\underbrace{p(z|x,{\theta }^t)}}dz=\underset{期望\text{最大化}}{\underbrace{\underset{\theta }{\mathrm{argmax}}\stackrel{\text{期望}}{\overbrace{{\mathbb{E}}_{z|x,{\theta }^t}[\log p(x,z|\theta )]}}}}$$

$p(z|x,{\theta }^t)$ 为隐变量概率分布（当$t=1$时，需要手动初始化其分布，与案例 b 中初始化 ${\hat{\theta }}_A^{(0)}$ 和 ${\hat{\theta }}_B^{(0)}$类似），上式具体可分为两步：

• E step：给定初始参数，计算 $\log p(x,z|\theta )$ 在概率分布 $p(z|x,{\theta }^t)$ 下的期望
• M step：计算得到使期望最大化的参数，新的参数可作为下一次迭代的输入

为什么迭代

怎么保证每次迭代都是朝着更好的方向？也就是需要证明：$\log p(x|{\theta }^t)\le \log p(x|{\theta }^{t+1})$
证明：$\log p(x|\theta )=\log p(z,x|\theta )-\log p(z|x,\theta )$，对左右两边求积分（关于概率分布$p(z|x,{\theta }^t)$ 的积分）： $$\begin{array}{rl}& Left:{\int }_{z}p(z|x,{\theta }^t)\log p(x|\theta )dz=\log p(x|\theta )\underset{1}{\underbrace{{\int }_{z}p(z|x,{\theta }^t)dz}}=\log p(x|\theta )\\ & Right:\underset{Q(\theta ,{\theta }^t)}{\underbrace{{\int }_{z}p(z|x,{\theta }^t)\log p(x,z|\theta )dz}}-\underset{H(\theta ,{\theta }^t)}{\underbrace{{\int }_{z}p(z|x,{\theta }^t)\log p(z|x,\theta )dz}}\\ & =Q(\theta ,{\theta }^t)-H(\theta ,{\theta }^t)\end{array}$$所以： $$\log p(x|\theta )=Q(\theta ,{\theta }^t)-H(\theta ,{\theta }^t)$$ 由于 $$Q(\theta ,{\theta }^t)={\int }_{z}p(z|x,{\theta }^t)\log p(x,z|\theta )dz$$而 ${\theta }^{t+1}$ 本身就是由最大化 $Q(\theta ,{\theta }^t)$ 得到：$${\theta }^{t+1}=\underset{\theta }{\mathrm{argmax}}{\int }_z\log [p(x,z|\theta )]p(z|x,{\theta }^t)dz$$所以自然有： $$Q({\theta }^{t+1},{\theta }^t)\ge Q(\theta ,{\theta }^t)$$$\theta$为随机变量，因此取 $\theta ={\theta }^t$ 也是成立的：$$Q({\theta }^{t+1},{\theta }^t)\ge Q({\theta }^t,{\theta }^t)$$要证 $\log p(x|{\theta }^t)\le \log p(x|{\theta }^{t+1})$，需证：$H({\theta }^t,{\theta }^t)\ge H({\theta }^{t+1},{\theta }^t)$： $$\begin{array}{rl}& H({\theta }^{t+1},{\theta }^t)-H({\theta }^t,{\theta }^t)\\ & ={\int }_{z}p(z|x,{\theta }^t)\log p(z|x,{\theta }^{t+1})dz-{\int }_{z}p(z|x,{\theta }^t)\log p(z|x,{\theta }^t)dz\\ & ={\int }_{z}p(z|x,{\theta }^t)\log \frac{p(z|x,{\theta }^{t+1})}{p(z|x,{\theta }^t)}\\ & =-KL(p(z|x,{\theta }^t),p(z|x,{\theta }^{t+1}))\le 0\end{array}$$ 综合上面的结果： $$\log p(x|{\theta }^t)\le \log p(x|{\theta }^{t+1})$$根据上面的证明可知似然函数每一步都会增大。

EM 公式推导

引入关于 $z$ 的分布 $q(z)\ne 0$：
$$\log p(x|\theta )=\log p(z,x|\theta )-\log p(z|x,\theta )=\log \frac{p(z,x|\theta )}{q(z)}-\log \frac{p(z|x,\theta )}{q(z)}$$

分别对两边求期望 $\mathbb{E}q(z)$： $$\begin{array}{rl}& Left:{\int }_{z}q(z)\log p(x|\theta )dz=\log p(x|\theta )\\ & Right:{\int }_{z}q(z)\log \frac{p(z,x|\theta )}{q(z)}dz-{\int }_{z}q(z)\log \frac{p(z|x,\theta )}{q(z)}dz\\ & =ELBO-KL(q(z),p(z|x,\theta ))\end{array}$$

上式中，Evidence Lower Bound(ELBO)为下界，所以 $\log p(x|\theta )\ge ELBO$，等于号取在 KL 散度为0是，即：$q(z)=p(z|x,\theta )$，EM 算法的目的是将 ELBO 最大化，根据上面的证明过程，在每一步 EM 后，求得了最大的ELBO，并根据这个使 ELBO 最大的参数代入下一步中： $$\hat{\theta }=\underset{\theta }{\mathrm{argmax}}ELBO=\underset{\theta }{\mathrm{argmax}}{\int }_{z}q(z)\log \frac{p(x,z|\theta )}{q(z)}dz$$

由于 $q(z)=p(z|x,{\theta }^t)$ 的时候，这一步的最大值才能取等号，所以：
$$\begin{array}{rl}\hat{\theta }& =\underset{\theta }{\mathrm{argmax}}ELBO=\underset{\theta }{\mathrm{argmax}}{\int }_{z}q(z)\log \frac{p(x,z|\theta )}{q(z)}dz\\ & =\underset{\theta }{\mathrm{argmax}}{\int }_{z}p(z|x,{\theta }^t)\log \frac{p(x,z|\theta )}{p(z|x,{\theta }^t)}dz\\ & =\underset{\theta }{\mathrm{argmax}}{\int }_{z}p(z|x,{\theta }^t)[\log p(x,z|\theta )-\underset{\theta \text{无关项,可看作常数}}{\underbrace{\log p(z|x,{\theta }^t)}}]dz\\ & =\underset{\theta }{\mathrm{argmax}}{\int }_{z}p(z|x,{\theta }^t)\log p(x,z|\theta )\end{array}$$

上式就是 EM 迭代公式

由 Jensen 不等式导出 EM 公式

从 Jensen 不等式出发，也可导出 EM 公式：
$$\begin{array}{rl}\log p(x|\theta )& =\log {\int }_{z}p(x,z|\theta )dz=\log {\int }_z\frac{p(x,z|\theta )q(z)}{q(z)}dz\\ & =\log \mathbb{E}q(z)\left[\frac{p(x,z|\theta )}{q(z)}\right]\ge \mathbb{E}q(z)\left[\log \frac{p(x,z|\theta )}{q(z)}\right]\end{array}$$

上式蓝色项就是 ELBO，等号在 $p(x,z|\theta )=Cq(z)$ 时成立，于是：
$$\begin{array}{r}{\int }_{z}q(z)dz=\frac{1}{C}{\int }_{z}p(x,z|\theta )dz=\frac{1}{C}p(x|\theta )=1\\ \Rightarrow q(z)=\frac{1}{p(x|\theta )}p(x,z|\theta )=p(z|x,\theta )\end{array}$$

其中 $,\theta$ 应该是上一次迭代得到的 ${\theta }^t$，于是就推导得到了 EM 迭代公式

广义 EM

EM 模型解决了概率生成模型的参数估计的问题，通过引入隐变量 $z$来学习模型参数$\theta$，不同的模型对 $z$ 可以有不同的假设。对学习任务 $p(x|\theta )$，就是学习任务 $\frac{p(x,z|\theta )}{p(z|x,\theta )}$。在这个式子中，假定了在 E 步骤中，$q(z)=p(z|x,\theta )$，如果$p(z|x,\theta )$ 无法求解，那么就需要使用采样（MCMC）或者变分推断等方法来近似推断这个后验。观察 KL 散度（衡量两个分布的距离）的表达式，为了最大化 ELBO，在固定的 $\theta$ 时，需要最小化 KL 散度，于是： $$\hat{q}(z)=\underset{q}{\mathrm{argmin}}KL(p,q)=\underset{q}{\mathrm{argmax}}ELBO$$

广义 EM 的基本思路：

E step： 固定模型参数的值，优化隐含数据的分布$${\hat{q}}^{t+1}(z)=\underset{q}{\mathrm{argmax}}{\int }_z{q}^t(z)\log \frac{p(x,z|\theta )}{q^t(z)}dz,fixed\theta$$

M step：固定隐含数据分布，优化模型参数的值 $$\hat{\theta }=\underset{\theta }{\mathrm{argmax}}{\int }_z{q}^{t+1}(z)\log \frac{p(x,z|\theta )}{q^{t+1}(z)}dz,fixed\hat{q}$$

对于上面的积分： $$ELBO={\int }_{z}q(z)\log \frac{p(x,z|\theta )}{q(z)}dz={\mathbb{E}}_{q(z)}[p(x,z|\theta )]+Entropy(q(z))$$

EM 的推广

EM 算法类似于坐标上升法，固定部分坐标，优化其他坐标，再一遍一遍的迭代。如果在 EM 框架中，无法求解 $z$ 后验概率，那么需要采用一些变种的 EM 来估算这个后验。

• 基于平均场的变分推断，VBEM/VEM
• 基于蒙特卡洛的EM，MCEM

参考文献

【1】期望最大
【2】EM算法原理总结
【3】从最大似然到EM算法浅解
【4】What is the expectation maximization algorithm?
【5】EM算法通俗实例（What is the expectation maximization algorithm?）
【6】EM算法求高斯混合模型參数预计——Python实现
【7】EM算法Python实战
【8】EM算法及其推广