【机器学习算法系列之一】EM算法实例分析

from:http://chenrudan.github.io/blog/2015/12/02/emexample.html

最近两天研究了一下EM算法，主要是基于《统计学习方法》和论文《What is the expectation maximization algorithm?》[1]，但是对两个文章里面给的实例求解过程都比较的困惑，搜索网上的一些博客也没有找到对应的求解过程，自己就仔细研究了一下，中间也遇到了一些坑，现在把解题思路给出来。因为书上和网上的博客[2]对EM算法的推导和证明解释的非常清楚，本文就不做解释了，如果对EM算法原理不清楚的建议先看看《统计学习方法》第9章或者博客[2][3]。本文只给出两个文章中的例子的求解过程。

(题目我会列出来，如果不是看这个两个文章而了解EM算法的也不要紧，题目是通用的)

本文中观测数据记为Y(因为两个例子都是输出是观测数据)，隐藏变量(未观测变量)记为z，模型参数记为θθ。

1.三硬币模型

假设有三枚硬币A、B、C，每个硬币正面出现的概率是π、p、qπ、p、q。进行如下的掷硬币实验：先掷硬币A，正面向上选B，反面选C；然后掷选择的硬币，正面记1，反面记0。独立的进行10次实验，结果如下：1，1，0，1，0，0，1，0，1，1。假设只能观察最终的结果(0 or 1)，而不能观测掷硬币的过程(不知道选的是B or C)，问如何估计三硬币的正面出现的概率？

首先针对某个输出y值，它在参数θ(θ=(π,p,q))θ(θ=(π,p,q))下的概率分布为

P(y|θ)=∑zP(y,z|θ)=∑zP(z|θ)P(y|z,θ)=πpy(1−p)1−y+(1−π)qy(1−q)1−y
P(y|θ)=∑zP(y,z|θ)=∑zP(z|θ)P(y|z,θ)=πpy(1−p)1−y+(1−π)qy(1−q)1−y
从而针对观测数据Y=(y1,y2,⋅⋅⋅,yn)TY=(y1,y2,⋅⋅⋅,yn)T的似然函数为

P(Y|θ)=∑zP(Y,z|θ)=∑zP(z|θ)P(Y|z,θ)=∏j=1nπpyj(1−p)1−yj+(1−π)qyj(1−q)1−yj
P(Y|θ)=∑zP(Y,z|θ)=∑zP(z|θ)P(Y|z,θ)=∏j=1nπpjy(1−p)1−yj+(1−π)qjy(1−q)1−yj
因此本题的目标是求解参数θθ的极大似然估计，即θ^=argmaxθlogP(Y|θ)θ^=argmaxθlogP(Y|θ)。直接对连乘的似然函数求导太复杂，所以一般用极大似然估计都会转化成对数似然函数，但是就算转化成了求和，如果这个式子对某个参数(例如ππ)求导，由于这个式子中有“和的对数”，求导非常复杂。因此这个问题需要用EM算法来求解。

E步:根据EM算法，在这一步需要计算的是未观测数据的条件概率分布，也就是每一个P(z|yj,θ)P(z|yj,θ)，μi+1μi+1表示在已知的模型参数θiθi下观测数据yjyj来自掷硬币B的概率，相应的来自掷C的概率就是1−μi+11−μi+1。

μi+1=πi(pi)yj(1−pi)1−yjπi(pi)yj(1−pi)1−yj+(1−πi)(qi)yj(1−qi)1−yj
μi+1=πi(pi)yj(1−pi)1−yjπi(pi)yj(1−pi)1−yj+(1−πi)(qi)yj(1−qi)1−yj
这里的分子就是z取掷硬币B和y的联合概率分布，需要注意的是，这里的μi+1μi+1通过E步的计算就已经是一个常数了，后面的求导不需要把这个式子代入。

M步:针对Q函数求导，Q函数的表达式是

Q(θ,θi)=∑j=1N∑zP(z|yj,θi)logP(yj,z|θ)=∑j=1Nμjlog(πpyj(1−p)1−yj)+(1−μj)log((1−π)qyj(1−q)1−yj)]
Q(θ,θi)=∑j=1N∑zP(z|yj,θi)logP(yj,z|θ)=∑j=1Nμjlog(πpyj(1−p)1−yj)+(1−μj)log((1−π)qyj(1−q)1−yj)]
最开始求导犯了一个大错，没有将表达式展开来求，这样就直接默认μjμj是一个系数，求导将它给约去了，这样就得不到最后的结果。

∂Q∂π=(μ1π−1−μ11−π)+⋅⋅⋅+(μNπ−1−μN1−π)=μ1−ππ(1−π)+⋅⋅⋅+μN−ππ(1−π)=∑Nj=1μj−Nππ(1−π)
∂Q∂π=(μ1π−1−μ11−π)+⋅⋅⋅+(μNπ−1−μN1−π)=μ1−ππ(1−π)+⋅⋅⋅+μN−ππ(1−π)=∑j=1Nμj−Nππ(1−π)
再令这个结果等于0，即获得πi+1=1N∑Nj=1μi+1jπi+1=1N∑j=1Nμji+1，其他两个也同理。

2.两硬币模型

