0、前言

EM算法，也叫最大期望算法，或者是期望最大化算法，是机器学习十大算法之一，它很简单，但是也同样很有深度。

• 简单是因为它就分两步求解问题：

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

• 深度在于它的数学推理涉及到比较繁杂的概率公式等。

所以本文会介绍很多概率方面的知识，不懂的同学可以先去了解一些知识，当然本文也会尽可能的讲解清楚这些知识，讲的不好的地方麻烦大家评论指出，后续不断改进完善。

1、EM算法引入

  概率模型有时候既含有观测变量，又含有隐变量或潜在变量，如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计方法估计模型参数，但是，当模型含有隐变量时，就不能简单的使用这些方法，EM算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验概率估计法，我们讨论极大似然估计，极大后验概率估计与其类似。

参考统计学习方法书中的一个例子来引入EM算法：

  假设有3枚硬币，分别记做A、B、C，这些硬币正面出现的概率分别是π\piπ、ppp、qqq，进行如下实验：

• 先掷硬币A，根据硬币A的结果选出硬币B和硬币C，正面选硬币B，反面选硬币C；
• 通过选择出的硬币（B或者C），记录掷硬币的结果出现正面为1，反面为0。

  也就是说一共需要掷两次硬币，第一次为硬币A，第二次硬币需要根据A的结果来选择，只记录一次结果，就是第二次掷的硬币的正反面。把这整个过程看作一次实验，如此独立地重复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

  假设现在我们只知道到记录掷硬币的结果，也就是上面的10次结果，不知道每次实验掷硬币的过程，也就是不知道每次实验第一次掷硬币A的结果；现在需要估计每次实验结果出现正面的概率，怎么求呢？

我们来构建这样一个三硬币模型：

硬币A的结果	第二次的硬币	第二次的结果	概率
正面	硬币B	正面	πp\pi pπp
正面	硬币B	反面	π(1−p)\pi (1-p)π(1−p)
反面	硬币C	正面	πq\pi qπq
反面	硬币C	反面	π(1−q)\pi(1-q)π(1−q)

  记每次实验记录的结果为yyy，则第二次硬币为B的概率为π\piπ，则B的结果概率为πpy(1−p)(1−y)\pi p^y(1-p)^{(1-y)}πpy(1−p)(1−y)，当出现正面，就令y=1y=1y=1，概率就为πp\pi pπp。对C同理。计算模型的概率P(y∣θ)P(y|\theta)P(y∣θ)，其中θ=(π,p,q)\theta=(\pi,p,q)θ=(π,p,q)是模型的参数，yyy受到硬币A的结果的影响，记硬币A的结果为zzz，则P(y∣θ)P(y|\theta)P(y∣θ)是联合概率P(y,z∣θ)P(y,z|\theta)P(y,z∣θ)的边际概率，即P(y∣θ)=∑zP(y,z∣θ)P(y|\theta)=\sum_zP(y,z|\theta)P(y∣θ)=∑z​P(y,z∣θ)，所以：

P(y∣θ)=∑zP(y,z∣θ)=∑zP(z∣θ)P(y∣z,θ)=P(z=正面∣θ)P(y∣z=正面,θ)+P(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) \\ &=P(z=正面|\theta)P(y|z=正面,\theta)+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,θ)=P(z=正面∣θ)P(y∣z=正面,θ)+P(z=反面∣θ)P(y∣z=反面,θ)=πpy(1−p)(1−y)+(1−π)qy(1−q)(1−y)​

• 若y=1y=1y=1，表示这此看到的是正面，这个正面有可能是B的正面，也可能是C的正面，则P(y=1∣θ)=πp+(1−π)qP(y=1|\theta)=\pi p+(1-\pi)qP(y=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)

  yyy是观测变量，表示一次观测结果是111或000，zzz是隐藏变量，表示掷硬币A的结果，这个是观测不到结果的，θ=(π,p,q)\theta=(\pi,p,q)θ=(π,p,q)表示模型参数。若考虑nnn次实验，将观测数据表示为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=0n[π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}^n[ \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∏n​[π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θmax​logP(Y∣θ)
  这个问题没有解析解，只有通过迭代方法来求解，EM算法就是可以用于求解这个问题的一种迭代算法，下面给出EM算法的迭代过程：

• 首先选取初始值，记做θ0=(π0,p0,q0)\theta^{0}=(\pi^{0},p^{0},q^{0})θ0=(π0,p0,q0)，第iii次的迭代参数的估计值为θ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的概率：
  P(zi=πi∣yi,θi)=μi+1=P(zi=πi,yi∣θi)P(yi∣θi)=P(zi=πi,yi∣θi)∑ziP(zi,yi∣θi)=P(zi=πi,yi∣θi)P(zi=πi,yi∣θi)+P(zi=(1−πi),yi∣θi)=πi(pi)yi(1−pi)1−yiπi(pi)yi(1−pi)1−yi+(1−πi)(qi)yi(1−pi)1−yi\begin{aligned} P(z_i=\pi^i|y_i,\theta^i)=\mu^{i+1}&=\frac{P(z_i=\pi^i,y_i|\theta^i)}{P(y_i|\theta^i)}\\ &=\frac{P(z_i=\pi^i,y_i|\theta^i)}{\sum\limits_{z_i}P(z_i,y_i|\theta^i)}\\ &=\frac{P(z_i=\pi^i,y_i|\theta^i)}{P(z_i=\pi^i,y_i|\theta^i)+P(z_i=(1-\pi^i),y_i|\theta^i)}\\ &=\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}} \end{aligned} P(zi​=πi∣yi​,θi)=μi+1​=P(yi​∣θi)P(zi​=πi,yi​∣θi)​=zi​∑​P(zi​,yi​∣θi)P(zi​=πi,yi​∣θi)​=P(zi​=πi,yi​∣θi)+P(zi​=(1−πi),yi​∣θi)P(zi​=πi,yi​∣θi)​=π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硬币的概率P(Y,Z∣θ)P(Y,Z|\theta)P(Y,Z∣θ)。

