EM算法

前言：
本来已经写过一遍了，但是很无奈公式的推导实在是看不懂，甚至觉得学习资料中的公式推导某些地方有些问题，或许是我跟不上思路。本来已经想放弃了，但是无意中看到一篇手写EM算法文章，文章地址，自觉惭愧，最后说到了我的心里。便决定从头再捋一遍。

• E步：求期望

• M步：求极大

EM算法的通俗理解

参考：https://www.jianshu.com/p/1121509ac1dc
如果使用基于最大似然估计的模型，模型中存在隐变量，就要用EM算法做参数估计。理解EM算法背后的idea，远比看懂它的数学推导重要。idea会让你有一个直观的感受，从而明白算法的合理性，数学推导只是将这种合理性用更加严谨的语言表达出来而已。
假设现在有两枚硬币1和2，,随机抛掷后正面朝上概率分别为P1，P2。为了估计这两个概率，做实验，每次取一枚硬币，连掷5下，记录下结果，如下：

硬币	结果	统计
1	正正反正反	3正-2反
2	反反正正反	2正-3反
1	正反反反反	1正-四反
2	正反反正正	3正-2反
1	反正正反反	2正-3反

可以很容易地估计出P1和P2，如下：

P1 = （3+1+2）/ 15 = 0.4
P2= （2+3）/10 = 0.5

到这里，一切似乎很美好，下面我们加大难度。

加入隐变量z

还是上面的问题，现在我们抹去每轮投掷时使用的硬币标记，如下：

硬币	结果	统计
Unknown	正正反正反	3正-2反
Unknown	反反正正反	2正-3反
Unknown	正反反反反	1正-四反
Unknown	正反反正正	3正-2反
Unknown	反正正反反	2正-3反

好了，现在我们的目标没变，还是估计P1和P2，要怎么做呢？

显然，此时我们多了一个隐变量z，可以把它认为是一个5维的向量（z1,z2,z3,z4,z5)，代表每次投掷时所使用的硬币，比如z1，就代表第一轮投掷时使用的硬币是1还是2。但是，这个变量z不知道，就无法去估计P1和P2，所以，我们必须先估计出z，然后才能进一步估计P1和P2。

但要估计z，我们又得知道P1和P2，这样我们才能用最大似然概率法则去估计z，这不是鸡生蛋和蛋生鸡的问题吗，如何破？

答案就是先随机初始化一个P1和P2，用它来估计z，然后基于z，还是按照最大似然概率法则去估计新的P1和P2，如果新的P1和P2和我们初始化的P1和P2一样，请问这说明了什么？（此处思考1分钟）

这说明我们初始化的P1和P2是一个相当靠谱的估计！

就是说，我们初始化的P1和P2，按照最大似然概率就可以估计出z，然后基于z，按照最大似然概率可以反过来估计出P1和P2，当与我们初始化的P1和P2一样时，说明是P1和P2很有可能就是真实的值。这里面包含了两个交互的最大似然估计。

如果新估计出来的P1和P2和我们初始化的值差别很大，怎么办呢？就是继续用新的P1和P2迭代，直至收敛。

这就是下面的EM初级版。

我们不妨这样，先随便给P1和P2赋一个值，比如：
P1 = 0.2
P2 = 0.7

然后，我们看看第一轮抛掷最可能是哪个硬币。
如果是硬币1，得出3正2反的概率为 0.20.20.20.80.8 = 0.00512
如果是硬币2，得出3正2反的概率为0.70.70.70.30.3=0.03087
然后依次求出其他4轮中的相应概率。做成表格如下：

轮数	若是硬币1	若是硬币2
1	0.00512	0.03087
2	0.02048	0.01323
3	0.08192	0.00567
4	0.00512	0.03087
5	0.02048	0.01323

按照最大似然法则：
第1轮中最有可能的是硬币2
第2轮中最有可能的是硬币1
第3轮中最有可能的是硬币1
第4轮中最有可能的是硬币2
第5轮中最有可能的是硬币1

我们就把上面的值作为z的估计值。然后按照最大似然概率法则来估计新的P1和P2。

