Datawhale-机器学习算法-Task3 EM算法

一.前言

EM算法是机器学习十大算法之一，它很简单，但是也同样很有深度，简单是因为它就分两步求解问题，

E步：求期望（expectation）
M步：求极大（maximization)

深度在于它的数学推理涉及到比较繁杂的概率公式等，所以本文会介绍很多概率方面的知识，不懂的同学可以先去了解一些知识，当然本文也会尽可能的讲解清楚这些知识，讲的不好的地方麻烦大家评论指出，后续不断改进完善。

二.EM算法引入

概率模型有时候既含有观测变量，又含有隐变量或潜在变量，如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计方法估计模型参数，但是当模型含有隐变量时，就不能简单的使用这些方法，EM算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验概率估计法，我们讨论极大似然估计，极大后验概率估计与其类似。
参考统计学习方法书中的一个例子来引入EM算法，假设有3枚硬币，分别记做A、B、C，这些硬币正面出现的概率分别是π\piπ、ppp、qqq，进行如下实验：

先掷硬币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,0,1,0,0,1,0,1,1
假设只通过观测到掷硬币的结果，不能观测掷硬币的过程，问如何估计三个硬币出现正面的概率？
我们来构建这样一个三硬币模型：
P(y∣θ)=∑zP(y,z∣θ)=∑zP(z∣θ)P(y∣z,θ)=πpy(1−p)1−y+(1−π)qy(1−q)1−y\begin{aligned} P(y|\theta) &=\sum_{z}P(y,z|\theta)=\sum_{z}P(z|\theta)P(y|z,\theta) \\ &=\pi p^{y}(1-p)^{1-y}+(1-\pi)q^{y}(1-q)^{1-y} \end{aligned} P(y∣θ)=z∑P(y,z∣θ)=z∑P(z∣θ)P(y∣z,θ)=πpy(1−p)1−y+(1−π)qy(1−q)1−y
若y=1y=1y=1，表示这此看到的是正面，这个正面有可能是B的正面，也可能是C的正面，则P(1∣θ)=πp+(1−π)qP(1|\theta)=\pi p+(1-\pi)qP(1∣θ)=πp+(1−π)q
若y=0y=0y=0，则P(0∣θ)=π(1−p)+(1−π)(1−q)P(0|\theta)=\pi (1-p)+(1-\pi)(1-q)P(0∣θ)=π(1−p)+(1−π)(1−q)

y是观测变量，表示一次观测结果是1或0，z是隐藏变量，表示掷硬币A的结果，这个是观测不到结果的，θ=(π,p,q)\theta=(\pi,p,q)θ=(π,p,q)表示模型参数，将观测数据表示为Y=(Y1,Y2,...,Yn)TY=(Y_1,Y_2,...,Y_n)^{T}Y=(Y1,Y2,...,Yn)T，未观测的数据表示为Z=(Z1,Z2,...,Zn)TZ=(Z_1,Z_2,...,Z_n)^{T}Z=(Z1,Z2,...,Zn)T，则观测函数的似然函数是：
P(Y∣θ)=∑ZP(Z∣θ)P(Y∣Z,θ)=∏i=0(πpyi(1−p)1−yi+(1−π)qyi(1−q)1−yi)\begin{aligned} P(Y|\theta)&=\sum_{Z}P(Z|\theta)P(Y|Z,\theta)\\ &=\prod_{i=0} ( \pi p^{y_i}(1-p)^{1-y_{i}}+(1-\pi)q^{y_{i}}(1-q)^{1-y_{i}}) \end{aligned} P(Y∣θ)=Z∑P(Z∣θ)P(Y∣Z,θ)=i=0∏(πpyi(1−p)1−yi+(1−π)qyi(1−q)1−yi)
考虑求模型参数θ=(π,p,q)\theta=(\pi,p,q)θ=(π,p,q)的极大似然估计，即：
θ^=argmax⁡θlogP(Y∣θ)\hat{\theta}=arg\max_{\theta}logP(Y|\theta) θ^=argθmaxlogP(Y∣θ)
这个问题没有解析解，只有通过迭代方法来求解，EM算法就是可以用于求解这个问题的一种迭代算法，下面给出EM算法的迭代过程：

首先选取初始值，记做θ0=(π0,p0,q0)\theta^{0}=(\pi^{0},p^{0},q^{0})θ0=(π0,p0,q0)，第i次的迭代参数的估计值为θi=(πi,pi,qi)\theta^{i}=(\pi^{i},p^{i},q^{i})θi=(πi,pi,qi)
E步：计算在模型参数πi，pi，qi\pi^{i}，p^{i}，q^{i}πi，pi，qi下观测变量yiy_iyi来源于硬币B的概率：
μi+1=πi(pi)yi(1−pi)1−yiπi(pi)yi(1−pi)1−yi+(1−πi)(qi)yi(1−pi)1−yi\mu^{i+1}=\frac{\pi^{i}(p^{i})^{y_i}(1-p^i)^{1-y_i}}{\pi^{i}(p^{i})^{y_i}(1-p^i)^{1-y_i}+(1-\pi^{i})(q^{i})^{y_i}(1-p^i)^{1-y_i}} μi+1=πi(pi)yi(1−pi)1−yi+(1−πi)(qi)yi(1−pi)1−yiπi(pi)yi(1−pi)1−yi
备注一下：这个公式的分母是P(Y∣θ)P(Y|\theta)P(Y∣θ)，分子表示是来源与B硬币的概率。
M步：计算模型参数的新估计值：
πi+1=1n∑j=1nμji+1\pi^{i+1}=\frac{1}{n}\sum_{j=1}^{n}\mu_{j}^{i+1} πi+1=n1j=1∑nμji+1
因为B硬币A硬币出现正面的结果，所以A硬币概率就是μj\mu_{j}μj的平均值。
pi+1=∑j=1nμji+1yj∑j=1nμji+1p^{i+1}=\frac{\sum_{j=1}^{n}\mu_{j}^{i+1}y_j}{\sum_{j=1}^{n}\mu_{j}^{i+1}} pi+1=∑j=1nμji+1∑j=1nμji+1yj
分子乘以yiy_{i}yi，所以其实是计算B硬币出现正面的概率。
qi+1=∑j=1n(1−μji+1)yj∑j=1n(1−μji+1)q^{i+1}=\frac{\sum_{j=1}^{n}(1-\mu_{j}^{i+1})y_j}{\sum_{j=1}^{n}(1-\mu_{j}^{i+1})} qi+1=∑j=1n(1−μji+1)∑j=1n(1−μji+1)yj
(1−μji+1)(1-\mu_{j}^{i+1})(1−μji+1)表示出现C硬币的概率。