• M步：计算模型参数的新估计值：

πi+1=1n∑j=1nμji+1\pi^{i+1}=\frac{1}{n}\sum_{j=1}^{n}\mu_{j}^{i+1} πi+1=n1​j=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\limits_{j=1}^{n}\mu_{j}^{i+1}y_j}{\sum\limits_{j=1}^{n}\mu_{j}^{i+1}} pi+1=j=1∑n​μji+1​j=1∑n​μji+1​yj​​

分子乘以yiy_{i}yi​，所以其实是计算B硬币出现正面的概率。

qi+1=∑j=1n(1−μji+1)yj∑j=1n(1−μji+1)q^{i+1}=\frac{\sum\limits_{j=1}^{n}(1-\mu_{j}^{i+1})y_j}{\sum\limits_{j=1}^{n}(1-\mu_{j}^{i+1})} qi+1=j=1∑n​(1−μji+1​)j=1∑n​(1−μji+1​)yj​​

其中，(1−μji+1)(1-\mu_{j}^{i+1})(1−μji+1​)表示出现C硬币的概率，即A出现反面的概率。

  闭环形成，从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硬币抛，这时候概率模型就存在着隐含变量！

2、具体的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为第iii次迭代参数θ\thetaθ的估计值，在第i+1i+1i+1次迭代的E步，计算Q(θ,θi)Q(\theta,\theta^i)Q(θ,θi)函数：

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)是在给定观测数据YYY和当前的参数估计θi\theta^iθi下隐变量数据ZZZ的条件概率分布；

• (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θmax​Q(θ,θi)

其中，Q(θ,θi)Q(\theta,\theta^{i})Q(θ,θi)是EM算法的核心，称为Q函数(Q function)，这个是需要自己构造的，可以看出这是一个期望，而M步是极大化这个期望值，这也是EM（Expectation Maximization ）算法名字的来源。

• (4) 重复第(2)步和第(3)步，直到收敛，收敛条件：

∣∣θi+1−θi∣∣<ε1或者：∣∣Q(θi+1,θi)−Q(θi,θi)∣∣<ε2|| \theta^{i+1}-\theta^{i} || < \varepsilon_1 \qquad或者：\qquad||Q(\theta^{i+1},\theta^{i})-Q(\theta^{i},\theta^{i})|| <\varepsilon_2 ∣∣θi+1−θi∣∣<ε1​或者：∣∣Q(θi+1,θi)−Q(θi,θi)∣∣<ε2​
收敛迭代就结束了。我们来拆解一下这个M步骤，

3、EM算法推导

3.1 Jensen不等式

  Jensen不等式，这个公式在算法推导和证明收敛都能用到，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))

3.2 EM推导

  如果模型不包含隐变量，我们在模型学习的时候，直接使用极大似然估计使P(Y∣θ)P(Y|\theta)P(Y∣θ)最大，其中Y是观测数据，θ\thetaθ是模型参数，意思就是在参数θ\thetaθ的情况下，观测数据Y出现的概率最大。然而当我们面对一个含有隐变量的概率模型时，我们也想使在参数θ\thetaθ的情况下，观测数据Y出现的概率最大，即最大化P(Y∣θ)P(Y|\theta)P(Y∣θ)。但是我们不能直接使用极大似然估计，因为参数θ\thetaθ包含有隐变量的信息，而Y与隐变量相关，估计Y需要知道Z，估计Z需要知道Y，形成一个死循环。可以将隐变量分离出来，对于每一个样本，我们需要确定Z来使联合概率P(Y,Z∣θ)P(Y,Z|\theta)P(Y,Z∣θ)最大，由条件概率公式：P(Y,Z∣θ)=P(Y∣Z,θ)P(Z∣θ)P(Y,Z|\theta)=P(Y|Z,\theta)P(Z|\theta)P(Y,Z∣θ)=P(Y∣Z,θ)P(Z∣θ)：

P(Y∣θ)=∑ZP(Y,Z∣θ)=∑ZP(Y∣Z,θ)P(Z∣θ)P(Y|\theta)=\sum\limits_ZP(Y,Z|\theta)=\sum\limits_ZP(Y|Z,\theta)P(Z|\theta)P(Y∣θ)=Z∑​P(Y,Z∣θ)=Z∑​P(Y∣Z,θ)P(Z∣θ)

  然后就可以使用极大似然估计使P(Y∣θ)P(Y|\theta)P(Y∣θ)出现的概率最大，为了方便计算，对似然函数取对数。所以我们的目标是极大化观测数据Y关于参数θ\thetaθ的对数似然函数L(θ)L(\theta)L(θ):

