本文摘自统计学习方法 李航著 清华大学出版社

EM算法介绍及总结

EM算法是一种迭代的算法，1977年由Dempster等人提出，用于含有隐变量（Hidden Variable）的概率模型参数的极大似然估计（不了解的可以参考我的另一篇博客，极大似然估计），或极大后验概率估计。EM算法的每次迭代由两步组成：E步，求期望（expectation）；M步，求极大（maximum）。所以这一算法称为期望极大算法（expectation maximum algorithm），简称EM算法。本文介绍EM算法，在下一篇博客中会介绍EM算法在高斯混合模型学习中的应用。

EM算法的引入

概率模型又是即含有观测变量（observe variable），又含有隐变量或潜在变量（latent variable）。如果概率模型的变量都是观测变量，那么给定的数据，可以直接用极大似然估计法，或贝叶斯估计法估计模型参数。但是，当模型含有隐变量时，就不能简单地使用这些估计方法。EM算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验估计法。我们仅讨论极大似然估计，极大后验概率与其相似。

EM算法

例子：三硬币模型

假设有3枚硬币，分别记做A，B，C。这些硬币正面出现的概率分别是π,p, 和 q。进行如下掷硬币试验：先掷硬币A，根据其结果选出硬币B或硬币C，正面选硬币B，反面选硬币C；然后掷选出的硬币，掷硬币的结果，出现正面记做1，出现反面记做0；独立地重复n次试验（这里， n=10），观测结果如下：

1, 1, 0, 1, 0, 0, 1, 0, 1, 1

假设只能观测到掷硬币的结果，不能观察到掷硬币的过程，问如何估计三枚硬币正面出现的概率，即三硬币模型的参数（π,p,q ）。

解答

三枚硬币的模型可以写作：

P(y|Θ)=∑zP(y,z|Θ)=∑zP(z|Θ)P(y|z,Θ) （贝叶斯公式，不懂的可以翻一番概率论）

P(y|Θ)=πpy(1−p)1−y+(1−π)qy(1−q)1−y

这里，随机变量y是观测变量（可以从题目中获得的数据），表示一次试验观测的结果是1或0；随机变量z是隐变量，表示未观测到的掷硬币A的结果；Θ=(π,p,q)是模型参数。这一模型是以上数据的生成模型。注意随机变量y的数据可以观测，随机变量z的数据不可观测（无法从题目中获得）。
将观测数据表示为Y=(Y1,Y2,...,Yn)T，未观测的数据表示为Z=(Z1,Z2,...,Zn)T，则观测的数据的似然函数为

P(Y|Θ)=∑zP(Z|Θ)P(Y|Z,Θ)

即

P(Y|Θ)=∏[πpyj(1−p)1−yj+(1−π)qyj(1−q)1−yj]

考虑求模型参数Θ=(π,p,q)的极大似然估计，即

Θ=argmaxΘlogP(Y|Θ)

这个问题没有解析解，只有通过迭代的方法求解。EM算法就是可以用于求解这个问题的一种迭代方法。下面给出针对以上问题的EM算法，其推到过程省略。
EM算法首先选取参数的初值，记做Θ(0)=(π0,p0,q0)，然后通过下面的步骤迭代计算参数的估计值，直到收敛为止。第i次迭代参数的估计值为Θ(i)=(πi,pi,qi)。EM算法的第i+1次迭代如下。

• E步：计算在模型参数πi,pi,qi下观测数据yj来自掷硬币B的概率
  

  
  μj(i+1)=π(i)(p(i))yi(1−p(i))1−yiπ(i)(p(i))yi(1−p(i))1−yi+(1−π(i))(q(i))yi(1−q(i))1−yi （1）

• M步：计算模型参数的新估计值：

π(i+1)=1n∑nj=1μ(i+1)j （2）

p(i+1)=∑nj=1μ(i+1)jyj∑nj=1μi+1j （3） （条件概率，在抽到硬币B的条件下翻到正面的概率）

q(i+1)=∑nj=1(1−μ(i+1)j)yj∑nj=1(1−μi+1j) （4）（条件概率，在抽到硬币C的条件下翻到正面的概率）

---

现在，如果你嫌公式太抽象，咋们来点数字计算
假设模型参数的初值取为

π(0)=0.5,p(0)=0.5,q(0)=0.5 （这样的取法比较实际，一枚硬币正面反面的概率大致相等）

