EM算法是一种迭代算法，是一种用于计算包含隐变量概率模型的最大似然估计方法，或极大后验概率。EM即expectation maximization，期望最大化算法。

1. 极大似然估计

  在概率模型中，若已知事件服从的分布或者其他概率模型的参数，那么我们可以通过计算得到某事件发生的概率。而在估计中，这些变成了方向过程：已知一组数据发生的结果，相当于获得了经验概率，通过这组数据假设模型服从什么分布，再通过经验概率求解模型参数。
  比如统计学校学生身高服从的概率分布，抽样1000人得到他们的身高数据，再假设身高服从正态分布，于是只需要求解\(u,\sigma\)即可。剩下的问题就是如何确定这两个参数，极大似然估计的思想就是求解在这两个参数下，得到的这组样本发生的概率最大的情况便是这两个参数最有可能的取值。而抽取每个人的身高都是独立的，让每个p最大化就是让其乘积最大化，于是得到最大似然函数
\[L(\theta) = \prod _xp(x | \theta)\]  再求解其关于\(\theta\)的极大化问题便能得到一个对\(\theta\)的估计
\[\widehat \theta = argmaxL(\theta)\]

2.EM算法

  上面是模型变量都是可观测的情况，这种情况可以直接根据观测到的样本数据，使用极大似然估计对模型参数进行估计。但若模型中含有隐变量的时候，就不能简单的使用极大似然估计进行计算。EM算法就是针对含有隐变量的概率模型参数的极大似然估计方法。
例子：
  假设三枚硬币，记为A、B、C，它们正面向上的概率分别为n、p、q。实验如下：先抛A，根据A的结果选择B或者C再次抛，将这次正面向上的结果记为1，反面记为0，其中A正面向上则选择B，反面则选择C。经过n次实验后得到n个观测值，根据这n个观测值估计n、p、q的参数值。实验过程中只能观测到最后结果，并不能观察实验过程，也就是A的结果是未知的。该模型可以表示如下
\[P(y|\theta) = \sum_zP(y,z|\theta) = \sum_zP(z|\theta)P(y|z,\theta)\]  其中y表示随机变量Y的观测值，表示该次结果是1或0。z代表隐变量，表示模型中未能观测到的A的结果，\(\theta=（n,p,q)\)是模型参数。其中Y是不完全数据，Y+Z是完全数据。
  若是没有隐变量Z，那么可直接使用对数极大似然函数估计
\[\widehat \theta=arg\underset \theta {max}L(\theta)=arg\underset \theta {max}\sum_y logP(y|\theta) \]  加入隐变量后，对数极大似然函数变为
\[L(\theta,z)=\sum_ylog\sum_zP(y,z|\theta) \]  求解
\[ \widehat \theta,z=arg\underset {\theta,z} {max}L(\theta,z) \]  按照极大似然估计中的方法，似乎应该对\(\theta,z\)求导然后令其为0解方程组，然后注意此时的\(L(\theta,z)\)函数，log内含有求和运算，若是直接求导计算十分困难，因此退而求其次，既然要极大化\(L(\theta,z)\)，就先找到一个它的下界，再想办法提高这个下界，同样能达到极大化L的效果。使用Jense不等式对似然函数进行放缩。

• Jense不等式（琴生不等式）

  凸函数：是定义在实数域上的函数，如果对于任意的实数，都有：
\[f''\geqslant 0\]   则该函数是凸函数。
  当函数是凸函数的时候，Jense不等式的含义是函数的期望大于等于期望的函数（凹函数相反）。图下图所示

  二维情况下可用凸函数定义来解释，当一个函数是凸函数的时候，它满足