L(θ)=log⁡P(Y∣θ)=log⁡∑ZP(Y∣Z,θ)P(Z∣θ)L(\theta)=\log P(Y|\theta)=\log\sum\limits_ZP(Y|Z,\theta)P(Z|\theta)L(θ)=logP(Y∣θ)=logZ∑​P(Y∣Z,θ)P(Z∣θ)

  想要直接极大化上面的式子，可以对Y，Z求偏导，但是log里面连加求导，会很麻烦，可以自己想象。可以看Jensen不等式正好将求和和函数连在一起，但是必须要使求和的那个系数之和为1，刚好和概率之和为1对上了呀。构造一个概率之和，而∑ZP(Z∣Y,θ)=1\sum\limits_ZP(Z|Y,\theta)=1Z∑​P(Z∣Y,θ)=1，刚好可以利用Jensen不等式。

  根据EM算法的定义，我们需要通过迭代的方式逐渐去极大化L(θ)L(\theta)L(θ)，每一次更新θ\thetaθ时，都能使L(θ)L(\theta)L(θ)增大，也就是假设当前第iii次迭代的参数估计值为θi\theta^iθi，则在这i+1i+1i+1次迭代中，我们需要重新估计参数值，使L(θ)L(\theta)L(θ)增大，即L(θ)>L(θi)L(\theta)>L(\theta^i)L(θ)>L(θi)。则L(θ)−L(θi)L(\theta)-L(\theta^{i})L(θ)−L(θi)：

L(θ)−L(θi)=log⁡∑ZP(Y∣Z,θ)P(Z∣θ)−log⁡P(Y∣θi)(构造Jensen不等式←)=log(∑ZP(Y∣Z,θ)P(Z∣θ)P(Z∣Y,θi)P(Z∣Y,θi))−log⁡P(Y∣θi)(∑ZP(Z∣Y,θi)=1←)=log(∑ZP(Z∣Y,θi)P(Y∣Zθ)P(Z∣θ)P(Z∣Y,θi))−log⁡P(Y∣θi)(log是凸函数←)≥∑ZP(Z∣Y,θi)log(P(Y∣Z,θ)P(Z∣θ)P(Z∣Y,θi))−log⁡P(Y∣θi)\begin{aligned} L(\theta)-L(\theta^{i})&=\log\sum\limits_ZP(Y|Z,\theta)P(Z|\theta)-\log P(Y|\theta^i)\\ (构造Jensen不等式\leftarrow)&=log(\sum_{Z} P(Y|Z,\theta)P(Z|\theta)\frac{P(Z|Y,\theta^i)}{P(Z|Y,\theta^i)})-\log P(Y|\theta^i)\\ (\sum\limits_ZP(Z|Y,\theta^i)=1\leftarrow)&=log(\sum_{Z}P(Z|Y,\theta^i)\frac{ P(Y|Z\theta)P(Z|\theta)}{P(Z|Y,\theta^i)})-\log P(Y|\theta^i)\\ (log是凸函数\leftarrow)&\ge \sum_{Z} P(Z|Y,\theta^i)log(\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^i)})-\log P(Y|\theta^i) \end{aligned} L(θ)−L(θi)(构造Jensen不等式←)(Z∑​P(Z∣Y,θi)=1←)(log是凸函数←)​=logZ∑​P(Y∣Z,θ)P(Z∣θ)−logP(Y∣θi)=log(Z∑​P(Y∣Z,θ)P(Z∣θ)P(Z∣Y,θi)P(Z∣Y,θi)​)−logP(Y∣θi)=log(Z∑​P(Z∣Y,θi)P(Z∣Y,θi)P(Y∣Zθ)P(Z∣θ)​)−logP(Y∣θi)≥Z∑​P(Z∣Y,θi)log(P(Z∣Y,θi)P(Y∣Z,θ)P(Z∣θ)​)−logP(Y∣θi)​

这一个步骤就是采用了凸函数的Jensen不等式做转换。因为ZZZ是隐藏变量，所以有：
∑ZP(Z∣Y,θi)=1，P(Z∣Y,θi)>0\sum_{Z} P(Z|Y,\theta^i)=1，P(Z|Y,\theta^i)>0Z∑​P(Z∣Y,θi)=1，P(Z∣Y,θi)>0

于是继续变：

L(θ)−L(θi)=log(∑ZP(Z∣Y,θi)P(Y∣Z,θ)P(Z∣θ)P(Z∣Y,θi))−log⁡P(Y∣θi)≥∑ZP(Z∣Y,θi)log(P(Y∣Z,θ)P(Z∣θ)P(Z∣Y,θi))−log⁡P(Y∣θi)=∑ZP(Z∣Y,θi)log(P(Y∣Z,θ)P(Z∣θ)P(Z∣Y,θi))−∑ZP(Z∣Y,θi)log⁡P(Y∣θi)=∑ZP(Z∣Y,θi)[logP(Y∣Z,θ)P(Z∣θ)P(Z∣Y,θi)−logP(Y∣θi)]=∑ZP(Z∣Y,θi)logP(Y∣Z,θ)P(Z∣θ)P(Z∣Y,θi)P(Y∣θi)\begin{aligned} L(\theta)-L(\theta^{i})&=log(\sum_{Z}P(Z|Y,\theta^i)\frac{ P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^i)})-\log P(Y|\theta^i)\\ &\ge \sum_{Z} P(Z|Y,\theta^i)log(\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^i)})-\log P(Y|\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)\log P(Y|\theta^i)\\ &= \sum_{Z} P(Z|Y,\theta^i)[log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^i)}-logP(Y|\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})} \end{aligned} L(θ)−L(θi)​=log(Z∑​P(Z∣Y,θi)P(Z∣Y,θi)P(Y∣Z,θ)P(Z∣θ)​)−logP(Y∣θi)≥Z∑​P(Z∣Y,θi)log(P(Z∣Y,θi)P(Y∣Z,θ)P(Z∣θ)​)−logP(Y∣θi)=Z∑​P(Z∣Y,θi)log(P(Z∣Y,θi)P(Y∣Z,θ)P(Z∣θ)​)−Z∑​P(Z∣Y,θi)logP(Y∣θi)=Z∑​P(Z∣Y,θi)[logP(Z∣Y,θi)P(Y∣Z,θ)P(Z∣θ)​−logP(Y∣θi)]=Z∑​P(Z∣Y,θi)logP(Z∣Y,θi)P(Y∣θi)P(Y∣Z,θ)P(Z∣θ)​​
也就是：L(θ)≥L(θi)+∑ZP(Z∣Y,θi)log(P(Y∣Z,θ)P(Z,θ)P(Y∣Z,θi)P(Y∣θi))L(\theta)\ge L(\theta^{i})+ \sum\limits_{Z} P(Z|Y,\theta^i)log(\frac{P(Y|Z,\theta)P(Z,\theta)}{P(Y|Z,\theta^i) P(Y|\theta^{i})})L(θ)≥L(θi)+Z∑​P(Z∣Y,θi)log(P(Y∣Z,θi)P(Y∣θi)P(Y∣Z,θ)P(Z,θ)​)，有下界（不等式右边就是下界），任何使下界增大的θ\thetaθ都可以使L(θ)L(\theta)L(θ)增大，因此我们需要最大化下界，来得到近似的极大值。

