「机器学习」EM算法

期望最大(EM)算法的大名想必大家都听过，说来惭愧，大概两年前就学过用过这个算法，但是到后来只记得EM这个名字，把算法的原理都给忘光光了，这几天重新捡起来，总结一下，也希望对想了解这个算法的人有点启发。

对我EM算法我总结为三点：为什么、是什么、怎么用。

为什么要用EM算法？

会看到这篇博客的大家想必都知道，EM算法是用来估计参数的，并且是以迭代的方式进行。看到估计参数这个字眼，你脑子里首先会想到什么方法？毫无疑问就是极大似然法吧，再帮大家回忆一下极大似然法的思想：什么参数对应的事件发生的概率最大，我们就认为这个参数就是未知参数的估计值。通俗一点讲就是你知道了事件的观测值 $y_1,…,y_N$，然后要寻找能够让事件发生的可能性 $p(y_1,…,y_N|\theta)$ 最大的参数 $\hat{\theta}$，也就是

或者通常为了计算方便采用对数的形式

其中 $\log p(y_1,…,y_N|\theta)$ 就称为对数似然函数，常常假设事件之间相互独立，那也就能进一步把联合概率写成独立概率乘积的形式：

这样一个无约束的优化问题利用极值条件(导数为0)也就能得到待估计参数的解析解了。

回到我们的问题为什么要用EM算法，那肯定是因为极大似然法不管用了呗，也就是当有隐变量(latent variable)存在的时候，原来的对数似然函数变为

相当于求完全数据 $Y,Z$ 的联合分布 $p(Y,Z|\theta)$ (把隐变量积分掉)的和，但是由于里面包含隐变量，直接用极大似然求解是不行滴。但是如果能够估计出隐变量 $Z$，那这个问题又回到了原来的极大似然估计。这样也就自然地与EM算法的思想挂上勾了：估计隐变量(E步)，再估计参数(M步)，不断迭代。

现在可以回答第一个问题了：为什么要用EM算法？EM算法是用来解决含有隐变量的参数估计问题的。

什么是EM算法？

那么这个所谓的迭代算法是怎么迭代的呢？我们还是先回到对数似然函数，可以将其表示为

第一行公式(完整数据减去后验)根据条件概率公式可得，第二行公式相当于减去一个 $\log q(Z)$ 再加上一个 $\log q(Z)$ 恒等变换。这里假设 $q(Z)$ 是 $Z$ 的一个概率分布 ($\int_Z q(Z) dZ = 1$)。

然后同时对两侧对$q(Z)$求期望

左边积分还是它本身：

右边积分则是

左边等于右边即

KL divergence的概念

$$ D_{\mathrm{KL}}(P \| Q)=-\sum_{i} P(i) \ln \frac{Q(i)}{P(i)} \ge 0 $$

由吉布斯不等式可知，当且仅当 $P = Q$ 时 $D_{KL}(P||Q)$为零。

因为相对熵大于等于0的性质，所以 ELBO 为 $\log p(Y|\theta)$ 提供了下界，所以EM算法的想法就是想让ELBO不断变大从而让对数似然变大，并且只有当 $q(Z) = p(Z|Y,\theta^{(t)})$ 时ELBO等于对数似然， 即

这也就解释了EM算法是怎么来的啦，当 $q(Z)$ 等于后验 $p(Z|Y,\theta^{{t}})$ 时，ELBO 等于log likelihood，先固定 $\theta^{(t)}$ 得到后验，然后求 ELBO (关于后验的期望)，滑动 $\theta$ 得到期望的最大值，最大期望对应的 $\theta$ 就是 $\theta^{(t+1)}$

从思路上是能说通了，你可能会疑惑为什么这么迭代求解就能得到最优解呢？或者说 $\theta_{t}$ 到最后会收敛吗？我们回到最初的目标最大化对数似然 $\log p(Y|\theta)$，那如果我们每次迭代求解都能使得对应的对数似然函数的值增加，最后是不是也就说明能向似然函数的最大值靠近了呢，也就是我们希望看到的结果是

事实也说明了这是成立的，我们再回到对数似然函数的表达式

两边求关于后验 $p(Z|Y,\theta^{(t)})$ 的期望，也就是求积分

左边积分还是它本身：

右边积分则是：

其中，$Q(\theta,\theta^{(t)}) = \int_Z p(Z|Y,\theta^{(t)}) \log p(Y,Z|\theta) dZ, H(\theta,\theta^{(t)}) = \int_Z p(Z|Y,\theta^{(t)}) \log p(Z|Y,\theta) dZ$

那要使得迭代收敛需要满足的条件就变为了

