机器学习算法基础：硬核高斯分布

提前放狠话：看不懂打si我吧，真的很细节了。
提前了解：
1、极大似然估计（MLE）、最大后验概率估计（MAP）、贝叶斯估计区别
2、正定矩阵与半正定矩阵定义性质与理解
3、马氏距离和欧式距离详解
4、参考Github笔记

高斯分布

正态分布（Normal distribution），也称“常态分布”，又名高斯分布（Gaussian distribution）。

下面将分一维和多维一一进行介绍。

拓展：
希腊字母表的第十八字母（$\Sigma$，$\sigma$）。$\Sigma$ 英语名称：Sigma 汉语名称：西格玛（大写$\Sigma$，小写$\sigma$），第十八个希腊字母。
大写$\Sigma$用于数学上的总和符号，比如：$\Sigma P_i$，其中$i=1,2,...,N$，即为求$P_1+P_2+...+P_N$的和。小写$\sigma$用于统计学上的标准差。

MLE参数估计

数据集$X$符合某分布，且是独立同分布的（$iid$），在 MLE 方法中，参数估计的方法如下：
$${\theta }_{MLE}=\underset{\theta }{\mathrm{argmax}}\log \prod _{i=1}^{N}p(x_i|\theta )=\underset{\theta }{\mathrm{argmax}}\sum _{i=1}^N\log p(x_i|\theta )$$

一维情况

数据符合一维高斯分布（独立同分布）：
$$X\sim \mathcal{N}(\mu ,{\sigma }^2)$$
一维高斯分布的概率密度函数PDF写为：
$$p(x|\mu ,{\sigma }^2)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{(x-\mu {)}^2}{2{\sigma }^2}}$$

我们要使用MLE进行估计的参数为：$\theta =(\mu ,{\sigma }^2)$

带入 MLE 中我们考虑一维的情况
$$\begin{array}{rl}\log p(X|\theta )=\sum _{i=1}^N\log p(x_i|\theta )& =\sum _{i=1}^N\log \frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{(x_i-\mu {)}^2}{2{\sigma }^2}}\\ & \\ & =\sum _{i=1}^N\log \frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{(x_i-\mu {)}^2}{2{\sigma }^2}}\\ & \\ & =\sum _{i=1}^N\log \frac{1}{\sqrt{2\pi }\sigma }+\sum _{i=1}^N\log e^{-\frac{(x_i-\mu {)}^2}{2{\sigma }^2}}\\ & \\ & =\sum _{i=1}^N[\log \frac{1}{\sqrt{2\pi }}+\log \frac{1}{\sigma }-\frac{(x_i-\mu {)}^2}{2{\sigma }^2}]\end{array}$$
首先对 $\mu$ 的极值可以得到 ：
$${\mu }_{MLE}=\underset{\mu }{\mathrm{argmax}}\log p(X|\theta )=\underset{\mu }{\mathrm{argmax}}-\sum _{i=1}^N(x_i-\mu {)}^2=\underset{\mu }{\mathrm{argmin}}\sum _{i=1}^N(x_i-\mu {)}^2$$
于是：
$$\frac{\partial }{\partial \mu }\sum _{i=1}^N(x_i-\mu {)}^2=0⟶{\mu }_{MLE}=\frac{1}{N}\sum _{i=1}^N{x}_i$$
其次对 $\theta$ 中的另一个参数 $\sigma$ ，有：
$$\begin{array}{rl}{\sigma }_{MLE}=\underset{\sigma }{\mathrm{argmax}}\log p(X|\theta )& =\underset{\sigma }{\mathrm{argmax}}\sum _{i=1}^N[-\log \sigma -\frac{(x_i-\mu {)}^2}{2{\sigma }^2}]\\ & =\underset{\sigma }{\mathrm{argmin}}\sum _{i=1}^N[\log \sigma +\frac{(x_i-\mu {)}^2}{2{\sigma }^2}]\end{array}$$
于是：
$$\frac{\partial }{\partial \sigma }\sum _{i=1}^N[\log \sigma +\frac{1}{2{\sigma }^2}(x_i-\mu {)}^2]=0⟶{\sigma }_{MLE}^2=\frac{1}{N}\sum _{i=1}^N(x_i-\mu {)}^2$$

有偏VS无偏