第i+1i+1i+1次迭代的参数估计值θi+1\theta^{i+1}θi+1，就是使下界最大化的θ\thetaθ：

θi+1=L(θi)+∑ZP(Z∣Y,θi)log(P(Y∣Z,θ)P(Z∣θ)P(Y∣Z,θi)P(Y∣θi))=L(θi)+∑ZP(Z∣Y,θi)logP(Y∣Z,θ)P(Z∣θ)P(Y∣Z,θi)−L(θi)=argmax⁡θ∑ZP(Z∣Y,θi)logP(Y∣Z,θ)P(Z∣θ)P(Z∣Y,θi)=argmax⁡θ∑ZP(Z∣Y,θi)logP(Y∣Z,θ)P(Z∣θ)−∑ZP(Z∣Y,θi)logP(Z∣Y,θi))=argmax⁡θ∑ZP(Z∣Y,θi)logP(Y∣Z,θ)P(Z∣θ)=argmax⁡θ∑ZP(Z∣Y,θi)logP(Y,Z∣θ)=argmax⁡θQ(θ,θi)\begin{aligned} \theta^{i+1}&=L(\theta^{i})+ \sum\limits_{Z} P(Z|Y,\theta^i)log(\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Y|Z,\theta^i) P(Y|\theta^{i})})\\ &=L(\theta^{i})+ \sum\limits_{Z} P(Z|Y,\theta^i)log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Y|Z,\theta^i)}-L(\theta^i)\\ &=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)logP(Y|Z,\theta)P(Z|\theta)-\sum\limits_ZP(Z|Y,\theta^i)log{P(Z|Y,\theta^i)})\\ &=arg \max_{\theta}\sum_{Z} P(Z|Y,\theta^i)logP(Y|Z,\theta)P(Z|\theta)\\ & =arg \max_{\theta}\sum_{Z} P(Z|Y,\theta^i)logP(Y,Z|\theta)\\ &=arg \max_{\theta}Q(\theta,\theta^i) \end{aligned} θi+1​=L(θi)+Z∑​P(Z∣Y,θi)log(P(Y∣Z,θi)P(Y∣θi)P(Y∣Z,θ)P(Z∣θ)​)=L(θi)+Z∑​P(Z∣Y,θi)logP(Y∣Z,θi)P(Y∣Z,θ)P(Z∣θ)​−L(θi)=argθmax​Z∑​P(Z∣Y,θi)logP(Z∣Y,θi)P(Y∣Z,θ)P(Z∣θ)​=argθmax​Z∑​P(Z∣Y,θi)logP(Y∣Z,θ)P(Z∣θ)−Z∑​P(Z∣Y,θi)logP(Z∣Y,θi))=argθmax​Z∑​P(Z∣Y,θi)logP(Y∣Z,θ)P(Z∣θ)=argθmax​Z∑​P(Z∣Y,θi)logP(Y,Z∣θ)=argθmax​Q(θ,θi)​

其中logloglog分母提出来是关于ZZZ的∑ZP(Z∣Y,θi)logP(Z∣Y,θi)\sum\limits_{Z} P(Z|Y,\theta^i)logP(Z|Y,\theta^i)Z∑​P(Z∣Y,θi)logP(Z∣Y,θi)，是一个常数可以去掉。其中Q(θ,θi)=∑ZP(Z∣Y,θi)logP(Y,Z∣θ)Q(\theta,\theta^i)=\sum\limits_{Z} P(Z|Y,\theta^i)logP(Y,Z|\theta)Q(θ,θi)=Z∑​P(Z∣Y,θi)logP(Y,Z∣θ)是EM算法中非常重要的Q函数，很多时候写成期望的形式：

Q(θ,θi)=∑ZP(Z∣Y,θi)logP(Y,Z∣θ)=EZ[logP(Y,Z∣θ)∣Y,θi]=EZ∣Y,θilogP(Y,Z∣θ)\begin{aligned}Q(\theta,\theta^i)&=\sum\limits_{Z} P(Z|Y,\theta^i)logP(Y,Z|\theta)\\ &=E_Z[logP(Y,Z|\theta)|Y,\theta^i]\\ &=E_{Z|Y,\theta^i}logP(Y,Z|\theta) \end{aligned}Q(θ,θi)​=Z∑​P(Z∣Y,θi)logP(Y,Z∣θ)=EZ​[logP(Y,Z∣θ)∣Y,θi]=EZ∣Y,θi​logP(Y,Z∣θ)​

