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

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

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

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

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

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

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

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

看回小球颜色的例子。

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

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

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

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

让 g g g 代表单高斯模型的概率密度函数，则混合高斯模型可以被写成： p ( x ) = ∑ k = 1 K w k g k ( x ∣ μ k , Σ k ) (1) p(\textbf{\textit{x}}) = \sum^{K}_{k=1} w_k g_k \left( \textbf{\textit{x}} | \boldsymbol{\mu}_k , \Sigma_k \right)\tag{1} p(x)=k=1∑K​wk​gk​(x∣μk​,Σk​)(1)上面表示 k k k 个不同参数的高斯模型的加权和，这些权重都是正的，其中：

g k g_k gk​ 表示均值为 μ k \boldsymbol{\mu}_k μk​，协方差为 Σ k \Sigma_k Σk​ 的单高斯模型。 w k w_k wk​ 表示混合系数(权值)。

权值的和为 1 1 1，保证了混合高斯模型的密度函数积分仍然是 1 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 K 本身也是一个参数，你需要确定它。

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

这里为了简单起见，我们只考虑给定的 K K K ，以及恒定权值 w w 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

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

我们不最大化联合概率，而是假设观测都是独立的，然后最大化这些概率的乘积： μ ^ , Σ ^ = arg ⁡ max ⁡ μ , Σ ∏ i = 1 N p ( x i ∣ μ , Σ ) \hat{\boldsymbol{\mu}},\hat{\Sigma} = \arg \max_{\boldsymbol{\mu},\Sigma} \prod^N_{i=1} p(\textbf{\textit{x}}_i |\boldsymbol{\mu},\Sigma ) μ^​,Σ^=argμ,Σmax​i=1∏N​p(xi​∣μ,Σ)

然后取对数，并且转为最大化对数概率的和： μ ^ , Σ ^ = arg ⁡ max ⁡ μ , Σ ∑ i = 1 N ln ⁡ p ( x i ∣ μ , Σ ) ‾ (2) \hat{\boldsymbol{\mu}},\hat{\Sigma} = \arg \max_{\boldsymbol{\mu},\Sigma} \sum^N_{i=1} \ln \underline{p(\textbf{\textit{x}}_i |\boldsymbol{\mu},\Sigma )}\tag{2} μ^​,Σ^=argμ,Σmax​i=1∑N​lnp(xi​∣μ,Σ)​(2)

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

把式 (3) 代入式 (2) 后，得到这个样子： μ ^ , Σ ^ = arg ⁡ max ⁡ μ , Σ ∑ i = 1 N ln ⁡ { 1 K ∑ k = 1 K g k ( x i ∣ μ k , Σ k ) } (4) \hat{\boldsymbol{\mu}},\hat{\Sigma} = \arg \max_{\boldsymbol{\mu},\Sigma} \sum^N_{i=1} \ln \left\{ \frac{1}{K} \sum^{K}_{k=1} g_k \left( \textbf{\textit{x}}_i | \boldsymbol{\mu}_k , \Sigma_k \right) \right\}\tag{4} μ^​,Σ^=argμ,Σmax​i=1∑N​ln{K1​k=1∑K​gk​(xi​∣μk​,Σk​)}(4)其中 g k ( x ) = 1 ( 2 π ) D / 2 ∣ Σ k ∣ 1 / 2 exp ⁡ { − 1 2 ( x − μ k ) T Σ k − 1 ( x − μ k ) } g_k(\textbf{\textit{x}})= \frac{1}{(2\pi)^{D/2} |\Sigma_k|^{1/2}} \exp\left\{ -\frac{1}{2}(\textbf{\textit{x}} - \boldsymbol{\mu}_k)^T \Sigma_k^{-1}(\textbf{\textit{x}} - \boldsymbol{\mu}_k)\right\} gk​(x)=(2π)D/2∣Σk​∣1/21​exp{−21​(x−μk​)TΣk−1​(x−μk​)}

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

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

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

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

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

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

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

现在来看看潜变量： z k i = g k ( x i ∣ μ k , Σ k ) ∑ k = 1 K ( x i ∣ μ k , Σ k ) z^{i}_k = \frac{g_k \left( \textbf{\textit{x}}_i | \boldsymbol{\mu}_k,\Sigma_k \right)} {\sum^K_{k=1} \left( \textbf{\textit{x}}_i | \boldsymbol{\mu}_k,\Sigma_k \right)} zki​=∑k=1K​(xi​∣μk​,Σk​)gk​(xi​∣μk​,Σk​)​

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

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

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

我们有数据点 x i \textbf{\textit{x}}_i xi​，还有两个高斯模型 g 1 g_1 g1​ 和 g 2 g_2 g2​。

如果： g 1 g_1 g1​ 产生 x i \textbf{\textit{x}}_i xi​ 的概率是 p 1 p_1 p1​； g 2 g_2 g2​ 产生 x i \textbf{\textit{x}}_i xi​ 的概率是 p 2 p_2 p2​；

则在数据点 x i \textbf{\textit{x}}_i xi​ 处，潜变量 z 1 i = p 1 p 1 + p 2 z^i_1 = \dfrac{p_1}{p_1+p_2} z1i​=p1​+p2​p1​​， z 2 i = p 2 p 1 + p 2 z^i_2 = \dfrac{p_2}{p_1+p_2} z2i​=p1​+p2​p2​​

