这是阅读周志华教授《机器学习》中第七章（贝叶斯分类器）中EM算法部分的笔记，书上总结得很精炼，有很多细枝末节的东西需要填充。
PDF版：http://pcoln8jiu.bkt.clouddn.com/em.pdf
简单总结
书中将EM算法总结为估计参数隐变量的方法，它是一种迭代式算法：
若参数$\Theta$已知，则根据训练数据$D$推断出最优隐变量$Z$；若$Z$的值已知，则可对$\Theta$求极大似然估计。
步骤简单地可以总结为下面两步：
E步：利用当前估计的参数值$\Theta$计算对数似然的期望。
M步：寻找可以使E步产生的似然期望最大的$\Theta$。
上面两个步骤交替执行，直至收敛。
补充
强烈建议观看徐亦达教授的EM算法视频：https://www.bilibili.com/video/av23901379
这部分之前写的不是很清楚，于是我决定重新看一下徐老师的视频和PPT，认真把EM算法框架推导和证明过一遍，再推导出高斯混合模型带入框架中的具体形式。
EM算法是一个非常好的解决计算困难问题的思路，巧妙地引入隐变量（latent variables）避开一些复杂计算。首先从高斯混合模型引入问题：
高斯混合模型（Gaussian Mixture Model, GMM）
基于贝叶斯估计通常首先假设数据中的特征之间独立同分布。我们遇到数据集的分布无法仅用一个Gaussian Distribution拟合时通常考虑使用混合模型，Gaussian mixture model就是将多个高斯模型加权求和达到对拟合数据的目的。使用GMM描述一个数据集的分布可以表示成：
$$
p(X)=\sum_{i=1}^k\alpha_iN(X|\mu_i,\Sigma_i)
$$
$$
\sum\alpha_i=1
$$
$k$为高斯分布个数，$\Theta$表示高斯混合模型，$\alpha$为每个高斯分布的权，$N$表示一个高斯分布的概率值，定义为$f(x;\theta_i)$，是参数为$\theta_i$的高斯分布概率密度函数$\dfrac{1}{\sqrt{2\pi}\Sigma_i}exp(-\dfrac{(x-\mu_i)^2}{2\Sigma_i^2})$。$\mu、\Sigma$分别为每个高斯分布的参数。
对于整个GMM $\Theta$，我们需要确定的参数有$2k-1$个，即$\mu_1,…,\mu_k,\Sigma_1,…,\Sigma_k,\alpha_1,…,\alpha_{k-1}$.
插一句关于$\alpha_i$的来历：因为数据分布（或者说宽度）并不一定非常平均，我随便搜了一张图来解释一下：
比如这个图，很明显两个高斯分布比例不一样，因此用权重更合适。
我思考为什么权值和是1，参考了问答[1]：因为GMM的定义本质上是一个概率密度函数，概率密度函数在负无穷到正无穷内的积分为1，因此我们首先要满足概率密度函数的性质。
我们希望GMM更好地拟合数据集，也就是希望$p(x)$尽可能地大，于是通常利用对数似然估计对高斯混合模型进行调参，使得拟合特定的数据集。
这里再插一句：为什么不用一般的极大似然估计(MAP)？
这是我和师兄同学在饭桌上讨论过的一个问题，总结一下，原因有两条：
- 因为极大似然估计是把所有概率乘在一起，每个概率值都是小数，数据量大了小数位精度会丢得非常厉害。
- 将连乘换成对数和的形式，求导方便。
针对高斯混合模型，我们构造损失函数$L(\Theta)$：
$$
L(\Theta|X)=\sum_{i=1}^n \log p(X)
$$
即：
$$
L(\Theta|X)=\sum_{i=1}^n \log \sum_{j=1}^k\alpha_jN(x_i|\mu_j,\Sigma_j)
$$
于是我们的优化目标就是：
$$
\Theta^{MLE}=\mathop{argmax}_\Theta L(\Theta|X)
$$
虽然上式可导，但是无法一步到位，计算过程过于复杂。
考虑是不是可以有一个迭代的方法，求得$\Theta^{(1)},\Theta^{(2)},…,\Theta^{(t)}$，使得我们最后收敛？
这是推导EM需要的前置知识，摆在这里。
Jensen不等式[2]
设$f$是定义域为实数的函数，如果对于所有的实数$x$，$f(x)$的二次导数大于等于0，那么$f是$凸函数。当$x$是向量时，如果其Hessian矩阵H是半正定的，那么$f$是凸函数。如果只大于0，不等于0，那么称$f$是严格凸函数。
Jensen不等式表述如下：
如果f是凸函数，$X$是随机变量，那么：$E[f(X)]\geq f(E[X])$，特别地，如果$f$是严格凸函数，当且仅当$X$是常量时，上式取等号。
图中，实线$f$是凸函数，$X$是随机变量，有0.5的概率是$a$，有0.5的概率是$b$。$X$的期望值就是$a$和$b$的中值了，图中可以看到$E[f(X)]\geq f(E[X])$成立。
顺便，中国数学界关于函数凹凸性定义和国外很多定义是反的。国内教材中的凹凸，是指曲线，而不是指函数，图像的凹凸与直观感受一致，却与函数的凹凸性相反[3].
EM算法
直接给出EM算法需要迭代解决的基本任务：
$$
\Theta^{(t+1)}=\mathop{argmax}_\Theta\int_z \log p(X,z|\Theta)p(z|X,\Theta^{(t)}) \mathop{dz}
$$
当然这是最普遍的定义，如果放在离散问题上，那么$z$的数量是有限的于是可以写成：
$$
\Theta^{(t+1)}=\mathop{argmax}\Theta\sum{i=1}^k\log p(X,z_i|\Theta)p(z_i|X,\Theta^{(t)})
$$
接下来会一步一步地推出上面这个迭代式。
我再摆一下原来的优化目标：
$$
\Theta^{MLE}=\mathop{argmax}\Theta\log \sum{j=1}^k\alpha_jN(X|\mu_j,\Sigma_j)
$$
我们发现有2个不一样的地方：
- 怎么多了个变量$z$？
- 这个$\sum$怎么跑$\log$左边去了？
怎么多了个变量？这对原目标不会有影响么？
首先，其实$z$是辅助变量，称为隐变量（latent variable），引入这个东西是EM算法的一个亮点。我记得在数据挖掘课上徐君老师点评SVM算法，说这是一个典型的“看似化简为繁”的工作，但是有了这一步，后面的计算往往会更加简便。
这和EM算法异曲同工，EM算法的提出者也看得很远，他对需要优化的损失函数引入辅助变量，而且必须保证加和不加的分布相同（保证加了$z$之后的边缘分布与原分布一致），也就是说，给出数据的概率分布，必须保证有：
$$
p(X)=\sum_{i=1}^kp(X|\Theta,z_i)p(\Theta,z_i)
$$
$p(X|\Theta,z_i)$表示对数据判断的似然（对应GMM中的概率密度函数$N$）
$p(\Theta,z_i)$表示先验（对应GMM中的权重$\alpha_i$）
我们这么定义，因此合适的$p$都不会有影响。
这个$\sum$怎么跑$\log$左边去了？
保留这个问题，首先分析我们的优化目标。由于是迭代的，那么最基本的要求就是我们每次迭代都要使当前结果更加接近最终结果，也就是：
$$
\log p(X|\Theta^{(t+1)})\geq \log p(X|\Theta^{(t)})
$$
我们的目标是推出这样一个式子，此时已经没有隐变量$z$了。我们希望通过优化目标下式来得到上式的结果：
$$
\Theta^{(t+1)}=\mathop{argmax}\Theta\sum{i=1}^k\log p(X,z_i|\Theta)p(z_i|X,\Theta^{(t)})
$$
这个假设非常大胆，我们很轻易地发现这里其实蕴含了一个矛盾，我们定义$\Theta$是最终的优化结果，而$\Theta^{(t)}$却是当前步的最优目标，我们在计算过程中实际上并不知道$\Theta$，这完全是个Chicken and eggs problem。
如何处理并不是很显然，我们尝试推一下。首先知道：
$$
p(X)=\dfrac{p(X)p(XZ)}{p(XZ)}=\dfrac{p(XZ)}{p(Z|X)}
$$
于是：
$$
\log p(X|\Theta)=\log \dfrac{p(XZ|\Theta)}{p(Z|X\Theta)}= \log p(XZ|\Theta)-\log p(Z|X\Theta)
$$
那么它们在$Z$下的期望也一定相等，即：
$$
E_{p(Z|X\Theta^{(t)})}[\log p(X|\Theta)]= E_{p(Z|X\Theta^{(t)})}[\log p(XZ|\Theta)-\log p(Z|X\Theta)]
$$
等式左边展开期望的定义：
$$
E_{p(Z|X\Theta^{(t)})}[\log p(X|\Theta)]=\int_Z[\log p(X|\Theta)]p(z|X\Theta^{(t)})dz
$$
$$
=\log p(X|\Theta)\int_zp(z|X\Theta^{(t)})dz
$$
$$
=\log p(X|\Theta)
$$
等式右边，由期望可加性拆成两项差，期望定义展开：
$$
=\int_Z\log p(XZ|\Theta)p(z|X\Theta^{(t)})dz-\int_Z\log p(Z|X\Theta)p(z|X\Theta^{(t)})dz
$$
这是整个EM算法里最tricky的部分（其次是引入latent variable），仔细观察就会发现EM算法迭代优化的只是前面那个积分$\int_Z\log p(XZ|\Theta)p(Z|X\Theta^{(t)})dz$，所以接下来就是消掉后面这个积分，这个被消掉的过程是有点厉害的。。
前半部分定义给个记号：$Q(\Theta,\Theta^{(t)})$，后半部分定义为：$H(\Theta,\Theta^{(t)})$.
重申下我们的目标，希望推出最终优化目标为：
$$
\mathop{argmax}_{\Theta}[\int_Z\log p(XZ|\Theta)p(z|X\Theta^{(t)})dz]=\Theta^{(t)}
$$
那么到$t+1$步的时候，$Q$和$H$会变大还是变小？
讨论$Q$：（因为我们的最终目标就只包含Q，所以Q肯定单调增加啦）显然$Q(\Theta^{(t+1)},\Theta^{(t)})Q>(\Theta,\Theta^{(t)})$.
讨论H：如果$H(\Theta^{(t)},\Theta^{(t)})\geq H(\Theta^{(t+1)},\Theta^{(t)})$，即$H$是单调递减的话就可以将优化目标的后半部分的积分约掉了，下面证明。
证明$H(\Theta^{(t)},\Theta^{(t)})\geq H(\Theta,\Theta^{(t)})$
$$
\int_z\log[p(z|X\Theta^{(t)})]p(z|X\Theta^{(t)})dz-\int_z\log[p(Z|X\Theta)]p(z|X\Theta^{(t)})dz
$$
即：
$$
\int_z\log[\dfrac{p(z|X\Theta^{(t)})}{p(Z|X\Theta)}]p(z|X\Theta^{(t)})dz
$$
希望上式$\geq0$，ummmm…似乎做不下去了。
冷静分析，上面实际上是个函数的期望$E_Z(log[\dfrac{p(z|X\Theta^{(t)})}{p(Z|X\Theta)}])$，里面有个$\log$实在太复杂，能不能把$\log$拿出来，比如写成$\log E(\dfrac{p(z|X\Theta^{(t)})}{p(Z|X\Theta)})$.
我们再仔细考虑，诶？这个分子是$p(z|X\Theta^{(t)})$，跟$Z$下的概率一样嘛。于是保证这个函数的单调性的同时做一下变换：
$$
E_Z(log[\dfrac{p(z|X\Theta^{(t)})}{p(Z|X\Theta)}])=E_Z(-log[\dfrac{p(Z|X\Theta)}{p(z|X\Theta^{(t)})}])
$$
我们要把$log$拿出来，于是作者想到了用Jensen不等式$E[f(X)]\geq f(E[X])$，函数的期望≥期望的函数，带进去就OK了：
$$
E_Z(-log[\dfrac{p(Z|X\Theta)}{p(z|X\Theta^{(t)})}])\geq -\log E_Z(\dfrac{p(Z|X\Theta)}{p(z|X\Theta^{(t)})})
$$
展开不等式右边：
$$
-\log E_Z(\dfrac{p(Z|X\Theta)}{p(z|X\Theta^{(t)})})=\int_Z \dfrac{p(Z|X\Theta)}{p(z|X\Theta^{(t)})}p(z|X\Theta^{(t)})dz=-log(1)=0
$$
于是：
$$
E_Z(-log[\dfrac{p(Z|X\Theta)}{p(z|X\Theta^{(t)})}])\geq -\log E_Z(\dfrac{p(Z|X\Theta)}{p(z|X\Theta^{(t)})})=0
$$
单调减的性质就证出来了。
所以我们的优化目标就是：
$$
\Theta^{(t+1)}=\mathop{argmax}\Theta\sum{i=1}^k\log p(X,z_i|\Theta)p(z_i|X,\Theta^{(t)})
$$
（$\sum$竟然跑到$\log$的左边去了）
Reference：
徐老师的视频链接：https://www.bilibili.com/video/av23901379
[1]. https://www.quora.com/Why-are-all-of-the-weights-in-Gaussian-mixture-models-supposed-to-be-positive