P1 = （2+1+2）/15 = 0.33
P2=（3+3）/10 = 0.6

设想我们是全知的神，知道每轮抛掷时的硬币就是如本文第001部分标示的那样，那么，P1和P2的最大似然估计就是0.4和0.5（下文中将这两个值称为P1和P2的真实值）。那么对比下我们初始化的P1和P2和新估计出的P1和P2：

初始化的P1	估计出的P1	真实的P1	初始化的P2	估计出的P2	真实的P2
0.2	0.33	0.4	0.7	0.6	0.5

看到没？我们估计的P1和P2相比于它们的初始值，更接近它们的真实值了！

可以期待，我们继续按照上面的思路，用估计出的P1和P2再来估计z，再用z来估计新的P1和P2，反复迭代下去，就可以最终得到P1 = 0.4，P2=0.5，此时无论怎样迭代，P1和P2的值都会保持0.4和0.5不变，于是乎，我们就找到了P1和P2的最大似然估计。

这里有两个问题：
1、新估计出的P1和P2一定会更接近真实的P1和P2？
答案是：没错，一定会更接近真实的P1和P2，数学可以证明，但这超出了本文的主题，请参阅其他书籍或文章。
2、迭代一定会收敛到真实的P1和P2吗？
答案是：不一定，取决于P1和P2的初始化值，上面我们之所以能收敛到P1和P2，是因为我们幸运地找到了好的初始化值。

EM进阶版

下面，我们思考下，上面的方法还有没有改进的余地？

我们是用最大似然概率法则估计出的z值，然后再用z值按照最大似然概率法则估计新的P1和P2。也就是说，我们使用了一个可能的z值，而不是所有可能的z值。

如果考虑所有可能的z值，对每一个z值都估计出一个新的P1和P2，将每一个z值概率大小作为权重，将所有新的P1和P2分别加权相加，这样的P1和P2应该会更好一些。

所有的z值有多少个呢？显然，有2^5=32种，需要我们进行32次估值？？

不需要，我们可以用期望来简化运算。

轮数	若是硬币1	若是硬币2
1	0.00512	0.03087
2	0.02048	0.01323
3	0.08192	0.00567
4	0.00512	0.03087
5	0.02048	0.01323

利用上面这个表，我们可以算出每轮抛掷中使用硬币1或者使用硬币2的概率。比如第1轮，使用硬币1的概率是：
0.00512/(0.00512+0.03087)=0.14
使用硬币2的概率是1-0.14=0.86
依次可以算出其他4轮的概率，如下：

轮数	z_i=硬币1	z_i=硬币2
1	0.14	0.86
2	0.61	0.39
3	0.94	0.06
4	0.14	0.86
5	0.61	0.39

上表中的右两列表示期望值。看第一行，0.86表示，从期望的角度看，这轮抛掷使用硬币2的概率是0.86。相比于前面的方法，我们按照最大似然概率，直接将第1轮估计为用的硬币2，此时的我们更加谨慎，我们只说，有0.14的概率是硬币1，有0.86的概率是硬币2，不再是非此即彼。这样我们在估计P1或者P2时，就可以用上全部的数据，而不是部分的数据，显然这样会更好一些。

这一步，我们实际上是估计出了z的概率分布，这步被称作E步。
结合下表：

硬币	结果	统计
Unknown	正正反正反	3正-2反
Unknown	反反正正反	2正-3反
Unknown	正反反反反	1正-四反
Unknown	正反反正正	3正-2反
Unknown	反正正反反	2正-3反

我们按照期望最大似然概率的法则来估计新的P1和P2：

以P1估计为例，第1轮的3正2反相当于
$0.14\ast 3=0.42$正
$0.14\ast 2=0.28$反
类似于将每一轮中是硬币1的概率作为权重进行计算。
依次算出其他四轮，列表如下：

轮数	正面	反面
1	0.42	0.28
2	1.22	1.83
3	0.94	3.76
4	0.42	0.28
5	1.22	1.83
总计	4.22	7.98

