EM算法(Expectation Maximization Algorithm)
EM算法(Expectation Maximization Algorithm)
1. 前言
这是本人写的第一篇博客(2013年4月5日发在cnblogs上,现在迁移过来)，是学习李航老师的《统计学习方法》书以及斯坦福机器学习课Andrew Ng的EM算法课后，对EM算法学习的介绍性笔记，如有写得不恰当或错误的地方，请指出，并多多包涵，谢谢。另外本人数学功底不是很好，有些数学公式我会说明的仔细点的，如果数学基础好，可直接略过。
2.基础数学知识
在正式介绍EM算法之前，先介绍推导EM算法用到的数学基础知识，包括凸函数，Jensen不等式。
2.1.凸函数
对于凸函数，凹函数，如果大家学过高等数学，都应该知道，需要注意的是国内教材如同济大学的《高等数学》的这两个概念跟国外刚好相反，为了能更好的区别，本文章把凹凸函数称之为上凸函数，下凸函数，具体定义如下：
上凸函数：函数 f(x) 满足对定义域上任意两个数 a , b 都有 f[(a+b)/2]≥[f(a)+f(b)]/2
下凸函数：函数 f(x) 满足对定义域上任意两个数 a , b 都有 f[(a+b)/2]≤[f(a)+f(b)]/2
更直观的可以看图2.1和2.2：
 |  |
