高维高斯分布基础

多位高斯分布的几何理解

多维高斯分布表达式为：
$p(x∣μ,Σ)=1(2π)p/2∣Σ∣1/2e−12(x−μ)TΣ−1(x−μ)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,μ∈Rp,Σ∈Rp×px,\mu\in\mathbb{R}^{p},\Sigma\in\mathbb{R}^{p\times p}$ ， $Σ\Sigma$ 为协方差矩阵，一般而言是半正定矩阵，这里我们强化一下条件，只考虑正定矩阵。

首先我们处理指数上的数字，指数上的数字可以记为 $x$ 和 $μ\mu$ 之间的马氏距离。

马氏距离(Mahalanobis Distance)是度量学习中一种常用的距离指标，同欧氏距离、曼哈顿距离、汉明距离等一样被用作评定数据之间的相似度指标。但却可以应对高维线性分布的数据中各维度间非独立同分布的问题。

两个向量 $x\bf{x}$ 和 $y\mathbf{y}{}$ 之间的马氏距离为：
$DM(x,y)=(x−y)TΣ−1(x−y))D_M(\bf{x},\bf{y})=\sqrt{(\bf{x}-\bf{y})^T\Sigma^{-1}(\bf{x}-\bf{y}))}$
其中 $Σ\Sigma$ 是多维随机变量的协方差矩阵， $μ\mu$ 为样本均值，如果协方差矩阵是单位向量（ $Σ=I\Sigma=I$ ），也就是各维度独立同分布，马氏距离就变成了欧氏距离。

关于马氏距离，详见：马氏距离(Mahalanobis Distance)。

对于对称的协方差矩阵可进行特征值分解， $Σ=UΛUT=(u1,u2,⋯,up)diag(λi)(u1,u2,⋯,up)T=∑i=1puiλiuiT\Sigma=U\Lambda U^{T}=(u_{1},u_{2},\cdots,u_{p})diag(\lambda_{i})(u_{1},u_{2},\cdots,u_{p})^{T}=\sum\limits _{i=1}^{p}u_{i}\lambda_{i}u_{i}^{T}$ ，于是：
$Σ−1=∑i=1pui1λiuiT\Sigma^{-1}=\sum\limits _{i=1}^{p}u_{i}\frac{1}{\lambda_{i}}u_{i}^{T}$

$Δ=(x−μ)TΣ−1(x−μ)=∑i=1p(x−μ)Tui1λiuiT(x−μ)=∑i=1pyi2λi\Delta=(x-\mu)^{T}\Sigma^{-1}(x-\mu)=\sum\limits _{i=1}^{p}(x-\mu)^{T}u_{i}\frac{1}{\lambda_{i}}u_{i}^{T}(x-\mu)=\sum\limits _{i=1}^{p}\frac{y_{i}^{2}}{\lambda_{i}}$

我们注意到 $y_{i}$ 是 $x−μx-\mu$ 在特征向量 $u_{i}$ 上的投影长度，因此上式子就是 $Δ\Delta$ 取不同值时的同心椭圆。例如，在维度 $P = 2$ 时，取 $Δ=1\Delta=1$ ，则有： $y12λ1+y22λ2=1\frac{y_1^2}{\lambda_1}+\frac{y_2^2}{\lambda_2}=1$ ，明显就是椭圆的曲线方程。

多维高斯模型的限制

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

参数 $Σ,μ\Sigma,\mu$ 的自由度为 $O(p^{2})$ 对于维度很高的数据其自由度太高。
- 解决方案：高自由度的来源是 $Σ\Sigma$ 有 $p(p+1)2\frac{p(p+1)}{2}$ 个自由参数，可以假设其是对角矩阵，甚至假设其对角线上的元素都相同，此时称为各向同性的高斯分布。前一种的算法有 Factor Analysis，后一种有概率 PCA (p-PCA) 。
第二个问题是单个高斯分布是单峰的，对有多个峰的数据分布不能得到好的结果。
- 解决方案：使用多个单高斯模型组合得到高斯混合模型 GMM。

多维高斯分布的边缘概率和条件概率

对于高斯模型的线性变换，有这样一个定理（暂不证明）：

定理：已知 $x∼N(μ,Σ),y∼Ax+bx\sim\mathcal{N}(\mu,\Sigma), y\sim Ax+b$ ， $x∈Rp,y∈Rpx\in\mathbb{R}^p,y\in\mathbb{R}^p$ ，那么 $y∼N(Aμ+b,AΣAT),Σ∈Rp×p,A∈R1×py\sim\mathcal{N}(A\mu+b, A\Sigma A^T),\ \ \Sigma \in \mathbb{R}^{p\times p },\ \ A\in\mathbb{R}^{1\times p}$ 。

我们将 $p$ 维样本数据拆分为 $m + n$ 维： $x=(x1,x2,⋯,xp)T=(xa,m×1,xb,n×1)Tx=(x_1, x_2,\cdots,x_p)^T=(x_{a,m\times 1}, x_{b,n\times1})^T$ 。