假设有两枚硬币A、B，以相同的概率随机选择一个硬币，进行如下的掷硬币实验：共做5次实验，每次实验独立的掷十次，结果如图中a所示，例如某次实验产生了H、T、T、T、H、H、T、H、T、H，H代表证明朝上。a是在知道每次选择的是A还是B的情况下进行，b是在不知道选择的硬币情况下进行，问如何估计两个硬币正面出现的概率？

1

图1 两硬币模型题目(图来源[1])
针对a情况，已知选择的A or B，重点是如何计算输出的概率分布，论文中直接统计了5次实验中A正面向上的次数再除以总次数作为A的θA^θA^，这其实也是极大似然求导求出来的。

argmaxθlogP(Y|θ)=log((θ5B(1−θB)5)(θ9A(1−θA))(θ8A(1−θA)2)(θ4B(1−θB)6)(θ7A(1−θA)3))=log((θ24A(1−θA)6)(θ9B(1−θB)11))
argmaxθlogP(Y|θ)=log((θB5(1−θB)5)(θA9(1−θA))(θA8(1−θA)2)(θB4(1−θB)6)(θA7(1−θA)3))=log((θA24(1−θA)6)(θB9(1−θB)11))
上面这个式子求导之后就能得出θA^=2424+6=0.80θA^=2424+6=0.80以及θB^=99+11=0.45θB^=99+11=0.45。

针对b情况，由于并不知道选择的是A还是B，因此采用EM算法。

E步:计算在给定的θ(0)A^θA(0)^和θ(0)B^θB(0)^下，选择的硬币可能是A or B的概率，例如第一个实验中选择A的概率为(由于选择A、B的过程是等概率的，这个系数被我省略掉了)

P(z=A|y1,θ)=P(z=A,y1|θ)P(z=A,y1|θ)+P(z=B,y1|θ)=(0.6)5∗(0.4)5(0.6)5∗(0.4)5+(0.5)10=0.45
P(z=A|y1,θ)=P(z=A,y1|θ)P(z=A,y1|θ)+P(z=B,y1|θ)=(0.6)5∗(0.4)5(0.6)5∗(0.4)5+(0.5)10=0.45
M步:针对Q函数求导，在本题中Q函数形式如下，参数设置参照例1，只是这里的yjyj代表的是每次正面朝上的个数。

Q(θ,θi)=∑j=1N∑zP(z|yj,θi)logP(yj,z|θ)=∑j=1Nμjlog(θyjA(1−θA)10−yj)+(1−μj)log(θyjB(1−θB)10−yj)]
Q(θ,θi)=∑j=1N∑zP(z|yj,θi)logP(yj,z|θ)=∑j=1Nμjlog(θAyj(1−θA)10−yj)+(1−μj)log(θByj(1−θB)10−yj)]
从而针对这个式子来对参数求导，例如对θAθA求导

∂Q∂θA=μ1(y1θA−10−y11−θA)+⋅⋅⋅+μ5(y5θA−10−y51−θA)=μ1(y1−10θAθA(1−θA))+⋅⋅⋅+μ5(y5−10θAθA(1−θA))=∑5j=1μjyj−∑5j=110μjθAθA(1−θA)
∂Q∂θA=μ1(y1θA−10−y11−θA)+⋅⋅⋅+μ5(y5θA−10−y51−θA)=μ1(y1−10θAθA(1−θA))+⋅⋅⋅+μ5(y5−10θAθA(1−θA))=∑j=15μjyj−∑j=1510μjθAθA(1−θA)
求导等于0之后就可得到图中的第一次迭代之后的参数值θ(1)A^=0.71θA(1)^=0.71和θ(1)B^=0.58θB(1)^=0.58。

这个例子可以非常直观的看出来，EM算法在求解M步是将每次实验硬币取A或B的情况都考虑进去了。

3.小结

EM算法将不完全数据补全成完全数据，而E步并不是只取最可能补全的未观测数据，而是将未观测的数据的所有补全可能都计算出对应的概率值，从而对这些所有可能的补全计算出它们的期望值，作为下一步的未观测数据。至于为什么取期望，一是因为这个未观测数据本身就是基于一组不完全正确的参数估计出来的，例如三硬币例子假如每次在进行maximization之前都只取某一个值(极端一点，每次结果都是认为B是最可能的观测数据，而不算C)，那么在更新参数时，也只有B的参数在更新。二是这种情况下JENSEN不等式不成立，那么对θθ的似然函数变换形式就不成立，收敛也不成立。

这两个例子想明白之后求解实际上非常简单，所以很多博主并没把它们列出来，但如果一开始思考的方向不对就会浪费很多时间，当我把上面的过程想清楚之后再去求解别的例子，发现很轻松就能解出来。当然EM算法的核心还是证明和推导，这点别的文章讲的非常清晰了我就不赘述了。这也是数学上常用的思路，当无法直接对某个含参式子求极大值时，考虑对它的下界求极大值，当确定下界取极大值的参数时也能让含参式子值变大，也就是不断求解下界的极大值逼近求解对数似然函数极大化(李航.《统计学习方法》)。

如果本文有错误，请一定要指出来，感谢～

4.参考：

[1] What is the expectation maximization
algorithm?

[2] （EM算法）The EM Algorithm

[3] 从最大似然到EM算法浅解