机器学习笔记补充——EM算法及其在GMM中的应用

  本文内容主体是基于李航老师的《统计学习方法》第二版中EM算法章节的内容，有兴趣的同学可以自行参考书籍资料。其中个人觉得不是很直观的部分进行了适当的解释，解释仅属于个人理解，若觉得解释更加难以理解的同学可以自行忽略。代码部分有时间再补充。若文中有编辑错误，烦请指正，谢谢。

1. EM算法的简介及案例介绍

  EM算法是一种迭代算法，用于含有隐变量的概率模型参数的极大似然估计或极大后验概率估计。注意，EM算法不是一种机器学习的模型，只是一种求参数的算法。EM算法的每次迭代都由两步组成：E步，求期望；M步，求极大值。所以该算法才被称为期望极大算法(Expectation Maximization algorithm)。
  从上述描述可知，EM算法求解极大似然估计和其他模型很不同的一点是，它可以求解含有隐变量的概率模型参数。通常，概率模型的变量都是可观测变量，这时就可以根据给定的数据，直接使用极大似然估计法或贝叶斯估计法求解模型参数。但是，如果概率模型涉及到的变量其中有不可观测的变量怎么办？这个不可观测的变量就叫做隐变量。所谓不可观测，直观理解就是对于目标问题并没有给出该变量对应的任何信息，一切都是隐于深处的。
  还是那句话，概念总是深奥且不易理解的，所以还是举个例子来加强对概念的理解。
  案例：
  假设有两个形状不规则的硬币A和B，我们希望求出每个硬币抛出正面的概率${\hat{\theta }}_A,{\hat{\theta }}_B$。目前，我们已知五组抛硬币的实验结果，每组抛10次，结果如下表。

实验组号	结果1(H表示正，T表示反)	结果2("A:"表示用A硬币做该组实验，H表示正，T表示反)
1	H T T T H H T H T H	B:H T T T H H T H T H
2	H H H H H H T H H H	A:H H H H H H T H H H
3	H T H H H H T H H H	A:H T H H H H T H H H
4	H T T T H T T H T H	B:H T T T H T T H T H
5	H H T H H H T H T H	A:H H T H H H T H T H

  如果问题的背景如表中结果2的类型，也就是每组实验都知道是用哪个硬币抛的，那么根据极大似然估计的方法，很容易就能求出：
$${\hat{\theta }}_A=\frac{24}{24+6}=0.8，{\hat{\theta }}_B=\frac{9}{9+11}=0.45$$
  但是问题如果稍微做一点改变，实验结果是如结果1所示，也就是这时候，你只有5组正反的结果，却不知道每组是用哪个硬币做的实验，这时候问题就来了，极大似然估计是没法用的，因为这个问题存在一个隐变量，也就是你在分析结果的时候需要思考，这组结果是用A硬币做的还是B硬币做的？而这个事件，就是一个隐随机变量，这时候，EM算法就能派上用场了。
  其实，现实中，包含隐变量的问题是很常见的，所以EM算法的应用是很广泛的。那么针对上述问题，EM算法是如何求解的呢？
  首先，我们按随便取初值${\hat{\theta }}_A=0.6,{\hat{\theta }}_B=0.5$，然后拿出一组结果数据(以第一组为例)，由于不知道该结果是用哪个硬币做的，所以我们需要计算每个硬币产生这种结果的概率：
$$P(resul{t}_{1A})=C_{10}^5{\hat{\theta }}_A^5(1-{\hat{\theta }}_A{)}^5=0.201,P(resul{t}_{1B})=C_{10}^5{\hat{\theta }}_B^5(1-{\hat{\theta }}_B{)}^5=0.246$$
  那么，结果1是由A硬币或B硬币产生的概率就是：
$$P(resul{t}_1|A)=\frac{0.201}{0.201+0.246}\approx 0.45，P(resul{t}_1|B)=\frac{0.246}{0.201+0.246}\approx 0.55$$
  同理，可以计算出其他结果的类似概率，如下表。

