EM算法：数学推导+实例演示

EM(Expectation Maximum)算法即期望最大化算法，是一种对不完全数据(因数据缺失或有未被观测等含有隐变量的数据)估计未知变量的迭代算法。

在隐变量这篇文章里用一个例子解释了什么是隐变量，本文会在此例基础上进行扩展来引出EM算法是什么、能解决什么问题以及是如何解决问题的。

EM算法实例演示

EM算法本质上是通过不断迭代去求出模型参数的算法。下面先通过一个具体例子来演示EM算法的流程，然后再从数学角度进行推理验证。

抛硬币游戏规则

有 $A、B、C$ 三枚硬币，记正面向上的概率分别为 $\pi$ 、$p$、$q$ 且均未知，记为模型参数 $\theta =(\pi ,p,q)$ 。现在想通过抛硬币实验来估计出模型参数。

先抛硬币 $A$ ，如果硬币 $A$ 正面向上，则选择硬币 $B$ 连续抛10次，记录这一轮10次的结果；如果硬币 $A$ 反面向上，则选择硬币 $C$ 连续抛10次，记录这一轮10次的结果。

表1是某次抛5轮的结果，其中H表示正面(Head)向上；T表示反面(Tail)向上。

对于这组观测数据，用极大似然估计就能求出 $p$ 和 $q$ 。

轮次	硬币A的结果(已知)	硬币B的结果(观测数据-24正6反)	硬币C的结果(观测数据-9正11反)
1	T		HTTTHHTHTH (5正5反)
2	H	HHHHTHHHHH (9正1反)
3	H	HTHHHHHTHH (8正2反)
4	T		HTHTTTHHTT (4正6反)
5	H	THHHTHHHTH (7正3反)

表1

极大似然估计

假设一枚硬币正面向上的概率为 $\theta$ ，在抛 $N$ 次硬币中有 $k$ 次正面向上且顺序已知的情况下对应的概率或似然函数可表示为：

$L(\theta )={\theta }^k(1-\theta {)}^{N-k}$

具体地，比如对于第二轮的结果 HTTTHHTHTH，$L(\theta )=P(HTTTHHTHTH|B)={\theta }_B^9(1-{\theta }_B{)}^1$

对数似然函数为：

$lnL(\theta )=ln{\theta }^k(1-\theta {)}^{N-k}=ln{\theta }^k+ln(1-\theta {)}^{N-k}=kln\theta +(N-k)ln(1-\theta )$

求能使得上式最大的 $\theta$ ，即通过最大似然估计出 $\theta$ 。为此，对 $\theta$ 求偏导：

$\frac{\partial L(\theta )}{\partial \theta }=\frac{\partial (kln\theta +(N-k)ln(1-\theta ))}{\partial \theta }=\frac{k}{\theta }-\frac{(N-k)}{1-\theta }$

$lnx$ 的导数为 $\frac{1}{x}$

$ln(1-x)$ 的导数为 $\frac{1}{(x-1)}$

令上式为0，则：

$\theta =\frac{k}{N}=\frac{硬币正面向上的次数}{硬币抛出的总次数}$

由上述结论可得 ：

$p=\frac{24}{24+6}=0.80$

$q=\frac{9}{9+11}=0.45$

也就是说，我们只要知道硬币 $B、C$ 正面向上次数，就可以估算出 $p$ 和 $q$ 。

需要注意的是，在上面的抛硬币游戏中，硬币 $A$ 的结果是已知的，也就是说不存在隐变量，$p、q$ 可以在一组完全数据(complete-data)下通过极大似然求得。

隐变量困境

如果如表2所示， 硬币 $A$ 的结果是未知的，这种情况下如何求 $p$ 和 $q$ ？

轮次	硬币A的结果(未知)	硬币B或C的结果(观测数据)
1	未知	HTTTHHTHTH (5正5反)
2	未知	HHHHTHHHHH (9正1反)
3	未知	HTHHHHHTHH (8正2反)
4	未知	HTHTTTHHTT (4正6反)
5	未知	THHHTHHHTH (7正3反)

表2

这里未知的硬币 $A$ 的结果，就是隐变量，这是一组包含了隐变量的不完全数据(incomplete-data)。