P1=4.22/(4.22+7.98)=0.35
可以看到，改变了z值的估计方法后，新估计出的P1要更加接近0.4。原因就是我们使用了所有抛掷的数据，而不是之前只使用了部分的数据。

这步中，我们根据E步中求出的z的概率分布，依据最大似然概率法则去估计P1和P2，被称作M步。

EM算法的原理（公式推导）

概率模型有时候既含有观测变量，又含有隐变量或潜在变量，如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计方法估计模型参数。
EM算法是针对模型中含有隐变量的情况，计算参数的极大似然估计法，或·者是极大后验概率估计法。我们讨论极大似然估计法，极大后验概率估计法与其类似。
假设有3枚硬币，分别记作A，B，C，正面出现的概率分别是$\pi ,p,q$。实验规则如下：

• 过程1：先抛掷硬币A，若是正面，选择硬币B，若是反面，选择硬币C。
• 过程2：根据选出的硬币进行抛掷，掷硬币的结果记为变量$y$，即为最终结果。
  则$y$为0-1变量，$y=1$代表最终抛掷硬币的结果为正面，$y=0$代表为反面。
  如此独立地重复n次实验，我们当前规定n=10，则10次的结果如下所示：
  $$1,1,0,1,0,0,1,0,1,1$$
  构建三硬币模型，来计算最终的结果。对于该模型，设过程1的结果为$z$，则$z$为不可观测变量。
  $$\begin{array}{rl}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{array}$$
  $\theta =(\pi ,p,q)$表示模型参数.
  将观测数据表示为$Y=(Y_1,Y_2,...,Y_n{)}^T$，未观测的数据表示为$Z=(Z_1,Z_2,...,Z_n{)}^T$。
  下面将试验次数的标记标在下面，将参数$\theta$的迭代次数标在上面，如${\theta }^i$。
  观测函数的似然函数是：
  $$\begin{array}{rl}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{array}$$
  考虑求模型参数$\theta =(\pi ,p,q)$的极大似然估计，即：
  $$\hat{\theta }=arg\underset{\theta }{\max }logP(Y|\theta )$$
  这个问题没有解析解，只有通过迭代方法来求解，EM算法就是可以用于求解这个问题的一种迭代算法，下面给出EM算法的迭代过程：
• 首先选取初始值，记做${\theta }^0=({\pi }^0,p^0,q^0)$，第i次的迭代参数的估计值为${\theta }^i=({\pi }^i,p^i,q^i)$
• E步：计算在模型参数${\pi }^i，p^i，q^i$下观测变量$y_i$来源于硬币B的概率：
  $${\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}}$$
  备注一下：这个公式的分母是$P(Y|\theta )$，分子表示是来源于B硬币的概率。
• M步：计算模型参数的新估计值：
  $${\pi }^{i+1}=\frac{1}{n}\sum _{j=1}^n{\mu }_j^{i+1}$$
  因为B硬币是A硬币出现正面的结果，$y$是0-1变量，所以A硬币出现正面的概率就是${\mu }_j$的平均值。
  $$p^{i+1}=\frac{{\sum }_{j=1}^n{\mu }_j^{i+1}{y}_j}{{\sum }_{j=1}^n{\mu }_j^{i+1}}$$
  分子乘以$y_i$，所以其实是计算B硬币出现正面的概率。
  $$q^{i+1}=\frac{{\sum }_{j=1}^n(1-{\mu }_j^{i+1})y_j}{{\sum }_{j=1}^n(1-{\mu }_j^{i+1})}$$
  $(1-{\mu }_j^{i+1})$表示出现C硬币的概率。
  闭环形成，从$P(Y|\theta )$ 到 $\pi 、p、q$一个闭环流程，接下来可以通过迭代法来做完成。针对上述例子，我们假设初始值为${\pi }^0=0.5，p^0=0.5，q^0=0.5$，因为对$y_i=1$和$y_i=0$均有${\mu }_j^1=0.5$，利用迭代公式计算得到${\pi }^1=0.5，p^1=0.6，q^1=0.6$，继续迭代得到最终的参数：$$\widehat{{\pi }^0}=0.5，\widehat{p^0}=0.6，\widehat{q^0}=0.6$$如果一开始初始值选择为：${\pi }^0=0.4，p^0=0.6，q^0=0.7$，那么得到的模型参数的极大似然估计是$$\hat{\pi }=0.4064，\hat{p}=0.5368，\hat{q}=0.6432$$，这说明EM算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。
  这个例子中你只观察到了硬币抛完的结果，并不了解A硬币抛完之后，是选择了B硬币抛还是C硬币抛，这时候概率模型就存在着隐含变量!