因为 $Q(\theta,\theta^{(t)})$ 对应的就是EM的迭代公式，那根据最大化问题可以得到 $Q(\theta,\theta^{(t)}) \le Q(\theta^{(t+1)}\theta^{(t)})$，那也就有 $Q(\theta^{(t)},\theta^{(t)}) \le Q(\theta^{(t+1)}\theta^{(t)})$，那再根据放缩的思想，只要我们能证明 $H(\theta^{(t)},\theta^{(t)}) \ge H(\theta^{(t+1)},\theta^{(t)})$，那一切就都顺理成章了。

直接对这两项作差

看出点什么没有，不然再回到上面 KL 散度的定义看看，这结果不就是 $-D_{\text{KL}}(p(Z|Y,\theta^{(t)})||p(Z|Y,\theta^{(t+1)}))$，根据散度大于等于0的性质，可以得到 $H(\theta^{(t+1)},\theta^{(t)}) - H(\theta^{(t)},\theta^{(t)}) \le 0$，这也就得证了EM迭代求解每次得到的似然函数值总是非减的( $\theta^{(t)}$ 会逐渐收敛)，那最终也就会取到最大似然咯，需要注意这并不能保证找到全局最优解的。

怎么用EM算法？

看到这里的话想必你对算法的原理也有了大概的认识，但是大家最关心的应该是在认识原理后如何去应用，这里根据李航老师的《统计学习方法》里的例子介绍一下算法的应用和实现。

三硬币模型：有A，B，C三枚硬币，它们正面向上的概率分别为 $a,b,c$，进行如下试验：先掷硬币 A，如果正面则再掷 B，如果反面则掷 C，掷硬币的结果正面记为1反面为0，独立重复10次试验，结果为1，1，0，1，0，0，1，0，1，1。我们只知道最终掷硬币的结果并不能观测掷硬币的过程，那么如何估计三枚硬币正面向上的概率呢？

这个例子里观测变量 $y_i$ 也就是掷硬币的结果0或1，要估计的参数为概率 $\theta=(a,b,c)$，由于我们不知道掷A的结果导致这一问题不能直接采用极大似然来估计，因此我们引入隐变量 $z$ 来表示掷 A的结果。

观测数据 $y_1,…,y_N$ 由以下模型生成

隐变量 $z_i = 1$ 表示第 $i$ 个观测来自掷硬币B的结果，$z_i = 0$ 表示第 $i$ 个观测来自掷硬币C的结果。

于是可以写出完全数据的似然函数：

那么，完全数据的对数似然函数为

按照我上面说的先固定 $\theta^{(t)}$ 计算后验：

有了隐变量的后验分布后，在对完全数据的对数似然函数求最大期望

又回到这个极大似然的问题，现在相当于隐变量已知咯，那直接按照传统的套路利用极值条件求导获得待参数的解析解。

对于参数 $a$，令 $\partial Q(\theta)/ \partial a = 0$

可得

同理可得

matlab 实现EM算法

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

% 给定观测结果
y = [1;1;0;1;0;0;1;0;1;1];
N = length(y);
% 先随便假设个初始参数 
a_old = 0.4;
b_old = 0.6;
c_old = 0.7;
theta0 = [a_old;b_old;c_old];
% 隐变量z=1的后验概率的初始值
p_z  = zeros(N,1);
% 用于存储更新后的参数
a_new = 0.0;
b_new = 0.0;
c_new = 0.0;
theta = [a_new;b_new;c_new];
% 给EM算法的迭代设置个终止条件
epsilon = 0.01;
% 迭代更新参数
while norm(theta - theta0) > epsilon
% E步，计算隐变量的后验
	for i = 1:N
		p_z(i) = (a_old*b_old^y(i)*(1-b_old)^(1-y(i)))/((a_old*(b_old)^y(i)*(1-b_old)^(1-y(i))) + ((1-a_old)*c_old^y(i)*(1-c_old)^(1-y(i))));
	end
	theta0 = [a_old;b_old;c_old];
% M步，极大似然更新待估计参数 
	a_new = 1/N*sum(p_z);
	b_new = p_z'*y/sum(p_z);
	c_new = (1-p_z)'*y/sum(1-p_z);
	theta = [a_new;b_new;c_new];
	a_old = a_new;
	b_old = b_new;
	c_old = c_new;
end

本文小结

首先从为什么需要用EM算法说起：含有隐变量的参数估计问题，然后从最大化ELBO的角度推导了EM算法迭代公式的由来并证明了算法的收敛性，最后以三硬币模型为例介绍了EM算法的使用以及matlab实现。