组号	1	2	3	4	5
A硬币	0.45	0.80	0.73	0.35	0.65
B硬币	0.55	0.20	0.27	0.65	0.35

  当得到每组结果是由A硬币或B硬币产生的概率后，就可以根据极大似然估计的方法来进行估计，首先我们可以求出A硬币和B硬币分别对每组结果中，出现正反次数的期望(E步)，如下表：

组号	1	2	3	4	5	总计
A硬币	0.45×5H=2.2H,0.45×5T=2.2T	7.2H,0.8T	5.9H,1.5T	1.4H,2.1T	4.5H,1.9T	21.3H,8.6T
B硬币	0.55×5H=2.8H,0.55×5T=2.8T	1.8H,0.2T	2.1H,0.5T	2.6H,3.9T	2.5H,1.1T	11.7H,8.4T

  然后按极大似然估计的方法就可以求出新一轮的参数：
$${\hat{\theta }}_A^{(1)}=\frac{21.3}{21.3+8.6}=0.71,{\hat{\theta }}_B^{(1)}=\frac{11.7}{11.7+8.4}=0.58$$
  最后，将新得到的${\hat{\theta }}_A^{(1)},{\hat{\theta }}_B^{(1)}$继续带入最开始的步骤，重新计算，又可以得到${\hat{\theta }}_A^{(2)},{\hat{\theta }}_B^{(2)}$，反复迭代，最终达到终止条件或者结果收敛为止。最后，可以得到${\hat{\theta }}_A^{(n)}=0.8,{\hat{\theta }}_B^{(n)}=0.52$。这就是EM算法得到的参数估计结果。

2. EM算法的推导

  一般地，用Y表示可观测变量，Z表示隐变量数据，两者的联合分布是$P(Y,Z|\theta )$，条件分布是$P(Z|Y,\theta )$。Y和Z连在一起称为完全数据，单独的Y又称为不完全数据。
  面对一个含有隐变量的概率模型，目标是极大化观测数据Y关于参数$\theta$的对数似然函数，即求下式的极大值：
$$L(\theta )=\log P(Y|\theta )=\log \sum _{Z}P(Y,Z|\theta )=\log (\sum _{Z}P(Y|Z,\theta )P(Z|\theta ))$$
  上式中很难求解的部分一是带有不可观测数据Z的变量；二是具有包含求和公式的对数。稍微解释一下上式中的各项：
  $P(Y|\theta )$表示已知参数$\theta$的情况下，可观测数据Y取某一值的概率，由于问题包含隐变量，且隐变量不可观测，也就导致了$P(Y|\theta )$其实是无法直接求解的，而是要将其进行分解；
  ${\sum }_{Z}P(Y,Z|\theta )$表示，在已知参数$\theta$的情况下，遍历所有隐变量Z的取值情况的变量Y取某一值的概率。很明显，Z是所有隐变量的集合，那么遍历所有Z的情况后，其概率之和自然就与$P(Y|\theta )$相同了，这一步也就是把$P(Y|\theta )$分解成了一个个子项，每个子项都与过程中的隐变量有关。
  $P(Y,Z|\theta )=P(Y|Z,\theta )P(Z|\theta )$这其实就是条件概率公式，数学公式很好理解，但是其实$P(Y|Z,\theta )、P(Z|\theta )$两项都无法求解，因为隐变量是无法确切得到的。

  事实上，EM算法的本质是通过迭代逐步接近$L(\theta )$的极大值的。假设在第i次迭代后$\theta$的估计值是${\theta }^{(i)}$。EM算法希望保证每一次新的估计值$\theta$都能使$L(\theta )$增加，即$L(\theta )>L({\theta }^{(i)})$，这样就能逐步达到极大值。因此，考虑两者之差：
$$L(\theta )-L({\theta }^{(i)})=\log (\sum _{Z}P(Y|Z,\theta )P(Z|\theta ))-\log P(Y|{\theta }^{(i)})$$
  利用Jensen不等式($\log {\sum }_j{\lambda }_j{y}_j\ge {\sum }_j{\lambda }_j\log y_j$,其中${\lambda }_j\ge 0,{\sum }_j{\lambda }_j=1$)得到其下界：