目前并不知道这5轮观测结果是来自硬币 $B$ 还是 $C$ ，因此不能通过统计正面向上的次数来估计 $p$ 和 $q$ 。

那么可以先用观测数据来推断每轮观测数据来自哪枚硬币吗？

观测结果的概率表示

实际上，表1每轮的结果可以用 来自于硬币 $B$ 或 $C$ 的概率 表示，由于硬币 $A$ 结果已知，所以概率比较绝对，非0即1：

轮次	硬币A的结果(已知)	硬币B或C的结果(观测数据)	来自硬币B的概率	来自硬币C的概率	硬币B正面向上的次数	硬币C正面向上的次数
1	T	HTTTHHTHTH (5正5反)	0	1	0*5=0	1*5=5
2	H	HHHHTHHHHH (9正1反)	1	0	1*9=9	0*9=0
3	H	HTHHHHHTHH (8正2反)	1	0	1*8=8	0*8=0
4	T	HTHTTTHHTT (4正6反)	0	1	0*4=0	1*4=4
5	H	THHHTHHHTH (7正3反)	1	0	1*7=7	0*7=0

表3

同理，也可以用表2已知的观测数据推断每轮抛硬币的结果来自硬币 $B$ 或 $C$ 的概率，然后就可以计算出硬币 $B$ 和 $C$ 正面向上的次数，继而求出 $p$ 和 $q$ 。

想象很丰满，现实很骨感。每轮的结果来自硬币 $B$ 或 $C$ 的概率怎么求呢？

可以借助贝叶斯定理！

贝叶斯定理公式：

$P(X|Y)=\frac{P(Y|X)P(X)}{P(Y)}$

以第1轮为例，要求第1轮抛硬币的结果来自硬币 $B$ 的概率，从概率的角度来看就是求 $P(B|HTTTHHTHTH)$ 。

由贝叶斯概率公式得： $P(B|HTTTHHTHTH)=\frac{P(HTTTHHTHTH|B)P(B)}{P(HTTTHHTHTH)}=\frac{P(HTTTHHTHTH|B)P(B)}{P(HTTTHHTHTH|B)P(B)+P(HTTTHHTHTH|C)P(C)}=\frac{\pi p^5(1-p{)}^5}{\pi p^5(1-p{)}^5+(1-\pi )q^5(1-q{)}^5}(1)$

$P(B)、P(C)$ 表示自然选取硬币的概率；

$P(HTTTHHTHTH)$ 表示由 $B$ 和 $C$ 划分空间的全概率，可以展开；

$P(HTTTHHTHTH|B)$ 表示在选取硬币 $B$ 的情况下抛出 $HTTTHHTHTH$ 这种结果的概率：

$P(HTTTHHTHTH|B)=p^5(1-p{)}^5$

同理可得

$P(HTTTHHTHTH|C)=q^5(1-q{)}^5$

EM算法流程

从式 $(1)$ 可以看出，想要计算 $P(B|HTTTHHTHTH)$，需要先知道 $p、q$ ，而要得到 $p、q$，又需要先知道 $P(B|HTTTHHTHTH)$。

这看起来是个矛盾的问题，前面我们说过，EM算法本质上是通过不断迭代求模型参数的算法。这个矛盾的问题就可以用迭代来解决。

事实上，迭代开始时，可以先随机初始化 $p、q$ 的值，记作 $p^i、q^i$ (第1次迭代的 $i=0$ )，由式 $(1)$ 可以求出 $P(B|HTTTHHTHTH)、P(C|HTTTHHTHTH)$，然后可以结合5轮观测数据计算出硬币 $B、C$ 正面向上的次数，继而利用极大似然估计得到一组新的 $p^{i+1}、q^{i+1}$ 。如果 $p^{i+1}、q^{i+1}$ 与 $p^i、q^i$ 一样，说明初始化的值是比较靠谱的；如果还有差距，那么使用估计出的 $p^{i+1}、q^{i+1}$ 替代 $p^i、q^i$ 进行下一次迭代，直到收敛。

问题是每次估计出的 $p^{i+1}、q^{i+1}$ 一定会越来越接近真实的 $p、q$ 吗？事实上确实可以，后面会用数学证明，暂且按下不表。