\[f(\frac {a+b} 2)\leqslant \frac {f(a)+f(b)} 2\]  左边其实相当于其变量x先求期望后带入函数，右边是先求函数值再计算函数值的期望，也就是
\[f(E(x))\leqslant E(f(x)) \]  再回到\(L(\theta,z)\)中来，目的是为了将对数函数中的和项去掉，便可利用jense不等式的性质，将期望的函数变为函数的期望。先进行放缩
\[ \begin{aligned} L(\theta,z)&=\sum_ylog\sum_zP(y,z|\theta)\\ &=\sum_ylog\sum_zQ_i(z)\frac {P(y,z|\theta)} {Q_i(z)}\\ &\geqslant \sum_y \sum_zQ_i(z)log\frac {P(y,z|\theta)} {Q_i(z)}\\ \end{aligned} \]  其中最后一步用到了Jense不等式，因为对数函数是凹函数，所以不等号反了过来，\(f(E(x))\geqslant E(f(x))\),此处\[f(x) = log(x), x=\frac {P(y,z|\theta)} {Q_i(z)} ，\]  并且所添加的\(Q_i(x)\)满足\[\sum_zQ_i(z) = 1\]  这是根据第三类Jense不等式的条件设定的，不同系数的加权求和期望只要满足系数之和为1就能使用Jense不等式。
  所以得到结论，\(log\frac {P(y,z|\theta)} {Q_i(z)}\)的加权平均就是\(L(\theta,z)\)的一个下界。这便是EM算法中E（期望）的来由。
  目前\(Q(z)\)还是未知的，需要根据一些条件来选择一个合适的函数，再次强调最终目的是极大化似然函数，现在我们得到了似然函数的一个下界，一个想法就是让这个下界去更好的逼近似然函数的真实值，下界是根据不等式放缩后得到的，那么当不等式取等号的时候便是最接近原函数的时候。回忆Jense不等式\(f(E(x)) \leqslant E(f(x))\)，显然当\(x\)为常数的时候不等式成立。即
\[x =\frac {P(y,z|\theta)} {Q_i(z)}=c\]  既然要确定的是\(Q_i(z)\)，不妨设\(\theta\)为常数，事实上，EM是一种迭代算法，也就是我们是先给出\(\theta\)的假设值后再进行迭代计算的，设上一次为第i次迭代，\(\theta ^{(i)}\)为第i次的参数值，由上式得
\[ \begin{aligned} P(y,z|\theta ^{(i)}) &=cQ_i(z) \\ \sum_zP(y,z|\theta ^{(i)})&=\sum_zcQ_i(z) \\ \sum_zP(y,z|\theta ^{(i)})&=c \end{aligned} \]
  将c带入x
\[ \begin{aligned} Q_i(z) &= \frac {P(y,z|\theta ^{(i)})} {\sum_zP(y,z|\theta ^{(i)})} \\ &=\frac {P(y,z|\theta ^{(i)})} {P(y|\theta ^{(i)})} \\ &= P(z|y,\theta ^{(i)}) \end{aligned} \]
  第二步用到边缘概率，第三步条件概率。至此，我们得出了Q(z)的表达式，Q(z)是已知样本观测和模型参数下的隐变量的分布。将Q_i(z)视作常数去掉后，得到下面表达式

