前言
EM算法是机器学习十大算法之一，它很简单，但是也同样很有深度，简单是因为它就分两步求解问题，
- E步：求期望（expectation）
- M步：求极大（maximization)
深度在于它的数学推理涉及到比较繁杂的概率公式等，所以本文会介绍很多概率方面的知识，不懂的同学可以先去了解一些知识，当然本文也会尽可能的讲解清楚这些知识，讲的不好的地方麻烦大家评论指出，后续不断改进完善。
EM算法引入
概率模型有时候既含有观测变量，又含有隐变量或潜在变量，如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计方法估计模型参数，但是当模型含有隐变量时，就不能简单的使用这些方法，EM算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验概率估计法，我们讨论极大似然估计，极大后验概率估计与其类似。 参考统计学习方法书中的一个例子来引入EM算法， 假设有3枚硬币，分别记做A、B、C，这些硬币正面出现的概率分别是𝜋π、𝑝p、𝑞q，进行如下实验：
-
先掷硬币A，根据结果选出硬币B和硬币C，正面选硬币B，反面选硬币C
-
通过选择出的硬币，掷硬币的结果出现正面为1，反面为0 如此独立地重复n次实验，我们当前规定n=10，则10次的结果如下所示：
1,1,0,1,0,0,1,0,1,11,1,0,1,0,0,1,0,1,1假设只通过观测到掷硬币的结果，不能观测掷硬币的过程，问如何估计三个硬币出现正面的概率？ 我们来构建这样一个三硬币模型：𝑃(𝑦|𝜃)=∑𝑧𝑃(𝑦,𝑧|𝜃)=∑𝑧𝑃(𝑧|𝜃)𝑃(𝑦|𝑧,𝜃)=𝜋𝑝𝑦(1−𝑝)1−𝑦+(1−𝜋)𝑞𝑦(1−𝑞)1−𝑦P(y|θ)=∑zP(y,z|θ)=∑zP(z|θ)P(y|z,θ)=πpy(1−p)1−y+(1−π)qy(1−q)1−y -
若𝑦=1y=1，表示这此看到的是正面，这个正面有可能是B的正面，也可能是C的正面，则𝑃(1|𝜃)=𝜋𝑝+(1−𝜋)𝑞P(1|θ)=πp+(1−π)q
-
若𝑦=0y=0，则𝑃(0|𝜃)=𝜋(1−𝑝)+(1−𝜋)(1−𝑞)P(0|θ)=π(1−p)+(1−π)(1−q)
y是观测变量，表示一次观测结果是1或0，z是隐藏变量，表示掷硬币A的结果，这个是观测不到结果的，𝜃=(𝜋,𝑝,𝑞)θ=(π,p,q)表示模型参数，将观测数据表示为𝑌=(𝑌1,𝑌2,...,𝑌𝑛)𝑇Y=(Y1,Y2,...,Yn)T，未观测的数据表示为𝑍=(𝑍1,𝑍2,...,𝑍𝑛)𝑇Z=(Z1,Z2,...,Zn)T，则观测函数的似然函数是：
考虑求模型参数𝜃=(𝜋,𝑝,𝑞)θ=(π,p,q)的极大似然估计，即：
这个问题没有解析解，只有通过迭代方法来求解，EM算法就是可以用于求解这个问题的一种迭代算法，下面给出EM算法的迭代过程：
-
首先选取初始值，记做𝜃0=(𝜋0,𝑝0,𝑞0)θ0=(π0,p0,q0)，第i次的迭代参数的估计值为𝜃𝑖=(𝜋𝑖,𝑝𝑖,𝑞𝑖)θi=(πi,pi,qi)
-
E步：计算在模型参数𝜋𝑖，𝑝𝑖，𝑞𝑖πi，pi，qi下观测变量𝑦𝑖yi来源于硬币B的概率：
𝜇𝑖+1=𝜋𝑖(𝑝𝑖)𝑦𝑖(1−𝑝𝑖)1−𝑦𝑖𝜋𝑖(𝑝𝑖)𝑦𝑖(1−𝑝𝑖)1−𝑦𝑖+(1−𝜋𝑖)(𝑞𝑖)𝑦𝑖(1−𝑝𝑖)1−𝑦𝑖μi+1=πi(pi)yi(1−pi)1−yiπi(pi)yi(1−pi)1−yi+(1−πi)(qi)yi(1−pi)1−yi备注一下：这个公式的分母是𝑃(𝑌|𝜃)P(Y|θ)，分子表示是来源与B硬币的概率。 -
M步：计算模型参数的新估计值：
𝜋𝑖+1=1𝑛∑𝑗=1𝑛𝜇𝑖+1𝑗πi+1=1n∑j=1nμji+1因为B硬币A硬币出现正面的结果，所以A硬币概率就是𝜇𝑗μj的平均值。𝑝𝑖+1=∑𝑛𝑗=1𝜇𝑖+1𝑗𝑦𝑗∑𝑛𝑗=1𝜇𝑖+1𝑗pi+1=∑j=1nμji+1yj∑j=1nμji+1分子乘以𝑦𝑖yi，所以其实是计算B硬币出现正面的概率。𝑞𝑖+1=∑𝑛𝑗=1(1−𝜇𝑖+1𝑗)𝑦𝑗∑𝑛𝑗=1(1−𝜇𝑖+1𝑗)qi+1=∑j=1n(1−μji+1)yj∑j=1n(1−μji+1)(1−𝜇𝑖+1𝑗)(1−μji+1)表示出现C硬币的概率。
闭环形成，从𝑃(𝑌|𝜃)P(Y|θ) 到 𝜋、𝑝、𝑞π、p、q一个闭环流程，接下来可以通过迭代法来做完成。针对上述例子，我们假设初始值为𝜋0=0.5，𝑝0=0.5，𝑞0=0.5π0=0.5，p0=0.5，q0=0.5，因为对𝑦𝑖=1yi=1和𝑦𝑖=0yi=0均有𝜇1𝑗=0.5μj1=0.5，利用迭代公式计算得到𝜋1=0.5，𝑝1=0.6，𝑞1=0.6π1=0.5，p1=0.6，q1=0.6，继续迭代得到最终的参数：
如果一开始初始值选择为：𝜋0=0.4，𝑝0=0.6，𝑞0=0.7π0=0.4，p0=0.6，q0=0.7，那么得到的模型参数的极大似然估计是𝜋ˆ=0.4064，𝑝ˆ=0.5368，𝑞ˆ=0.6432π^=0.4064，p^=0.5368，q^=0.6432，这说明EM算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。
这个例子中你只观察到了硬币抛完的结果，并不了解A硬币抛完之后，是选择了B硬币抛还是C硬币抛，这时候概率模型就存在着隐含变量！
EM算法
输入：观测变量数据Y，隐变量数据Z，联合分布𝑃(𝑌,𝑍|𝜃)P(Y,Z|θ)，条件分布𝑃(𝑍|𝑌,𝜃)P(Z|Y,θ)； 输出：模型参数𝜃θ
-
(1)选择参数的初值𝜃0θ0，开始迭代
-
(2) E步：记𝜃𝑖θi为第i次迭代参数𝜃θ的估计值，在第i+1次迭代的E步，计算
𝑄(𝜃,𝜃𝑖)=𝐸𝑍[𝑙𝑜𝑔𝑃(𝑌,𝑍|𝜃)|𝑌,𝜃𝑖]=∑𝑍𝑙𝑜𝑔𝑃(𝑌,𝑍|𝜃)𝑃(𝑍|𝑌,𝜃𝑖)Q(θ,θi)=EZ[logP(Y,Z|θ)|Y,θi]=∑ZlogP(Y,Z|θ)P(Z|Y,θi)这里，𝑃(𝑍|𝑌,𝜃𝑖)P(Z|Y,θi)是在给定观测数据Y和当前的参数估计𝜃𝑖θi下隐变量数据Z的条件概率分布； -
(3) M步：求使𝑄(𝜃,𝜃𝑖)Q(θ,θi)极大化的𝜃θ，确定第i+1次迭代的参数的估计值𝜃𝑖+1θi+1，
𝜃𝑖+1=𝑎𝑟𝑔max𝜃𝑄(𝜃,𝜃𝑖)θi+1=argmaxθQ(θ,θi)𝑄(𝜃,𝜃𝑖)Q(θ,θi)是EM算法的核心，称为Q函数(Q function)，这个是需要自己构造的。 -
(4) 重复第(2)步和第(3)步，直到收敛，收敛条件：
||𝜃𝑖+1−𝜃𝑖||<𝜀1||θi+1−θi||<ε1或者：||𝑄(𝜃𝑖+1,𝜃𝑖)−𝑄(𝜃𝑖,𝜃𝑖)||<𝜀2||Q(θi+1,θi)−Q(θi,θi)||<ε2收敛迭代就结束了。我们来拆解一下这个M步骤，
推导逼近
主要讲解Jensen不等式，这个公式在推导和收敛都用到，主要是如下的结论：
- 𝑓(𝑥)f(x)是凸函数
𝑓(𝐸(𝑋))≤𝐸(𝑓(𝑥))f(E(X))≤E(f(x))
- 𝑓(𝑥)f(x)是凹函数
𝑓(𝐸(𝑋))≥𝐸(𝑓(𝑥))f(E(X))≥E(f(x))
推导出Em算法可以近似实现对观测数据的极大似然估计的办法是找到E步骤的下界，让下届最大，通过逼近的方式实现对观测数据的最大似然估计。统计学习基础中采用的是相减方式，我们来看下具体的步骤。
- 增加隐藏变量
𝐿(𝜃)=∑𝑍𝑙𝑜𝑔𝑃(𝑌|𝑍,𝜃)𝑃(𝑍,𝜃)L(θ)=∑ZlogP(Y|Z,θ)P(Z,θ)则𝐿(𝜃)−𝐿(𝜃𝑖)L(θ)−L(θi)为：𝐿(𝜃)−𝐿(𝜃𝑖)=𝑙𝑜𝑔(∑𝑍𝑃(𝑌|𝑍,𝜃𝑖)𝑃(𝑌|𝑍,𝜃)𝑃(𝑍,𝜃)𝑃(𝑌|𝑍,𝜃𝑖))−𝐿(𝜃𝑖)≥∑𝑍𝑃(𝑌|𝑍,𝜃𝑖)𝑙𝑜𝑔(𝑃(𝑌|𝑍,𝜃)𝑃(𝑍,𝜃)𝑃(𝑌|𝑍,𝜃𝑖))−𝐿(𝜃𝑖)L(θ)−L(θi)=log(∑ZP(Y|Z,θi)P(Y|Z,θ)P(Z,θ)P(Y|Z,θi))−L(θi)≥∑ZP(Y|Z,θi)log(P(Y|Z,θ)P(Z,θ)P(Y|Z,θi))−L(θi)≥≥这一个步骤就是采用了凹函数的Jensen不等式做转换。因为𝑍Z是隐藏变量，所以有∑𝑍𝑃(𝑌|𝑍,𝜃𝑖)==1，𝑃(𝑌|𝑍,𝜃𝑖)>0∑ZP(Y|Z,θi)==1，P(Y|Z,θi)>0，于是继续变：
也就是：𝐿(𝜃)≥𝐿(𝜃𝑖)+∑𝑍𝑃(𝑍|𝑌,𝜃𝑖)𝑙𝑜𝑔(𝑃(𝑌|𝑍,𝜃)𝑃(𝑍,𝜃)𝑃(𝑌|𝑍,𝜃𝑖)𝐿(𝜃𝑖))L(θ)≥L(θi)+∑ZP(Z|Y,θi)log(P(Y|Z,θ)P(Z,θ)P(Y|Z,θi)L(θi))，有下界，最大化下界，来得到近似值。这里有一个细节：𝑃(𝑌|𝑍,𝜃𝑖)P(Y|Z,θi) 变为𝑃(𝑍|𝑌,𝜃𝑖)P(Z|Y,θi)？如果要满足Jensen不等式的等号，则有：
c为一个常数，而∑𝑍𝑃(𝑌|𝑍,𝜃𝑖)=1∑ZP(Y|Z,θi)=1则：
大家是不是很奇怪𝑃(𝑌|𝑍,𝜃)𝑃(𝑍,𝜃)P(Y|Z,θ)P(Z,θ)加上∑∑之后等于什么，其实有的博客这里使用𝑃(𝑍,𝜃)=𝑃(𝑌𝑖,𝑍𝑖,𝜃𝑖)P(Z,θ)=P(Yi,Zi,θi)来替代𝑃(𝑌|𝑍,𝜃)P(Y|Z,θ)参与计算，这样∑𝑍𝑃(𝑌𝑖,𝑍𝑖,𝜃𝑖)∑ZP(Yi,Zi,θi)，这样就方便理解来了。
于是最大化如下：
其中𝑙𝑜𝑔log分母提出来是关于𝑍Z的∑𝑍𝑃(𝑍|𝑌,𝜃𝑖)𝑙𝑜𝑔𝑃(𝑍|𝑌,𝜃𝑖)∑ZP(Z|Y,θi)logP(Z|Y,θi)，可以去掉。当然也有博客写的形式是：
形式其实一样，表示的不一样而已。
证明收敛
我们知道已知观测数据的似然函数是𝑃(𝑌,𝜃)P(Y,θ)，对数似然函数为：
要证明收敛，就证明单调递增，∑𝑀𝑖=1𝑙𝑜𝑔𝑃(𝑦𝑖,𝜃𝑗+1)>∑𝑀𝑖=1𝑙𝑜𝑔𝑃(𝑦𝑖,𝜃𝑗)∑i=1MlogP(yi,θj+1)>∑i=1MlogP(yi,θj) 由上文知道：
我们构造一个函数𝐻H，让他等于：
让𝑄(𝜃,𝜃𝑖)−𝐻(𝜃,𝜃𝑖)Q(θ,θi)−H(θ,θi)：
所以：
该公式左边已经被证明是大于0，证明右边：𝐻(𝜃𝑖+1,𝜃𝑖)−𝐻(𝜃𝑖,𝜃𝑖)<0H(θi+1,θi)−H(θi,θi)<0：
其中不等式是由于Jensen不等式，由此证明了∑𝑀𝑖=1𝑙𝑜𝑔𝑃(𝑦𝑖,𝜃𝑗+1)>∑𝑀𝑖=1𝑙𝑜𝑔𝑃(𝑦𝑖,𝜃𝑗)∑i=1MlogP(yi,θj+1)>∑i=1MlogP(yi,θj)，证明了EM算法的收敛性。但不能保证是全局最优，只能保证局部最优。
高斯混合分布
EM算法的一个重要应用场景就是高斯混合模型的参数估计。高斯混合模型就是由多个高斯模型组合在一起的混合模型（可以理解为多个高斯分布函数的线性组合，理论上高斯混合模型是可以拟合任意类型的分布），例如对于下图中的数据集如果用一个高斯模型来描述的话显然是不合理的：
两个高斯模型可以拟合数据集，如图所示：
如果有多个高斯模型，公式表示为：
𝜙(𝑦|𝜃𝑘)ϕ(y|θk)表示为第k个高斯分布密度模型，定义如上，其中𝑎𝑘ak表示被选中的概率。在本次模型𝑃(𝑦|𝜃)P(y|θ)中，观测数据是已知的，而观测数据具体来自哪个模型是未知的，有点像之前提过的三硬币模型，我们来对比一下，A硬币就像是概率𝑎𝑘ak，用来表明具体的模型，而B、C硬币就是具体的模型，只不过这里有很多个模型，不仅仅是B、C这两个模型。我们用𝛾𝑗𝑘γjk来表示，则：
所以一个观测数据𝑦𝑗yj的隐藏数据(𝛾𝑗1,𝛾𝑗2,...,𝛾𝑗𝑘)(γj1,γj2,...,γjk)，那么完全似然函数就是：
取对数之后等于：
- E 步 ：
𝑄(𝜃.𝜃𝑖)=𝐸[𝑙𝑜𝑔(𝑃(𝑦,𝛾|𝜃))]=∑𝐾𝑘=1(∑𝑗=1𝑁(𝐸𝛾𝑗𝑘)𝑙𝑜𝑔(𝑎𝑘)+∑𝑗=1𝑁(𝐸𝛾𝑗𝑘)[𝑙𝑜𝑔(12𝜋⎯⎯⎯⎯√)−𝑙𝑜𝑔(𝛿𝑘)−(𝑦𝑖−𝜇𝑘)22𝛿2𝑘])Q(θ.θi)=E[log(P(y,γ|θ))]=∑Kk=1(∑j=1N(Eγjk)log(ak)+∑j=1N(Eγjk)[log(12π)−log(δk)−(yi−μk)22δk2])其中我们定义𝛾𝑗𝑘^γjk^：𝛾𝑗𝑘^=𝐸(𝛾𝑗𝑘|𝑦,𝜃)=𝑎𝑘𝜙(𝑦𝑖|𝜃𝑘)∑𝐾𝑘=1𝑎𝑘𝜙(𝑦𝑖|𝜃𝑘)𝑗=1,2,..,𝑁；𝑘=1,2,...,𝐾𝑛𝑘=∑𝑗=𝑖𝑁𝐸𝛾𝑗𝑘γjk^=E(γjk|y,θ)=akϕ(yi|θk)∑k=1Kakϕ(yi|θk)j=1,2,..,N；k=1,2,...,Knk=∑j=iNEγjk于是化简得到：𝑄(𝜃.𝜃𝑖)=∑𝐾𝑘=1(𝑛𝑘𝑙𝑜𝑔(𝑎𝑘)+∑𝑗=1𝑁(𝐸𝛾𝑗𝑘)[𝑙𝑜𝑔(12𝜋⎯⎯⎯⎯√)−𝑙𝑜𝑔(𝛿𝑘)−(𝑦𝑖−𝜇𝑘)22𝛿2𝑘])Q(θ.θi)=∑Kk=1(nklog(ak)+∑j=1N(Eγjk)[log(12π)−log(δk)−(yi−μk)22δk2])
E 步 在代码设计上只有𝛾𝑗𝑘^γjk^有用，用于M步的计算。
- M步，
𝜃𝑖+1=𝑎𝑟𝑔max𝜃𝑄(𝜃,𝜃𝑖)θi+1=argmaxθQ(θ,θi)对𝑄(𝜃,𝜃𝑖)Q(θ,θi)求导，得到每个未知量的偏导，使其偏导等于0，求解得到：𝜇𝑘^=∑𝑁𝑗=1𝛾𝑗𝑘^𝑦𝑖∑𝑁𝑗=1𝛾𝑗𝑘^𝛿𝑘^=∑𝑁𝑗=1𝛾𝑗𝑘^(𝑦𝑖−𝜇𝑘)2∑𝑁𝑗=1𝛾𝑗𝑘^𝑎𝑘^=∑𝑁𝑗=1𝛾𝑗𝑘^𝑁μk^=∑j=1Nγjk^yi∑j=1Nγjk^δk^=∑j=1Nγjk^(yi−μk)2∑j=1Nγjk^ak^=∑j=1Nγjk^N给一个初始值，来回迭代就可以求得值内容。这一块主要用到了𝑄(𝜃.𝜃𝑖)Q(θ.θi)的导数，并且用到了E步的𝛾𝑗𝑘^γjk^。
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