$$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)})$$
$$\because \log P(Y|{\theta }^{(i)})=1\times \log P(Y|{\theta }^{(i)})=\sum _{Z}P(Z|Y,{\theta }^{(i)})\log \times P(Y|{\theta }^{(i)})=\sum _Z(P(Z|Y,{\theta }^{(i)})\log P(Y|{\theta }^{(i)}))$$
$$\therefore 原式\ge \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)})}$$
$$令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 }^{(i)})=B({\theta }^{(i)},{\theta }^{(i)})$。所以，理论上，任何可以让$B(\theta ,{\theta }^{(i)})$增大的$\theta$，也可以让$L(\theta )$增大。为了让$L(\theta )$尽可能的增大，所以选择新一轮迭代的参数值为：
$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }B(\theta ,{\theta }^{(i)})$$
  因此，为了求${\theta }^{(i+1)}$,忽略$B(\theta ,{\theta }^{(i)})$中所有与$\theta$无关的常数项(留下$P(Z|Y,{\theta }^{(i)})$)：
$${\theta }^{(i+1)}=\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 ))$$
  记$Q(\theta ,{\theta }^{(i)})={\sum }_{Z}P(Z|Y,{\theta }^{(i)})\log P(Y,Z|\theta )$，称为Q函数，是可以求解为关于$\theta$的表达式，并根据表达式求解极大值的。
  下图是$L(\theta )$和$B(\theta ,{\theta }^{(i)})$的曲线(来自李航《统计学习方法》)，根据上述推导，在$\theta ={\theta }^{(i)}$处，两者是相同的，$B(\theta ,{\theta }^{(i)})$取极大值后，$L(\theta )$也会上升。但是从图中也可以看出，EM算法并不能保证找到全局最优解。

3. EM算法

3.1 算法步骤：

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

1. 选择参数的初值${\theta }^{(0)}$，开始EM迭代；
2. E步：记${\theta }^{(i)}$为第i次迭代参数$\theta$的估计值，在第i+1次迭代的E步，计算Q函数：
   $$Q(\theta ,{\theta }^{(i)})=E_Z[\log P(Y,Z|\theta )|Y,{\theta }^{(i)}]=\sum _Z\log P(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})$$