对应的高斯模型的参数也进行拆分：均值 $μ=(μa,m×1,μb,n×1)T\mu=(\mu_{a,m\times1}, \mu_{b,n\times1})^T$ ，协方差矩阵 $Σ=(ΣaaΣabΣbaΣbb)\Sigma=\begin{pmatrix}\Sigma_{aa}&\Sigma_{ab}\\\Sigma_{ba}&\Sigma_{bb}\end{pmatrix}$ ，已知 $x∼N(μ,Σ)x\sim\mathcal{N}(\mu,\Sigma)$ 。

下面介绍如何求多维高斯分布的边缘概率和条件概率 $p(x_a),p(x_b),p(x_a|x_b),p(x_b|x_a)$ 。

求边缘概率

构造 $xa=(Im×mOm×n)(xaxb)x_a=\begin{pmatrix}{I}_{m\times m}&{O}_{m\times n}\end{pmatrix}\begin{pmatrix}x_a\\x_b\end{pmatrix}$ ，其中 $I, O$ 分别是单位矩阵和零矩阵，代入上述定理中得到：
$E[xa]=(IO)(μaμb)=μaVar[xa]=(IO)(ΣaaΣabΣbaΣbb)(IO)=Σaa\mathbb{E}[x_a]=\begin{pmatrix}{I}&{O}\end{pmatrix}\begin{pmatrix}\mu_a\\\mu_b\end{pmatrix}=\mu_a\\ Var[x_a]=\begin{pmatrix}{I}&{O}\end{pmatrix}\begin{pmatrix}\Sigma_{aa}&\Sigma_{ab}\\\Sigma_{ba}&\Sigma_{bb}\end{pmatrix}\begin{pmatrix}{I}\\{O}\end{pmatrix}=\Sigma_{aa}$
所以 $xa∼N(μa,Σaa)x_a\sim\mathcal{N}(\mu_a,\Sigma_{aa})$ ，边缘概率 $p(x_a)$ 就得到了。类似的， $xb∼N(μb,Σbb)x_b\sim\mathcal{N}(\mu_b,\Sigma_{bb})$ 。

求条件概率

对于两个条件概率，通常是用配方法（如 PRML 的证明），这里我们用一种构造法。首先引入三个量，令：
$xb⋅a=xb−ΣbaΣaa−1xaμb⋅a=μb−ΣbaΣaa−1μaΣbb⋅a=Σbb−ΣbaΣaa−1Σabx_{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}$
特别的，最后一个式子叫做 $Σaa\Sigma_{aa}$ 的 Schur Complementary。可以看到：
$xb⋅a=(−ΣbaΣaa−1In×n)(xaxb)x_{b\cdot a}=\begin{pmatrix}-\Sigma_{ba}\Sigma_{aa}^{-1}&{I}_{n\times n}\end{pmatrix}\begin{pmatrix}x_a\\x_b\end{pmatrix}$
再有定理，有：
$E[xb⋅a]=(−ΣbaΣaa−1In×n)(μaμb)=μb⋅aVar[xb⋅a]=(−ΣbaΣaa−1In×n)(ΣaaΣabΣbaΣbb)(−Σaa−1ΣbaTIn×n)=Σbb⋅a\mathbb{E}[x_{b\cdot a}]=\begin{pmatrix}-\Sigma_{ba}\Sigma_{aa}^{-1}&\mathbb{I}_{n\times n}\end{pmatrix}\begin{pmatrix}\mu_a\\\mu_b\end{pmatrix}=\mu_{b\cdot a}\\ Var[x_{b\cdot a}]=\begin{pmatrix}-\Sigma_{ba}\Sigma_{aa}^{-1}&{I}_{n\times n}\end{pmatrix}\begin{pmatrix}\Sigma_{aa}&\Sigma_{ab}\\\Sigma_{ba}&\Sigma_{bb}\end{pmatrix}\begin{pmatrix}-\Sigma_{aa}^{-1}\Sigma_{ba}^T\\{I}_{n\times n}\end{pmatrix}=\Sigma_{bb\cdot a}$
则对于我们构造的 $xb⋅ax_{b\cdot a}$ ，有 $xb⋅a∼N(μb⋅a,Σbb⋅a)x_{b\cdot a}\sim\mathcal{N}(\mu_{b\cdot a},\Sigma_{bb\cdot a})$ ，这里可以看到最初这个构造的设计中，核心的构造就是 $xb⋅ax_{b\cdot a}$ ，而 $μb⋅a,Σbb⋅a\mu_{b\cdot a},\ \Sigma_{bb\cdot a}$ 只是两个记号，在这种构造的推导下来表示一下均值和方差。

由我们最初的构造，有 $xb=xb⋅a+ΣbaΣaa−1xax_b=x_{b\cdot a}+\Sigma_{ba}\Sigma_{aa}^{-1}x_a$ ，再由定理：
$E[xb∣xa]=μb⋅a+ΣbaΣaa−1xa\mathbb{E}[x_b|x_a]=\mu_{b\cdot a}+\Sigma_{ba}\Sigma_{aa}^{-1}x_a$

