高斯分布的极大似然估计

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

前面通过一个例子简单地介绍了极大似然估计的意思，现在来对高斯分布做极大似然估计。

一维高斯分布

概率密度函数

一维高斯分布（Gaussian Distribution）的概率密度函数如下：
$$p(x)=\frac{1}{\sqrt{2\pi }\sigma }\exp \left\{-\frac{(x-\mu {)}^2}{2{\sigma }^2}\right\}$$

高斯分布非常有用，而且非常重要：

$\bullet$ 描述高斯分布只需要 2 个参数，均值 $\mu$ 和 方差 ${\sigma }^2$，它们就是该分布的本质信息，同也容易计算和解释。

$\bullet$ 高斯分布具有一些很好的数学性质，比如多个高斯分布的乘积可以形成另一个高斯分布，所以你不需要担心遇到其它形式的分布。

$\bullet$ 中心极限定理告诉我们，任何随机变量的样本均值的期望，都收敛于高斯分布，这说明高斯分布是一个为噪声和不确定性建模的合适选择。

下面举一个例子，图片处理的问题，看看在目标颜色的建模中如何运用高斯分布。

小球颜色

下面的图片里有一个黄色小球。人眼一下就能看出，但是电脑怎么知道它的颜色？

统计它每个像素的色相值，得到下面的直方图：

纵轴可以认为是像素点的数量，横轴是像素点的色相值。
大部分像素点的色相值在 $50\sim 60$ ，这部分刚好是黄色所在色相区间。（黑色部分先忽略）
按照经验可以认为它是黄色小球。

如果这个图片的像素是 $1920\times 1080$ 呢？
有两百多万像素，统计这样的直方图需要巨大内存，非常消耗资源。

---

极大似然估计

为了减少工作量，我们从中抽取部分像素点，记它们的色调值为 $(x_1,x_2,x_3,\dots ,x_i)$。
通过这一部分样本，来估算这张图片所有像素点的色相的分布，推算小球的颜色。

我们剧透一下结果：
由中心极限定理或者上面的直方图，可以认为该样本满足高斯分布。
求出概率密度函数中的 $\mu$ 和 $\sigma$，即可得到像素点色调值的分布，大概会像这样：

根据 $\mu$ 可以推断小球颜色。

我们把抽取的像素点成为观测数据，如何通过观测数据，求出其概率模型的参数？
方法就是用极大似然估计（maximum likelihood estimation）。
上一篇文章通过一个例子简单介绍过它的意思。

现在我们运用高斯分布模型，它的参数是 $\mu$ 和 $\sigma$。

对于观测数据 $(x_1,x_2,x_3,\dots ,x_i)$，我们关心的是能将似然函数最大化的参数，$\hat{\mu }$ 和 $\hat{\sigma }$。

用数学的方式来表达： μ ^ , σ ^ = arg ⁡ max ⁡ μ , σ p ( { x i } ∣ μ , σ ) \hat{\mu},\hat{\sigma} = \arg\max_{\mu,\sigma} p(\{x_i\}|\mu,\sigma) μ^​,σ^=argμ,σmax​p({xi​}∣μ,σ)

上标 $\hat{}$ 表示估计值。

我们需要最大化的似然函数是所有样本数据的联合概率，假设每一个观测之间都是相互独立的，联合似然函数就可以被简单地写成：关于每个样本的似然函数的乘积 。于是： p ( { x i } ∣ μ , σ ) = ∏ i = 1 N p ( x i ∣ μ , σ ) p\left(\{x_i\}|\mu,\sigma \right)= \prod^{N}_{i=1}p(x_i|\mu,\sigma) p({xi​}∣μ,σ)=i=1∏N​p(xi​∣μ,σ)

所以似然函数是：$$\hat{\mu },\hat{\sigma }=\arg \underset{\mu ,\sigma }{\max }\prod _{i=1}^{N}p(x_i|\mu ,\sigma )$$

---

要求解这个似然函数，我们用到一些对数函数的特性。先画一个对数函数的样子：

它是关于 $x$ 单调递增的，原函数取对数后，单调性不变。

利用这个特性，我们改变一下求解目标：最大化对数似然函数。

$$\arg \underset{\mu ,\sigma }{\max }\prod _{i=1}^{N}p(x_i|\mu ,\sigma )=\arg \underset{\mu ,\sigma }{\max }\ln \left\{\prod _{i=1}^{N}p(x_i|\mu ,\sigma )\right\}\text{\hspace{1em}(1)}$$
虽然改变了目标函数，但是使目标函数达到最大值的参数值是一样的，最大化对数似然函数往往更简单。

