1 混合模型与高斯混合模型
混合模型在当前科学研究以及实际应用中都有广泛使用，作为分析复杂现象的一个灵活有力的模型,随着数据量的不断增大，混合模型相比单一模型有着明显的优势，所以它被广泛应用于医学、人脸识别、模式识别、图像处理等多个领域。本文用 EM 算法应用到混合模型的参数估计问题中。
混合模型是指随机变量 $Y$ 的概率密度函数为如下形式：
$$ P(y|\theta) = \sum_{k=1}^K \alpha_k\phi(y|\theta_k),\ \alpha_k\ge0\text{且}\sum_{k=1}^K \alpha_k=1 $$
此处的 $\theta = (\alpha_1,\cdots, \alpha_K,\theta_1,\cdots,\theta_K)$ ，即混合模型有 $K$ 个分支组成，每个分支的权值为 $\alpha_k$ 。特别地，当每个分支服从的分布都是高斯分布时，则称混合分布为有 $K$ 个分支的高斯混合模型（Gaussian miixture model，GMM）。
若为 GMM 模型，则第 $k$ 个分支为：
$$ \phi(y|\theta_k) = \frac1{\sqrt{2\pi}\sigma_k}\exp{(-\frac{(y-\mu_k)^2}{2\sigma_k^2})},\ k=1,2,\cdots,K,\ \theta_k=(\mu_k,\sigma_k^2) $$
一般混合模型可以由任意概率密度组合，在此仅介绍最常用的高斯混合模型。
2 高斯混合模型参数估计的 EM 算法
假设观测数据 $y_1,y_2, \cdots, y_N$ 由高斯混合模型生产，
$$ P(y|\theta) = \sum_{k=1}^K \alpha_k\phi(y|\theta_k) $$
其中， $\theta = (\alpha_1,\cdots, \alpha_K,\theta_1,\cdots,\theta_K)$ 。我们用 EM 算法估计高斯混合模型的参数 $\theta$。
2.1 明确隐变量，写出完全数据的对数似然函数
我们的观测数据为 $y_1,y_2, \cdots, y_N$，这些数据的产生方法如下：
- 依据概率 $\alpha_k$ 选择第 $k$ 个高斯分布分模型 $\phi(y|\theta_k)$；
- 依据第 $k$ 个高斯分布分模型 $\phi(y|\theta_k)$ 生成观测数据 $y_j,j=1,2,\cdots,N$
因此，我们的观测数据 $y_1,y_2, \cdots, y_N$ 是已知的。但是我们不知道观测数据 $y_j$ 来自哪一个分模型，我们将 $y_j$ 来自第 $k$ 个分模型这一未能观测到的数据用隐变量 $\gamma_{jk}$ 来表示，其定义如下：
$$ \gamma_{jk} = \begin{cases} 1,\text{第 $j$ 个观测来自第 $k$ 个分模型}\\ 0,\text{其他} \end{cases} $$
其中，$j=1,2,\cdots,M;\ k=1,2,\cdots,K$ 。
有了观测数据 $y_j$ 及未观测数据 $\gamma_{jk}$ ，则完全数据为：
$$ (y_j, \gamma_{j1}, \gamma_{j2}, \cdots, \gamma_{jK}),\ j=1,2,\cdots,N $$
于是，可以写出完全数据的似然函数：
$$ \begin{aligned} P(y,\gamma|\theta) &= \prod_{j=1}^N P(y_j, \gamma_{j1}, \gamma_{j2}, \cdots, \gamma_{jK} | \theta)\\ &= \prod_{k=1}^K \prod_{j=1}^N [\alpha_k\phi(y|\theta_k)]^{\gamma_{jk}}\\ &= \prod_{k=1}^K \alpha_k^{n_k} \prod_{j=1}^N [\phi(y|\theta_k)]^{\gamma_{jk}}\\ &= \prod_{k=1}^K \alpha_k^{n_k} \prod_{j=1}^N \Big[\frac1{\sqrt{2\pi}\sigma_k}\exp{(-\frac{(y_j-\mu_k)^2}{2\sigma_k^2})}\Big]^{\gamma_{jk}} \end{aligned}其中，$n_k = \sum_{j=1}^N\gamma_{jk}$ ，er $$
其中，$n_k = \sum_{j=1}^N\gamma_{jk}$ ，而 $ \sum_{k=1}^Kn_k = N$ 。
那么，完全数据的对数似然函数为
$$ \log P(y,\gamma|\theta) = \sum_{k=1}^K\Bigg\{n_k\log \alpha_k + \sum_{j=1}^N\gamma_{jk}\Big[ \log (\frac1{\sqrt{2\pi}})-\log\sigma_k - \frac1{2\sigma_k^2}(y_j-\mu_k)^2\Big] \Bigg\} $$
2.2 EM 算法的 E 步：确定 Q 函数
$$ \begin{aligned} Q(\theta|\theta^{(i)}) &= E\{\log P(y,\gamma|\theta)|y,\theta^{(i)}\}\\ &= E\Bigg\{ \sum_{k=1}^K\Bigg\{n_k\log \alpha_k + \sum_{j=1}^N\gamma_{jk}\Big[ \log (\frac1{\sqrt{2\pi}})-\log\sigma_k - \frac1{2\sigma_k^2}(y_j-\mu_k)^2\Big] \Bigg\} \Bigg\} \\ &= \sum_{k=1}^K \Bigg\{ \sum_{j=1}^N (E\gamma_{jk})\log\alpha_k + \sum_{j=1}^N(E\gamma_{jk}) \Big[ \log (\frac1{\sqrt{2\pi}})-\log\sigma_k - \frac1{2\sigma_k^2}(y_j-\mu_k)^2\Big] \Bigg\} \end{aligned} $$
这里最后一步用到了 $n_k = \sum_{j=1}^N\gamma_{jk}$
上式需要计算 $E(\gamma_{jk}|y,\theta)$，记为 $\hat \gamma_{jk}$ 。
$$ \begin{aligned} \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)} \\ &= \frac{P(\gamma_{jk}=1,y_j|\theta)}{\sum_{k=1}^K P(\gamma_{jk}=1,y_j|\theta)} \\ &= \frac{P(y_j|\gamma_{jk}=1,\theta)P(\gamma_{jk}=1|\theta)}{\sum_{k=1}^KP(y_j|\gamma_{jk}=1,\theta)P(\gamma_{jk}=1|\theta)}\\ &= \frac{\alpha_k\phi(y_j|\theta_k)}{\sum_{k=1}^K\alpha_k\phi(y_j|\theta_k)},\ j=1,2,\cdots,N;\ k=1,2,\cdots,K \end{aligned} $$
由于 $\gamma_{jk}$ 是一个 0-1 随机变量，因此其期望就等于 $P(\gamma_{jk}=1|y,\theta)$，故第二个等式成立；
$\gamma_{jk}$ 只跟观测数据 $y_j$ 有关，但是观测数据 $y_j$ 跟 $\gamma_{j1},\gamma_{j2},\cdots,\gamma_{jK}$ 都有关，结合贝叶斯公式和边际分布可成立第三、四、五个等式。
$\hat \gamma_{jk}$ 是在当前模型参数下第 $j$ 个观测数据来自第 $k$ 个分类模型的概率，称为分模型 $k$ 对观测数据 $y_j$ 的响应度。
将 $\hat \gamma_{jk} = E\gamma_{jk}$ 及 $n_k = \sum_{j=1}^N\gamma_{jk}$ 代入 $Q$ 函数中，即得：
$$ Q(\theta|\theta^{(i)}) = \sum_{k=1}^K\Bigg\{ n_k\log\alpha_k+\sum_{j=1}^N\hat\gamma_{jk}\Big[ \log (\frac1{\sqrt{2\pi}})-\log\sigma_k - \frac1{2\sigma_k^2}(y_j-\mu_k)^2 \Big] \Bigg\} $$
2.3 确定 EM 算法的 M 步
迭代的 M 步是求函数 $Q(\theta,\theta^{(i)})$ 对 $\theta$ 的极大值，即求新一轮迭代的模型参数：
$$ \theta^{(i+1)} = arg\max_\theta Q(\theta,\theta^{(i)}) $$
用 $\hat\mu_k$，$\hat\sigma_k^2$ 及 $\hat\alpha_k$ $k=1,2,\cdots,K$ ，表示 $\theta^{(i+1)}$ 的各参数。求 $\hat\mu_k$，$\hat\sigma_k^2$ 只需将 2.2 节最后得到的 $Q$ 函数 分别对 $\hat\mu_k$，$\hat\sigma_k^2$ 求偏导数并令其为 0，即可得到；求 $\hat\alpha_k$ 是在 $\sum_{k=1}^N\alpha_k=1$ 条件下求偏导数并令其为 0 得到的。结果如下：
$$ \hat\mu_k = \frac{\sum_{j=1}^N\hat\gamma_{jk}y_j}{\sum_{j=1}^N\hat\gamma_{jk}},\ k=1,2,\cdots,K $$
$$ \hat\sigma_k^2 = \frac{\sum_{j=1}^N\hat\gamma_{jk}(y_j-\mu_k)^2}{\sum_{j=1}^N\hat\gamma_{jk}},\ k=1,2,\cdots,K $$
$$ \hat\alpha_k = \frac{n_k}{N} = \frac{\sum_{j=1}^N\hat\gamma_{jk}}{N},\ k=1,2,\cdots,K $$
重复以上计算，知道对数似然函数值不再有明显的变化为止。
2.4 高斯混合模型参数估计的 EM 算法步骤总结
输入：观测数据 $y_1,y_2,\cdots,y_N$，高斯混合模型；
输出：高斯混合模型参数。
(1) 取参数的初始值开始迭代；
(2) E 步：依据当前模型参数，计算分模型 $k$ 对观测数据 $y_j$ 的响应度
$$ \hat\gamma_{jk}= \frac{\alpha_k\phi(y_j|\theta_k)}{\sum_{k=1}^K\alpha_k\phi(y_j|\theta_k)},\ j=1,2,\cdots,N;\ k=1,2,\cdots,K $$
(3) M 步：计算新一轮迭代的模型参数
$$ \hat\mu_k = \frac{\sum_{j=1}^N\hat\gamma_{jk}y_j}{\sum_{j=1}^N\hat\gamma_{jk}},\ k=1,2,\cdots,K $$
$$ \hat\sigma_k^2 = \frac{\sum_{j=1}^N\hat\gamma_{jk}(y_j-\mu_k)^2}{\sum_{j=1}^N\hat\gamma_{jk}},\ k=1,2,\cdots,K $$
$$ \hat\alpha_k = \frac{n_k}{N} = \frac{\sum_{j=1}^N\hat\gamma_{jk}}{N},\ k=1,2,\cdots,K $$
(4) 重复第 2 和第 3 步，直至收敛。