3. M步：求使$Q(\theta ,{\theta }^{(i)})$达到极大值的$\theta$，确定第i+1次迭代的参数估计值${\theta }^{(i+1)}$
   $${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$$
4. 重复第2,3步，直到收敛，或达到设置迭代停止条件(通常是$||{\theta }^{(i+1)}-{\theta }^{(i)}||<{ϵ}_1,{ϵ}_1$是极小的正数)

  EM算法也可以用在无监督学习中，当训练数据只有输入，没有对应标签时，需要学习模型的任务称为无监督学习。EM算法可以把无监督学习的训练数据看作是联合概率分布P(X,Y)产生的数据，其中X是可观测数据，Y是不可观测数据。

3.2 EM算法的收敛性

  定理1：设$P(Y|\theta )$为观测数据的似然函数，${\theta }^{(i)}(i=1,2,...)$为EM算法得到的参数估计序列，$P(Y|{\theta }^{(i)})(i=1,2,...)$为对应的似然函数序列，则$P(Y|{\theta }^{(i)})$是单调递增的，即：$P(Y|{\theta }^{(i+1)})\ge P(Y|{\theta }^{(i)})$。
  证明：
$$P(Y|\theta )=\frac{P(Z|Y,\theta )P(Y|\theta )}{P(Z|Y,\theta )}=\frac{P(Y,Z|\theta )}{P(Z|Y,\theta )}$$
$$取对数得：\log P(Y|\theta )=\log P(Y,Z|\theta )-\log P(Z|Y,\theta )$$
$$\{\begin{array}{l}已知：Q(\theta ,{\theta }^{(i)})={\sum }_Z\log P(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})\\ 设定：H(\theta ,{\theta }^{(i)})={\sum }_Z\log P(Z|Y,\theta )P(Z|Y,{\theta }^{(i)})\end{array}$$
$$所以对数似然函数可以写成：\log P(Y|\theta )=Q(\theta ,{\theta }^{(i)})-H(\theta ,{\theta }^{(i)})$$
  令$\theta ={\theta }^{(i)}$和${\theta }^{(i+1)}$代入上式可得：
$$\log P(Y|{\theta }^{(i+1)})-\log P(Y|{\theta }^{(i)})=[Q({\theta }^{(i+1)},{\theta }^{(i)})-Q({\theta }^{(i+1)},{\theta }^{(i)})]-[H({\theta }^{(i+1)},{\theta }^{(i)})-H({\theta }^{(i+1)},{\theta }^{(i)})]$$
  要证明$P(Y|{\theta }^{(i+1)})\ge P(Y|{\theta }^{(i)})$，则根据log函数单调性可知，上式右端是非负的即可。根据Q函数的定义，新选的${\theta }^{(i+1)}$需要使Q函数达到极大值，因此显然有：
$$Q({\theta }^{(i+1)},{\theta }^{(i)})-Q({\theta }^{(i+1)},{\theta }^{(i)})\ge 0$$
$$又\because H({\theta }^{(i+1)},{\theta }^{(i)})-H({\theta }^{(i+1)},{\theta }^{(i)})=\sum _Z\log \frac{P(Z|Y,{\theta }^{(i+1)})}{P(Z|Y,{\theta }^{(i)})}P(Z|Y,{\theta }^{(i)})\le \log \sum _Z(\frac{P(Z|Y,{\theta }^{(i+1)})}{P(Z|Y,{\theta }^{(i)})}P(Z|Y,{\theta }^{(i)}))$$
$$=\log \sum _{Z}P(Z|Y,{\theta }^{(i+1)})=\log 1=0$$
  所以，定理1得证。

  定理2：设$L(\theta )=\log P(Y|\theta )$为观测数据的对数似然函数，${\theta }^{(i)}(i=1,2,...)$为EM算法得到的参数估计序列，$L({\theta }^{(i)})(i=1,2,...)$为对应的对数似然函数序列。则：
  1).如果$P(Y|\theta )$有上界，则$L({\theta }^{(i)})=\log P(Y|{\theta }^{(i)})$收敛到某一值$L^{\ast }$
  2).在函数$Q(\theta ,{\theta }^{\prime })$与$L(\theta )$满足一定条件下，由EM算法得到的参数估计序列${\theta }^{(i)}$的收敛值${\theta }^{\ast }$是$L(\theta )$的稳定点。
  第一条的证明显然$P(Y|\theta )\le 1$，再根据$L(\theta )$的单调性就可以得到上述结论；
  第二条的证明较为复杂，可参见文献“On the convergence properties of the EM algorithm”。

4. EM算法应用——求解高斯混合模型(GMM)的参数

4.1 高斯混合模型(Gaussian mixture model,GMM)

  高斯混合模型是指某随机变量的概率分布模型是多个不同高斯模型的加法模型：
$$P(y|\theta )=\sum _{k=1}^K{\alpha }_k\varphi (y|{\theta }_k)$$
  其中，${\alpha }_k$是系数，${\alpha }_k\ge 0,{\sum }_{k=1}^K{\alpha }_k=1$；$\varphi (y|{\theta }_k)$是高斯分布密度，${\theta }_k=({\mu }_k,{\sigma }_k^2)$
$$第k个分模型为：\varphi (y|{\theta }_k)=\frac{1}{\sqrt{2\pi }{\sigma }_k}\exp (-\frac{(y-{\mu }_k{)}^2}{2{\sigma }_k^2})$$