EM算法的具体实现

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

• (1)选择参数的初值${\theta }^0$，开始迭代
• (2) E步：记${\theta }^i$为第i次迭代参数$\theta$的估计值，在第i+1次迭代的E步，计算$$\begin{array}{rl}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{array}$$以z为主变量求期望，对应通俗解释里的利用期望的最大似然概率估计，为M步做铺垫。
  通俗解释：在已知结果Y以及参数$\theta$下，Y,Z在参数$\theta$下概率的对数的期望，${\sum }_{z}P(Z|Y,{\theta }^i)=1$。
  这里，$P(Z|Y,{\theta }^i)$（对应通俗解释中的0.14）是在给定观测数据Y和当前的参数估计${\theta }^i$下隐变量数据Z的条件概率分布；
• (3) M步：求使$Q(\theta ,{\theta }^i)$极大化的$\theta$，确定第i+1次迭代的参数的估计值${\theta }^{i+1}$，$${\theta }^{i+1}=arg\underset{\theta }{\max }Q(\theta ,{\theta }^i)$$$Q(\theta ,{\theta }^i)$是EM算法的核心，称为Q函数(Q function)，这个是需要自己构造的。
• (4) 重复第(2)步和第(3)步，直到收敛，收敛条件：$$||{\theta }^{i+1}-{\theta }^i||<{\epsilon }_1$$或者：$$||Q({\theta }^{i+1},{\theta }^i)-Q({\theta }^i,{\theta }^i)||<{\epsilon }_2$$收敛迭代就结束了。

证明部分

Jensen不等式：

• $f(x)$是凸函数
  $$f(E(X))\le E(f(x))$$
• $f(x)$是凹函数
  $$f(E(X))\ge E(f(x))$$
  推到EM算法近似实现对观测数据的极大似然估计的办法是找到似然函数的下界，让下界最大，通过逼近的方式实现对观测数据的最大似然估计。具体证明过程如下：
  对数似然函数：
  $$L(\theta )=L(Y|\theta )=logP(Y|\theta )=log\sum _{z}P(Y|Z,\theta )P(Z|\theta )$$
  $$L(\theta )-L({\theta }^i)=log(\sum _{z}P(Y|Z,\theta )P(Z|\theta ))-L({\theta }^i)$$
  $$\begin{array}{r}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)})-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)\end{array}$$
  $\ge$这一个步骤就是采用了凹函数的Jensen不等式做转换。
  
  