\[ \begin{aligned} L(\theta,z)&=\sum_y \sum_zQ_i(z)log\frac {P(y,z|\theta)} {Q(z)}\\ &=\sum_y \sum_zQ_i(z)logP(y,z|\theta)\\ &=\sum_y\sum_zP(z| y,\theta ^{(i)} )logP(y,z|\theta) \end{aligned} \]  到这儿已经完成了E步，对对数似然函数下界的期望的表示，这个期望是\(logP(y,z|\theta)\)关于\(Q_i(z)=P(z|y,\theta^{(i)})\)的期望，接下来需要做的便是极大化这个期望，也就是M步。极大化的方法就是我们常见的优化问题，使用拉格朗日乘数法即可，添加约束条件\(\sum_zP(z|y,\theta ^{(i)})=1\)后对\(\theta,z,\lambda\)求偏导令为0。求解
\[\widehat \theta = arg\underset \theta {max}L(\theta,z)\]  迭代的停止条件常常为\(||\theta^{(i+1)}-\theta^{(i)}|| \leqslant \varepsilon ,or,||L(\theta^{(i+1)})-L(\theta^{(i)}|| \leqslant \varepsilon\)

3.EM算法的收敛性

  观测数据的似然函数是\(L(\theta)=P(y|\theta)\)，我们先证明似然函数关于每次迭代后的\(\theta\)是单调递增的
proof：
\[ P(Y|\theta) = \frac {P(Y,Z|\theta)} {P(Z|Y,\theta)} \]取对数
\[ logP(Y|\theta) = logP(Y,Z|\theta)-logP(Z|Y,\theta) \]由前面结论
\[ L(\theta,z) =\sum_zP(Z|Y,\theta^{(i)})logP(Y,Z|\theta) \]令
\[ H(\theta,z) = \sum_zP(Z|Y,\theta^{(i)})logP(Z|Y,\theta) \]于是，对数似然函数可写作
\[logP(Y|\theta) = L(\theta,z) - H(\theta,z)\]那么
\[logP(Y|\theta^{(i+1)})-logP(Y|\theta^{(i)}=[L(\theta^{(i+1)},z)-L(\theta^{(i)},z)]-[H(\theta^{(i+1)},z)-H(\theta^{(i)},z)] \]只需证右边非负，因\(\theta^{(i+1)}\)使原似然函数\(L(\theta^{(i+1)},z)\)达到极大，所以
\[L(\theta^{(i+1)},z)-L(\theta^{(i)},z) \geqslant 0 \]对于第二项
\[ \begin{aligned} &H(\theta^{(i+1)},z)-H(\theta^{(i)},z)\\ &=\sum_zP(Z|y,\theta^{(i)})log \frac {P(Z|Y,\theta^{(i+1)})} {P(Z|Y,\theta^{(i)})}\\ (Jensen 不等式) &\leqslant log\left ( \sum_zP(Z|Y,\theta^{(i)})\frac {P(Z|Y,\theta^{(i+1)})} {P(Z|Y,\theta^{(i)})} \right )\\ &=log\sum_zP(Z|Y,\theta^{(i+1)})\\ &=log1=0 \end{aligned} \]  因此\(logP(Y|\theta)\)关于迭代的\(\theta\)单调递增，那么如果\(P(Y|\theta)\)有上界，经过\(\theta\)的单增迭代一定会收敛到某一值，还有一个定理是EM算法得到的\(\theta\)收敛值是似然函数的稳定点，人生苦短，证明就免了吧。

4.EM算法在高斯混合模型中的参数估计

  通过前面的介绍已经知道EM算法用于解决含隐变量的概率模型，EM算法可以认为是一种解决问题的思想，在不知道隐变量发生情况的条件下，我们通过计算已知观测数据，求出不同数据在隐变量的某个事件下发生的期望，再通过这个期望去最大化该期望下的似然函数，如此迭代下去便能得到一个局部极大的结果。
  高斯混合模型也是这样的一个概率模型，它由多个正态分布加权组合而成。也就是可以观测到混合模型产生的样本的分布，但是不知道这每一个样本是由其中哪个高斯模型产生的。因此，隐变量就是哪个样本选择了哪个模型。

• 高斯混合模型

  假设观测数据\(y_1,y_2,...,y_n\)由高斯混合模型产生
\[ P(y|\theta) = \sum_{k=1} ^K \alpha_k \Phi(y|\theta_k) \]  其中\(\theta = (\alpha_1,...,\alpha_k;\theta_1,...,\theta_k)\)，由前分析得隐变量是每个样本来自哪个模型，j个样本和k个模型，隐变量\(\gamma_{jk}\)表示如下
\[\gamma_{jk}＝ \left\{\begin{aligned} &1,&第j个观测来自第k个分模型 \\& 0,&否则 \end{aligned}\right.\]  于是有\(j\)个观测数据\(j*k\)个位观测数据，总共\((j+1)*k\)个数据便组成了完全数据，完全数据的似然函数为
\[ \begin{aligned} P(y,\gamma|\theta)&=\prod_j^NP(y_j,\gamma_{j1},\gamma_{j2},...,\gamma_{jk}|\theta)\\ &=\prod_k \prod_j[\alpha_k \Phi(y_j|\theta_k)]^{\gamma_{jk}}\\ &=\prod_k \alpha_k^{n_k} \prod_j[\Phi(y_j|\theta_k)]^{\gamma_{jk}}\\ \end{aligned}\]  其中\(n_k = \sum_{j=1}^N\gamma_{jk},\sum_{k=1}^Kn_k = N\),
\[\Phi(y_j|\theta_k) = \frac {1} {\sqrt{2\pi}\delta_k}exp\left ( \frac {(y_i-u_k)^2} {2{\delta_k} ^2}\right )\]  显然\(n_k\)表示的是第k个模型出现的样本数量，因此所有模型出现的样本数量和就是N，完全数据的对数似然函数为
\[ logP(y,\gamma|\theta) = \sum_k \left \{ n_klog\alpha_k +\sum_j \gamma_{jk}[log(\frac 1 {\sqrt{2\pi} } -log\delta_k-\frac {(y_j-u_k)^2} {2{\delta_k}^2})] \right \} \]  确定好完全数据的对数似然函数后就可以直接套用EM算法的框架了，首先是E-step，求解对数似然函数关于隐变量的后验概率的期望
\[ \begin{aligned} L(\theta,\gamma) &=E[logP(y,\gamma|\theta);y,\theta^{(i)}]\\ &=E\left \{ \sum_k \left \{ n_klog\alpha_k +\sum_j \gamma_{jk}[log(\frac 1 {\sqrt{2\pi} } -log\delta_k-\frac {(y_j-u_k)^2} {2{\delta_k}^2})] \right \} \right \}\\ &=\sum_k \left \{ \sum_jE (\gamma _{jk})log\alpha_k +\sum_j E\gamma_{jk}[log(\frac 1 {\sqrt{2\pi} } -log\delta_k-\frac {(y_j-u_k)^2} {2{\delta_k}^2})] \right \} \end{aligned} \]  注意其中的\(n_k\)是关于\(\gamma_{jk}\)的函数\(n_k=\sum_j\gamma_{jk}\)，本次期望相当于\(E(P(\gamma;y,\theta))\),而对于离散随机变量\[E(g(X))=\sum_ig(x_i)*p_i\]  计算\(E(\gamma_{jk})\),记为\(\widehat \gamma\)
\[ \begin{aligned} \widehat \gamma &= E(\gamma_{jk}|y,\theta)\\ &= \frac {P(\gamma_{jk} = 1,y_j | \theta)} {\sum_kP(\gamma_{jk}=1,y_j |\theta)}\\ &= \frac {P(\gamma_{jk} = 1,y_j | \theta) P(\gamma_{jk}=1|\theta)} {\sum_kP(\gamma_{jk}=1,y_j |\theta)P(\gamma_{jk}=1|\theta)}\\ & = \frac {\alpha_k \Phi(y_j|\theta_k)} {\sum_k \alpha_k \Phi(y_j|\theta_k)} \end{aligned} \]  带入\(L(\theta,\gamma)\)得
\[L(\theta,\gamma) = \sum_k \left \{ n_klog\alpha_k +\sum_j \widehat \gamma_{jk}[log(\frac 1 {\sqrt{2\pi} } -log\delta_k-\frac {(y_j-u_k)^2} {2{\delta_k}^2})] \right \}\]  接下就是M步，只需最大化这个L即可
\[\theta^{(i+1)} = arg \underset \theta {max} L(\theta,\gamma)\]  求解方式同样是拉格朗日，约束条件\(\sum_k\alpha_k = 1\)