4.2 EM算法估计高斯混合模型的参数

  根据高斯混合模型的介绍可知，需要估计的参数为$\theta =({\alpha }_1,{\alpha }_2,...,{\alpha }_K;{\theta }_1,{\theta }_2,...,{\theta }_K)$，其中${\theta }_k=({\mu }_k,{\sigma }_k^2)$。
  1).首先要明确隐变量，写出完全数据的对数似然函数。
  可以这么理解高斯混合模型得到的数据，第一步是根据概率${\alpha }_k$选择第k个高斯分布分模型$\varphi (y|{\theta }_k)$,然后根据该分模型得到观测数据。因此，已知的是观测数据，隐参数是该数据是来自第几个分模型。所以设定隐变量如下：
$${\gamma }_{jk}=\{\begin{array}{l}1,第j个观测来自第k个分模型\\ 0,其他\end{array},j=1,2,...,N;k=1,2,...,K$$
  综上，完全数据就是$(y_j,{\gamma }_{j1},{\gamma }_{j2},...,{\gamma }_{jK}),j=1,2,...,N$
  所以，完全数据的对数似然函数为：
$$P(y,\gamma |\theta )=\prod _{j=1}^{N}P(y_j,{\gamma }_{j1},{\gamma }_{j2},...,{\gamma }_{jK}|\theta )$$
  理解可知，$({\gamma }_{j1},{\gamma }_{j2},...,{\gamma }_{jK})$中，只有一个为1，其他全为0，而$P(y_j,{\gamma }_{j1},{\gamma }_{j2},...,{\gamma }_{jK}|\theta )$是在${\gamma }_{jk}=1$时，用对应的分模型求出的，所以上式可以写成：
$$P(y,\gamma |\theta )=\prod _{k=1}^K\prod _{j=1}^N[{\alpha }_k\varphi (y_i|{\theta }_k){]}^{{\gamma }_{jk}}=\prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N[\varphi (y_i|{\theta }_k){]}^{{\gamma }_{jk}}=\prod _{k=1}^K{\alpha }_k^{n_k}\prod _{j=1}^N[\frac{1}{\sqrt{2\pi }{\sigma }_k}\exp (-\frac{(y_j-{\mu }_k{)}^2}{2{\sigma }_k^2}){]}^{{\gamma }_{jk}}$$
  其中，$n_k={\sum }_{j=1}^N{\gamma }_{jk},{\sum }_{k=1}^K{n}_k=N$，注意，$n_k$中是遍历下标j，而不是下标k。
  所以，完全数据的对数似然函数即为：
$$\log P(y,\gamma |\theta )=\sum _{k=1}^K\{n_k\log {\alpha }_k+\sum _{j=1}^N{\gamma }_{jk}[\log \frac{1}{\sqrt{2\pi }}-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}(y_j-{\mu }_k{)}^2\left]\right\}$$

  2).EM算法E步，确定Q函数
$$Q(\theta ,{\theta }^{(i)})=E[\log P(y,\gamma |\theta )|y,{\theta }^{(i)}]=E\left\{\sum _{k=1}^K\right\{n_k\log {\alpha }_k+\sum _{j=1}^N{\gamma }_{jk}[\log \frac{1}{\sqrt{2\pi }}-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}(y_j-{\mu }_k{)}^2\left]\right\}\}$$
$$=\sum _{k=1}^K\left\{\sum _{j=1}^N\right(E_{{\gamma }_{jk}})\log {\alpha }_k+\sum _{j=1}^N{E}_{{\gamma }_{jk}}[\log \frac{1}{\sqrt{2\pi }}-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}(y_j-{\mu }_k{)}^2]\}$$
  其中，记${\hat{\gamma }}_{jk}=E_{{\gamma }_{jk}}=E({\gamma }_{jk}|y,\theta )={\sum }_{k=1}^K{\gamma }_{jk}P({\gamma }_{jk}|y,\theta )=P({\gamma }_{jk}=1|y,\theta )$
