高斯混合模型（Gaussian Mixture Model，GMM）和期望最大化（Expectation Maximization，EM）算法

本文是关于 coursera 上 《Robotics: Estimation and Learning》 课程的笔记。

前面讲了一维和多维高斯分布的相关知识。

但是在某些情况下，使用 单高斯模型(Gaussian single model, GSM） 会有一些局限。

在现实世界中我们需要学习的目标可能符合这样的分布 ：

如上图所示，当你用单高斯模型去拟合它时，得到这样的曲线。

显然它不能很好地表征目标。这样的目标有多种模式，或者缺乏对称性。

---

你将看到混合高斯模型的表现力则很好，好到可以建模任意的分布。
简单地说，混合高斯模型是一些高斯模型的和。

下面这个图中，那些彩色的线是十个随机的高斯模型曲线，黑线则是它们的和。
一个混合后的曲线可以有非常怪异的形状，这是无法用一个简单的方程描述出来的。
所以，如果我们选择了合适的高斯模型，我们就可以表征任意一个分布 。

---

看回小球颜色的例子。

右边的 $2$ 维图展示了红色和绿色通道的颜色分布，我们可以用一个 $2\text{D}$ 高斯模型去建模这个分布：

或者我们可以用两个高斯模型的混合去建模：

用混合的方式，看起来更好滴描述了这两个颜色的分布情况。

现在用我们最喜爱的工具 ---- 数学😂，来更加正规地描述混合高斯模型 。

高斯混合模型（Gaussian Mixture Model，GMM）

让 $g$ 代表单高斯模型的概率密度函数，则混合高斯模型可以被写成：
$$p(\text{x})=\sum _{k=1}^K{w}_k{g}_k\left(\text{x}|{\mathbf{\mu }}_k,{\Sigma }_k\right)\text{\hspace{1em}(1)}$$上面表示 $k$ 个不同参数的高斯模型的加权和，这些权重都是正的，其中：

$g_k$ 表示均值为 ${\mathbf{\mu }}_k$，协方差为 ${\Sigma }_k$ 的单高斯模型。
$w_k$ 表示混合系数(权值)。

权值的和为 $1$，保证了混合高斯模型的密度函数积分仍然是 $1$。

如果我们可以使用任意多的单模型，并且这些模型有任意小的方差，我们在理论上就可以表示任意形状的分布，这正是混合高斯模型强大的地方 。

虽然混合高斯模型很强，但是它的代价是需要很多参数： μ = { μ 1 , μ 2 , … , μ K } Σ = { Σ 1 , Σ 2 , … , Σ K } w = { w 1 , w 2 , … , w K } K : Number of Components \boldsymbol{\mu} = \left\{ \boldsymbol{\mu}_1, \boldsymbol{\mu}_2, \dots , \boldsymbol{\mu}_K \right\} \\ \Sigma = \left\{ \Sigma_1, \Sigma_2, \dots, \Sigma_K \right\} \\ \textbf{\textit{w}} = \left\{ w_1, w_2, \dots, w_K \right\} \\ K: \text{Number of Components} μ={μ1​,μ2​,…,μK​}Σ={Σ1​,Σ2​,…,ΣK​}w={w1​,w2​,…,wK​}K:Number of Components

在同样的维数下，混合高斯模型比单高斯模型有更多的参数。
当混合的模型数量增加，均值和协方差矩阵的数量也随之增加 。
而且还多了一个新的参数：权值。
还有高斯模型的数量 $K$ 本身也是一个参数，你需要确定它。

这么多的参数带来一些副作用：
首先，参数估计变得更困难，后面会看到，我们没有混合高斯模型的解析解。
第二，犯错的概率变大了，特别是我们可能遇到过拟合的问题。谨记，我们需要对模型的复杂程度特别小心。

这里为了简单起见，我们只考虑给定的 $K$ ，以及恒定权值 $w$ 的情况：
μ = { μ 1 , μ 2 , … , μ K } Σ = { Σ 1 , Σ 2 , … , Σ K } w = w = 1 / K K : 指定的数字 \boldsymbol{\mu} = \left\{ \boldsymbol{\mu}_1, \boldsymbol{\mu}_2, \dots , \boldsymbol{\mu}_K \right\} \\ \Sigma = \left\{ \Sigma_1, \Sigma_2, \dots, \Sigma_K \right\} \\ \textbf{\textit{w}} = w = 1/K \\ K: \text{指定的数字} μ={μ1​,μ2​,…,μK​}Σ={Σ1​,Σ2​,…,ΣK​}w=w=1/KK:指定的数字在此基础上探讨估计模型的均值和协方差矩阵。

将参数量简化后，我们来尝试做似然估计，求解混合高斯模型的参数。

---

简单回顾一下，和前面一样，我们的似然函数为：
μ ^ , Σ ^ = arg ⁡ max ⁡ μ , Σ p ( { x i } ∣ μ , Σ ) μ = { μ k } Σ = { Σ k } k = 1 , 2 , … , K \hat{\boldsymbol{\mu}},\hat{\Sigma} = \arg \max_{\boldsymbol{\mu},\Sigma} p(\{\textbf{\textit{x}}_i\} |\boldsymbol{\mu},\Sigma ) \\[1em] \boldsymbol{\mu} = \left\{ \boldsymbol{\mu}_k\right\} \\ \Sigma = \left\{ \Sigma_k \right\} \\ k = 1,2,\dots,K μ^​,Σ^=argμ,Σmax​p({xi​}∣μ,Σ)μ={μk​}Σ={Σk​}k=1,2,…,K

最大似然估计意味着，我们希望找到模型的一组参数，使得产生观测数据的可能性最大。

我们不最大化联合概率，而是假设观测都是独立的，然后最大化这些概率的乘积：
$$\hat{\mathbf{\mu }},\hat{\Sigma }=\arg \underset{\mathbf{\mu },\Sigma }{\max }\prod _{i=1}^{N}p({\text{x}}_i|\mathbf{\mu },\Sigma )$$

然后取对数，并且转为最大化对数概率的和：$$\hat{\mathbf{\mu }},\hat{\Sigma }=\arg \underset{\mathbf{\mu },\Sigma }{\max }\sum _{i=1}^N\ln \underline{p({\text{x}}_i|\mathbf{\mu },\Sigma )}\text{\hspace{1em}(2)}$$

到目前为止，结果看起来和单高斯模型一样。
但是混合高斯模型的概率密度函数是多个单高斯分布的加权和。
对于式 (2) 下划线的部分，混合高斯模型的情况为：$$\underline{p(\text{x})}=\sum _{k=1}^K{w}_k{g}_k\left(\text{x}|{\mathbf{\mu }}_k,{\Sigma }_k\right)\text{\hspace{1em}(3)}$$ $K$ 表示单高斯模型的数量；
$g_k$ 是第 $k$ 个单高斯模型的概率密度函数，其均值为 ${\mathbf{\mu }}_k$，协方差为 ${\Sigma }_k$。
$w_k$ 是第 $k$ 个单高斯模型的权值。（这里 $w_k$ 都是 $1/K$）

把式 (3) 代入式 (2) 后，得到这个样子：$$\hat{\mathbf{\mu }},\hat{\Sigma }=\arg \underset{\mathbf{\mu },\Sigma }{\max }\sum _{i=1}^N\ln \left\{\frac{1}{K}\sum _{k=1}^K{g}_k\left({\text{x}}_i|{\mathbf{\mu }}_k,{\Sigma }_k\right)\right\}\text{\hspace{1em}(4)}$$其中$$g_k(\text{x})=\frac{1}{(2\pi {)}^{D/2}|{\Sigma }_k{|}^{1/2}}\exp \left\{-\frac{1}{2}(\text{x}-{\mathbf{\mu }}_k{)}^T{\Sigma }_k^{-1}(\text{x}-{\mathbf{\mu }}_k)\right\}$$

到这里再也无法解析地化简这个方程了，因为对数方程里出现了高斯模型的和。这意味着我们只能通过迭代计算，估计参数。并且得到的结果可能不是全局最优的。

Don’t worry，下面来简要地讲讲EM算法，这种算法能在特定条件下给出合理的解。

期望最大化算法（Expectation-Maximization，EM）

这里先用一个特例来看一下如何计算混合高斯模型的参数，然后扩展到 EM 的一般情况。

首先来看看用 EM 算法计算 GMM 的参数，我们需要两个新的东西：

$\bullet$ 初始的 $\mathbf{\mu }$ 和 $\Sigma$。（可以通过猜测来设定）
$\bullet$ 潜变量 $z$（Latent variable）。

利用一个猜测值去求解我们要估计的参数，听起来可能有点奇怪。

在这类称为非凸的优化问题里，有许多次优解，称为局部最小值，而初始的猜测将影响我们找到的解。猜测初始值也是一门技术活，本教程暂时不研究这个。

现在来看看潜变量：
$$z_k^i=\frac{g_k\left({\text{x}}_i|{\mathbf{\mu }}_k,{\Sigma }_k\right)}{{\sum }_{k=1}^K\left({\text{x}}_i|{\mathbf{\mu }}_k,{\Sigma }_k\right)}$$

$z_k^i$ 是关于第 $k$ 个高斯模型的第 $i$ 个数据点的潜变量。它定义为第 $k$ 个高斯模型在该数据点处的相对比，本质上表征了第 $i$ 个数据点由第 $k$ 个高斯模型生成的概率。

（每个数据点上都有 $k$ 个潜变量，所以共有 $k\times i$ 个潜变量）

来看一个 $K$ 等于 $2$ 时的 $1$ 维高斯模型的例子：

我们有数据点 ${\text{x}}_i$，还有两个高斯模型 $g_1$ 和 $g_2$。

如果：
$g_1$ 产生 ${\text{x}}_i$ 的概率是 $p_1$；
$g_2$ 产生 ${\text{x}}_i$ 的概率是 $p_2$；

则在数据点 ${\text{x}}_i$ 处，潜变量 $z_1^i=\frac{p_1}{p_1+p_2}$，$z_2^i=\frac{p_2}{p_1+p_2}$

这个例子中 $p_1$ 比 $p_2$ 大，所以 ${\text{x}}_i$ 更可能由 $g_1$ 产生。

现在给定所有数据点关于所有模型的这个潜变量，我们重新定义带权重的均值和协方差矩阵。以 $z$ 为权值，如果一个数据点属于第 $k$ 个高斯模型的概率很小，则这个点对计算第 $k$ 个高斯模型的参数的贡献就很小。
$$\begin{array}{rl}{\hat{\mathbf{\mu }}}_k& =\frac{1}{z_k}\sum _{i=1}^N{z}_k^i{\text{x}}_i\\ {\hat{\Sigma }}_k& =\frac{1}{z_k}\sum _{i=1}^N{z}_k^i({\text{x}}_i-\hat{\mathbf{\mu }})({\text{x}}_i-\hat{\mathbf{\mu }}{)}^T\\ z_k& =\sum _{i=1}^N{z}_k^i\end{array}$$
这个求解的形式和多维高斯分布很像，只是多了一个权重项 $z$，可以回顾一下前面作为比较。

综合之前讨论的东西，我们可以迭代计算模型参数和潜变量。

首先，我们初始化参数 $\mathbf{\mu }$ 和 $\Sigma$：

然后用初始的 $\mathbf{\mu }$ 和 $\Sigma$ 来计算潜变量 $z$：

当 $z$ 更新后，我们就固定 $z$，再用它来更新 $\mathbf{\mu }$ 和 $\Sigma$：

往复地进行这两个步骤，直到参数的变化量变得非常小，得到估计值 ${\hat{\mathbf{\mu }}}_k$ 和 ${\hat{\Sigma }}_k$：

我们将不会在这里证明这个算法，但是迭代最终会收敛到一个局部最优值上。

至此已经介绍了如何利用迭代的 EM 算法来估计 GMM 的参数，通过一个初始设定的均值和协方差矩阵，及引入潜变量，我们神奇地将问题变成了在局部最优意义下可解的问题。

下面介绍更一般化的 EM 算法。

---

回顾上面公式 (4)，这个极大似然估计不好求解，求和里面又有对数，对数里面又有求和。

如果我们能找到似然函数的下边界，并且使下边界最大化，这就相当于将似然函数最大化。

现在要做的是，把 EM 算法看作是对目标函数下界的最大化过程。我们需要找到一个表达式，来描述目标函数的下界。

为了找到这个表达式，这里用到了 Jensen不等式 。它是这个算法的基础，EM算法的主要思想可以用Jensen不等式解释。

Jensen不等式

考虑图所示的一个凸函数，在定义域中取两个点 $x_1$ 和 $x_2$，以及他们的中点 $\frac{x_1+x_2}{2}$。

现在来比较两个值，首先是中点的函数值 $f(\frac{x_1+x_2}{2})$（红色的点），它是曲线上的一个点。

第二个是两点函数值的平均值 $\frac{f(x_1)+f(x_2)}{2}$（蓝色的点），它落在 $f(x_1)$ 和 $f(x_2)$ 连线的直线上。

容易看到在这个区间里，直线上的点永远大于曲线上的点。于是有图中那条不等式：$$f(\frac{x_1+x_2}{2})\le \frac{f(x_1)+f(x_2)}{2}$$

我们把这个想法一般化，下面这个不等式在任意大于 $0$ 的权值下都成立：$$\begin{gathered} f(a_1{x}_1+a_2{x}_2)\le a_{1}f(x_1)+a_{2}f(x_2) \\ \begin{array}{rl}{a}_1+a_2& =1\\ a_1& \ge 0\\ a_2& \ge 0\end{array} \end{gathered}$$

这是 $2$ 个点的情况，事实上对于多个点，权值都为正的情况下同样适用。
甚至是在多维的情况下，只要我们的函数是凸的：$$\begin{gathered} f\left(\sum a_i{x}_i\right)\le \sum a_{i}f(x_i) \\ \begin{array}{rl}\sum a_i& =1\\ a_i& \ge 0\end{array} \end{gathered}$$

如果函数是凹的，我们只需要翻转这个不等式。
我们将要处理的函数是对数函数，它总是凹的：$$\begin{gathered} \ln \left(\sum a_i{p}_i\right)\ge \sum a_i\ln p_i \\ \begin{array}{rl}\sum a_i& =1\\ a_i& \ge 0\end{array} \end{gathered}$$

（这里要正确区分凹凸函数的概念，请查询相关资料。为了方便解释仅从形状上描述凹凸性）

利用Jensen不等式，我们能获得目标函数的下界。例如这里 $\ln \left(\sum a_i{p}_i\right)$ 的下界就是 $\sum a_i\ln p_i$。

---

我们引入我们永远也无法确切知道的隐变量，隐变量的含义在一些特定的应用中是清楚的，例如前面提到的混合高斯模型的参数估计。

我们将对隐变量取边缘概率，根据边缘概率（marginal probability）的定义：$$p(X|\theta )=\sum _{Z}p(X,Z|\theta )$$          （可以理解为 $X$ 是样本集合， $\theta$ 是所有参数集合？）

像之前一样，我们考虑概率的对数：$$\ln p(X|\theta )=\ln \sum _{Z}p(X,Z|\theta )$$

利用Jensen不等式，我们能获得这样的下界：$$\ln \sum _{Z}p(X,Z|\theta )=\ln \sum _{Z}q(Z)\frac{p(X,Z|\theta )}{q(Z)}\ge \sum _{Z}q(Z)\ln \frac{p(X,Z|\theta )}{q(Z)}$$          注：$q(Z)$ 是 $Z$ 上的有效概率分布。

我们无法直接最大化概率和的对数，但是我们可以利用函数的下界，对于它我们是可以计算的。

获得下界意味着找到了隐变量 $z$ 的一个分布 $q(Z)$，这个下界又恰好是一个关于旧的参数猜测的函数。

现在来看一下怎样具体实施这套算法：

首先有一个初始的猜测 ${\theta }_0$，并且通过利用猜测，得到的 $q$，进而得到了下界 $G$。

如我说过的，对数形式的概率可以被看作是参数的函数。

假设 $F$ 和 $G$ 是如图所示的 $\theta$ 的函数，$F$ 和 $G$ 并不完全一致 ，但是他们在 ${\theta }_0$ 附近有一些局部的相似性。

记得我们的任务是，希望求一个参数 $\theta$ 的最大似然估计。

给定当前的下界 $G$，我们可以找到一个更好的参数，记作 ${\theta }^{\ast }$：

这仍然不是 $F$ 的最大值，但它比之前的估计更优。

由于我们有了一个更好的参数的猜测，我们可以计算改进后的下界 $G$：

同样的，新的 $G$ 也不那么像 $F$，但它在 ${\theta }_1$ 附近和 $F$ 很接近。

我们又可以根据改进后的 $G$ 再来计算一个更好的 $\theta$。

如果我们重复这个过程，参数最终会收敛到一个局部最优值：

总结一下，在 E 步骤，我们通过给定的参数计算对数概率的下界：

再次说明，这意味着我们获得了在给定参数 ($\hat{\theta }$) 上的隐变量的分布 ($q(Z)$)

在 M 步骤，我们在当前的下界上最大化参数：

重复这个过程，直到产生的变化小于我们的阈值。

我们已经学习了如何利用 EM 算法寻找复杂问题的最大似然解 。
前面课程中提到的混合高斯模型的参数估计，是应用这个算法的一个例子。
在后面探讨定位（localization）时你将看到EM算法的另一个应用。