闭环形成，从P(Y∣θ)P(Y|\theta)P(Y∣θ) 到 π、p、q\pi、p、qπ、p、q一个闭环流程，接下来可以通过迭代法来做完成。针对上述例子，我们假设初始值为π0=0.5，p0=0.5，q0=0.5\pi^{0}=0.5，p^{0}=0.5，q^{0}=0.5π0=0.5，p0=0.5，q0=0.5，因为对yi=1y_i=1yi=1和yi=0y_i=0yi=0均有μj1=0.5\mu_j^{1}=0.5μj1=0.5，利用迭代公式计算得到π1=0.5，p1=0.6，q1=0.6\pi^{1}=0.5，p^{1}=0.6，q^{1}=0.6π1=0.5，p1=0.6，q1=0.6，继续迭代得到最终的参数：
π0^=0.5，p0^=0.6，q0^=0.6\widehat{\pi^{0}}=0.5，\widehat{p^{0}}=0.6，\widehat{q^{0}}=0.6π0=0.5，p0=0.6，q0=0.6
如果一开始初始值选择为：π0=0.4，p0=0.6，q0=0.7\pi^{0}=0.4，p^{0}=0.6，q^{0}=0.7π0=0.4，p0=0.6，q0=0.7，那么得到的模型参数的极大似然估计是π^=0.4064，p^=0.5368，q^=0.6432\widehat{\pi}=0.4064，\widehat{p}=0.5368，\widehat{q}=0.6432π=0.4064，p=0.5368，q=0.6432，这说明EM算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。

这个例子中你只观察到了硬币抛完的结果，并不了解A硬币抛完之后，是选择了B硬币抛还是C硬币抛，这时候概率模型就存在着隐含变量！

三.EM算法

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

(1)选择参数的初值θ0\theta^0θ0，开始迭代
(2) E步：记θi\theta^iθi为第i次迭代参数θ\thetaθ的估计值，在第i+1次迭代的E步，计算
Q(θ,θi)=EZ[logP(Y,Z∣θ)∣Y,θi]=∑ZlogP(Y,Z∣θ)P(Z∣Y,θi)\begin{aligned} Q(\theta,\theta^i)&=E_{Z}[logP(Y,Z|\theta)|Y,\theta^i]\\ &=\sum_{Z}logP(Y,Z|\theta)P(Z|Y,\theta^i) \end{aligned} Q(θ,θi)=EZ[logP(Y,Z∣θ)∣Y,θi]=Z∑logP(Y,Z∣θ)P(Z∣Y,θi)
这里，P(Z∣Y,θi)P(Z|Y,\theta^i)P(Z∣Y,θi)是在给定观测数据Y和当前的参数估计θi\theta^iθi下隐变量数据Z的条件概率分布；
(3) M步：求使Q(θ,θi)Q(\theta,\theta^i)Q(θ,θi)极大化的θ\thetaθ，确定第i+1次迭代的参数的估计值θi+1\theta^{i+1}θi+1，
θi+1=argmax⁡θQ(θ,θi)\theta^{i+1}=arg \max \limits_{\theta}Q(\theta,\theta^{i}) θi+1=argθmaxQ(θ,θi)
Q(θ,θi)Q(\theta,\theta^{i})Q(θ,θi)是EM算法的核心，称为Q函数(Q function)，这个是需要自己构造的。
(4) 重复第(2)步和第(3)步，直到收敛，收敛条件：
∣∣θi+1−θi∣∣<ε1|| \theta^{i+1}-\theta^{i} || < \varepsilon_1 ∣∣θi+1−θi∣∣<ε1
或者：
∣∣Q(θi+1,θi)−Q(θi,θi)∣∣<ε2||Q(\theta^{i+1},\theta^{i})-Q(\theta^{i},\theta^{i})|| <\varepsilon_2 ∣∣Q(θi+1,θi)−Q(θi,θi)∣∣<ε2
收敛迭代就结束了。我们来拆解一下这个M步骤，

四.推导逼近

主要讲解Jensen不等式，这个公式在推导和收敛都用到，主要是如下的结论：

f(x)f(x)f(x)是凸函数
f(E(X))≤E(f(x))f(E(X)) \le E(f(x)) f(E(X))≤E(f(x))
f(x)f(x)f(x)是凹函数
f(E(X))≥E(f(x))f(E(X)) \ge E(f(x)) f(E(X))≥E(f(x))

推导出Em算法可以近似实现对观测数据的极大似然估计的办法是找到E步骤的下界，让下届最大，通过逼近的方式实现对观测数据的最大似然估计。统计学习基础中采用的是相减方式，我们来看下具体的步骤。