|---|---|
| 图2.1. 上凸函数 | 图2.2. 下凸函数 |
可以清楚地看到图2.1上凸函数中，f[(a+b)/2]≥[f(a)+f(b)]/2，而且不难发现，如果f(x)是上凸函数，那么 −f(x) 是下凸函数。
当 a≠b 时，f[(a+b)/2]>[f(a)+f(b)]/2 成立，那么称 f(x) 为严格的上凸函数，等号成立的条件当且仅当 a=b,下凸函数与其类似。
2.2.Jensen不等式
有了上述凸函数的定义后，我们就能很清楚的Jensen不等式的含义，它的定义如下：
如果f是上凸函数，X 是随机变量，那么 f(E[X])≥E[f(X)]
特别地，如果f是严格上凸函数，那么 E[f(X)]=f(E[X]) 当且仅当 p(X=E[X])=1，也就是说 X 是常量。
那么很明显 logx 函数是上凸函数，可以利用这个性质。
有了上述的数学基础知识后，我们就可以看具体的EM算法了。
3.EM算法所解决问题的例子
在推导EM算法之前，先引用《统计学习方法》中EM算法的例子：
例1. (三硬币模型)假设有3枚硬币,分别记作 A,B,C 。这些硬币正面出现的概率分别为 π，p 和 q。投币实验如下，先投 A，如果 A 是正面，即 A=1，那么选择投 B；A=0，投 C。最后，如果 B 或者 C 是正面，那么 y=1；是反面，那么 y=0；独立重复 n 次试验 (n=10)，观测结果如下： 1,1,0,1,0,0,1,0,1,1 假设只能观测到投掷硬币的结果，不能观测投掷硬币的过程。问如何估计三硬币正面出现的概率，即 π，p 和 q 的值。
解：设随机变量 y 是观测变量，则投掷一次的概率模型为:
P(y|θ)=πpy(1−p)1−y+(1−π)qy(1−q)1−y有 n 次观测数据 Y，那么观测数据 Y 的似然函数为:P(Y|θ)=∏nj=1[πpyj(1−p)1−yj+(1−π)qyj(1−q)1−yj]那么利用最大似然估计求解模型解，即
θˆ=argmaxθlogP(Y|θ)=argmaxθ∑j=110logP(yj|θ)=argmaxθ∑j=110log[πpyj(1−p)1−yj+(1−π)qyj(1−q)1−yj](1)(2)(3)这里将概率模型公式和似然函数代入 (1) 式中，可以很轻松地推出 (1)⇒(2)⇒(3)，然后选取 θ(π,p,q)，使得 (3) 式值最大，即最大似然。然后，我们会发现因为 (3) 中右边多项式 + 符号的存在，使得 (3) 直接求偏导等于 θ 或者用梯度下降法都很难求得 θ值。
这部分的难点是因为 (3) 多项式中 + 符号的存在，而这是因为这个三硬币模型中，我们无法得知最后得结果是硬币 B 还是硬币 C抛出的这个隐藏参数。那么我们把这个latent 随机变量加入到 log-likelihood 函数中，得
l(θ)=∑j=110log∑i=12P(yj,zi∣θ)=∑j=110log∑i=12Qj(zi)P(yj,zi∣θ)Qj(zi)≥∑j=110∑i=12Qj(zi)logP(yj,zi∣θ)Qj(zi)(4)(5)(6)
略看一下，好像很复杂，其实很简单，首先是公式 (4)，这里将 zi 做为隐藏变量，当 z1 为结果由硬币 B 抛出，z2 为结果由硬币C抛出，不难发现:∑i=12P(yj,zi∣θ)=P(yj∣θ)=πpyj(1−p)1−yj+πqyj(1−q)1−yj接下来公式说明 (4)⇒(5) (其中 (5) 中 Q(z) 表示的是关于 z 的某种分布，∑iQj(zi)=1，很直接，在 P 的分子分母同乘以 Q(zi)。最后是 (5)⇒(6)，到了这里终于用到了第二节介绍的Jensen不等式，数学好的人可以很快发现，∑2i=1Qj(zi)P(yj,zi|θ)Qj(zi) 就是 [P(yj,zi|θ)Qj(zi)] 的期望值（如果不懂，可google期望公式并理解），且log是上凸函数，所以就可以利用Jensen不等式得出这个结论。因为我们要让log似然函数 l(θ)最大，那么这里就要使等号成立。根据Jensen不等式可得，要使等号成立，则要使 P(yj,zi|θ)Qj(zi)=c 成立。
再因为∑iQj(zi)=1，所以得∑iP(yj,zi|θ)=c，c 为常数，那么(这里感谢网友@无影随想指出错误)Q(zi)=P(yj,zi∣θ)/∑iP(yj,zi∣θ)=P(yj,zi)/P(yj∣θ)=P(zi∣yj,θ)这里可以发现Qj(z1)=πpyj(1−p)1−yjπpyj(1−p)1−yj+(1−π)qyj(1−q)1−yjQj(z2)=(1−π)qyj(1−q)1−yjπpyj(1−p)1−yj+(1−π)qyj(1−q)1−yjOK,到这里，可以发现公式 (6) 中右边多项式已经不含有“+”符号了，如果知道 Q(z) 的所有值，那么可以容易地进行最大似然估计计算，但是 Q 的计算需要知道 θ 的值。这样的话，我们是不是可以先对θ进行人为的初始化 θ0，然后计算出 Q 的所有值 Q1(在 θ0 固定的情况下，可在 Q1 取到公式 (6) 的极大值)，然后在对公式 (6) 最大似然估计，得出新的 θ1 值（在固定Q1的情况下，取到公式 (6) 的极大值），这样又可以计算新的 Q 值 Q1,然后依次迭代下去。答案当然是可以。因为 Q1 是在 θ0 的情况下产生的，可以调节公式 (6) 中 θ 值，使公式 (6) 的值再次变大，而 θ 值变了之后有需要调节 Q 使 (6) 等号成立，结果又变大，直到收敛（单调有界必收敛），如果到现在还不是很清楚，具体清晰更广义的证明可以见下部分EM算法说明。
另外对公式 (6) 进行求偏导等于 θ，求最大值，大家可以自己练习试试，应该很简单的，这里不做过多陈述。
在《统计学习方法》书中，进行两组具体值的计算
- π0=0.5, p0=0.5, q0=0.5，迭代结果为 π=0.5, p=0.6, q=0.5
- π0=0.4, p0=0.6, q0=0.7，迭代结果为 π=0.4064, p=0.5368, q=0.6432
两组值的最后结果不相同，这说明EM算法对初始值敏感，选择不同的初值可能会有不同的结果，只能保证参数估计收敛到稳定点。因此实际应用中常用的办法就是选取多组初始值进行迭代计算，然后取结果最好的值。
在进行下部分内容之前，还需说明下一个东西。在上面的举例说明后，其实可以发现上述的解决方法跟一个简单的聚类方法很像，没错，它就是K-means聚类。K-means算法先假定k个中心，然后进行最短距离聚类，之后根据聚类结果重新计算各个聚类的中心点，一次迭代，是不是很像，而且K-means也是初始值敏感，因此其实K-means算法也包含了EM算法思想，只是这边EM算法中用P概率计算，而K-means直接用最短距离计算。所以EM算法可以用于无监督学习。在下一篇文章，我准备写下典型的用EM算法的例子，高斯混合模型(GMM,Gaussian Mixture Model)。
4.EM算法
4.1.模型说明
考虑一个参数估计问题，现有 y1,y2,…,yn 共 n 个训练样本，需有多个参数 π 去拟合数据，那么这个 log 似然函数是：
可能因为 θ 中多个参数的某种关系（如上述例子中以及高斯混合模型中的3类参数），导致上面的 log 似然函数无法直接或者用梯度下降法求出最大值时的 θ 值，那么这时我们需要加入一个隐藏变量 z，以达到简化 l(θ),迭代求解 l(θ) 极大似然估计的目的。
4.2.EM算法推导
这小节会对EM算法进行具体推导，许多跟上面例子的解法推导是相同的，如果已经懂了，可以加速阅读。首先跟“三硬币模型”一样，加入隐变量 z 后，假设 Q(z) 是关于隐变量 z 的某种分布，那么有如下公式：
公式 (7) 是加入隐变量，(7)⇒(8) 是在基础上分子分母同乘以，(8)⇒(9) 用到Jensen不等式（跟“三硬币模型”一样），等号成立的条件是，c 是常数。再因为，则有如下 Q 的推导：
再一次重复说明，要使 (9) 等式成立，则 Qj(zi) 为 yj, z的后验概率。算出 Qj(zi) 后对 (9) 就可以进行求偏导，以剃度下降法求得 θ 值，那么又可以计算新 Qj(zi) 的值，依次迭代，EM算法就实现了。
EM 算法(1)
选取初始值 θ0 初始化 θ，t=0
Repeat {
E步：Qtj(zi)=P(yj,zi∣θt)/∑iP(yj,zi∣θt)=P(yj,zi∣θt)/P(yj∣θt)=P(zi∣yj,θt)M步：
θt+1t=argmaxθ∑j=1n∑iQtj(zi)logP(yj,zi∣θ)Qtj(zi)=t+1}直到收敛
4.3.EM算法收敛性证明
当 θ 取到 θt 值时，求得
那么可得如下不等式：