由公式（1）可得，无论yj=1或者yj=0，都可以算得μ(1)j=0.5。当yj=1时，计算过程如下图，同理可计算yj=0。

利用迭代公式（2~4），得到

π(1)=0.5,p(1)=0.6,q(1)=0.6

再根据公式（1）计算μ(2)j=0.5,j=1,2,3,...,10
继续迭代，得：

π(2)=0.5,p(2)=0.6,q(2)=0.6

于是得到模型参数Θ的极大似然估计：

π=0.5,p=0.6,q=0.6

π=0.5表示硬币A是均匀的，这一结果容易理解。

计算的结果依赖于初值，如果取初值π(0)=0.4,p(0)=0.6,q(0)=0.7，那么得到的模型参数的极大似然估计是π=0.4064,p=0.5368,q=0.6432。计算出来的结果非常接近于实际值，误差仅为1.6%。这说明，EM算法与初始值的选择有关，选择不同的初值可能得到不同的参数估计。所以初值选择一定要符合实际情况，这样计算出来的结果才会接近于现实。
一般的，用Y表示观测随机变量的数据，Z表示隐随机变量的数据。Y和Z连在一起成为完全数据（complete-data），观测数据Y又称为不完全数据（incomplete-data）。假设给定观测数据Y，其概率分布是P(Y|Θ)，其中Θ是需要估计的模型参数，那么不完全数据Y的似然函数是L(Θ)=logP(Y|Θ)；假设Y和Z的联合概率分布是P(Y,Z|Θ)，那么完全数据的对数似然函数是logP(Y,Z|Θ)。
EM算法通过迭代求L(Θ)=logP(Y|Θ)的极大似然估计。每次迭代包含两步：E步，求期望；M步，求极大化。下面来介绍EM算法。

EM算法

输入： 观测变量数据Y，隐变量数据Z，联合分布P(Y,Z|Θ)，条件分布P(Z|Y,Θ)；
输出： 模型参数Θ。

（1） 选择参数的初始值Θ(0)，开始迭代。
（2） E步：记Θ(i)为第i次迭代参数Θ的估计值，在第i+1次迭代的E步，计算

Q(Θ,Θ(i))=Ez[logP(Y,Z|Θ)|Y,Θ(i)]=∑zlogP(Y,Z|Θ)P(Z|Y,Θ(i)) 公式（5）

---

注解：其中，Θ(i)代表使用其数据用来迭代下一次计算，在观测模型Y的条件下计算期望，在选定的模型参数和本次观测数据y(i+1)条件下，隐变量Z发生的概率，再乘以对应的概率，即得到期望，对应公式（1）。P(Z|Y,Θ(i))是在给定观测数据Y和当前的参数估计Θ(i)下的隐变量数据Z的条件概率分布；

---

（3）M步：求使Q(Θ,Θ(i))极大化的Θ，确定第i+1次迭代的参数估计Θ(i+1)

Θ(i+1)=argmaxΘQ(Θ,Θ(i)) 公式（6）

（4）重复第（2）步和第（3）步，直到收敛。
其中，公式（5）是EM算法的核心，称为Q函数（Q function）。

Q函数

完全数据的对数似然函数logP(Y,Z|Θ)关于在给定的观测数据Y和当前参数Θ(i)下对未观测数据Z的条件概率分布P(Z|Y,Θ(i))的期望称为Q函数，即

Q(Θ,Θ(i))=Ez[logP(Y,Z|Θ)|Y,Θ(i)] 公式（7）

---

关于EM算法的几点说明：

1. 步骤（1）参数的初值可以任意选择，但需要注意EM算法对初值是敏感的。
2. 步骤（2）E步求Q(Θ,Θ(i))。Q函数式中Z是未观测数据，Y是观测数据。注意，Q(Θ,Θ(i))的第一个变元表示要极大化的参数，第二个变元表示参数的当前估计值。每次迭代实际在求Q函数及其极大。
3. 步骤（3）M步求Q(Θ,Θ(i))的极大化，得到Θ(i+1)，完成一次迭代Θ(i)⇒Θ(i+1)。后面证明每次迭代是似然函数增大达到局部极值。
4. 步骤（4）给出停止迭代条件，一般是对较小的正数Σ1,Σ2，若满足
   
   ∥∥Θ(i+1)−Θ(i)∥∥<Σ1 或者 ∥∥Q(Θ(i+1),Θ(i))−Q(Θ(i),Θ(i−1))∥∥

则停止迭代。

---