增加隐藏变量
L(θ)=∑ZlogP(Y∣Z,θ)P(Z,θ)L(\theta)=\sum_{Z}logP(Y|Z,\theta)P(Z,\theta) L(θ)=Z∑logP(Y∣Z,θ)P(Z,θ)
则L(θ)−L(θi)L(\theta)-L(\theta^{i})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)\begin{aligned} L(\theta)-L(\theta^{i})=log(\sum_{Z} P(Y|Z,\theta^i)\frac{P(Y|Z,\theta)P(Z,\theta)}{P(Y|Z,\theta^i)})-L(\theta^{i})\\ \ge \sum_{Z} P(Y|Z,\theta^i)log(\frac{P(Y|Z,\theta)P(Z,\theta)}{P(Y|Z,\theta^i)})-L(\theta^{i}) \end{aligned} L(θ)−L(θi)=log(Z∑P(Y∣Z,θi)P(Y∣Z,θi)P(Y∣Z,θ)P(Z,θ))−L(θi)≥Z∑P(Y∣Z,θi)log(P(Y∣Z,θi)P(Y∣Z,θ)P(Z,θ))−L(θi)
≥\ge≥这一个步骤就是采用了凹函数的Jensen不等式做转换。因为ZZZ是隐藏变量，所以有∑ZP(Y∣Z,θi)==1，P(Y∣Z,θi)>0\sum_{Z} P(Y|Z,\theta^i)==1，P(Y|Z,\theta^i)>0∑ZP(Y∣Z,θi)==1，P(Y∣Z,θi)>0，于是继续变：