记B(θ,θi)=L(θi)+∑ZP(Z∣Y,θi)log(P(Y∣Z,θ)P(Z,θ)P(Y∣Z,θi)P(Y∣θi))B(\theta,\theta^i)=L(\theta^i)+ \sum\limits_{Z} P(Z|Y,\theta^i)log(\frac{P(Y|Z,\theta)P(Z,\theta)}{P(Y|Z,\theta^i) P(Y|\theta^{i})})B(θ,θi)=L(θi)+Z∑​P(Z∣Y,θi)log(P(Y∣Z,θi)P(Y∣θi)P(Y∣Z,θ)P(Z,θ)​)，也就是L(θ)L(\theta)L(θ)的下界，可以看出L(θ)≥B(θ,θi)L(\theta)\geq B(\theta,\theta^i)L(θ)≥B(θ,θi)，当θ=θi\theta=\theta^iθ=θi时，等号成立。函数B(θ,θi)B(\theta,\theta^i)B(θ,θi)的增加会导致L(θ)L(\theta)L(θ)在每次迭代都是增加的，我们想尽快迭代到最优值，就要求在每一次迭代的时候都取到函数B(θ,θi)B(\theta,\theta^i)B(θ,θi)的最大值，也就是θi+1\theta^{i+1}θi+1，也就是Q函数的最大值。然后基于θi+1\theta^{i+1}θi+1，重新计算B(θ,θi+1)B(\theta,\theta^{i+1})B(θ,θi+1)，会得到类似的图，继续寻找最大值θi+2\theta^{i+2}θi+2，以此类推下去。在这个过程中，对数似然函数L(θ)L(\theta)L(θ)不断增大，直到达到全局最优值。

3.3 EM算法的收敛性

  EM算法的最大优点是简单性和普适性，但是EM算法得到的估计序列是否收敛呢？是收敛到全局最大值还是局部最大值？

定理一：设P(Y∣θ)P(Y|\theta)P(Y∣θ)是观测数据的似然函数，θi(i=1,2,...)\theta^i(i=1,2,...)θi(i=1,2,...)是EM算法得到的参数估计序列，P(Y∣θi)(i=1,2,...)P(Y|\theta^i)(i=1,2,...)P(Y∣θi)(i=1,2,...)是对应的似然函数序列，则P(Y∣θi)P(Y|\theta^i)P(Y∣θi)是单调递增的，即：

P(Y∣θi+1)≥P(Y∣θi)P(Y|\theta^{i+1})\geq P(Y|\theta^i)P(Y∣θi+1)≥P(Y∣θi)

  证明过程：
由条件概率P(Y,Z∣θ)=P(Y∣θ)P(Z∣Y,θ)P(Y,Z|\theta)=P(Y|\theta)P(Z|Y,\theta)P(Y,Z∣θ)=P(Y∣θ)P(Z∣Y,θ)可得：

P(Y∣θ)=P(Y,Z∣θ)P(Z∣Y,θ)P(Y|\theta)=\frac{P(Y,Z|\theta)}{P(Z|Y,\theta)}P(Y∣θ)=P(Z∣Y,θ)P(Y,Z∣θ)​

取对数有：

logP(Y∣θ)=logP(Y,Z∣θ)−logP(Z∣Y,θ)logP(Y|\theta)=logP(Y,Z|\theta)-logP(Z|Y,\theta)logP(Y∣θ)=logP(Y,Z∣θ)−logP(Z∣Y,θ)

由上文可知:

Q(θ,θi)=∑ZP(Z∣Y,θi)logP(Y,Z∣θ)Q(\theta,\theta^i)=\sum\limits_ZP(Z|Y,\theta^i)logP(Y,Z|\theta)Q(θ,θi)=Z∑​P(Z∣Y,θi)logP(Y,Z∣θ)

令:

H(θ,θi)=∑ZP(Z∣Y,θi)logP(Z∣Y,θ)H(\theta,\theta^i)=\sum\limits_ZP(Z|Y,\theta^i)logP(Z|Y,\theta)H(θ,θi)=Z∑​P(Z∣Y,θi)logP(Z∣Y,θ)

Q(θ,θi)−H(θ,θi)=∑ZP(Z∣Y,θi)logP(Y,Z∣θ)−∑ZP(Z∣Y,θi)logP(Z∣Y,θ)=∑ZP(Z∣Y,θi)logP(Y,Z∣θ)P(Z∣Y,θ)=∑ZP(Z∣Y,θi)logP(Y∣θ)=logP(Y∣θ)\begin{aligned} Q(\theta,\theta^i)-H(\theta,\theta^i)&=\sum\limits_ZP(Z|Y,\theta^i)logP(Y,Z|\theta)-\sum\limits_ZP(Z|Y,\theta^i)logP(Z|Y,\theta)\\ &=\sum_ZP(Z|Y,\theta^i)log\frac{P(Y,Z|\theta)}{P(Z|Y,\theta)}\\ &=\sum_ZP(Z|Y,\theta^i)logP(Y|\theta)\\ &=logP(Y|\theta) \end{aligned}Q(θ,θi)−H(θ,θi)​=Z∑​P(Z∣Y,θi)logP(Y,Z∣θ)−Z∑​P(Z∣Y,θi)logP(Z∣Y,θ)=Z∑​P(Z∣Y,θi)logP(Z∣Y,θ)P(Y,Z∣θ)​=Z∑​P(Z∣Y,θi)logP(Y∣θ)=logP(Y∣θ)​

则：

logP(Y∣θi+1)−logP(Y∣θi)=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)]\begin{aligned} logP(Y|\theta^{i+1})-logP(Y|\theta^i)&=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)] \end{aligned}logP(Y∣θi+1)−logP(Y∣θi)​=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)]​