$Var[xb∣xa]=Σbb⋅aVar[x_b|x_a]=\Sigma_{bb\cdot a}$

所以 $xb∣xa∼N(μb⋅a+ΣbaΣaa−1xa,Σbb⋅a)x_b|x_a\sim \mathcal{N}(\mu_{b\cdot a}+\Sigma_{ba}\Sigma_{aa}^{-1}x_a,\Sigma_{bb\cdot a})$ 。类似的， $xa∣xb∼N(μa⋅b+ΣabΣbb−1xb,Σaa⋅b)x_a|x_b\sim \mathcal{N}(\mu_{a\cdot b}+\Sigma_{ab}\Sigma_{bb}^{-1}x_b,\Sigma_{aa\cdot b})$ 。

根据边缘概率和条件概率求联合概率

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

这种类型的问题在线性高斯模型、PCA降维等机器学习模型中经常出现。
这里的 $Λ,L\Lambda, L$ 称为精度矩阵，它们是协方差矩阵的逆。

解：令 $y=Ax+b+ϵ,ϵ∼N(0,L−1)y=Ax+b+\epsilon,\epsilon\sim\mathcal{N}(0,L^{-1})$ ，且 $ϵ\epsilon$ 与 $x$ 相互独立，还是根据上节的定理，有
$E[y]=E[Ax+b+ϵ]=Aμ+bVar[y]=AΛ−1AT+L−1\mathbb{E}[y]=\mathbb{E}[Ax+b+\epsilon]=A\mu+b\\ Var[y]=A \Lambda^{-1}A^T+L^{-1}$
此时，就已经得到 $y∼N(Aμ+b,L−1+AΛ−1AT)y\sim\mathcal{N}(A\mu+b,L^{-1}+A\Lambda^{-1}A^T)$ ，即 $p(y)=N(Aμ+b,L−1+AΛ−1AT)p(y)=\mathcal{N}(A\mu+b,L^{-1}+A\Lambda^{-1}A^T)$ 。

因此：
$Var[x∣y]=Λ−1−Λ−1AT(L−1+AΛ−1AT)−1AΛ−1Var[x|y]=\Lambda^{-1}-\Lambda^{-1}A^T(L^{-1}+A\Lambda^{-1}A^T)^{-1}A\Lambda^{-1}$
引入 $z=(xy)z=\begin{pmatrix}x\\y\end{pmatrix}$ ，我们可以得到 $Cov[x,y]=E[(x−E[x])(y−E[y])T]Cov[x,y]=\mathbb{E}[(x-\mathbb{E}[x])(y-\mathbb{E}[y])^T]$ 。对于这个协方差可以直接计算：
$Cov(x,y)=E[(x−E[x])(y−E[y])T]=E[(x−μ)(Ax+b−Aμ−b+ϵ)T]=E[(x−μ)(Ax−Aμ+ϵ)T]=E[(x−μ)(Ax−Aμ)T+(x−μ)ϵT]=E[(x−μ)(Ax−Aμ)T]=E[(x−μ)(x−μ)T]AT=Var[x]AT=Λ−1AT\begin{align} Cov(x,y)&=\mathbb{E}[(x-\mathbb{E}[x])(y-\mathbb{E}[y])^T]\\ &=\mathbb{E}[(x-\mu)(Ax+b-A\mu-b+\epsilon)^T]\\ &=\mathbb{E}[(x-\mu)(Ax-A\mu+\epsilon)^T]\\ &=\mathbb{E}[(x-\mu)(Ax-A\mu)^T+(x-\mu)\epsilon^T]\\ &=\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{align}$
注意到协方差矩阵的对称性，所以 $p(z)=N(μAμ+b),(Λ−1Λ−1ATAΛ−1L−1+AΛ−1AT))p(z)=\mathcal{N}\begin{pmatrix}\mu\\A\mu+b\end{pmatrix},\begin{pmatrix}\Lambda^{-1}&\Lambda^{-1}A^T\\A\Lambda^{-1}&L^{-1}+A\Lambda^{-1}A^T\end{pmatrix})$ 。根据上一节的公式，我们可以得到：
$E[x∣y]=μ+Λ−1AT(L−1+AΛ−1AT)−1(y−Aμ−b)\mathbb{E}[x|y]=\mu+\Lambda^{-1}A^T(L^{-1}+A\Lambda^{-1}A^T)^{-1}(y-A\mu-b)$

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

故得到： $p(x∣y)=N(μ+Λ−1AT(L−1+AΛ−1AT)−1(y−Aμ−b),Λ−1−Λ−1AT(L−1+AΛ−1AT)−1AΛ−1)p(x|y)=\mathcal{N}(\mu+\Lambda^{-1}A^T(L^{-1}+A\Lambda^{-1}A^T)^{-1}(y-A\mu-b),\Lambda^{-1}-\Lambda^{-1}A^T(L^{-1}+A\Lambda^{-1}A^T)^{-1}A\Lambda^{-1})$ 。