$${\hat{\gamma }}_{jk}=E({\gamma }_{jk}|y,\theta )=P({\gamma }_{jk}=1|y,\theta )=\frac{P({\gamma }_{jk}=1,y_j|\theta )P(y_j|\theta )}{{\sum }_{t=1}^{K}P({\gamma }_{jt}=1,y_j|\theta )P(y_j|\theta )}=\frac{P({\gamma }_{jk}=1,y_j|\theta )}{{\sum }_{t=1}^{K}P({\gamma }_{jt}=1,y_j|\theta )}$$
  此处j是固定量，由于$({\gamma }_{j1},{\gamma }_{j2},...,{\gamma }_{jK})$中，只有一个为1，其他全为0，所以期望$E({\gamma }_{jk}|y,\theta )$是已知观测值和参数时，选中特定分模型的概率。
$$\therefore {\hat{\gamma }}_{jk}=\frac{P(y_j|{\gamma }_{jk}=1,\theta )P({\gamma }_{jk}=1|\theta )}{{\sum }_{t=1}^{K}P(y_j|{\gamma }_{jt}=1,\theta )P({\gamma }_{jt}=1|\theta )}=\frac{{\alpha }_k\varphi (y_j|{\theta }_k)}{{\sum }_{t=1}^K{\alpha }_t\varphi (y_j|{\theta }_t)}$$
  注意，上述推导中分母求和的变量原书里使用的是’k’，很容易和分子里的k混起来，不建议按原书的写法，容易引起歧义。另外，求解${\hat{\gamma }}_{jk}$所用到的参数${\alpha }_t,{\mu }_t,{\sigma }_t$等都是当前模型下已知的参数，因此${\hat{\gamma }}_{jk}$是可求的，同理下文中的$n_k$也是可求的。
  将${\hat{\gamma }}_{jk}=E_{{\gamma }_{jk}}$和$n_k={\sum }_{j=1}^N{E}_{{\gamma }_{jk}}$代入Q函数中可得
$$Q(\theta ,{\theta }^{(i)})=\sum _{k=1}^K\{n_k\log {\alpha }_k+\sum _{j=1}^N{\hat{\gamma }}_{jk}[\log \frac{1}{\sqrt{2\pi }}-\log {\sigma }_k-\frac{1}{2{\sigma }_k^2}(y_j-{\mu }_k{)}^2]\}$$
  3). 确定EM算法的M步
  求函数$Q(\theta ,{\theta }^{(i)})$对参数$\theta$的极大值，从而得到新一轮迭代的模型参数：
$${\theta }^{(i+1)}=\arg \underset{\theta }{\max }Q(\theta ,{\theta }^{(i)})$$
  用${\hat{\mu }}_t,{\hat{\sigma }}_t^2$和${\hat{\alpha }}_t,t=1,2,...,K$表示参数${\theta }^{(i+1)}$中的具体参数。求${\hat{\mu }}_t,{\hat{\sigma }}_t^2$只需将$Q(\theta ,{\theta }^{(i)})$对求${\hat{\mu }}_t,{\hat{\sigma }}_t^2$偏导并令其为0即可。