因为$Z$是隐藏变量，所以有${\sum }_{Z}P(Z|Y,{\theta }^i)==1，P(Z|Y,{\theta }^i)>0$，于是继续变：
$$\begin{array}{rl}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)})-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)log(\frac{P(Y|Z,{\theta }^i)P(Z|{\theta }^i)}{P(Z|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)})\\ & \ge 0\end{array}$$
$$L({\theta }^i)=L(Y|{\theta }^i)=logP(Y|{\theta }^i)=log\sum _{z}P(Y|Z,{\theta }^i)P(Z|{\theta }^i)$$
令$B(\theta ,{\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)})$
则$L(\theta )\ge B(\theta ,{\theta }^i)$,所以$B(\theta ,{\theta }^i)$为$L(\theta )$的一个下界，因此为了使$L(\theta )$尽可能的大，所应该选择使$B(\theta ,{\theta }^i)$最大的那个${\theta }^{i+1}$。
$$\begin{array}{rl}{\theta }^{i+1}& =arg\underset{\theta }{\max }[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)})]\\ & =arg\underset{\theta }{\max }\sum _{Z}P(Z|Y,{\theta }^i)log(P(Y|Z,\theta )P(Z,\theta ))\\ & =arg\underset{\theta }{\max }\sum _{Z}P(Z|Y,{\theta }^i)log(P(Y,Z|\theta ))\\ & =arg\underset{\theta }{\max }Q(\theta ,{\theta }^i)\end{array}$$
上式中的等号是等价的意思，并不是绝对的等式转换关系。
下面证明收敛
我们知道已知观测数据的似然函数是$P(Y|\theta )$，对数似然函数为：
$$\begin{array}{rl}L(Y|\theta )=\sum _{i=1}^{n}logP(y_i,\theta )& =\sum _{i=1}^{n}log(\frac{P(y_i,Z|\theta )}{P(Z|y_i,\theta )})\\ & =\sum _{i=1}^{n}logP(y_i,Z|\theta )-\sum _{i=1}^{n}logP(Z|y_i,\theta )\end{array}$$
要证明收敛，就证明单调递增，$$\sum _{i=1}^{n}logP(y_i,Z|{\theta }^{j+1})>\sum _{i=1}^{n}logP(Z|y_i,{\theta }^j)$$即$L(Y|{\theta }^{j+1})\ge L(Y|{\theta }^j)$
由上文知道：
$$\begin{array}{rl}Q(\theta ,{\theta }^i)& =\sum _{Z}logP(Y,Z|\theta )P(Z|Y,{\theta }^i)\\ & =\sum _{i=1}^n\sum _{Z}logP(y_i,Z|\theta )P(Z|y_i,{\theta }^i)\end{array}$$
我们构造一个函数$G({\theta }^i,\theta )$，
$$\begin{array}{rl}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)})-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)\end{array}$$$$G({\theta }^i,\theta )=\sum _{Z}P(Z|Y,{\theta }^i)log(\frac{P(Y|Z,\theta )P(Z|\theta )}{P(Z|Y,{\theta }^i)})$$$$\begin{array}{rl}G({\theta }^i,{\theta }^i)-L({\theta }^i)& =\sum _{Z}P(Z|Y,{\theta }^i)log(\frac{P(Y,Z|{\theta }^i)}{P(Z|Y,{\theta }^i)})-logP(Y|{\theta }^i)\\ & =\sum _{Z}P(Z|Y,{\theta }^i)log(\frac{P(Y|Z,{\theta }^i)P(Z|{\theta }^i)}{P(Z|Y,{\theta }^i)})-\sum _{z}P(Z|Y,{\theta }^i)logP(Y|{\theta }^i)\\ & =\sum _{Z}P(Z|Y,{\theta }^i)log(\frac{P(Y|Z,{\theta }^i)P(Z|{\theta }^i)}{P(Z|Y,{\theta }^i)P(Y|{\theta }^i)})\\ & =\sum _{z}P(Z|Y,{\theta }^i)log1\\ & =0\end{array}$$
所以似然函数$L(\theta )$与$G({\theta }^i,\theta )$在${\theta }^i$这一点处相等，而对于$G({\theta }^i,\theta )$,$$G({\theta }^i,{\theta }^i)\le G({\theta }^i,{\theta }^{i+1})$$又因为$G({\theta }^i,\theta )$是$L(\theta )$的下限函数，所以：$$G({\theta }^i,{\theta }^{i+1})\le L({\theta }^{i+1})$$
综上：$L({\theta }^{i+1})\ge L({\theta }^i)$

高斯混合分布

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