要证明P(Y∣θi+1)≥P(Y∣θi)P(Y|\theta^{i+1})\geq P(Y|\theta^i)P(Y∣θi+1)≥P(Y∣θi)，只需证明上式的右边是非负的。其中Q(θi+1,θi)−Q(θi,θi)Q(\theta^{i+1},\theta^i)-Q(\theta^{i},\theta^i)Q(θi+1,θi)−Q(θi,θi)是必大于0的，因为θi+1\theta^{i+1}θi+1是Q(θ,θi)Q(\theta,\theta^i)Q(θ,θi)的最大值。则：

H(θi+1,θi)−H(θi,θi)=∑ZP(Z∣Y,θi+1)logP(Z∣Y,θ)−∑ZP(Z∣Y,θi)logP(Z∣Y,θ)(利用Jensen不等式←)=∑ZP(Z∣Y,θi)logP(Z∣Y,θi+1)P(Z∣Y,θi)=log(∑ZP(Z∣Y,θi)P(Z∣Y,θi+1)P(Z∣Y,θi))=log(∑ZP(Z∣Y,θi+1))=log1=0\begin{aligned} H(\theta^{i+1},\theta^{i}) - H(\theta^{i},\theta^{i}) &=\sum\limits_ZP(Z|Y,\theta^{i+1})logP(Z|Y,\theta)-\sum\limits_ZP(Z|Y,\theta^i)logP(Z|Y,\theta)\\ (利用Jensen不等式\leftarrow)&=\sum_{Z}P(Z|Y,\theta^i)log\frac{P(Z|Y,\theta^{i+1})}{P(Z|Y,\theta^i)}\\ &=log\bigg(\sum_{Z}P(Z|Y,\theta^i)\frac{P(Z|Y,\theta^{i+1})}{P(Z|Y,\theta^i)}\bigg)\\ &=log\bigg(\sum_{Z}P(Z|Y,\theta^{i+1})\bigg)\\ &=log1\\ &=0 \end{aligned} H(θi+1,θi)−H(θi,θi)(利用Jensen不等式←)​=Z∑​P(Z∣Y,θi+1)logP(Z∣Y,θ)−Z∑​P(Z∣Y,θi)logP(Z∣Y,θ)=Z∑​P(Z∣Y,θi)logP(Z∣Y,θi)P(Z∣Y,θi+1)​=log(Z∑​P(Z∣Y,θi)P(Z∣Y,θi)P(Z∣Y,θi+1)​)=log(Z∑​P(Z∣Y,θi+1))=log1=0​

由此可以证明P(Y∣θi+1)≥P(Y∣θi)P(Y|\theta^{i+1})\geq P(Y|\theta^i)P(Y∣θi+1)≥P(Y∣θi)。

定理二：设L(θ)=logP(Y∣θ)L(\theta)=logP(Y|\theta)L(θ)=logP(Y∣θ)为观测数据的对数似然函数，θi(i=1,2,...)\theta^i(i=1,2,...)θi(i=1,2,...)是EM算法得到的参数估计序列，L(θi)(i=1,2,...)L(\theta^i)(i=1,2,...)L(θi)(i=1,2,...)是对应的对数似然函数序列。