$$\frac{\partial Q(\theta ,{\theta }^{(i)})}{\partial {\hat{\mu }}_t}=\frac{\partial {\sum }_{j=1}^N{\hat{\gamma }}_{jt}[\log \frac{1}{\sqrt{2\pi }}-\log {\sigma }_t-\frac{1}{2{\sigma }_t^2}(y_j-{\mu }_t{)}^2]}{\partial {\hat{\mu }}_t}=\sum _{j=1}^N\frac{{\hat{\gamma }}_{jt}}{{\sigma }_t^2}(y_j-{\mu }_t)$$
$$令\frac{\partial Q(\theta ,{\theta }^{(i)})}{\partial {\hat{\mu }}_t}=0，得{\hat{\mu }}_t=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jt}{y}_j}{{\sum }_{j=1}^N{\hat{\gamma }}_{jt}},t=1,2,...K$$
$$\frac{\partial Q(\theta ,{\theta }^{(i)})}{\partial {\hat{\sigma }}_t^2}=\frac{\partial {\sum }_{j=1}^N{\hat{\gamma }}_{jt}[\log \frac{1}{\sqrt{2\pi }}-\log {\sigma }_t-\frac{1}{2{\sigma }_t^2}(y_j-{\mu }_t{)}^2]}{\partial {\hat{\sigma }}_t^2}=\sum _{j=1}^N-{\hat{\gamma }}_{jt}(\frac{1}{{\sigma }_t}\cdot \frac{1}{2}\cdot \frac{1}{{\sigma }_t}-\frac{1}{2{\sigma }^4}(y_j-{\mu }_t{)}^2)$$
$$令\frac{\partial Q(\theta ,{\theta }^{(i)})}{\partial {\hat{\sigma }}_t^2}=0,得{\hat{\sigma }}_t^2=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jt}(y_j-{\mu }_t{)}^2}{{\sum }_{j=1}^N{\hat{\gamma }}_{jt}},t=1,2,...K$$
  注意，上述推导中参数变量下标原书里使用的是’k’，求导时很容易Q函数里的求和下标k混起来，不建议按原书的写法，容易引起歧义。所以此处用t下标代替，下述推导也同样如此。
  对${\hat{\alpha }}_t$的求导略微复杂，因为该求导是有约束条件的求导，约束条件为${\sum }_{k=1}^K{\alpha }_k=1$。因此，需要使用拉格朗日乘子法求解极值。首先，去除$Q(\theta ,{\theta }^{(i)})$中与${\alpha }_k$无关的项：
$$Q(\theta ,{\theta }^{(i)})\to \sum _{k=1}^K\left\{n_k\log {\alpha }_k\right\}$$
  然后，用拉格朗日乘子法构造：
$$L({\alpha }_t)=\sum _{k=1}^K\left\{n_k\log {\alpha }_k\right\}+\beta (\sum _{k=1}^K{\alpha }_k-1)$$
  $L({\alpha }_t)$对${\alpha }_t$求导可得：
$$\frac{\partial L({\alpha }_t)}{\partial {\alpha }_t}=\frac{n_t}{{\alpha }_t}+\beta =0\Rightarrow \beta {\alpha }_t=-n_t\Rightarrow \sum _{t=1}^K\beta {\alpha }_t=-\sum _{t=1}^K\sum _{j=1}^N{\hat{\gamma }}_{jt}\Rightarrow \beta =-N$$
  上述推导中，有${\sum }_{t=1}^K{\sum }_{j=1}^N{\hat{\gamma }}_{jt}={\sum }_{j=1}^N{\sum }_{t=1}^K{\hat{\gamma }}_{jt}={\sum }_{j=1}^{N}1=N$。将$\beta =-N$代入偏导方程$\frac{\partial L({\alpha }_t)}{\partial {\alpha }_t}$可得：
$${\hat{\alpha }}_t=\frac{n_t}{N}=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jt}}{N},t=1,2,...,K$$
  4).重复上述2，3步，直到参数收敛即可。

  算法步骤总结：

• 输入：观测数据$y_1,y_2,..,y_N$，高斯混合模型；
• 输出：高斯混合模型参数；
  1). 设置参数的初始值，准备开始迭代；
  2).E步：根据当前模型参数，计算分模型$k$对观测数据$y_j$的响应度，即观测数据$y_j$是从分模型$k$中产生的概率；
  $${\hat{\gamma }}_{jk}=\frac{{\alpha }_k\varphi (y_j|{\theta }_k)}{{\sum }_{t=1}^K{\alpha }_t\varphi (y_j|{\theta }_t)},j=1,2,...,N;t=1,2,...,K$$
  3).M步：计算新一轮迭代的模型参数；
  $${\hat{\mu }}_t=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jt}{y}_j}{{\sum }_{j=1}^N{\hat{\gamma }}_{jt}},t=1,2,...K$$
  $${\hat{\sigma }}_t^2=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jt}(y_j-{\mu }_t{)}^2}{{\sum }_{j=1}^N{\hat{\gamma }}_{jt}},t=1,2,...K$$
  $${\hat{\alpha }}_t=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jt}}{N},t=1,2,...K$$
  4).重复上述2，3步，直到参数收敛即可。

5.EM算法的推广——广义期望极大算法(GEM)

  待补充