如果有多个高斯模型，公式表示为：
$$\begin{gathered} P(y|\theta )=\sum _{k=1}^K{a}_k\varphi (y|{\theta }_k) \\ \varphi (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 \end{gathered}$$
$\varphi (y|{\theta }_k)$表示为第k个高斯分布密度模型，定义如上，其中$a_k$表示被选中的概率。在本次模型$P(y|\theta )$中，观测数据是已知的，而观测数据具体来自哪个模型是未知的，有点像之前提过的三硬币模型，我们来对比一下，A硬币就像是概率$a_k$，用来表明具体的模型，而B、C硬币就是具体的模型，只不过这里有很多个模型，不仅仅是B、C这两个模型。我们用${\gamma }_{jk}$来表示，则：
$${\gamma }_{jk}=\{\begin{array}{ll}1& \text{第j个观测数据来源于第k个模型}\\ 0& \text{否则}\end{array}$$
所以一个观测数据$y_j$的隐藏数据$({\gamma }_{j1},{\gamma }_{j2},...,{\gamma }_{jk})$，那么完全似然函数就是：

$$P(y,\gamma |\theta )=\prod _{k=1}^K\prod _{j=1}^N[a_k\varphi (y|{\theta }_k){]}^{{\gamma }_{jk}}$$

取对数之后等于：

$$\begin{array}{rl}log(P(y,\gamma |\theta ))& =log(\prod _{k=1}^K\prod _{j=1}^N[a_k\varphi (y|{\theta }_k){]}^{{\gamma }_{jk}})\\ & =\sum _K^{k=1}(\sum _{j=1}^N({\gamma }_{jk})log(a_k)+\sum _{j=1}^N({\gamma }_{jk})[log(\frac{1}{\sqrt{2\pi }})-log({\delta }_k)-\frac{(y_i-{\mu }_k{)}^2}{2{\delta }_k^2}\left]\right)\end{array}$$

• E 步 ：
  $$\begin{array}{rl}Q(\theta .{\theta }^i)& =E[log(P(y,\gamma |\theta ))]\\ & =\sum _K^{k=1}(\sum _{j=1}^N(E{\gamma }_{jk})log(a_k)+\sum _{j=1}^N(E{\gamma }_{jk})[log(\frac{1}{\sqrt{2\pi }})-log({\delta }_k)-\frac{(y_i-{\mu }_k{)}^2}{2{\delta }_k^2}\left]\right)\end{array}$$
  其中我们定义$\widehat{{\gamma }_{jk}}$：
  $$\begin{gathered} \widehat{{\gamma }_{jk}}=E({\gamma }_{jk}|y,\theta )=\frac{a_k\varphi (y_i|{\theta }_k)}{{\sum }_{k=1}^K{a}_k\varphi (y_i|{\theta }_k)} \\ j=1,2,..,N；k=1,2,...,K \\ {n}_k=\sum _{j=i}^{N}E{\gamma }_{jk} \end{gathered}$$
  于是化简得到：
  $$\begin{array}{rl}Q(\theta .{\theta }^i)& =\sum _K^{k=1}(n_{k}log(a_k)+\sum _{j=1}^N(E{\gamma }_{jk})[log(\frac{1}{\sqrt{2\pi }})-log({\delta }_k)-\frac{(y_i-{\mu }_k{)}^2}{2{\delta }_k^2}\left]\right)\end{array}$$

E 步 在代码设计上只有$\widehat{{\gamma }_{jk}}$有用，用于M步的计算。

• M步，
  $${\theta }^{i+1}=arg\underset{\theta }{\max }Q(\theta ,{\theta }^i)$$
  对$Q(\theta ,{\theta }^i)$求导，得到每个未知量的偏导，使其偏导等于0，求解得到：
  $$\begin{gathered} \widehat{{\mu }_k}=\frac{{\sum }_{j=1}^N\widehat{{\gamma }_{jk}}{y}_i}{{\sum }_{j=1}^N\widehat{{\gamma }_{jk}}} \\ \widehat{{\delta }_k}=\frac{{\sum }_{j=1}^N\widehat{{\gamma }_{jk}}(y_i-{\mu }_k{)}^2}{{\sum }_{j=1}^N\widehat{{\gamma }_{jk}}} \\ \widehat{a_k}=\frac{{\sum }_{j=1}^N\widehat{{\gamma }_{jk}}}{N} \end{gathered}$$
  给一个初始值，来回迭代就可以求得值内容。这一块主要用到了$Q(\theta .{\theta }^i)$的导数，并且用到了E步的$\widehat{{\gamma }_{jk}}$。