EM 算法（Expectation-Maximization algorithm）是一种最大期望算法。主要包括 E-step 和 M-step 两部分，通过反复迭代来获得相关参数 $\theta$ 的估计。
贝叶斯估计
对于 EM 算法中需要利用到贝叶斯估计公式，总结如下：
$$ \mathrm{P}(A\mid B) = \frac{\mathrm{P}(B \mid A)\mathrm{P}(A)}{\mathrm{P}(B)} $$
$$ \mathrm{P}(A, B) = \mathrm{P}(A\mid B)\mathrm{P}(B) = \mathrm{P}(B\mid A)\mathrm{P}(A) $$
E-step 与完全数据对数似然 （L 函数）
E-step 描述求期望的过程，是 EM 算法的核心，主要包括两个函数，一个被称为完全数据的对数似然，公式表达为：
$$ \log{L(X, Z, \theta)} = \log{\mathrm{P}(X, Z|\theta)} $$
其中 $X$ 为观察数据，$Z$ 为隐参数而 $\theta$ 为模型参数。描述在 $\theta$ 模型参数下，实验数据 $X$ 与 隐参数 $Z$ 之间的联合分布。
另一个参数用来描述完全数据的对数似然期望，公式表达为：
$$ \begin{aligned} Q(\theta | \theta_{old}) & = \mathrm{P}(Z|X, \theta)\cdot \log{L(X,Z,\theta)} \ \\ &= \sum_{Z}{\mathrm{P}(Z|X, \theta_{old}}) \cdot \log{L(X, Z, \theta)} \end{aligned} $$
为解释 $Q(\theta\mid \theta_{old})$ 从何而来，需要利用到 Jenson 不等式。
Jenson 不等式
Jenson 不等式可对凹凸函数进行研究，假设 $\mathrm{g}(x)$ 为凹函数，权重 $\lambda_i \geq 0$，且 $\sum{\lambda_i} = 1$ ，则有如下关系：
$$ \mathrm{g}(\sum_{i}\lambda_i x_i) \leq \sum_i \lambda_i \mathrm{g}(x_i) $$
Jenson 不等式证明
引入支撑超平面概念：对于凸函数 $f(x)$ ，在函数曲线上任一点可以知道支撑超平面（切线）。引入一点 $z$ ，对于权重 $0 \leq \lambda \leq 1$ ，假设点 $z$ 在点 $x$ 与点 $y$ 之间，则有
$$ z = \lambda x + (1-\lambda)y $$
给出曲线上 $z$ 点处的切线函数式：
$$ t(z) = f(x) + \nabla f(x)^T (z-x) $$
代表函数梯度向量 $\nabla f(x)$ 在从点 x 到点 z 之间的向量投影，描述点 x 到点 z 的向量方向上函数值的变化程度。 同样也会有,
$$ f(z) \leq f(x) + \nabla f(x)^T(z-x) $$
这是凸函数上任一点的切线的直观理解，而假设这一点 $z$ 来到凸函数 $x$ 与 $y$ 之间线段上，则与切线情况不同，即
$$ f(z) \geq \lambda f(x) + (1-\lambda)f(y) $$
带入 $z$ 的定义值，有
$$ f(\lambda x + (1-\lambda)y) \geq \lambda f(x) + (1-\lambda)f(y) $$
考虑如果存在多个点 $x_i$，对应权重分别为 $\lambda_i$ ，则类似于 $x$ 、$y$ 两点的关系，可以得到 Jenson 不等式的结论：
$$ \mathrm{g}(\sum_{i}\lambda_i x_i) \leq \sum_i \lambda_i \mathrm{g}(x_i) $$
需要注意的是，这里的切线的情况并不是证明 Jenson 不等式所必须的，只是另一种临界条件帮助理解。
E-step 与 Q 函数
在证明 Jenson 不等式后，我们可以通过完全数据对数似然 $L(X, Z, \theta )$ 来研究其在隐参数变量 $Z$ 上的期望，引入一个新的函数，$Q(\theta \mid \theta_{old})$ 即在 $\theta$ 参数下，观察数据 $X$ 以及 隐参数 $Z$ 联合概率密度分布的期望。
$$ \begin{aligned} Q(\theta \mid \theta_{old}) &= \text{E}_{Z}\left[L(X,Z,\theta)\right] \\ \\ & = \sum_{Z}{\mathrm{P}(Z\mid X,\theta_{old})}\log{L(X,Z,\theta)} \end{aligned} $$
这里会有一些奇怪的地方，即为何突然出现 $Q(\theta \mid \theta_{old})$ ，这个函数的意义为何？
为了解释清楚这个问题，我们需要从通理解 EM 算法的算法目的，即求解似然概率：
$$ EM \rightarrow \arg\max_{\theta} P(X\mid \theta) $$
以图像重构为例，即在给出相关参数（ $\theta$ ）的前提下，不断迭代出概率最大的观察数据。
而 $Q$ 函数是一种人为设立的辅助函数，通过利用 $\theta_{old}$ 这一参数的不断迭代使得似然函数不断增加。而之所以这么折腾一下，是因为在复杂模型中，存在很多的隐参数 $Z$ ，会导致对于似然函数的直接求解存在很大的困难。
那么进一步分析为什么使用 $Q(\theta, \theta_{old})$ 函数，首先之前已经给出了 EM 的目的，即
$$ \arg\max_{\theta} P(X\mid \theta) $$
所以先对 $P(X\mid \theta)$ 开始着手：
$$ P(X\mid \theta) = \sum_{Z}P(X,Z\mid \theta) $$
代表 $X$ 和 $Z$ 之间联合概率在所有隐函数 $Z$ 取值求和，即似然概率。等式两端同时取对数有，
$$ \begin{aligned} \log{P(X\mid \theta)} &= \log \sum_{Z}P(X,Z\mid \theta) \\ \\ & = \log{\sum_{Z}}{P(Z\mid X, \theta)}\cdot \frac{P(X,Z\mid \theta)}{P(Z\mid X, \theta)} \end{aligned} $$
由 Jenson 不等式，可以的出
$$ \log{P(X\mid \theta)} \geq \sum_{Z}P(Z\mid X, \theta)\log{\frac{P(X,Z \mid \theta)}{P(Z\mid X, \theta_{old})}} $$
这样我们令
$$ Q(\theta, \theta_{old}) = \sum_{Z}P(Z\mid X, \theta_{old})\log{P(X,Z\mid \theta)} $$
而实际上，完整的带有 $Q$ 函数的 Jenson 不等式应该为：
$$ \log{P(X\mid \theta)} \geq \sum_{Z}P(Z\mid X, \theta_{old})\log{P(X,Z \mid \theta)} - \sum_{Z}P(Z\mid X, \theta_{old})\log{P(Z\mid X, \theta_{old})} $$
因为这里我们主要想关注的是如何最大化 $\theta^{(t+1)}$ ，所以这里后面部分是无关项（求偏导后会因为常数被约去）。这样就有了完整的不等式形式：
$$ \log{P(X\mid \theta)} \geq Q(\theta, \theta_{old}) + H\left(P(Z\mid X, \theta_{old})\right) $$
其中 $H$ 表示隐变量 $Z$ 的后验概率熵，是一个衡量随机变量不确定性的度量，数学上定义即为
$$ H(Z) = -\sum_{Z}P(Z)\log{P(Z)} $$
M-step
在构建 $Q$ 函数后，就可以通过 $\theta_{old}$ 最大化 $Q$ 函数来对 $\theta$ 参数进行迭代，为了更好的帮助理解迭代，用 $\theta^{t}$ 代替 $\theta_{old}$ :
$$ \theta^{(t+1)} = \arg\max_{\theta}{Q(\theta, \theta^{t})} $$
简单的说，例如对于较为常见的高斯混合模型，可以通过求偏导为零的方法，分别更新 $\theta$ 参数中包括的期望 $\mu$ 、标准差 $\sigma$ 以及混合系数 $\pi$ 。接下来分别展开详细介绍：
混合系数 $\pi_k$ 的求解
在讨论 $\pi_k$ 的详细求解之前，首先引入两个和 $Z_i$ 相关的两个后验概率，分别是 $\mathrm{P}(Z_i=k\mid X_i, \theta^{(t)})$ ，以及 $\mathrm{P}(Z_i, X_i = k \mid \theta)$ 分别表示 $Z_i = k$ 在 $X_i$ 观测数据以及 $\theta^{(t)}$ 参数下的后验概率以及 $Z_i, X_i$ 在 $\theta$ 参数下的联合后验概率。
首先先谈一下 $Z_i = k$ 后验概率，以 GMM 为例，代入 GMM 中有，
$$ \mathrm{P}(Z_i = k\mid X_i, \theta^{(t)}) = \frac{\pi_k\mathcal{N}(X_i\mid \mu_k, \sigma_k)}{\sum_{j=1}^{K} \pi_j\mathcal{N}(X_i\mid \mu_j, \sigma_j)} $$
这个表达式显而易见，描述全部混合系数下所有正太分布模型情况之中，使用 $\pi_k$ 作为当前正太分布混合系数的概率。其中比较值得注意的地方在于：
- 首先这个后验概率中，被采用的先验 $\theta^{(t)}$ 是为在 M-step 中求解后验 $\theta^{(t+1)}$ 而来；
- 这里 $k$ 指的某一个观察样本属于第 $k$ 个高斯分布，其混合系数为 $\pi_k$。
然后再分析一下关于 $Z_i$ 与 $X_i$ 的联合分布概率，同样代入 GMM，根据 bayes theorem，
$$ \begin{aligned} \mathrm{P}(X_i, z_i=k\mid \theta) & = \mathrm{P}(X_i\mid Z_i=k, \theta)\mathrm{P}(Z_i = k\mid \theta) \\ \\ & = \mathrm{P}(X_i\mid Z_i=k, \theta)\pi_k \end{aligned} $$
回顾 Q 函数:
$$ Q(\theta, \theta^{(t)}) = \sum_{Z}\mathrm{P}(Z\mid X, \theta_{(t)})\log{\mathrm{P}(X,Z\mid \theta)} $$
在此基础上，求解以下函数式：
$$ \frac{\partial Q(\theta\mid \theta^{(t)})}{\partial{\pi_k}} = \frac{\sum_{Z}\mathrm{P}(Z_i = k\mid X_i, \theta^{(t)})}{\pi_k}\tag{1} $$
但是因为 $Q$ 函数基于 Jenson 不等式的结论推导而来，所以需要规定权重 $\pi_k$ 一个约束，即
$$ \sum_{k=1}^{K}\pi_k = 1 \tag{2} $$
为了将这个月入代入 (1)，我们需要利用拉格朗日乘数法处理这个约束条件：
$$ \mathcal{L}(\theta, \lambda) = Q(\theta, \theta^{(t)}) + \lambda(1-\sum_{k=1}^{K}\pi_k) $$
对 $\mathcal{L}(\theta, \lambda)$ 求关于 $\pi_k$ 的偏导，有
$$ \frac{\partial \mathcal{L}}{\partial\pi_k}=\sum_{i=1}^{N}P(Z_i = k | X_i, \theta^{(t)})\cdot\frac{1}{\pi_k} - \lambda\tag{3} $$
这里需要注意在求算关于 $\pi_k$ 偏导时，$Z_i$ 后验概率项前是没有 $\sum_{k=1}^{K}$ 项的，因为已经确定了 $\pi_k$。
为研究
$$ \arg\max_{\pi_k, \lambda} \mathcal{L(\theta, \lambda)} $$
令 (3) 式为零，
$$ \pi_k^{(t+1)} = \frac{1}{\lambda}\sum_{i=1}^{N}\mathrm{P}(Z_i =k\mid X_i, \theta^{(t)})\tag{4} $$
这里需要带入 (2) 中这一限制，即在等式两端同时求和有
$$ 1 = \sum_{k=1}^{K}\pi_k^{(t+1)} = \sum_{k=1}^{K}\frac{1}{\lambda}\sum_{i=1}^{N}\mathrm{P}(Z_i =k\mid X_i, \theta^{(t)}) $$
可以求出
$$ \lambda = \sum_{k=1}^{K}\sum_{i=1}^{N}\mathrm{P}(Z_i = k\mid X_i, \theta^{(t)}) $$
又由
$$ \sum_{k=1}^{K}P(Z_i=k \mid X_i, \theta^{(t)}) = 1 $$
所以可以得到，
$$ \lambda =\sum_{i=1}^{N} 1 = N $$
带入 (4)，就可以求出更新后的 $\pi_k$
$$ \pi_k^{(t+1)} = \frac{1}{N}\sum_{i=1}^{N}\mathrm{P}(Z_i=k\mid X_i, \theta) $$
相似的其他迭代参数则更容易出，因为不需要考虑例如 （2）中的限制，
可通过对应的偏导求出：
$$ \mu_k^{(t+1)} = \frac{\sum_{i=1}^{N}\mathrm{P}(Z_i = k \mid X_i, \theta^{(t)})X_i}{\sum_{i=1}^{N}\mathrm{P}(Z_i=k\mid X_i, \theta^{(t)})} $$
$$ \sigma_k^{(t+1)}=\frac{\sum_{i=1}^N P\left(Z_i=k \mid X_i, \theta^{(t)}\right)\left(X_i-\mu_k^{(t+1)}\right)\left(X_i-\mu_k^{(t+1)}\right)^T}{\sum_{i=1}^N P\left(Z_i=k \mid X_i, \theta^{(t)}\right)} $$