假设我们已经知道了真实的 $p、q$ ，即上面计算的：

$p=\frac{24}{24+6}=0.80$

$q=\frac{9}{9+11}=0.45$

现在来按以上迭代流程，看看是否能求出 $p、q$ 。

① 初始化概率

令：

${\pi }^0=0.5$

$p^0=0.6$

$q^0=0.5$

② 观测结果的归属

代入式 $(1)$ 可得第1轮抛硬币的结果来自硬币 $B$ 的概率：

$P(B|HTTTHHTHTH)=\frac{p^5(1-p{)}^5}{p^5(1-p{)}^5+q^5(1-q{)}^5}=\frac{0.6^5(1-0.6{)}^5}{0.6^5(1-0.6{)}^5+0.5^5(1-0.5{)}^5}\approx 0.45$

第1轮抛硬币的结果来自硬币 $C$ 的概率：

$P(C|HTTTHHTHTH)=1-P(B|HTTTHHTHTH)\approx 0.55$

③ 期望最大化

依次计算第2-5轮，并将表2整理如下：

轮次	硬币A的结果(未知)	硬币B或C的结果(观测数据)	来自硬币B的概率$(p^0=0.6、q^0=0.5)$	来自硬币C的概率$(p^0=0.6、q^0=0.5)$	硬币B正面(H)、反面(T)向上的次数	硬币C正面(H)、反面(T)向上的次数
1	未知	HTTTHHTHTH (5正5反)	0.45	0.55	0.45*5=2.2 H 0.45*5=2.2 T	0.55*5=2.8 H 0.55*5=2.8 T
2	未知	HHHHTHHHHH (9正1反)	0.80	0.20	0.80*9=7.2 H 0.80*1=0.8 T	0.20*9=1.8 H 0.20*1=0.2 T
3	未知	HTHHHHHTHH (8正2反)	0.73	0.27	0.73*8=5.8 H 0.73*2=1.5 T	0.27*8=2.2 H 0.27*2=0.5 T
4	未知	HTHTTTHHTT (4正6反)	0.35	0.65	0.35*4=1.4 H 0.35*6=2.1 T	0.65*4=2.6 H 0.65*6=3.9 T
5	未知	THHHTHHHTH (7正3反)	0.65	0.35	0.65*7=4.5 H 0.65*3=1.9 T	0.35*7=2.5 H 0.35*3=1.1 T

表4

有了硬币 $B$ 和 $C$ 正面向上的次数，可得：

$p^1=\frac{硬币正面向上的次数}{硬币抛出的总次数}=\frac{2.2+7.2+5.8+1.4+4.5}{(2.2+7.2+5.8+1.4+4.5)+(2.2+0.8+1.5+2.1+1.9)}=\frac{21.3}{21.3+8.6}\approx 0.71$

$q^1=\frac{硬币正面向上的次数}{硬币抛出的总次数}=\frac{2.8+1.8+2.2+2.6+2.5}{(2.8+1.8+2.2+2.6+2.5)+(2.8+0.2+0.5+3.9+1.1)}=\frac{11.7}{11.7+8.4}\approx 0.58$

④ 迭代

以本次迭代中得到的 $p^1、q^1$ 替代 $p^0、q^0$ ，重复前三步可在第2次迭代后得到 $p^2、q^2$ ，以此类推，经过10次迭代后:

$p^{10}=0.80$

$q^{10}=0.45$

数次迭代后得到了与真实值相同的参数，以上几个步骤就是EM算法的推演过程，

EM算法实例的数学解释

假设还是以上面的抛硬币游戏规则进行了 $n$ 次抛硬币实验，其中观测结果记为 $Y=(y_1,y_2,...,y_j,...,y_n)$ ， $y_j$ 表示其中某次抛硬币的结果(正面或反面向上)，记正面向上时 $y_j=1$ 、反面向上时 $y_j=0$ 。其中未观测的隐变量记为 $Z=(z_1,z_2,...,z_j,...,z_n)$ ， $z_j$ 表示其中某次抛硬币的隐变量，即某次抛硬币 $A$ 是正面向上还是反面向上。

$Y$ 中出现 $y_j$ 的可能情况如下表所示：