1. 如果P(Y∣θ)P(Y|\theta)P(Y∣θ)有上界，则$L(θi)=logP(Y∣θi)L(\theta^i)=logP(Y|\theta^i)L(θi)=logP(Y∣θi)收敛到某一值L∗L^*L∗；
2. 在函数Q(θ,θ′)Q(\theta,\theta^{'})Q(θ,θ′)与L(θ)L(\theta)L(θ)满足一定条件下，由EM算法得到的参数估计序列θi\theta^iθi的收敛值θ∗\theta^*θ∗是L(θ)L(\theta)L(θ)的稳定点。

这里就不证明了

4、EM算法在高斯混合模型中的应用

4.1 高斯混合模型

  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∑K​ak​ϕ(y∣θk​)ϕ(y∣θk​)=2π​δk​1​exp(−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∏K​j=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∏K​j=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=1K​ak​ϕ(yi​∣θk​)ak​ϕ(yi​∣θk​)​j=1,2,..,N；k=1,2,...,Knk​=j=i∑N​Eγ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​(nk​log(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θmax​Q(θ,θ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​)2​ak​^​=N∑j=1N​γjk​^​​
  给一个初始值，来回迭代就可以求得值内容。这一块主要用到了Q(θ.θi)Q(\theta.\theta^i)Q(θ.θi)的导数，并且用到了E步的γjk^\hat{\gamma_{jk}}γjk​^​。

4.2 混合高斯分布模型python实现

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.4, 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)

E步主要计算内容：
其中我们定义γ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=1K​ak​ϕ(yi​∣θk​)ak​ϕ(yi​∣θk​)​j=1,2,..,N；k=1,2,...,Knk​=j=i∑N​Eγjk​

M步 主要计算内容
这一步骤主要是对Q函数求导后的数据进行计算，利用了 E 步 的γjk^\hat{\gamma_{jk}}γjk​^​

import math
import copy
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D#生成随机数据，4个高斯模型
def generate_data(sigma,N,mu1,mu2,mu3,mu4,alpha):global X                  #可观测数据集X = np.zeros((N, 2))       # 初始化X，2行N列。2维数据，N个样本X=np.matrix(X)global mu                 #随机初始化mu1，mu2，mu3，mu4mu = np.random.random((4,2))mu=np.matrix(mu)global excep              #期望第i个样本属于第j个模型的概率的期望excep=np.zeros((N,4))global alpha_             #初始化混合项系数alpha_=[0.25,0.25,0.25,0.25]for i in range(N):if np.random.random(1) < 0.1:  # 生成0-1之间随机数X[i,:]  = np.random.multivariate_normal(mu1, sigma, 1)     #用第一个高斯模型生成2维数据elif 0.1 <= np.random.random(1) < 0.3:X[i,:] = np.random.multivariate_normal(mu2, sigma, 1)      #用第二个高斯模型生成2维数据elif 0.3 <= np.random.random(1) < 0.6:X[i,:] = np.random.multivariate_normal(mu3, sigma, 1)      #用第三个高斯模型生成2维数据else:X[i,:] = np.random.multivariate_normal(mu4, sigma, 1)      #用第四个高斯模型生成2维数据print("可观测数据：\n",X)       #输出可观测样本print("初始化的mu1，mu2，mu3，mu4：",mu)      #输出初始化的mu# E 期望
#  \hat{\gamma_{jk}}
def e_step(sigma,k,N):global Xglobal muglobal excepglobal alpha_for i in range(N):denom=0for j in range(0,k):#  sigma.I 表示矩阵的逆矩阵# np.transpose ：矩阵转置   np.linalg.det():矩阵求行列式denom += alpha_[j]*  math.exp(-(X[i,:]-mu[j,:])*sigma.I*np.transpose(X[i,:]-mu[j,:]))  /np.sqrt(np.linalg.det(sigma))       #分母for j in range(0,k):numer = math.exp(-(X[i,:]-mu[j,:])*sigma.I*np.transpose(X[i,:]-mu[j,:]))/np.sqrt(np.linalg.det(sigma))        #分子excep[i,j]=alpha_[j]*numer/denom      #求期望print("隐藏变量：\n",excep)def m_step(k,N):global excepglobal Xglobal alpha_for j in range(0,k):denom=0   #分母numer=0   #分子for i in range(N):numer += excep[i,j]*X[i,:]denom += excep[i,j]mu[j,:] = numer/denom    #求均值alpha_[j]=denom/N        #求混合项系数#     #可视化结果
def plotShow():# 画生成的原始数据plt.subplot(221)plt.scatter(X[:,0].tolist(), X[:,1].tolist(),c='b',s=25,alpha=0.4,marker='o')    #T散点颜色，s散点大小，alpha透明度，marker散点形状plt.title('random generated data')#画分类好的数据plt.subplot(222)plt.title('classified data through EM')order=np.zeros(N)color=['b','r','k','y']for i in range(N):for j in range(k):if excep[i,j]==max(excep[i,:]):order[i]=j     #选出X[i,:]属于第几个高斯模型probility[i] += alpha_[int(order[i])]*math.exp(-(X[i,:]-mu[j,:])*sigma.I*np.transpose(X[i,:]-mu[j,:]))/(np.sqrt(np.linalg.det(sigma))*2*np.pi)    #计算混合高斯分布plt.scatter(X[i, 0], X[i, 1], c=color[int(order[i])], s=25, alpha=0.4, marker='o')      #绘制分类后的散点图#绘制三维图像ax = plt.subplot(223, projection='3d')plt.title('3d view')for i in range(N):ax.scatter(X[i, 0], X[i, 1], probility[i], c=color[int(order[i])])plt.show()

if __name__ == '__main__':iter_num=1000  #迭代次数N=500         #样本数目k=4            #高斯模型数probility = np.zeros(N)    #混合高斯分布u1=[5,35]u2=[30,40]u3=[20,20]u4=[45,15]sigma=np.matrix([[30, 0], [0, 30]])               #协方差矩阵alpha=[0.1,0.2,0.3,0.4]         #混合项系数generate_data(sigma,N,u1,u2,u3,u4,alpha)     #生成数据#迭代计算for i in range(iter_num):err=0     #均值误差err_alpha=0    #混合项系数误差Old_mu = copy.deepcopy(mu)Old_alpha = copy.deepcopy(alpha_)e_step(sigma,k,N)     # E步m_step(k,N)           # M步print("迭代次数:",i+1)print("估计的均值:",mu)print("估计的混合项系数:",alpha_)for z in range(k):err += (abs(Old_mu[z,0]-mu[z,0])+abs(Old_mu[z,1]-mu[z,1]))      #计算误差err_alpha += abs(Old_alpha[z]-alpha_[z])if (err<=0.001) and (err_alpha<0.001):     #达到精度退出迭代print(err,err_alpha)break

# 画图
plotShow()

5 EM算法在HMM模型中的应用

5.1 HMM模型

  HMM模型也是一个含有隐变量的模型，具体参考HMM模型详解

5.2 EM算法在HMM中应用

  HMM模型在学习参数的时候，可以使用极大似然估计，也可以使用EM算法（在HMM称为bamu-welch算法）

class HMM(object):def __init__(self, n, m, a=None, b=None, pi=None):# 可能的隐藏状态数self.N = n# 可能的观测数self.M = m# 状态转移概率矩阵self.A = a# 观测概率矩阵self.B = b# 初始状态概率矩阵self.Pi = pi# 观测序列self.X = None# 状态序列self.Y = None# 序列长度self.T = 0# 定义前向算法self.alpha = None# 定义后向算法self.beta = Nonedef baum_welch(self, x, criterion=0.001):self.T = len(x)self.X = xwhile True:# 为了得到alpha和beta的矩阵_ = self.forward(self.X)_ = self.backward(self.X)xi = np.zeros((self.T - 1, self.N, self.N), dtype=float)for t in range(self.T - 1):# 笨办法# for i in range(self.N):# gamma[t][i] = self.calculate_gamma(t, i)#         for j in range(self.N):#             xi[t][i][j] = self.calculate_psi(t, i, j)# for i in range(self.N):#     gamma[self.T-1][i] = self.calculate_gamma(self.T-1, i)# numpy矩阵的办法denominator = np.sum(np.dot(self.alpha[t, :], self.A) *self.B[:, self.X[t + 1]] * self.beta[t + 1, :])for i in range(self.N):molecular = self.alpha[t, i] * self.A[i, :] * self.B[:, self.X[t+1]]*self.beta[t+1, :]xi[t, i, :] = molecular / denominatorgamma = np.sum(xi, axis=2)prod = (self.alpha[self.T-1, :]*self.beta[self.T-1, :])gamma = np.vstack((gamma, prod / np.sum(prod)))new_pi = gamma[0, :]new_a = np.sum(xi, axis=0) / np.sum(gamma[:-1, :], axis=0).reshape(-1, 1)new_b = np.zeros(self.B.shape, dtype=float)for k in range(self.B.shape[1]):mask = self.X == knew_b[:, k] = np.sum(gamma[mask, :], axis=0) / np.sum(gamma, axis=0)if np.max(abs(self.Pi - new_pi)) < criterion and \np.max(abs(self.A - new_a)) < criterion and \np.max(abs(self.B - new_b)) < criterion:breakself.A, self.B, self.Pi = new_a, new_b, new_pi

参考

主要参考统计学习方法

统计学习方法-代码解读

EM算法 - 期望极大算法

机器学习最易懂之EM算法详解与python实现相关推荐

1. 【机器学习基础】EM算法详解及其收敛性证明

   EM算法详解 (一)单高斯模型 1.1 一维高斯分布: 1.2 多维高斯分布: (二)最大似然估计 2.1 最大似然估计的数学概念: 2.2 最大似然估计的基本步骤: 2.2.1 构造似然函数: 2. ...

2. AdaBoost算法详解与python实现

   AdaBoost算法详解与python实现 https://tangshusen.me/2018/11/18/adaboost/

3. RRT（Rapidly-Exploring Random Trees）算法详解及python实现

   RRT(Rapidly-Exploring Random Trees)算法详解及python实现 前言 一.原理 二.伪代码 三.代码详解 总结 前言 快速探索随机树(RRT):作为一种随机数据结构, ...

4. 排序算法(五)——堆排序算法详解及Python实现

   本文目录 一.简介 二.算法介绍 三.代码实现 排序算法系列--相关文章 一.简介 堆排序(Heap Sort)算法,属于选择排序类,不稳定排序,时间复杂度O(nlogn). 堆排序由Floyd和Wi ...

5. Simhash算法详解及python实现

   Simhash算法详解及python实现 GoogleMoses Charikar发表的一篇论文"detecting near-duplicates for web crawling&quo ...

6. kmeans算法详解和python代码实现

   kmeans算法详解和python代码实现 kmeans算法 无监督学习和监督学习 监督学习: 是通过已知类别的样本分类器的参数,来达到所要求性能的过程 简单来说,就是让计算机去学习我们已经创建好了的 ...

7. 编辑距离算法详解和python代码

   编辑距离(Levenshtein Distance)算法详解和python代码 最近做NLP用到了编辑距离,网上学习了很多,看到很多博客写的有问题,这里做一个编辑距离的算法介绍,步骤和多种python ...

8. python直线拟合_RANSAC算法详解(附Python拟合直线模型代码)

   之前只是简单了解RANSAC模型,知道它是干什么的.然后今天有个课程设计的报告,上去讲了一下RANSAC,感觉这个东西也没那么复杂,所以今天就总结一些RASAC并用Python实现一下直线拟合. RA ...

9. 【机器学习】集成学习及算法详解

   集成学习及算法详解 前言 一.随机森林算法原理 二.随机森林的优势与特征重要性指标 1.随机森林的优势 2.特征重要性指标 三.提升算法概述 四.堆叠模型简述 五.硬投票和软投票 1.概念介绍 2.硬 ...

最新文章

1. ST公司STM32F4与STM32F1的区别
2. 全球及中国车载扫地机行业销售前景态势与运营盈利分析报告2022版
3. 我的DWR学习（一）
4. 每日一题（53）—— 评价代码片段
5. 知方可补不足~数据库名称和数据库别名不同了怎么办
6. 两万字详细爬虫知识储备，数据采集与清洗基础习题（一）头歌参考答案
7. R变量索引 - 什么时候使用 @或$
8. 类字面常量和静态代码执行顺序
9. CentOS7安装Portainer实现docker可视化操作
10. PN结的形成及PN结工作原理（单向导电）讲解
11. HTML+5.2+新特性,HTML 5中的新特性
12. /etc/issue和/etc/motd
13. 数据分析真题日刷 | 网易2018校园招聘数据分析工程师笔试卷
14. nginx中配置gzip_static on提示nginx: [emerg] unknown directive “gzip_static“ in
15. 根据业绩确定提成比例并计算提成
16. 8通道温度采集器工作特性介绍
17. Bot console WeChat机器人
18. 酒店管理系统Python#qt
19. SP10707 COT2 - Count on a tree II【树上莫队】
20. 盛世昊通山东省运营中心盛大启动，迈向新的篇章

热门文章

1. STM32单片机RGB红蓝调光植物补光系统红光蓝光PWM调色调节亮度
2. Android网络加载框架Glide的使用
3. mybatis使用foreach进行批量操作 The error may involve defaultParameterMap
4. 项目经理如何面对困境
5. java图片加气泡文字,动态图片加气泡文字 微信动态图片加文字教程
6. 基于奇异值分解的图片压缩
7. 测试开发工程师mac电脑常用软件推荐
8. python基金比较上机题_手把手教你用python选基金
9. AI识别PS篡改图像
10. Android 集成百度身份证识别