对数函数的还一个特点是：$$\log (x_1{x}_2\dots x_k)=\log (x_1)+\log (x_2)+\cdots +\log (x_k)$$

所以将式子 (1) 展开得到：$$\begin{array}{rl}\arg \underset{\mu ,\sigma }{\max }\prod _{i=1}^{N}p(x_i|\mu ,\sigma )& =\arg \underset{\mu ,\sigma }{\max }\ln \left\{\prod _{i=1}^{N}p(x_i|\mu ,\sigma )\right\}\\ & =\arg \underset{\mu ,\sigma }{\max }\sum _{i=1}^N\ln p(x_i|\mu ,\sigma )\end{array}$$

将高斯分布的概率密度函数代进去：$$\begin{array}{rl}p(x_i|\mu ,\sigma )& =\frac{1}{\sqrt{2\pi }\sigma }\exp \left\{-\frac{(x_i-\mu {)}^2}{2{\sigma }^2}\right\}\end{array}$$

所以：$$\begin{array}{rl}\ln p(x_i|\mu ,\sigma )& =\ln \frac{1}{\sqrt{2\pi }\sigma }\exp \left\{-\frac{(x_i-\mu {)}^2}{2{\sigma }^2}\right\}\\ & =\left\{-\frac{(x_i-\mu {)}^2}{2{\sigma }^2}-\ln \sigma -\ln \sqrt{2\pi }\right\}\end{array}$$

常数项 $\ln \sqrt{2\pi }$ 可以去掉，因为它不影响参数的估计。

再把负号去掉，把 $\max$ 变成 $\min$，从而把这个方程变成最小化问题。
两个方程是等价的，而且最小化形式是优化问题的标准形式。

所以我们的求解问题变成了这样：$$\hat{\mu },\hat{\sigma }=\arg \underset{\mu ,\sigma }{\min }\sum _{i=1}^N\left\{\frac{(x_i-\mu {)}^2}{2{\sigma }^2}+\ln \sigma \right\}$$

把上面 橙色 部分整体用 $J(\mu ,\sigma )$ 表示，这是对最小化问题代价函数的一个常用记号 。

利用凸优化问题的判据，令 $J$ 的一阶偏导为 $0$ 可以算出使似然函数最大的 $\mu$ 和 $\sigma$：
$$\frac{\partial J}{\partial \mu }=0⟶\hat{\mu }$$$$\frac{\partial J(\hat{\mu },\sigma )}{\partial \sigma }=0⟶\hat{\sigma }$$

最后得到求解公式是：$$\begin{array}{rl}\hat{\mu }& =\frac{1}{N}\sum _{i=1}^N{x}_i\\ {\hat{\sigma }}^2& =\frac{1}{N}\sum _{i=1}^N(x_i-\hat{\mu }{)}^2\end{array}$$

代入所有观测数据 $x_i$，得到 $\mu =52.3$ ， $\sigma =1.5$。

---

推导

偏导的过程是这样求得：$$\begin{gathered} \begin{array}{rl}\frac{\partial J}{\partial \mu }& =\frac{\partial }{\partial \mu }\sum _{i=1}^N\left\{\frac{(x_i-\mu {)}^2}{2{\sigma }^2}+\ln \sigma \right\}\\ & =\sum _{i=1}^N\left\{\frac{\partial }{\partial \mu }\frac{(x_i-\mu {)}^2}{2{\sigma }^2}\right\}\\ & =\frac{1}{{\sigma }^2}\sum _{i=1}^N(x_i-\mu )=0\end{array} \\ \Downarrow \\ \hat{\mu }=\frac{1}{N}\sum _{i=1}^N{x}_i \end{gathered}$$

$$\begin{gathered} \begin{array}{rl}\frac{\partial J}{\partial \sigma }& =\frac{\partial }{\partial \sigma }\sum _{i=1}^N\left\{\frac{(x_i-\mu {)}^2}{2{\sigma }^2}+\ln \sigma \right\}\\ & =\left(\frac{\partial }{\partial \sigma }\frac{1}{2{\sigma }^2}\right)\left(\sum _{i=1}^N(x_i-\hat{\mu }{)}^2\right)+\frac{N}{\sigma }\\ & =\frac{1}{\sigma }\left(N-\frac{1}{{\sigma }^2}\sum _{i=1}^N(x_i-\hat{\mu }{)}^2\right)=0\end{array} \\ \Downarrow \\ {\hat{\sigma }}^2=\frac{1}{N}\sum _{i=1}^N(x_i-\hat{\mu }{)}^2 \end{gathered}$$

---

多维高斯分布

前面的小球颜色，是通过HSV颜色空间的色相（H）来确定的，而RGB是一种更常用的色彩模式：