硬币A的结果	硬币B或C的结果	概率
硬币A为 H (概率为 $\pi$)	硬币B为 H (概率为 $p$)	$\pi p$
硬币A为 H (概率为 $\pi$)	硬币B为 T (概率为 $1-p$)	$\pi (1-p)$
硬币A为 T (概率为 $1-\pi$)	硬币C为 H (概率为 $q$)	$(1-\pi )q$
硬币A为 T (概率为 $1-\pi$)	硬币C为 T (概率为 $1-q$)	$(1-\pi )(1-q)$

表5

其生成模型的概率表示为：

$P(y_j|\theta )=\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}(2)$

$y_j$ 来自硬币 $B$ 的概率记为 ${\mu }_j$ ：

${\mu }_j=\frac{\pi p^{y_j}(1-p{)}^{1-y_j}}{\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}}(3)$

由式 $(2)$ 可知这 $n$ 次抛硬币结果的对数似然函数为：
$$\begin{array}{rl}L(\theta )& =logP(Y|\theta )\\ & =log\sum _{z}P(Y,Z|\theta )\\ & =log\sum _{z}P(Z|\theta )P(Y|Z,\theta )\\ & =log\prod _{j_{=}1}^n[\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}]\\ & =\sum _{j_{=}1}^{n}log[\pi p^{y_j}(1-p{)}^{1-y_j}+(1-\pi )q^{y_j}(1-q{)}^{1-y_j}](4)\end{array}$$

公式解释：

Z：隐变量数据，对应抛硬币中硬币 $A$ 正面向上还是反面向上；

Y：观测变量数据，对应每轮次抛硬币中硬币 $B、C$ 的正反面结果；

θ：模型参数，对应抛硬币中硬币 $B、C$ 正面向上的概率 ${\theta }_B$、${\theta }_C$ ；

P(Y|θ)：已知模型参数情况下出现已知观测值的概率，即似然函数，能使之达到最大的参数即为待估模型参数；

P(Z|Y,θ)：已知模型参数和观测值情况下隐变量的概率，即硬币 $A$ 正面/反面向上的概率；

P(Y,Z|θ)：已知模型参数情况下观测值与隐变量的联合概率(两个事件的联合概率 $P(A,B)$ 是指两个事件同时发生的概率，如果两个事件 $P(A)、P(B)$ 独立，则$P(A,B)=P(A)P(B)$ )。

求式 $(4)$ 的极大似然得到的参数即为要求的参数：

$\arg \underset{\theta }{\max }log{\sum }_{z}P(Y,Z|\theta )(5)$

由于存在隐变量与和的对数，不能直接求解式 $(5)$ ，回想一下，在上面实例演示中，我们先得到了观测数据来自硬币 $B$ 或 $C$ 的概率，然后利用观测数据估计出的 $p、q$ 。那么是否也可以将式 $(5)$ 中的似然函数改造成包含 观测数据来自硬币 $B$ 或 $C$ 的概率 并且一样通过迭代来求解模型的参数呢？

为此，将 $log{\sum }_{z}P(Y,Z|\theta )$ 做个变换：

$log{\sum }_{z}P(Y,Z|\theta )=log(\frac{{\sum }_{z}P(Y,Z|\theta )}{P(Z|Y,{\theta }^i)}P(Z|Y,{\theta }^i))(6)$

其中 $P(Z|Y,{\theta }^i)$ 表示观测结果来自硬币B/C的概率，参数 ${\theta }^i$ 已知(第1次迭代时是随机初始化的) ，且满足 ${\sum }_{j=1}^{n}P(Z_j|Y,{\theta }^i)=1$ 、$P(Z_j|Y,{\theta }^i)\ge 0$ ，可以看作一个分布函数。