L(θ)−L(θi)=log(∑ZP(Y∣Z,θi)P(Y∣Z,θ)P(Z,θ)P(Y∣Z,θi))−L(θi)≥∑ZP(Z∣Y,θi)log(P(Y∣Z,θ)P(Z,θ)P(Z∣Y,θi))−L(θi)=∑ZP(Z∣Y,θi)log(P(Y∣Z,θ)P(Z,θ)P(Z∣Y,θi))−∑ZP(Z∣Y,θi)L(θi)=∑ZP(Z∣Y,θi)log(P(Y∣Z,θ)P(Z,θ)P(Z∣Y,θi)(P(Y∣θi))≥0\begin{aligned} L(\theta)-L(\theta^{i})&=log(\sum_{Z} P(Y|Z,\theta^i)\frac{P(Y|Z,\theta)P(Z,\theta)}{P(Y|Z,\theta^i)})-L(\theta^{i})\\ &\ge \sum_{Z} P(Z|Y,\theta^i)log(\frac{P(Y|Z,\theta)P(Z,\theta)}{P(Z|Y,\theta^i)})-L(\theta^{i})\\ &=\sum_{Z} P(Z|Y,\theta^i)log(\frac{P(Y|Z,\theta)P(Z,\theta)}{P(Z|Y,\theta^i)})-\sum_{Z} P(Z|Y,\theta^i)L(\theta^{i})\\ &= \sum_{Z} P(Z|Y,\theta^i)log(\frac{P(Y|Z,\theta)P(Z,\theta)}{P(Z|Y,\theta^i) (P(Y|\theta^{i})}) \\ & \ge0 \end{aligned} L(θ)−L(θi)=log(Z∑P(Y∣Z,θi)P(Y∣Z,θi)P(Y∣Z,θ)P(Z,θ))−L(θi)≥Z∑P(Z∣Y,θi)log(P(Z∣Y,θi)P(Y∣Z,θ)P(Z,θ))−L(θi)=Z∑P(Z∣Y,θi)log(P(Z∣Y,θi)P(Y∣Z,θ)P(Z,θ))−Z∑P(Z∣Y,θi)L(θi)=Z∑P(Z∣Y,θi)log(P(Z∣Y,θi)(P(Y∣θi)P(Y∣Z,θ)P(Z,θ))≥0
也就是：L(θ)≥L(θi)+∑ZP(Z∣Y,θi)log(P(Y∣Z,θ)P(Z,θ)P(Y∣Z,θi)L(θi))L(\theta)\ge L(\theta^{i})+ \sum_{Z} P(Z|Y,\theta^i)log(\frac{P(Y|Z,\theta)P(Z,\theta)}{P(Y|Z,\theta^i) L(\theta^{i})})L(θ)≥L(θi)+∑ZP(Z∣Y,θi)log(P(Y∣Z,θi)L(θi)P(Y∣Z,θ)P(Z,θ))，有下界，最大化下界，来得到近似值。这里有一个细节：P(Y∣Z,θi)P(Y|Z,\theta^i)P(Y∣Z,θi) 变为P(Z∣Y,θi)P(Z|Y,\theta^i)P(Z∣Y,θi)？如果要满足Jensen不等式的等号，则有：
P(Y∣Z,θ)P(Z,θ)P(Y∣Z,θi)=c\frac{P(Y|Z,\theta)P(Z,\theta)}{P(Y|Z,\theta^i)} = c P(Y∣Z,θi)P(Y∣Z,θ)P(Z,θ)=c
c为一个常数，而∑ZP(Y∣Z,θi)=1\sum_{Z}P(Y|Z,\theta^i)=1∑ZP(Y∣Z,θi)=1则：
∑ZP(Y∣Z,θ)P(Z,θ)=c∑ZP(Y∣Z,θi)=c=P(Y∣Z,θ)P(Z,θ)P(Y∣Z,θi)P(Y∣Z,θ)=P(Y∣Z,θ)P(Z,θ)∑ZP(Y∣Z,θ)P(Z,θ)=P(Y,Z,θ)P(Y,θ)=P(Z∣Y,θ)\begin{aligned} \sum_{Z}P(Y|Z,\theta)P(Z,\theta)= c\sum_{Z}P(Y|Z,\theta^i)&=c\\ &=\frac{P(Y|Z,\theta)P(Z,\theta)}{P(Y|Z,\theta^i)}\\ P(Y|Z,\theta)=\frac{P(Y|Z,\theta)P(Z,\theta)}{\sum_{Z}P(Y|Z,\theta)P(Z,\theta)}=\frac{P(Y,Z,\theta)}{P(Y,\theta)}=P(Z|Y,\theta) \end{aligned} Z∑P(Y∣Z,θ)P(Z,θ)=cZ∑P(Y∣Z,θi)P(Y∣Z,θ)=∑ZP(Y∣Z,θ)P(Z,θ)P(Y∣Z,θ)P(Z,θ)=P(Y,θ)P(Y,Z,θ)=P(Z∣Y,θ)=c=P(Y∣Z,θi)P(Y∣Z,θ)P(Z,θ)

大家是不是很奇怪P(Y∣Z,θ)P(Z,θ)P(Y|Z,\theta)P(Z,\theta)P(Y∣Z,θ)P(Z,θ)加上∑\sum∑之后等于什么，其实有的博客这里使用P(Z,θ)=P(Yi,Zi,θi)P(Z,\theta) = P(Y^i,Z^i,\theta^i)P(Z,θ)=P(Yi,Zi,θi)来替代P(Y∣Z,θ)P(Y|Z,\theta)P(Y∣Z,θ)参与计算，这样∑ZP(Yi,Zi,θi)\sum_{Z}P(Y^i,Z^i,\theta^i)∑ZP(Yi,Zi,θi)，这样就方便理解来了。

于是最大化如下：

θi+1=argmax⁡θ∑ZP(Z∣Y,θi)log(P(Y∣Z,θ)P(Z,θ)P(Z∣Y,θi))=argmax⁡θ∑ZP(Z∣Y,θi)log(P(Y∣Z,θ)P(Z,θ))=argmax⁡θ∑ZP(Z∣Y,θi)log(P(Y,Z∣θ))=argmax⁡θQ(θ,θi)\begin{aligned} \theta^{i+1}&=arg \max_{\theta}\sum_{Z} P(Z|Y,\theta^i)log(\frac{P(Y|Z,\theta)P(Z,\theta)}{P(Z|Y,\theta^i)})\\ &=arg \max_{\theta}\sum_{Z} P(Z|Y,\theta^i)log(P(Y|Z,\theta)P(Z,\theta))\\ & =arg \max_{\theta}\sum_{Z} P(Z|Y,\theta^i)log(P(Y,Z|\theta))\\ &=arg \max_{\theta}Q(\theta,\theta^i) \end{aligned} θi+1=argθmaxZ∑P(Z∣Y,θi)log(P(Z∣Y,θi)P(Y∣Z,θ)P(Z,θ))=argθmaxZ∑P(Z∣Y,θi)log(P(Y∣Z,θ)P(Z,θ))=argθmaxZ∑P(Z∣Y,θi)log(P(Y,Z∣θ))=argθmaxQ(θ,θi)
其中logloglog分母提出来是关于ZZZ的∑ZP(Z∣Y,θi)logP(Z∣Y,θi)\sum_{Z} P(Z|Y,\theta^i)logP(Z|Y,\theta^i)∑ZP(Z∣Y,θi)logP(Z∣Y,θi)，可以去掉。当然也有博客写的形式是：
argmax⁡θ∑i=1M∑ZiP(Zi∣Yi,θi)log(P(Yi,Zi;θ))arg \max_{\theta}\sum_{i=1}^{M}\sum_{Z^{i}} P(Z^{i}|Y^{i},\theta^i)log(P(Y^{i},Z^{i};\theta))\\ argθmaxi=1∑MZi∑P(Zi∣Yi,θi)log(P(Yi,Zi;θ))
形式其实一样，表示的不一样而已。

五.证明收敛

我们知道已知观测数据的似然函数是P(Y,θ)P(Y,\theta)P(Y,θ)，对数似然函数为：
L()=∑i=1MlogP(yi,θ)=∑i=1Mlog(P(yi,Z∣θ)P(Z∣yi,θ))=∑i=1MlogP(yi,Z∣θ)−∑i=1MlogP(Z∣yi,θ)\begin{aligned} L()=\sum_{i=1}^{M}logP(y^{i},\theta) &=\sum_{i=1}^{M}log(\frac{P(y^i,Z|\theta)}{P(Z|y^i,\theta)})\\ &=\sum_{i=1}^{M}logP(y^i,Z|\theta) - \sum_{i=1}^{M}logP(Z|y^i,\theta) \end{aligned} L()=i=1∑MlogP(yi,θ)=i=1∑Mlog(P(Z∣yi,θ)P(yi,Z∣θ))=i=1∑MlogP(yi,Z∣θ)−i=1∑MlogP(Z∣yi,θ)
要证明收敛，就证明单调递增，∑i=1MlogP(yi,θj+1)>∑i=1MlogP(yi,θj)\sum_{i=1}^{M}logP(y^{i},\theta^{j+1})>\sum_{i=1}^{M}logP(y^{i},\theta^{j})∑i=1MlogP(yi,θj+1)>∑i=1MlogP(yi,θj)
由上文知道：
Q(θ,θi)=∑ZlogP(Y,Z∣θ)P(Z∣Y,θi)=∑i=1M∑ZjlogP(yi,Zj∣θ)P(Zj∣yi,θi)\begin{aligned} Q(\theta,\theta^i)&=\sum_{Z}logP(Y,Z|\theta)P(Z|Y,\theta^i)\\ &=\sum_{i=1}^{M}\sum_{Z^j}logP(y^i,Z^j|\theta)P(Z^j|y^i,\theta^i) \end{aligned} Q(θ,θi)=Z∑logP(Y,Z∣θ)P(Z∣Y,θi)=i=1∑MZj∑logP(yi,Zj∣θ)P(Zj∣yi,θi)
我们构造一个函数HHH，让他等于：
H(θ,θi)=∑i=1M∑Zjlog(P(Z∣yi,θ)P(Z∣yi,θi))H(\theta,\theta^{i})=\sum_{i=1}^{M}\sum_{Z^j}log(P(Z|y^i,\theta)P(Z|y^i,\theta^i)) H(θ,θi)=i=1∑MZj∑log(P(Z∣yi,θ)P(Z∣yi,θi))
让Q(θ,θi)−H(θ,θi)Q(\theta,\theta^i)-H(\theta,\theta^{i})Q(θ,θi)−H(θ,θi)：
Q(θ,θi)−H(θ,θi)=∑i=1M∑ZjlogP(yi,Zj∣θ)P(Zj∣yi,θi)−∑i=1M∑Zjlog(P(Zj∣yi,θ)P(Zj∣yi,θi))=∑i=1M∑Zjlog(P(yi,Zj∣θ)−P(Zj∣yi,θ))=∑i=1MlogP(yi,θ)\begin{aligned} Q(\theta,\theta^i)-H(\theta,\theta^{i})&=\sum_{i=1}^{M}\sum_{Z^j}logP(y^i,Z^j|\theta)P(Z^j|y^i,\theta^i) - \sum_{i=1}^{M}\sum_{Z^j}log(P(Z^j|y^i,\theta)P(Z^j|y^i,\theta^i)) \\ &=\sum_{i=1}^{M}\sum_{Z^j}log\bigg(P(y^i,Z^j|\theta)-P(Z^j|y^i,\theta)\bigg) \\ &=\sum_{i=1}^{M}logP(y^{i},\theta) \end{aligned} Q(θ,θi)−H(θ,θi)=i=1∑MZj∑logP(yi,Zj∣θ)P(Zj∣yi,θi)−i=1∑MZj∑log(P(Zj∣yi,θ)P(Zj∣yi,θi))=i=1∑MZj∑log(P(yi,Zj∣θ)−P(Zj∣yi,θ))=i=1∑MlogP(yi,θ)所以：
∑i=1MlogP(yi,θj+1)−∑i=1MlogP(yi,θj)=Q(θi+1,θi)−H(θi+1,θi)−(Q(θi,θi)−H(θi,θi))=Q(θi+1,θi)−Q(θi,θi)−(H(θi+1,θi)−H(θi,θi))\sum_{i=1}^{M}logP(y^{i},\theta^{j+1})-\sum_{i=1}^{M}logP(y^{i},\theta^{j}) \\ = Q(\theta^{i+1},\theta^i)-H(\theta^{i+1},\theta^{i}) - (Q(\theta^{i},\theta^{i})-H(\theta^{i},\theta^{i}))\\ = Q(\theta^{i+1},\theta^i)- Q(\theta^{i},\theta^{i}) -( H(\theta^{i+1},\theta^{i}) - H(\theta^{i},\theta^{i})) i=1∑MlogP(yi,θj+1)−i=1∑MlogP(yi,θj)=Q(θi+1,θi)−H(θi+1,θi)−(Q(θi,θi)−H(θi,θi))=Q(θi+1,θi)−Q(θi,θi)−(H(θi+1,θi)−H(θi,θi))
该公式左边已经被证明是大于0，证明右边：H(θi+1,θi)−H(θi,θi)<0H(\theta^{i+1},\theta^{i}) - H(\theta^{i},\theta^{i})<0H(θi+1,θi)−H(θi,θi)<0：
H(θi+1,θi)−H(θi,θi)=∑Zj(log(P(Zj∣Y,θi+1)P(Zj∣Y,θi)))P(Zj∣Y,θi)=log(∑ZjP(Zj∣Y,θi+1)P(Zj∣Y,θi)P(Zj∣Y,θi))=logP(Z∣Y,θi+1)=log1=0\begin{aligned} H(\theta^{i+1},\theta^{i}) - H(\theta^{i},\theta^{i}) &=\sum_{Z^j}\bigg(log(\frac{P(Z^j|Y,\theta^{i+1})}{P(Z^j|Y,\theta^i)}) \bigg)P(Z^j|Y,\theta^i) \\ &=log\bigg(\sum_{Z^j}\frac{P(Z^j|Y,\theta^{i+1})}{P(Z^j|Y,\theta^i)}P(Z^j|Y,\theta^i) \bigg)\\ &=logP(Z|Y,\theta^{i+1})=log1=0 \end{aligned} H(θi+1,θi)−H(θi,θi)=Zj∑(log(P(Zj∣Y,θi)P(Zj∣Y,θi+1)))P(Zj∣Y,θi)=log(Zj∑P(Zj∣Y,θi)P(Zj∣Y,θi+1)P(Zj∣Y,θi))=logP(Z∣Y,θi+1)=log1=0
其中不等式是由于Jensen不等式，由此证明了∑i=1MlogP(yi,θj+1)>∑i=1MlogP(yi,θj)\sum_{i=1}^{M}logP(y^{i},\theta^{j+1})>\sum_{i=1}^{M}logP(y^{i},\theta^{j})∑i=1MlogP(yi,θj+1)>∑i=1MlogP(yi,θj)，证明了EM算法的收敛性。但不能保证是全局最优，只能保证局部最优。

六.高斯混合分布

EM算法的一个重要应用场景就是高斯混合模型的参数估计。高斯混合模型就是由多个高斯模型组合在一起的混合模型（可以理解为多个高斯分布函数的线性组合，理论上高斯混合模型是可以拟合任意类型的分布），例如对于下图中的数据集如果用一个高斯模型来描述的话显然是不合理的：

两个高斯模型可以拟合数据集，如图所示：

如果有多个高斯模型，公式表示为：
P(y∣θ)=∑k=1Kakϕ(y∣θk)ϕ(y∣θk)=12πδkexp(−(y−μk)22δk2)ak>0,∑ak=1P(y|\theta)=\sum_{k=1}^{K}a_k\phi(y|\theta_{k}) \\ \phi(y|\theta_{k})=\frac{1}{\sqrt{2\pi}\delta_{k}}exp(-\frac{(y-\mu_{k})^2}{2 \delta_{k}^{2}}) \\ a_k>0,\sum a_k =1 P(y∣θ)=k=1∑Kakϕ(y∣θk)ϕ(y∣θk)=2πδk1exp(−2δk2(y−μk)2)ak>0,∑ak=1
ϕ(y∣θk)\phi(y|\theta_{k})ϕ(y∣θk)表示为第k个高斯分布密度模型，定义如上，其中aka_kak表示被选中的概率。在本次模型P(y∣θ)P(y|\theta)P(y∣θ)中，观测数据是已知的，而观测数据具体来自哪个模型是未知的，有点像之前提过的三硬币模型，我们来对比一下，A硬币就像是概率aka_kak，用来表明具体的模型，而B、C硬币就是具体的模型，只不过这里有很多个模型，不仅仅是B、C这两个模型。我们用γjk\gamma_{jk}γjk来表示，则：
γjk={1第j个观测数据来源于第k个模型0否则\gamma_{jk} = \begin{cases} 1& \text{第j个观测数据来源于第k个模型}\\ 0& \text{否则} \end{cases} γjk={10第j个观测数据来源于第k个模型否则
所以一个观测数据yjy_jyj的隐藏数据(γj1,γj2,...,γjk)(\gamma_{j1},\gamma_{j2},...,\gamma_{jk})(γj1,γj2,...,γjk)，那么完全似然函数就是：

P(y,γ∣θ)=∏k=1K∏j=1N[akϕ(y∣θk)]γjkP(y,\gamma|\theta)= \prod_{k=1}^{K}\prod_{j=1}^{N}[a_{k}\phi(y|\theta_{k})]^{\gamma_{jk}} P(y,γ∣θ)=k=1∏Kj=1∏N[akϕ(y∣θk)]γjk

取对数之后等于：

log(P(y,γ∣θ))=log(∏k=1K∏j=1N[akϕ(y∣θk)]γjk)=∑Kk=1(∑j=1N(γjk)log(ak)+∑j=1N(γjk)[log(12π)−log(δk)−(yi−μk)22δk2])\begin{aligned} log(P(y,\gamma|\theta))&=log( \prod_{k=1}^{K}\prod_{j=1}^{N}[a_{k}\phi(y|\theta_{k})]^{\gamma_{jk}})\\ &=\sum_{K}^{k=1}\bigg(\sum_{j=1}^{N}(\gamma_{jk}) log(a_k)+\sum_{j=1}^{N}( \gamma_{jk})\bigg[log(\frac{1}{\sqrt{2\pi}})-log(\delta_{k})-\frac{(y_i-\mu_{k})^2}{2 \delta_{k}^{2}}\bigg]\bigg) \end{aligned} log(P(y,γ∣θ))=log(k=1∏Kj=1∏N[akϕ(y∣θk)]γjk)=K∑k=1(j=1∑N(γjk)log(ak)+j=1∑N(γjk)[log(2π1)−log(δk)−2δk2(yi−μk)2])

E 步：
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])\begin{aligned} Q(\theta.\theta^i) &= E[log(P(y,\gamma|\theta))]\\ &=\sum_{K}^{k=1}\bigg(\sum_{j=1}^{N}(E\gamma_{jk}) log(a_k)+\sum_{j=1}^{N}(E\gamma_{jk})\bigg[log(\frac{1}{\sqrt{2\pi}})-log(\delta_{k})-\frac{(y_i-\mu_{k})^2}{2 \delta_{k}^{2}}\bigg]\bigg) \end{aligned} Q(θ.θi)=E[log(P(y,γ∣θ))]=K∑k=1(j=1∑N(Eγjk)log(ak)+j=1∑N(Eγjk)[log(2π1)−log(δk)−2δk2(yi−μk)2])
其中我们定义γjk^\hat{\gamma_{jk}}γjk^：
γjk^=E(γjk∣y,θ)=akϕ(yi∣θk)∑k=1Kakϕ(yi∣θk)j=1,2,..,N；k=1,2,...,Knk=∑j=iNEγjk\hat{\gamma_{jk}} = E(\gamma_{jk}|y,\theta)=\frac{a_k\phi(y_i|\theta_{k})}{\sum_{k=1}^{K}a_k\phi(y_i|\theta_{k}) }\\ j=1,2,..,N；k=1,2,...,K\\ n_k=\sum_{j=i}^{N}E\gamma_{jk} γjk^=E(γjk∣y,θ)=∑k=1Kakϕ(yi∣θk)akϕ(yi∣θk)j=1,2,..,N；k=1,2,...,Knk=j=i∑NEγjk
于是化简得到：
Q(θ.θi)=∑Kk=1(nklog(ak)+∑j=1N(Eγjk)[log(12π)−log(δk)−(yi−μk)22δk2])\begin{aligned} Q(\theta.\theta^i) &= \sum_{K}^{k=1}\bigg(n_k log(a_k)+\sum_{j=1}^{N}(E\gamma_{jk})\bigg[log(\frac{1}{\sqrt{2\pi}})-log(\delta_{k})-\frac{(y_i-\mu_{k})^2}{2 \delta_{k}^{2}}\bigg]\bigg) \end{aligned} Q(θ.θi)=K∑k=1(nklog(ak)+j=1∑N(Eγjk)[log(2π1)−log(δk)−2δk2(yi−μk)2])

E 步在代码设计上只有γjk^\hat{\gamma_{jk}}γjk^有用，用于M步的计算。

M步，
θi+1=argmax⁡θQ(θ,θi)\theta^{i+1}=arg \max_{\theta}Q(\theta,\theta^i) θi+1=argθmaxQ(θ,θi)
对Q(θ,θi)Q(\theta,\theta^i)Q(θ,θi)求导，得到每个未知量的偏导，使其偏导等于0，求解得到：
μk^=∑j=1Nγjk^yi∑j=1Nγjk^δk^=∑j=1Nγjk^(yi−μk)2∑j=1Nγjk^ak^=∑j=1Nγjk^N\hat{\mu_k}=\frac{\sum_{j=1}^{N}\hat{\gamma_{jk}}y_i}{\sum_{j=1}^{N}\hat{\gamma_{jk}}} \\ \\ \hat{\delta_k}=\frac{\sum_{j=1}^{N}\hat{\gamma_{jk}}(y_i-\mu_k)^2}{\sum_{j=1}^{N}\hat{\gamma_{jk}}} \\ \\ \\ \hat{a_k}=\frac{\sum_{j=1}^{N}\hat{\gamma_{jk}} }{N} μk^=∑j=1Nγjk^∑j=1Nγjk^yiδk^=∑j=1Nγjk^∑j=1Nγjk^(yi−μk)2ak^=N∑j=1Nγjk^
给一个初始值，来回迭代就可以求得值内容。这一块主要用到了Q(θ.θi)Q(\theta.\theta^i)Q(θ.θi)的导数，并且用到了E步的γjk^\hat{\gamma_{jk}}γjk^。

七.总结

这里其实还有很多问题没讲，大家想了解的可以去学习统计学习方法这本书，讲解的还是挺全的。

八.高斯混合模型代码实验

EM算法更多是一种思想，用概率来解决问题的一种方法，具体的代码看自己选用模型，所以并没有通用的模型，本此代码主要是讲解混合高斯分布模型的

这其中的M步完全按照了公式来计算。

import numpy as np
import random
import math
import time

'''
数据集：伪造数据集（两个高斯分布混合）
数据集长度：1000
------------------------------
运行结果：
----------------------------
the Parameters set is:
alpha0:0.3, mu0:0.7, sigmod0:-2.0, alpha1:0.5, mu1:0.5, sigmod1:1.0
----------------------------
the Parameters predict is:
alpha0:0.7, mu0:0.6, sigmod0:-1.7, alpha1:0.7, mu1:0.7, sigmod1:0.9
----------------------------
'''def loadData(mu0, sigma0, mu1, sigma1, alpha0, alpha1):'''初始化数据集这里通过服从高斯分布的随机函数来伪造数据集:param mu0: 高斯0的均值:param sigma0: 高斯0的方差:param mu1: 高斯1的均值:param sigma1: 高斯1的方差:param alpha0: 高斯0的系数:param alpha1: 高斯1的系数:return: 混合了两个高斯分布的数据'''# 定义数据集长度为1000length = 1000# 初始化第一个高斯分布，生成数据，数据长度为length * alpha系数，以此来# 满足alpha的作用data0 = np.random.normal(mu0, sigma0, int(length * alpha0))# 第二个高斯分布的数据data1 = np.random.normal(mu1, sigma1, int(length * alpha1))# 初始化总数据集# 两个高斯分布的数据混合后会放在该数据集中返回dataSet = []# 将第一个数据集的内容添加进去dataSet.extend(data0)# 添加第二个数据集的数据dataSet.extend(data1)# 对总的数据集进行打乱（其实不打乱也没事，只不过打乱一下直观上让人感觉已经混合了# 读者可以将下面这句话屏蔽以后看看效果是否有差别）random.shuffle(dataSet)#返回伪造好的数据集return dataSet


# 高斯分布公式，没有什么特殊的
def calcGauss(dataSetArr, mu, sigmod):'''根据高斯密度函数计算值依据：“9.3.1 高斯混合模型” 式9.25注：在公式中y是一个实数，但是在EM算法中(见算法9.2的E步)，需要对每个j都求一次yjk，在本实例中有1000个可观测数据，因此需要计算1000次。考虑到在E步时进行1000次高斯计算，程序上比较不简洁，因此这里的y是向量，在numpy的exp中如果exp内部值为向量，则对向量中每个值进行exp，输出仍是向量的形式。所以使用向量的形式1次计算即可将所有计算结果得出，程序上较为简洁:param dataSetArr: 可观测数据集:param mu: 均值:param sigmod: 方差:return: 整个可观测数据集的高斯分布密度（向量形式）'''# 计算过程就是依据式9.25写的，没有别的花样result = (1 / (math.sqrt(2*math.pi)*sigmod**2)) * np.exp(-1 * (dataSetArr-mu) * (dataSetArr-mu) / (2*sigmod**2))# 返回结果return resultdef E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1):'''EM算法中的E步依据当前模型参数，计算分模型k对观数据y的响应度:param dataSetArr: 可观测数据y:param alpha0: 高斯模型0的系数:param mu0: 高斯模型0的均值:param sigmod0: 高斯模型0的方差:param alpha1: 高斯模型1的系数:param mu1: 高斯模型1的均值:param sigmod1: 高斯模型1的方差:return: 两个模型各自的响应度'''# 计算y0的响应度# 先计算模型0的响应度的分子gamma0 = alpha0 * calcGauss(dataSetArr, mu0, sigmod0)#print("gamma0=",gamma0.shape) # 1000, 维向量# 模型1响应度的分子gamma1 = alpha1 * calcGauss(dataSetArr, mu1, sigmod1)# 两者相加为E步中的分布sum = gamma0 + gamma1# 各自相除，得到两个模型的响应度gamma0 = gamma0 / sumgamma1 = gamma1 / sum# 返回两个模型响应度return gamma0, gamma1def M_step(muo, mu1, gamma0, gamma1, dataSetArr):# 依据算法9.2计算各个值# 这里没什么花样，对照书本公式看看这里就好了# np.dot 点积：[1,2] [2,3] = [2,6]mu0_new = np.dot(gamma0, dataSetArr) / np.sum(gamma0)mu1_new = np.dot(gamma1, dataSetArr) / np.sum(gamma1)# math.sqrt  平方根 sigmod0_new = math.sqrt(np.dot(gamma0, (dataSetArr - muo)**2) / np.sum(gamma0))sigmod1_new = math.sqrt(np.dot(gamma1, (dataSetArr - mu1)**2) / np.sum(gamma1))alpha0_new = np.sum(gamma0) / len(gamma0)alpha1_new = np.sum(gamma1) / len(gamma1)# 将更新的值返回return mu0_new, mu1_new, sigmod0_new, sigmod1_new, alpha0_new, alpha1_new## 训练主函数
def EM_Train(dataSetList, iter=500):'''根据EM算法进行参数估计算法依据“9.3.2 高斯混合模型参数估计的EM算法” 算法9.2:param dataSetList:数据集（可观测数据）:param iter: 迭代次数:return: 估计的参数'''# 将可观测数据y转换为数组形式，主要是为了方便后续运算dataSetArr = np.array(dataSetList)# 步骤1：对参数取初值，开始迭代alpha0 = 0.5mu0 = 0sigmod0 = 1alpha1 = 0.5mu1 = 1sigmod1 = 1# 开始迭代step = 0while (step < iter):# 每次进入一次迭代后迭代次数加1step += 1# 步骤2：E步：依据当前模型参数，计算分模型k对观测数据y的响应度gamma0, gamma1 = E_step(dataSetArr, alpha0, mu0, sigmod0, alpha1, mu1, sigmod1)# 步骤3：M步mu0, mu1, sigmod0, sigmod1, alpha0, alpha1 = M_step(mu0, mu1, gamma0, gamma1, dataSetArr)# 迭代结束后将更新后的各参数返回return alpha0, mu0, sigmod0, alpha1, mu1, sigmod1

if __name__ == '__main__':start = time.time()# 设置两个高斯模型进行混合，这里是初始化两个模型各自的参数# 见“9.3 EM算法在高斯混合模型学习中的应用”# alpha是“9.3.1 高斯混合模型” 定义9.2中的系数α# mu0是均值μ# sigmod是方差σ# 在设置上两个alpha的和必须为1，其他没有什么具体要求，符合高斯定义就可以alpha0 = 0.3  # 系数αmu0 = -2  # 均值μsigmod0 = 0.5  # 方差σalpha1 = 0.7  # 系数αmu1 = 0.5  # 均值μsigmod1 = 1  # 方差σ# 初始化数据集dataSetList = loadData(mu0, sigmod0, mu1, sigmod1, alpha0, alpha1)#打印设置的参数print('---------------------------')print('the Parameters set is:')print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f' % (alpha0, alpha1, mu0, mu1, sigmod0, sigmod1))# 开始EM算法，进行参数估计alpha0, mu0, sigmod0, alpha1, mu1, sigmod1 = EM_Train(dataSetList)# 打印参数预测结果print('----------------------------')print('the Parameters predict is:')print('alpha0:%.1f, mu0:%.1f, sigmod0:%.1f, alpha1:%.1f, mu1:%.1f, sigmod1:%.1f' % (alpha0, alpha1, mu0, mu1, sigmod0, sigmod1))# 打印时间print('----------------------------')print('time span:', time.time() - start)

#the Parameters predict is:
#alpha0:0.4, mu0:0.6, sigmod0:-1.9, alpha1:0.7, mu1:0.6, #sigmod1:0.8

九.关于EM应用到聚类原理

EM算法一般用于聚类，也就是无监督模型里面，因为无监督学习没有标签（即y值)，EM算法可以先给无监督学习估计一个隐状态（即标签），有了标签，算法模型就可以转换成有监督学习，这时就可以用极大似然估计法求解出模型最优参数。其中估计隐状态流程应为EM算法的E步，后面用极大似然估计为M步。
相比于 K-Means 算法，EM 聚类更加灵活，比如下面这两种情况，K-Means 会得到下面的聚类结果。

因为 K-Means 是通过距离来区分样本之间的差别的，且每个样本在计算的时候只能属于一个分类，称之为是硬聚类算法。而 EM 聚类在求解的过程中，实际上每个样本都有一定的概率和每个聚类相关，叫做软聚类算法。
你可以把 EM 算法理解成为是一个框架，在这个框架中可以采用不同的模型来用 EM 进行求解。常用的 EM 聚类有 GMM 高斯混合模型和 HMM 隐马尔科夫模型。GMM（高斯混合模型）聚类就是 EM 聚类的一种。比如上面这两个图，可以采用 GMM 来进行聚类。
和 K-Means 一样，我们事先知道聚类的个数，但是不知道每个样本分别属于哪一类。通常，我们可以假设样本是符合高斯分布的（也就是正态分布）。每个高斯分布都属于这个模型的组成部分（component），要分成 K 类就相当于是 K 个组成部分。这样我们可以先初始化每个组成部分的高斯分布的参数，然后再看来每个样本是属于哪个组成部分。这也就是 E 步骤。

再通过得到的这些隐含变量结果，反过来求每个组成部分高斯分布的参数，即 M 步骤。反复 EM 步骤，直到每个组成部分的高斯分布参数不变为止。

这样也就相当于将样本按照 GMM 模型进行了 EM 聚类。

所以说很多KMeans解决不了的问题，EM聚类是可以解决的。在 EM 框架中，我们将潜在类别当做隐藏变量，样本看做观察值，把聚类问题转化为参数估计问题，最终把样本进行聚类。
最后再多啰嗦一句，EM 算法相当于一个框架，你可以采用不同的模型来进行聚类，比如 GMM（高斯混合模型），或者 HMM（隐马尔科夫模型）来进行聚类。
参考：公众号https://mp.weixin.qq.com/s/oidooJb3PbbS-FSGj_aFrw

十.参考

datawhale 学习内容
主要参考统计学习方法这本书，书籍地址：http://www.dgt-factory.com/uploads/2018/07/0725/%E7%BB%9F%E8%AE%A1%E5%AD%A6%E4%B9%A0%E6%96%B9%E6%B3%95.pdf。

统计学习方法-代码解读

EM算法 - 期望极大算法