在3D图中画出所有像素的RGB值，就能得到如下的色彩分布：

现在考虑用RGB三通道模型对红球的色彩建模，看看这里如何运用高斯分布。

概率密度函数

多维高斯分布的表达式如下：$$p(x\text{x}x)=\frac{1}{(2\pi {)}^{D/2}|\Sigma {|}^{1/2}}\exp \left\{-\frac{1}{2}(x\text{x}x-\mu \text{μ}\mu {)}^T{\Sigma }^{-1}(x\text{x}x-\mu \text{μ}\mu )\right\}$$其中：
$D$ 表示维度数。
$x\text{x}x$ 是变量，注意它是加粗的，是个向量。
$\mu \text{μ}\mu$ 是均值，注意它是加粗的，是个向量。
$\Sigma$ 是协方差矩阵。(这个符号读做 Sigma)

和一维的情况相比，这里的变量 $x\text{x}x$ 是一个向量，均值 $\mu \text{μ}\mu$ 也是对应的向量形式，协方差 $\Sigma$ 现在是一个方阵。

在协方差矩阵中，有两个重要部分：对角项和非对角项。
下面是一个 $2$ 维协方差矩阵的例子：$$\Sigma =\left[\begin{array}{cc}{\sigma }_{x_1}^2& {\sigma }_{x_1}{\sigma }_{x_2}\\ {\sigma }_{x_2}{\sigma }_{x_1}& {\sigma }_{x_2}^2\end{array}\right]$$
对角项是两个变量 $x_1$ 和 $x_2$ 各自的方差；非对角项表示两个变量的相关性，即一个变量和另一个变量的相关程度。分母中 $\Sigma$ 两边的竖线代表 $\Sigma$ 矩阵的行列式。

下面是一个蓝色小球的例子：

这是对于处理3维变量RGB，变量向量包含我们的样本像素，即每个像素的红色，绿色，蓝色的值。
均值 $\mu \text{μ}\mu$ 是一个 $3\times 1$ 的向量，协方差矩阵 $\Sigma$ 是 $3\times 3$ 的矩阵；
$p(x\text{x}x)$ 是这个像素由小球产生的概率：在假设已知小球RGB模型的均值和方差的条件下，该样本像素从小球中抽出的概率。

下面是一个二维情况的例子：
$$p(x)=\frac{1}{2\pi }\exp \left\{-\frac{x^2+y^2}{2}\right\}$$其中： $D=2$ ，$x\text{x}x={\left[\begin{array}{cc}x& y\end{array}\right]}^T$， $\mu \text{μ}\mu ={\left[\begin{array}{cc}0& 0\end{array}\right]}^T$，$\Sigma =\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$

这是在给定参数下，简化的概率密度函数。它的图看起来像这样：

从图上看，$0$ 均值的二维高斯分布，像个只有一个峰的小山。如果你从中间切开这个曲面，它正是一个 $1$ 维的高斯分布。有时在 $2$ 维空间而非 $3$ 维空间画出分布更有利于理解。等高线上的所有 $x$ 具有相同的概率值，由于例子的协方差矩阵是个单位阵，等高线都是圆形，最中心的点表示了山峰，越外围的圆环代表概率越低。

当只有均值发生变化时，分布会发生平移：

当方差变大时，分布变胖，山峰的值变小；当方差变小，山峰的 $p(x)$ 变大，同时分布变瘦：

但是多变量高斯分布的协方差矩阵有一些特性，在 $1$ 维高斯分布中是没有的。
正如上面提到，$\Sigma$ 协方差矩阵的非对角元素包含相关项。
如果 $\Sigma$ 包含非零的对角元，则高斯分布的形状会被歪斜，这个现象在单变量情形下是不会发生的：

关于 $\Sigma$ 还有一些性质。

首先，协方差矩阵必须是对称正定的，这意味着 $\Sigma$ 的元素是关于对角线对称的，而且矩阵的特征值必须都是正的。

第二，就算 $\Sigma$ 具有非零的相关项，我们也能找到一种坐标变换，让分布的形状变成对称的。可以使用特征值分解算法将协方差矩阵分解，进而得到这些变换。（$\Sigma$ 可以被分解为 ${\text{UDU}}^T$ 的形式，其中 $\text{D}$ 是一个对角阵）

极大似然估计

前面已经介绍过一维高斯分布的情况，多维的情况也是一样，下面直接给出表达式： μ ^ , Σ ^ = arg ⁡ max ⁡ μ , Σ p ( { x i } ∣ μ , Σ ) \hat{\boldsymbol{\mu}},\hat{\Sigma} = \arg \max_{\boldsymbol{\mu},\Sigma}p(\{\textbf{\textit{x}}_i\} |\boldsymbol{\mu},\Sigma ) μ^​,Σ^=argμ,Σmax​p({xi​}∣μ,Σ)