这个例子中 p 1 p_1 p1​ 比 p 2 p_2 p2​ 大，所以 x i \textbf{\textit{x}}_i xi​ 更可能由 g 1 g_1 g1​ 产生。

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

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

首先，我们初始化参数 μ \boldsymbol{\mu} μ 和 Σ \Sigma Σ：

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

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

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

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

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

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

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

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

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

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

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

现在来比较两个值，首先是中点的函数值 f ( x 1 + x 2 2 ) \textcolor{e05d89}{f(\dfrac{x_1 + x_2}{2})} f(2x1​+x2​​)（红色的点），它是曲线上的一个点。

第二个是两点函数值的平均值 f ( x 1 ) + f ( x 2 ) 2 \textcolor{02a0da}{\dfrac{f(x_1) + f(x_2)}{2}} 2f(x1​)+f(x2​)​（蓝色的点），它落在 f ( x 1 ) f(x_1) f(x1​) 和 f ( x 2 ) f(x_2) f(x2​) 连线的直线上。

容易看到在这个区间里，直线上的点永远大于曲线上的点。于是有图中那条不等式： f ( x 1 + x 2 2 ) ≤ f ( x 1 ) + f ( x 2 ) 2 f(\dfrac{x_1 + x_2}{2}) \leq \dfrac{f(x_1) + f(x_2)}{2} f(2x1​+x2​​)≤2f(x1​)+f(x2​)​

我们把这个想法一般化，下面这个不等式在任意大于 0 0 0 的权值下都成立： f ( a 1 x 1 + a 2 x 2 ) ≤ a 1 f ( x 1 ) + a 2 f ( x 2 ) a 1 + a 2 = 1 a 1 ≥ 0 a 2 ≥ 0 f(a_1 x_1 + a_2 x_2) \leq a_1 f(x_1) + a_2 f(x_2) \\[1.2em] \begin{aligned} a_1+a_2 & = 1 \\ a_1 & \geq 0 \\ a_2 & \geq 0 \end{aligned} f(a1​x1​+a2​x2​)≤a1​f(x1​)+a2​f(x2​)a1​+a2​a1​a2​​=1≥0≥0​

这是 2 2 2 个点的情况，事实上对于多个点，权值都为正的情况下同样适用。 甚至是在多维的情况下，只要我们的函数是凸的： f ( ∑ a i x i ) ≤ ∑ a i f ( x i ) ∑ a i = 1 a i ≥ 0 f\left( \sum a_i x_i \right) \leq \sum a_i f(x_i)\\[1.2em] \begin{aligned} \sum a_i &= 1 \\ a_i & \geq 0 \end{aligned} f(∑ai​xi​)≤∑ai​f(xi​)∑ai​ai​​=1≥0​

如果函数是凹的，我们只需要翻转这个不等式。 我们将要处理的函数是对数函数，它总是凹的： ln ⁡ ( ∑ a i p i ) ≥ ∑ a i ln ⁡ p i ∑ a i = 1 a i ≥ 0 \textcolor{blue}{\ln \left( \sum a_i p_i \right)} \textcolor{red}{\geq} \textcolor{green}{\sum a_i \ln p_i} \\[1.2em] \begin{aligned} \sum a_i &= 1 \\ a_i & \geq 0 \end{aligned} ln(∑ai​pi​)≥∑ai​lnpi​∑ai​ai​​=1≥0​

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

利用Jensen不等式，我们能获得目标函数的下界。例如这里 ln ⁡ ( ∑ a i p i ) \textcolor{blue}{\ln \left( \sum a_i p_i \right)} ln(∑ai​pi​) 的下界就是 ∑ a i ln ⁡ p i \textcolor{green}{\sum a_i \ln p_i} ∑ai​lnpi​。

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

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

像之前一样，我们考虑概率的对数： ln ⁡ p ( X ∣ θ ) = ln ⁡ ∑ Z p ( X , Z ∣ θ ) \ln p(X|\theta) = \ln \sum_Z p(X,Z|\theta) lnp(X∣θ)=lnZ∑​p(X,Z∣θ)

利用Jensen不等式，我们能获得这样的下界： ln ⁡ ∑ Z p ( X , Z ∣ θ ) = ln ⁡ ∑ Z q ( Z ) p ( X , Z ∣ θ ) q ( Z ) ≥ ∑ Z q ( Z ) ln ⁡ p ( X , Z ∣ θ ) q ( Z ) \ln \sum_Z p(X,Z|\theta) = \textcolor{blue}{\ln \sum_Z q(Z) \frac{ p(X,Z|\theta)}{q(Z)}} \geq \textcolor{green}{\sum_Z q(Z)\ln \frac{ p(X,Z|\theta)}{q(Z)}} lnZ∑​p(X,Z∣θ)=lnZ∑​q(Z)q(Z)p(X,Z∣θ)​≥Z∑​q(Z)lnq(Z)p(X,Z∣θ)​          注： q ( Z ) q(Z) q(Z) 是 Z Z Z 上的有效概率分布。

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

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

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

首先有一个初始的猜测 θ 0 \theta_0 θ0​，并且通过利用猜测，得到的 q q q，进而得到了下界 G G G。

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

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

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

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

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

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

同样的，新的 G G G 也不那么像 F F F，但它在 θ 1 \theta_1 θ1​ 附近和 F F F 很接近。

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

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

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

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

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

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

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