值得注意的是，上面的推导中，首先对 $\mu$ 求 MLE， 然后利用这个结果求 ${\sigma }_{MLE}$ ，因此可以预期的是对 ${\mu }_{MLE}$求期望时， 可以发现${\mu }_{MLE}$是无偏差的：
$${\mathbb{E}}_{\mathcal{D}}[{\mu }_{MLE}]={\mathbb{E}}_{\mathcal{D}}[\frac{1}{N}\sum _{i=1}^N{x}_i]=\frac{1}{N}\sum _{i=1}^N{\mathbb{E}}_{\mathcal{D}}[x_i]=\mu$$
但是当对 ${\sigma }_{MLE}$ 求 期望的时候由于使用了单个数据集的 ${\mu }_{MLE}$，因为${\mu }_{MLE}$不一定与所有数据的$\mu$相等，因此对${\sigma }_{MLE}$ 求期望的时候我们会发现 ${\sigma }_{MLE}$ 是有偏差的：

$$\begin{array}{rl}{\mathbb{E}}_{\mathcal{D}}[{\sigma }_{MLE}^2]& ={\mathbb{E}}_{\mathcal{D}}[\frac{1}{N}\sum _{i=1}^N(x_i-{\mu }_{MLE}{)}^2]={\mathbb{E}}_{\mathcal{D}}[\frac{1}{N}\sum _{i=1}^N(x_i^2-2x_i{\mu }_{MLE}+{\mu }_{MLE}^2)\\ & \\ & ={\mathbb{E}}_{\mathcal{D}}[\frac{1}{N}\sum _{i=1}^N{x}_i^2-{\mu }_{MLE}^2]={\mathbb{E}}_{\mathcal{D}}[\frac{1}{N}\sum _{i=1}^N{x}_i^2-{\mu }^2+{\mu }^2-{\mu }_{MLE}^2]\\ & \\ & ={\mathbb{E}}_{\mathcal{D}}[\frac{1}{N}\sum _{i=1}^N{x}_i^2-{\mu }^2]-{\mathbb{E}}_{\mathcal{D}}[{\mu }_{MLE}^2-{\mu }^2]={\sigma }^2-({\mathbb{E}}_{\mathcal{D}}[{\mu }_{MLE}^2]-{\mu }^2)\\ & \\ & ={\sigma }^2-({\mathbb{E}}_{\mathcal{D}}[{\mu }_{MLE}^2]-{\mathbb{E}}_{\mathcal{D}}^2[{\mu }_{MLE}])={\sigma }^2-Var[{\mu }_{MLE}]\\ & \\ & ={\sigma }^2-Var[\frac{1}{N}\sum _{i=1}^N{x}_i]={\sigma }^2-\frac{1}{N^2}\sum _{i=1}^{N}Var[x_i]=\frac{N-1}{N}{\sigma }^2\end{array}$$

所以，无偏估计量${\hat{\sigma }}^2$如下式所示：
$${\hat{\sigma }}^2=\frac{1}{N-1}\sum _{i=1}^N(x_i-\mu {)}^2$$

多维情况

多维高斯分布的概率密度函数PDF写为：
$$p(x|\mu ,\Sigma )=\frac{1}{(2\pi {)}^{p/2}|\Sigma {|}^{1/2}}{e}^{-\frac{1}{2}(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )}$$
其中 $x,\mu \in {\mathbb{R}}^p,\Sigma \in {\mathbb{R}}^{p\times p}$ ，$\Sigma$ 为协方差矩阵，一般也是半正定矩阵且是对称矩阵。这里我们只考虑正定矩阵。
补充:
1、设$A$是$n$阶方阵，如果对任何非零向量$x$，都有$x^{T}Ax>0$，其中$x^T$表示$x$的转置，就称$A$为正定矩阵。
2、设$A$是实对称矩阵，如果对任意的实非零列向量$x$有$x^{T}Ax\ge 0$，就称$A$为半正定矩阵。
3、马氏距离公式：$\sqrt{(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )}$。我们可以发现如果${\Sigma }^{-1}$是单位矩阵的时候，马氏距离简化为欧氏距离。马氏距离的这一项${\Sigma }^{-1}$起到了方差归一化的作用，可以拿一维的高斯分布举例，$-\frac{1}{2}\frac{(x-\mu {)}^2}{{\sigma }^2}$，${\sigma }^2$等价于$\Sigma$。

首先我们处理指数上的数字，指数上的数字可以记为 $x$ 和 $\mu$ 之间的马氏距离。对于对称的协方差矩阵可进行特征值分解，$\Sigma =U\Lambda U^T=(u_1,u_2,\cdots ,u_p)diag({\lambda }_i)(u_1,u_2,\cdots ,u_p{)}^T=\sum _{i=1}^p{u}_i{\lambda }_i{u}_i^T$ ，其中$u_i\in {\mathbb{R}}^p$，$U,diag({\lambda }_i)\in {\mathbb{R}}^{p\times p}$，$U{U}^T=U^{T}U=I$，于是：

$${\Sigma }^{-1}=\sum _{i=1}^p{u}_i\frac{1}{{\lambda }_i}{u}_i^T$$

$$\Delta =(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )=\sum _{i=1}^p(x-\mu {)}^T{u}_i\frac{1}{{\lambda }_i}{u}_i^T(x-\mu )=\sum _{i=1}^p\frac{y_i^2}{{\lambda }_i}$$
其中$y_i=(x-\mu {)}^T{u}_i=y_i^T=u_i^T(x-\mu )$，$y_i$是一个实数。

我们注意到 $y_i$ 是 $x-\mu$ 在特征向量 $u_i$ 上的投影长度，
因此上式子在二维情况下，$\sum _{i=1}^2\frac{y_i^2}{{\lambda }_i}=\frac{y_1^2}{{\lambda }_1}+\frac{y_2^2}{{\lambda }_2}=\Delta$就是 $\Delta$ 取不同值时的同心椭圆。

我来换一个通俗易懂的讲法：
使用二维空间更直观

• X = { x 1 , x 2 , . . . , x N } X=\{x_1,x_2,...,x_N\} X={x1​,x2​,...,xN​}，其中$x_i=(x_i^1;x_i^2),i=1,2,...,N$，这些数据集$X$构成了一个二维空间。
• 使用$X$分别在两个维度上确定均值，${\mu }^1$和${\mu }^2$，即可得到$\mu =({\mu }^1;{\mu }^2)$。
• 如果使用的是欧式距离，$u_i$ 其实一一对应坐标轴$x_i$所指方向，因此$(x-\mu {)}^T{u}_i$两个向量相乘，$y_i$大小等于$(x-\mu {)}^T{u}_i$，且${\lambda }_1={\lambda }_2$，在二维空间上即为一个圆。
• 在上式为马氏距离，$u_i$ 不对应坐标轴$x_i$所指方向，因此$y_i$大小等于$x-\mu$ 在特征向量 $u_i$ 上的投影长度，在原二维空间上即为一个斜的椭圆。
• 因此马氏距离涵盖了欧式距离，它使用了各变量各自的方差去进行计算，起到了方差归一化的作用。
• 我们也可以直观地看到高斯二维分布其实为一个椭圆型。
  

化简了指数部分，我们重新看一下PDF：
$$p(x|\mu ,\Sigma )=\frac{1}{(2\pi {)}^{p/2}|\Sigma {|}^{1/2}}{e}^{-\frac{1}{2}\Delta }$$
当我们确定了$p(x|\mu ,\Sigma )$的值，那么我们也确定了$\Delta$的值，就可以得到满足该条件的特征向量形成的空间，其实在二维空间，构成了一座山的等高线，椭圆形的线。

下面我们看多维高斯模型在实际应用时的两个局限性

1. 参数过多。参数 $\Sigma ,\mu$ 的自由度为 $O(p^2)$ 对于维度很高的数据其自由度太高。高自由度的来源是 $\Sigma$ 有 $\frac{p(p+1)}{2}$ 个自由参数。
   解决方法：
   可以假设其是对角矩阵，这样就不需要进行特征值分解，$u_i$的形式就满足欧式距离出现的第一种情况，然后直观展示为是正的椭圆而不是斜的椭圆。
   甚至在各向同性假设中假设其对角线上的元素都相同。$u_i$的形式就满足欧式距离出现的所有情况，然后直观展示为一个正圆形。
   前一种的算法有 Factor Analysis，后一种有概率 PCA(p-PCA) 。

2. 第二个问题是单个高斯分布是单峰的，对有多个峰的数据分布不能得到好的模型结果。
   解决方法：
   高斯混合GMM 模型。

边缘概率以及条件概率

我们记 $x=(x_1,x_2,\cdots ,x_p{)}^T=(x_{a,m\times 1},x_{b,n\times 1}),\mu =({\mu }_{a,m\times 1},{\mu }_{b,n\times 1}),\Sigma =\left(\begin{array}{cc}{\Sigma }_{aa}& {\Sigma }_{ab}\\ {\Sigma }_{ba}& {\Sigma }_{bb}\end{array}\right)$，已知 $x\sim \mathcal{N}(\mu ,\Sigma )$。

首先是一个高斯分布的定理：

定理：已知 $x\sim \mathcal{N}(\mu ,\Sigma ),y\sim Ax+b$，那么 $y\sim \mathcal{N}(A\mu +b,A\Sigma A^T)$。

证明：$\mathbb{E}[y]=\mathbb{E}[Ax+b]=A\mathbb{E}[x]+b=A\mu +b$，$Var[y]=Var[Ax+b]=Var[Ax]=A\cdot Var[x]\cdot A^T$。

下面利用这个定理得到 $p(x_a),p(x_b),p(x_a|x_b),p(x_b|x_a)$ 这四个量。

1. $x_a=\left(\begin{array}{cc}{\mathbb{I}}_{m\times m}& {\mathbb{O}}_{m\times n}\end{array}\right)\left(\begin{array}{c}{x}_a\\ x_b\end{array}\right)$，代入上面的定理中可以得到：
   $$\begin{gathered} \mathbb{E}[x_a]=\left(\begin{array}{cc}\mathbb{I}& \mathbb{O}\end{array}\right)\left(\begin{array}{c}{\mu }_a\\ {\mu }_b\end{array}\right)={\mu }_a \\ Var[x_a]=\left(\begin{array}{cc}\mathbb{I}& \mathbb{O}\end{array}\right)\left(\begin{array}{cc}{\Sigma }_{aa}& {\Sigma }_{ab}\\ {\Sigma }_{ba}& {\Sigma }_{bb}\end{array}\right)\left(\begin{array}{c}\mathbb{I}\\ \mathbb{O}\end{array}\right)={\Sigma }_{aa} \end{gathered}$$
   所以 $x_a\sim \mathcal{N}({\mu }_a,{\Sigma }_{aa})$。

2. 同样的，$x_b\sim \mathcal{N}({\mu }_b,{\Sigma }_{bb})$。

3. 对于两个条件概率，我们引入三个量（第一个式子是一个构造性证明，记住就行）：
   $$\begin{gathered} x_{b\cdot a}=x_b-{\Sigma }_{ba}{\Sigma }_{aa}^{-1}{x}_a \\ {\mu }_{b\cdot a}={\mu }_b-{\Sigma }_{ba}{\Sigma }_{aa}^{-1}{\mu }_a \\ {\Sigma }_{bb\cdot a}={\Sigma }_{bb}-{\Sigma }_{ba}{\Sigma }_{aa}^{-1}{\Sigma }_{ab} \end{gathered}$$
   特别的，最后一个式子叫做 ${\Sigma }_{aa}$ 的 Schur Complementary。可以看到：
   $$x_{b\cdot a}=\left(\begin{array}{cc}-{\Sigma }_{ba}{\Sigma }_{aa}^{-1}& {\mathbb{I}}_{n\times n}\end{array}\right)\left(\begin{array}{c}{x}_a\\ x_b\end{array}\right)$$
   我们重点使用第一条式子推导第二、三条式子，所以：
   $$\begin{gathered} \mathbb{E}[x_{b\cdot a}]=\left(\begin{array}{cc}-{\Sigma }_{ba}{\Sigma }_{aa}^{-1}& {\mathbb{I}}_{n\times n}\end{array}\right)\left(\begin{array}{c}{\mu }_a\\ {\mu }_b\end{array}\right)={\mu }_{b\cdot a} \\ Var[x_{b\cdot a}]=\left(\begin{array}{cc}-{\Sigma }_{ba}{\Sigma }_{aa}^{-1}& {\mathbb{I}}_{n\times n}\end{array}\right)\left(\begin{array}{cc}{\Sigma }_{aa}& {\Sigma }_{ab}\\ {\Sigma }_{ba}& {\Sigma }_{bb}\end{array}\right)\left(\begin{array}{c}-{\Sigma }_{aa}^{-1}{\Sigma }_{ba}^T\\ {\mathbb{I}}_{n\times n}\end{array}\right)={\Sigma }_{bb\cdot a} \end{gathered}$$
   利用这三个量可以得到 $x_b=x_{b\cdot a}+{\Sigma }_{ba}{\Sigma }_{aa}^{-1}{x}_a$。

   我们先推导一下下面的公式
   若$x\sim \mathcal{N}(\mu ,{\Sigma }_{})$，则$Mx\perp Nx\iff M\Sigma N^T=0$
   如果两个随机变量均服从高斯分布，那么“不相关”等价于“独立”。

   证明：$$\begin{array}{rl}& \because x\sim \mathcal{N}(\mu ,{\Sigma }_{})\\ & \therefore Mx\sim \mathcal{N}(M\mu ,M{\Sigma }_{}{M}^T),Nx\sim \mathcal{N}(N\mu ,N{\Sigma }_{}{N}^T)\\ & \therefore \mathbb{C}ov(Mx,Nx)\\ & =\mathbb{E}[(Mx-M\mu )(Nx-N\mu {)}^T]\\ & =\mathbb{E}[M(x-\mu )(x-\mu {)}^{T}N]\\ & =M\mathbb{E}[(x-\mu )(x-\mu {)}^T]N\\ & =M\Sigma N^T\\ & \because Mx\perp N且均为高斯\\ & \therefore \mathbb{C}ov(Mx,Nx)=M\Sigma N^T=0\end{array}$$

   我们回到原来的式子：
   $$x_b=x_{b\cdot a}+{\Sigma }_{ba}{\Sigma }_{aa}^{-1}{x}_a$$
   我们辛苦构造的$x_{b\cdot a}$有什么作用呢，其实它和$x_a$是相互独立的，可以使用上面的公式以及前面的构造进行推导。
   $$\begin{gathered} \because x_{b\cdot a}\perp x_a\Rightarrow x_{b\cdot a}|x_a=x_{b\cdot a} \\ \therefore x_b|x_a=x_{b\cdot a}|x_a+{\Sigma }_{ba}{\Sigma }_{aa}^{-1}{x}_a|x_a=x_{b\cdot a}+{\Sigma }_{ba}{\Sigma }_{aa}^{-1}{x}_a \end{gathered}$$
   所以
   $$\mathbb{E}[x_b|x_a]={\mu }_{b\cdot a}+{\Sigma }_{ba}{\Sigma }_{aa}^{-1}{x}_a$$

   $$Var[x_b|x_a]={\Sigma }_{bb\cdot a}$$

   这里同样用到了定理。

4. 同样：
   $$\begin{gathered} x_{a\cdot b}=x_a-{\Sigma }_{ab}{\Sigma }_{bb}^{-1}{x}_b \\ {\mu }_{a\cdot b}={\mu }_a-{\Sigma }_{ab}{\Sigma }_{bb}^{-1}{\mu }_b \\ {\Sigma }_{aa\cdot b}={\Sigma }_{aa}-{\Sigma }_{ab}{\Sigma }_{bb}^{-1}{\Sigma }_{ba} \end{gathered}$$
   所以：
   $$\mathbb{E}[x_a|x_b]={\mu }_{a\cdot b}+{\Sigma }_{ab}{\Sigma }_{bb}^{-1}{x}_b$$

   $$Var[x_a|x_b]={\Sigma }_{aa\cdot b}$$

求$p(y),p(x|y)$

线性高斯模型常见求解

已知：$p(x)=\mathcal{N}(\mu ,{\Lambda }^{-1}),p(y|x)=\mathcal{N}(Ax+b,L^{-1})$，求解：$p(y),p(x|y)$。

解：令 $y=Ax+b+ϵ,ϵ\sim \mathcal{N}(0,L^{-1})$，所以 $\mathbb{E}[y]=\mathbb{E}[Ax+b+ϵ]=A\mu +b$，$Var[y]=A{\Lambda }^{-1}{A}^T+L^{-1}$，因此：
$$p(y)=\mathcal{N}(A\mu +b,L^{-1}+A{\Lambda }^{-1}{A}^T)$$
引入 $z=\left(\begin{array}{c}x\\ y\end{array}\right)$，我们可以得到 $Cov[x,y]=\mathbb{E}[(x-\mathbb{E}[x])(y-\mathbb{E}[y]{)}^T]$。对于这个协方差可以直接计算：
$$\begin{array}{rl}Cov(x,y)& =\mathbb{E}[(x-\mu )(Ax-A\mu +ϵ{)}^T]\\ & =\mathbb{E}[(x-\mu )(x-\mu {)}^T{A}^T]\\ & =Var[x]A^T\\ & ={\Lambda }^{-1}{A}^T\end{array}$$
注意到协方差矩阵的对称性，
所以 $p(z)=\mathcal{N}(\left(\begin{array}{c}\mu \\ A\mu +b\end{array}\right),\left(\begin{array}{cc}{\Lambda }^{-1}& {\Lambda }^{-1}{A}^T\\ A{\Lambda }^{-1}& L^{-1}+A{\Lambda }^{-1}{A}^T\end{array}\right))$。
根据之前的公式，我们可以得到：
$$\mathbb{E}[x|y]=\mu +{\Lambda }^{-1}{A}^T(L^{-1}+A{\Lambda }^{-1}{A}^T{)}^{-1}(y-A\mu -b)$$

$$Var[x|y]={\Lambda }^{-1}-{\Lambda }^{-1}{A}^T(L^{-1}+A{\Lambda }^{-1}{A}^T{)}^{-1}A{\Lambda }^{-1}$$