再回顾一下似然的定义，似然是一个给定的包含未知参数的模型，产生这个观测数据的概率。

这里的模型是多维高斯分布模型，参数是均值向量 $\mathbf{\mu }$ 和协方差矩阵 $\Sigma$。

我们感兴趣的是获得均值和协方差矩阵，来使得这个似然在给定的观测数据集上最大，这是我们目标的数学解释。

似然函数是所有数据概率的联合，通常是难以处理的。当我们假设各观测是独立的，总概率函数可以被写成每个独立似然函数的乘积，在这个假设条件下，求解问题变成这样： μ ^ , Σ ^ = 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​}∣μ,Σ)

然后取对数，变成求和的问题：$$\hat{\mathbf{\mu }},\hat{\Sigma }=\arg \underset{\mathbf{\mu },\Sigma }{\max }\sum _{i=1}^N\ln p({\text{x}}_i|\mathbf{\mu },\Sigma )$$

再把概率密度函数代进去：$$\hat{\mathbf{\mu }},\hat{\Sigma }=\arg \underset{\mathbf{\mu },\Sigma }{\max }\sum _{i=1}^N\left\{-\frac{1}{2}({\text{x}}_i-\mathbf{\mu }{)}^T{\Sigma }^{-1}({\text{x}}_i-\mathbf{\mu })-\frac{1}{2}\ln |\Sigma |-\frac{D}{2}\ln (2\pi )\right\}$$

常数项不影响求解，忽略它，然后再取负数，变成最小化问题：$$\hat{\mathbf{\mu }},\hat{\Sigma }=\arg \underset{\mathbf{\mu },\Sigma }{\min }\sum _{i=1}^N\left\{\frac{1}{2}({\text{x}}_i-\mathbf{\mu }{)}^T{\Sigma }^{-1}({\text{x}}_i-\mathbf{\mu })+\frac{1}{2}\ln |\Sigma |\right\}$$

最后解得最小化代价函数 $J$ 的最优解，这就得到了对均值和协方差矩阵的极大似然估计：$$J(\mathbf{\mu },\Sigma )=\sum _{i=1}^N\left\{\frac{1}{2}({\text{x}}_i-\mathbf{\mu }{)}^T{\Sigma }^{-1}({\text{x}}_i-\mathbf{\mu })+\frac{1}{2}\ln |\Sigma |\right\}$$$$\frac{\partial J}{\partial \mathbf{\mu }}=0⟶\hat{\mathbf{\mu }}$$$$\frac{\partial J(\hat{\mathbf{\mu }},\Sigma )}{\partial \Sigma }=0⟶\hat{\Sigma }$$

最后得到求解公式是：$$\begin{array}{rl}\hat{\mathbf{\mu }}& =\frac{1}{N}\sum _{i=1}^N{\text{x}}_i\\ \hat{\Sigma }& =\frac{1}{N}\sum _{i=1}^N({\text{x}}_i-\hat{\mathbf{\mu }})({\text{x}}_i-\hat{\mathbf{\mu }}{)}^T\end{array}$$

我们得到的最终解是向量形式的样本均值，和样本的协方差矩阵。

回到小球颜色的例子，图中展现了小球颜色在蓝色和红色两个维度上的分布。

使用观测数据和我们获得的求解方程，我们可以计算参数的极大似然估计。
从图中的等高线上可以看出，在我们的模型中红色和蓝色通道是负相关的。
目前为止我们的例子是完美对称且只有一个峰，然而，一些其它奇怪的分布是可能存在的。
后面会讲如何利用多维高斯分布，去创建一个混合模型，用它可以表示各种不同的分布。

---

推导

这个多维高斯分布的求解(偏导)是这样来的：

$\mathbf{\mu }$ 的估计：

$$\frac{\partial J}{\partial \mathbf{\mu }}=\frac{\partial }{\partial \mathbf{\mu }}\sum _{i=1}^N\left\{\frac{1}{2}({\text{x}}_i-\mathbf{\mu }{)}^T{\Sigma }^{-1}({\text{x}}_i-\mathbf{\mu })+\frac{1}{2}\ln |\Sigma |\right\}\text{\hspace{1em}(2)}$$