再利用 Jensen不等式 ，令 $f(x)=log(x)、{\lambda }_j=P(Z_j|Y,{\theta }^i)、x_j=\frac{P(Y,Z|\theta )}{P(Z_j|Y,{\theta }^i)}$ ，可得：
$$\begin{array}{rl}L(\theta )& =log\sum _{z}P(Y,Z|\theta )\\ & =log\sum _z\frac{P(Y,Z|\theta )}{P(Z|Y,{\theta }^i)}P(Z|Y,{\theta }^i)\\ & \ge \sum _{z}P(Z|Y,{\theta }^i)log\frac{P(Y,Z|\theta )}{P(Z|Y,{\theta }^i)}=\sum _{z}P(Z|Y,{\theta }^i)logP(Y,Z|\theta )-\sum _{z}P(Z|Y,{\theta }^i)logP(Z|Y,{\theta }^i)(7)\end{array}$$

Jensen不等式：

若函数 $f(x)$ 是定义在区间 [a,b] 上的 上凸函数(比如ln函数) ，则对任意的 $x_1,x_2,...,x_n\in [a,b]$，有不等式：

$f({\sum }_{j=1}^n{\lambda }_j{x}_j)\ge {\sum }_{j=1}^n{\lambda }_{j}f(x_j)(8)$

其中，${\sum }_{j=1}^n{\lambda }_j=1$ 且 ${\lambda }_1,{\lambda }_2,...,{\lambda }_n\ge 0$ ，当且仅当 $x_1=x_2=...=x_n$时等号成立。

由期望的定义，式 $(8)$ 相当于：

$f(E(x))\ge E(f(x))$

因为 $P(Z|Y,{\theta }^i)$ 已知，因此式 $(7)$ 可进一步化简为：

$log{\sum }_{z}P(Y,Z|\theta )\ge {\sum }_{z}P(Z|Y,{\theta }^i)logP(Y,Z|\theta )(9)$

式 ${\sum }_{z}P(Z|Y,{\theta }^i)logP(Y,Z|\theta )$ 就是EM算法中E步的期望公式，然后再对其最大化得到模型参数，这是EM算法里的M步。

因此式 $(4)$ 可以改写为：

$L(\theta )=logP(Y|\theta )=log{\sum }_{z}P(Y,Z|\theta )\ge {\sum }_{j_{=}1}^n[{\mu }_j^{i+1}log[p^{y_j}(1-p{)}^{1-y_j}]+(1-{\mu }_j^{i+1})log[q^{y_j}(1-q{)}^{1-y_j}]](10)$

其中 ${\mu }_j^{i+1}$ 为在模型参数 ${\theta }^i=({\pi }^i,p^i,q^i)$ 时观测数据来自硬币 $B$ 的概率。

将式 $(3)$ 代入式 $(10)$ 中并求偏导，可得：

${\pi }^{i+1}=\frac{1}{n}{\sum }_{j=1}^n{\mu }_j^{i+1}$

$p^{i+1}=\frac{{\sum }_{j=1}^n{\mu }_j^{i+1}{y}_i}{{\sum }_{j=1}^n{\mu }_j^{i+1}}$

$q^{i+1}=\frac{{\sum }_{j=1}^n(1-{\mu }_j{)}^{i+1}{y}_i}{{\sum }_{j=1}^n(1-{\mu }_j{)}^{i+1}}$

看情况继续迭代直到收敛…

从式 $(9)$ 可以看出，因为 $L(\theta )=log{\sum }_{z}P(Y,Z|\theta )$ 不好直接求解，因此先构建了一个下界(E步) ，然后再优化这个下界(M步)，逐步迭代最终得到模型参数。

那么每次迭代后都会使 $L(\theta )$ 比上次迭代时 $L({\theta }^i)$ 的大吗？下面我们从数学上来论证。

EM算法的推导

开始时，我们想通过极大似然估计来求解 $L(\theta )=log{\sum }_{z}P(Y,Z|\theta )$ ，然而隐变量的存在使其行不通。

那么我们转变一下思路，通过迭代来逐步逼近 $L(\theta )$ 。

假设第 $i$ 次迭代后估计出的模型参数为 ${\theta }^i$ ，对应似然函数 $L({\theta }^i)=logP(Y|{\theta }^i)$ 。我们希望下一次估计出的模型参数 $\theta$ 能使 $L(\theta )$ 相比上次迭代的更大，即 $L(\theta )-L({\theta }^i)\ge 0$ 。

