理解EM算法

介绍EM算法的材料里，我目前看过且觉得比较好的就是NG老师的CS229讲义和李航老师的统计学习方法。
我也提不出什么新东西，就结合混合高斯分布，在这两位牛人的基础上，谈一点自己觉得看待EM算法很重要的2个必须弄清楚的问题：为什么要有EM算法，为什么叫E步和M步，还解释了一些介绍EM算法时免不了要提到的公式。如果不把这些问题和公式解释清楚真的能理解em吗？我想可能不能

1. 为什么要有EM算法

我把EM算法当做最大似然估计的拓展，解决难以给出解析解的最大似然估计（MLE）问题。
考虑高斯分布，它的最大似然估计是这样的：

θ∗=argmaxθ∑XlogL(θ|X)(1)

其中，θ=(μ,σ), θ∗=(μ∗,σ∗), logL(θ|X)=logP(X;θ)是对数似然函数，分号左边是随机变量，右边是模型参数。P(X;θ)表示X的概率值函数，它是一个以θ为参数的函数（很多人看不懂EM算法就是因为符号问题）。这里对θ求导很容易解出θ∗。

但如果这是个含有隐量Y的模型比如混合高斯模型，

P(X;θ)=∑k=1KπkN(x;μk,σk)=∑YP(Y;π)P(X|Y;μ,σ)(2)

上面假设共有K个高斯模型混合.每个高斯模型的参数为θk=(μk,σk),每个高斯模型占总模型的比重为πk。隐变量Y∈{y1,y2,...,yK}表示样本xi来自于哪一个高斯分布。分布列为：

Y y1 y2 y3 … p(y) π1 π2 π3 …

可以认为，混合高斯分布的观测值是这样产生的：先以概率πk抽取一个高斯分布yk，再以该高斯分布N(x;μk,σk)去生成观测x。其实这里的πk 就是Y的先验分布P(Y;π) (这里特地加上； π 表示P(Y)的参数是 π ,你需要求出 π 才能表示这个先验分布),而 N(x;μk,σk) 就是给定Y下的条件概率 P(X|Y;μ,σ)
这时，令θ=(μ,σ,π), θ∗=(μ∗,σ∗,π∗), 最大似然估计成了

θ∗=argmaxθ∑XlogP(X;θ)=argmaxθ∑Xlog∑YP(Y;π)P(X|Y;μ,σ)=argmaxθ∑Xlog∑YP(X,Y;θ)(3)

据群众反映，求和、取对数、再求和，这种形式求偏导较为费劲（到底有多费劲。。。其实放在混合高斯这里也不是那么费劲，有的情形远比混合高斯复杂）要是能把\log 拿到求和的最里层就好了，直接对最里层的式子求偏导岂不快哉？于是就有了EM算法

2. 为什么要分E步和M步

为了解决这个问题，有人想到了Jensen（琴生）不等式. \log 是个凹函数，以隐变量Y的任一函数f(Y)举个例子：

logE[f(Y)]=log∑YP(Y)f(Y)≥∑YP(Y)logf(Y)=E[logf(Y)](4)

根据琴生不等式的性质，当随机变量函数 f(Y) 为常数时，不等式取等号。上式中的期望换成条件期望，分布 P(Y) 换成条件分布也是一样的。

注意(3)中的联合分布P(X,Y;θ)在执行∑Y时可以把X看做是定值，此时我们可以把这个联合分布当做Y的随机变量函数（它除以P(Y)当然还是Y的随机变量函数）来考虑，并且引入一个关于Y的分布Q(Y)，具体是啥分布还不清楚,可能是给定某某的条件分布，只知道它是一个关于θ的函数：

max=maxθ∑Xlog∑YP(X,Y;θ)=maxθ∑Xlog∑YQ(Y;θ)⋅P(X,Y;θ)Q(Y;θ)=maxθ∑XlogEQ[P(X,Y;θ)P(Y;θ)]≥maxθ∑XEQ[logP(X,Y;θ)Q(Y;θ)]=maxθ∑X∑YQ(Y;θ)logP(X,Y;θ)Q(Y;θ)(5)

只有当

P(X,Y;θ)Q(Y;θ)=c(6)

式(5)才能取等号，注意到Q是Y的某一分布，有∑YQ(Y;θ)=1这个性质，因此

Q(Y;θ)=P(X,Y;θ)c=P(X,Y;θ)c⋅∑YQ(Y;θ)=P(X,Y;θ)∑Yc⋅Q(Y;θ)=P(X,Y;θ)∑YP(X,Y;θ)=P(X,Y;θ)P(X;θ)=P(Y|X;θ)(7)

所以只需要把Q取为给定X下，Y的后验分布，就能使式(5)取等号，下一步只需要最大化就行了.这时(5)为

θ∗=argmaxθ∑X∑YP(Y|X;θ)logP(X,Y;θ)P(Y|X;θ)(8)

其中：

P(X,Y;θ)=P(Y;π)P(X|Y;μ,σ)=πkN(xi;μk,σk)(9)

P(Y|X;θ)=P(X,Y;θ)∑YP(X,Y;θ)=πkN(xi;μk,σk)∑Kk=1πkN(xi;μk,σk)(10)

好吧，直接对(μ,σ,π)求导还是很麻烦，不过已经可以用迭代来最大化啦。

1）先根据式(10)，由(μ(j),σ(j),π(j))求后验概率

Q(j)=P(Y|X;θ(j))

2）再把Q(j)带入(8)中，

θ(j+1)=argmaxθ∑X∑YQ(j)logP(X,Y;θ)Q(j)=argmaxθ∑X∑Y(Q(j)logP(X,Y;θ)−Q(j)logQ(j))=argmaxθ∑X∑YQ(j)logP(X,Y;θ)(11)

就只需要最大化联合分布P(X,Y;θ)了，最大化求出(μ(j+1),σ(j+1),π(j+1))后重复这2步。

M步很显然，就是最大化那一步，E步又从何谈起呢？式(11)可以写成

θ(j+1)=argmaxθ∑X∑YQ(j)logP(X,Y;θ)=argmaxθ∑XEQ(j)[logP(X,Y;θ)]=argmaxθ∑XEY|X;θ(j)[logP(X,Y;θ)]=argmaxθ∑XEY[logP(X,Y;θ)|X;θ(j)](12)

其实，E步就是求给定X下的条件期望，也就是后验期望，使得式(5)的琴生不等式能够取等号，是对琴声不等式中,小的那一端进行放大，使其等于大的那一端，这是一次放大；M步最大化联合分布，通过0梯度，拉格朗日法等方法求极值点，又是一次放大。只要似然函数是有界的，只要M步中的0梯度点是极大值点，一直放大下去就能找到最终所求。

我的知乎回答

3.另一个EM算法的简单实现

Kmeans是一种原型聚类算法（prototype-based clustering）, 需要给数据集指定类簇的数量k, 然后找到k个聚类中心（原型），样本被归类为距离最近的原型所代表的类簇。考虑数据集X, 每个样本所属的类簇为隐变量Y, 将k个原型看做模型的参数θ=(μ1,μ2,...,μk)，……未完待续