首先：$${\text{x}}^T{\Sigma }^{-1}\mathbf{\mu }=\left[\begin{array}{cccc}{x}_1& x_2& \cdots & x_n\end{array}\right]{\Sigma }^{-1}\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\\ ⋮\\ {\mu }_n\end{array}\right]=\left[\begin{array}{cccc}{\mu }_1& {\mu }_2& \cdots & {\mu }_n\end{array}\right]{\Sigma }^{-1}\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]={\mathbf{\mu }}^T{\Sigma }^{-1}\text{x}$$

然后：$$\begin{array}{rl}(\text{x}-\mathbf{\mu }{)}^T{\Sigma }^{-1}(\text{x}-\mathbf{\mu })& =({\text{x}}^T-{\mathbf{\mu }}^T){\Sigma }^{-1}(\text{x}-\mathbf{\mu })\\ & =({\text{x}}^T-{\mathbf{\mu }}^T){\Sigma }^{-1}\text{x}-({\text{x}}^T-{\mathbf{\mu }}^T){\Sigma }^{-1}\mathbf{\mu }\\ & ={\text{x}}^T{\Sigma }^{-1}\text{x}-{\mathbf{\mu }}^T{\Sigma }^{-1}\text{x}-{\text{x}}^T{\Sigma }^{-1}\mathbf{\mu }+{\mathbf{\mu }}^T{\Sigma }^{-1}\mathbf{\mu }\\ & ={\text{x}}^T{\Sigma }^{-1}\text{x}-2{\mathbf{\mu }}^T{\Sigma }^{-1}\text{x}+{\mathbf{\mu }}^T{\Sigma }^{-1}\mathbf{\mu }\end{array}$$

由于是对 $\mathbf{\mu }$ 求偏导，忽略掉常数和其它变量，因此公式 (2) 变为：
$$\frac{\partial J}{\partial \mathbf{\mu }}=\frac{\partial }{\partial \mathbf{\mu }}\sum _{i=1}^N\left\{\frac{1}{2}{\mathbf{\mu }}^T{\Sigma }^{-1}\mathbf{\mu }-{\mathbf{\mu }}^T{\Sigma }^{-1}{\text{x}}_i\right\}\text{\hspace{1em}(3)}$$

公式 (3) 有两项关于指数的求导。

首先，对于 ${\mathbf{\mu }}^T{\Sigma }^{-1}\text{x}$，记 ${\Sigma }^{-1}\text{x}$ 为 $\text{A}$，则 ${\mathbf{\mu }}^T{\Sigma }^{-1}\text{x}={\mathbf{\mu }}^T\text{A}$。

由求导公式 $\frac{\text{d}(\text{A}\text{X}{)}^T}{\text{d}\text{X}}={\text{A}}^T$ 得：

$$\frac{\partial ({\mathbf{\mu }}^T{\Sigma }^{-1}\text{x})}{\partial \mathbf{\mu }}=\frac{\partial ({\mathbf{\mu }}^T\text{A})}{\partial \mathbf{\mu }}=\frac{\partial ({\text{A}}^T\mathbf{\mu }{)}^T}{\partial \mathbf{\mu }}=\text{A}={\Sigma }^{-1}\text{x}\text{\hspace{1em}(4)}$$

然后，对于${\mathbf{\mu }}^T{\Sigma }^{-1}\mathbf{\mu }$，根据求导公式 $\frac{\text{d}({\text{X}}^T\text{A}\text{X})}{\text{d}\text{X}}=\text{A}\text{X}+{\text{A}}^T\text{X}$ 得：

$$\frac{\partial ({\mathbf{\mu }}^T{\Sigma }^{-1}\mathbf{\mu })}{\partial \mathbf{\mu }}={\Sigma }^{-1}\mathbf{\mu }+{{\Sigma }^{-1}}^T\mathbf{\mu }\text{\hspace{1em}(5)}$$

由于 ${\Sigma }^{-1}$ 是一个对称矩阵，所以 ${\Sigma }^{-1}={{\Sigma }^{-1}}^T$，把 (4) 和 (5) 代入 (3) 得：

∂ J ∂ μ = Σ − 1 ∑ i = 1 N { μ − x i } \frac{\partial J}{\partial \boldsymbol{\mu}} = \Sigma^{-1} \sum^N_{i=1} \left\{ \boldsymbol{\mu} - \textbf{\textit{x}}_i \right\} ∂μ∂J​=Σ−1i=1∑N​{μ−xi​}

$$\frac{\partial J}{\partial \mathbf{\mu }}=0⟶\hat{\mathbf{\mu }}=\frac{1}{N}\sum _{i=1}^N{\text{x}}_i$$

$\Sigma$ 的估计：

有点麻烦，推理可以看这里：https://www.cnblogs.com/bigmonkey/p/11379144.html

用的时候直接记答案好了