现在来推导一下(同样利用Jensen不等式)。
$$\begin{array}{rl}L(\theta )-L({\theta }^i)& =log\sum _{z}P(Y,Z|\theta )-logP(Y|{\theta }^i)\\ & =log\sum _z\frac{P(Y,Z|\theta )}{P(Z|Y,{\theta }^i)}P(Z|Y,{\theta }^i)-logP(Y|{\theta }^i)\\ & \ge \sum _{z}P(Z|Y,{\theta }^i)log\frac{P(Y,Z|\theta )}{P(Z|Y,{\theta }^i)}-\sum _{z}P(Z|Y,{\theta }^i)logP(Y|{\theta }^i)\\ & =\sum _{z}P(Z|Y,{\theta }^i)log\frac{P(Y,Z|\theta )}{P(Z|Y,{\theta }^i)logP(Y|{\theta }^i)}\\ & =\sum _{z}P(Z|Y,{\theta }^i)log\frac{P(Y,Z|\theta )}{P(Y,Z|{\theta }^i)}(11)\end{array}$$

上式中，因为 ${\sum }_{z}P(Z|Y,{\theta }^i)=1$ ，故 $logP(Y|{\theta }^i)={\sum }_{z}P(Z|Y,{\theta }^i)logP(Y|{\theta }^i)$ ；

$P(Y,Z|{\theta }^i)=P(Z|Y,{\theta }^i)logP(Y|{\theta }^i)$

由式 $(11)$ 可得：

$L(\theta )\ge L({\theta }^i)+{\sum }_{z}P(Z|Y,{\theta }^i)log\frac{P(Y,Z|\theta )}{P(Y,Z|{\theta }^i)}(12)$

令 $B(\theta ,{\theta }^i)=L({\theta }^i)+{\sum }_{z}P(Z|Y,{\theta }^i)log\frac{P(Y,Z|\theta )}{P(Y,Z|{\theta }^i)}(13)$

则 $L(\theta )\ge B(\theta ,{\theta }^i)$

函数 $B(\theta ,{\theta }^i)$ 在 $L(\theta )$ 的下方。

对于函数 $B(\theta ,{\theta }^i)$ ，当 $\theta$ 取值 ${\theta }^i$ 时：

$B({\theta }^i,{\theta }^i)=L({\theta }^i)+{\sum }_{z}P(Z|Y,{\theta }^i)log\frac{P(Y,Z|{\theta }^i)}{P(Y,Z|{\theta }^i)}=L({\theta }^i)+{\sum }_{z}P(Z|Y,{\theta }^i)log1=L({\theta }^i)\le L(\theta )(14)$

从下图可以看出，函数 $B(\theta ,{\theta }^i)$ 是 $L(\theta )$ 的下界。求下界 $B(\theta ,{\theta }^i)$ 的极大值就可以逼近似然函数的极大值。假设 ${\theta }^{i+1}$ 可以使 $B(\theta ,{\theta }^i)$ 取得极大值，即：
$$\begin{array}{rl}{\theta }^{i+1}& =\arg \underset{\theta }{\max }B(\theta ,{\theta }^i)\\ & =\arg \underset{\theta }{\max }[L({\theta }^i)+\sum _{z}P(Z|Y,{\theta }^i)log\frac{P(Y,Z|\theta )}{P(Y,Z|{\theta }^i)}]\\ & =\arg \underset{\theta }{\max }[\sum _{z}P(Z|Y,{\theta }^i)logP(Y,Z|\theta )-\sum _{z}P(Z|Y,{\theta }^i)P(Y,Z|{\theta }^i)+L({\theta }^i)]\\ & =\arg \underset{\theta }{\max }[\sum _{z}P(Z|Y,{\theta }^i)logP(Y,Z|\theta )](15)\end{array}$$

${\sum }_{z}P(Z|Y,{\theta }^i)P(Y,Z|{\theta }^i)、L({\theta }^i)$ 均为已知，对求极大值没有影响，故可以省略。

仔细看式 $(15)$ ，是不是有点面熟，是的，我们在上面式 $(7)$ 处已经推导过一次了。

式 $(15)$ 中的 ${\sum }_{z}P(Z|Y,{\theta }^i)logP(Y,Z|\theta )$ 就是EM算法E步得到的期望，对之最大化就是EM算法的一次迭代。

以上就是EM算法的推导过程。