随机过程（2.2）—— 多维高斯分布

• 随机过程（2.1）—— 多维随机变量（随机向量）
• 多维高斯分布是一种特殊的多维随机分布，应用非常广泛。本文参考【机器学习】【白板推导系列】

---

1. 定义

• 若 $n$ 维 r.v. $X$ 的概率密度函数为
  $$f(x)=f(x_1,...,x_n)=\frac{1}{(2\pi {)}^{n/2}|B{|}^{1/2}}{e}^{-\frac{1}{2}(\mathbf{x}-\mathbf{a}{)}^{\top }{\mathbf{B}}^{-1}(\mathbf{x}-\mathbf{a})}$$ 其中 $B$ 为 $n$ 阶实对称正定矩阵（所有特征值 > 0），则称随机向量 $X$ 服从期望为 $a$，协方差矩阵为 $B$ 的多维正态分布，记为 $X\sim N(a,B)$
• 注意这里 $B$ 是一个协方差矩阵，展开为
  $$\begin{array}{rl}B& =\left[\begin{array}{cccc}D{X}_1& \text{Cov}(X_1,X_2)& \dots & \text{Cov}(X_1,X_n)\\ \text{Cov}(X_2,X_1)& D{X}_2& \dots & \text{Cov}(X_2,X_n)\\ ⋮& ⋮& & ⋮\\ \text{Cov}(X_n,X_1)& \text{Cov}(X_n,X_2)& \dots & D{X}_n\end{array}\right]\\ & =\left[\begin{array}{cccc}{\sigma }_1^2& {\rho }_{12}{\sigma }_1{\sigma }_2& \dots & {\rho }_{1n}{\sigma }_1{\sigma }_n\\ {\rho }_{21}{\sigma }_2{\sigma }_1& {\sigma }_2^2& \dots & {\rho }_{2n}{\sigma }_2{\sigma }_n\\ ⋮& ⋮& & ⋮\\ {\rho }_{n,1}{\sigma }_n{\sigma }_1& {\rho }_{n,2}{\sigma }_n{\sigma }_2& \dots & {\sigma }_n^2\end{array}\right]\end{array}$$ 显然它 $B$ 是实对称矩阵，通常情况下（只要 $X_i$ 不为常数）其任意阶顺序主子式 > 0，故正定。这种实对称正定矩阵是一种正定 hermitian 矩阵，可以做乔里斯基分解 $B=L{L}^{\top }$（$L$ 是某下三角矩阵），这是处理多维高斯分布问题的常用技巧之一
• 显然 $f(x)$ 在任意取值下都是非负的，下面证明 ${\int }_{R^n}f(x)dx=1$ 以说明它是一个合法的概率密度函数，这个证明中用到两个技巧
  1. 上面提到的 $B=L{L}^{\top }$
  2. 向量对向量求导会得到一个 Jacobi 行列式，令 $x=Ly+a$，则 $\frac{\partial \mathbf{x}}{\partial \mathbf{y}}=|L|$
     

2. 理解

2.1 几何意义

• $X\sim N(a,B)$ 的概率密度函数 $f(x)=\frac{1}{(2\pi {)}^{n/2}|B{|}^{1/2}}{e}^{-\frac{1}{2}(\mathbf{x}-\mathbf{a}{)}^{\top }{\mathbf{B}}^{-1}(\mathbf{x}-\mathbf{a})}$，关注其中的指数部分，这是一个二次型，设为 $$△=(\mathbf{x}-\mathbf{a}{)}^{\top }{\mathbf{B}}^{-1}(\mathbf{x}-\mathbf{a})$$

  1. 首先明确 $(\mathbf{x}-\mathbf{a})$ 尺寸为 $n\times 1$，$\mathbf{B}$ 尺寸为 $n\times n$，$△$ 是一个数，可以看作 $\mathbf{x}$ 和 $\mathbf{a}$ 间的马氏距离

  2. 注意 $B$ 是实对称矩阵，一般情况下（$X_i$ 不为常数）认为它是正定矩阵，其具有以下性质

     1. 正定对称 $\Rightarrow$ 乔里斯基分解：$B=L{L}^{\top }$，其中 $L$ 是下三角矩阵
     2. 正定对称 $\Rightarrow$ LDL分解：$B=LD{L}^{\top }$，其中 $D$ 是一个对角阵，$L$ 是主对角线全为 1 的下三角矩阵
     3. 实对称 $\Rightarrow$ $B_{n\times n}$ 正交相似于其特征值组成的对角阵，即 $B=U\Lambda U^{\top }$，其中 $U$ 是正交矩阵（有 $U{U}^{\top }=U^{\top }U=E$），其列向量为 $B$ 的特征向量， Λ = diag { λ 1 , λ 2 , . . . , λ p } \Lambda = \text{diag} \{\lambda_1,\lambda_2,...,\lambda_p\} Λ=diag{λ1​,λ2​,...,λp​}（${\lambda }_i$ 是 $B$ 的特征值）。这是因为 $n$ 阶实对称矩阵的所有特征值都是实数，各个特征值的代数重数和几何重数相等（有 $n$ 个线性无关特征向量），且所有特征向量相互正交（参考此处）。此结论也可理解为 Schur 定理在实数域上的推论
     4. 可相似对角化 $\Rightarrow$ $B_{n\times n}$ 可以 谱分解：$B=P\Lambda P^{-1}$，其中 $P$ 的列向量是 $B$ 的（右）特征向量，$P^{-1}$ 的行向量是 $B$ 的左特征向量（参考此处）， Λ = diag { λ 1 , λ 2 , . . . , λ p } \Lambda = \text{diag} \{\lambda_1,\lambda_2,...,\lambda_p\} Λ=diag{λ1​,λ2​,...,λp​}（${\lambda }_i$ 是 $B$ 的特征值）。注意性质 3 中的 $U$ 是正交矩阵，有 $U^{-1}=U^{\top }$，因此实对称矩阵 $B$ 的谱分解和 3 中的相似对角化是一样的

     利用上述性质 3/4 分析 $B^{-1}$，可如下展开（其中 $u_i$ 是列向量，尺寸 $n\times 1$）
     $$\begin{array}{rl}{B}^{-1}& =(U\Lambda U^{\top }{)}^{-1}\\ & =(U^{\top }{)}^{-1}{\Lambda }^{-1}{U}^{-1}\\ & =U{\Lambda }^{-1}{U}^{\top }\\ & =\sum _{i=1}^n{u}_i\frac{1}{{\lambda }_i}{u}_i^{\top }\end{array}$$

• 进一步考虑多元高斯分布指数部分的二次型
  $$\begin{array}{rl}△& =(\mathbf{x}-\mathbf{a}{)}^{\top }{\mathbf{B}}^{-1}(\mathbf{x}-\mathbf{a})\\ & =(\mathbf{x}-\mathbf{a}{)}^{\top }\sum _{i=1}^n{u}_i\frac{1}{{\lambda }_i}{u}_i^{\top }(\mathbf{x}-\mathbf{a})\\ & =\sum _{i=1}^n(\mathbf{x}-\mathbf{a}{)}^{\top }{u}_i\frac{1}{{\lambda }_i}{u}_i^{\top }(\mathbf{x}-\mathbf{a})\\ & \stackrel{设{\mathbf{y}}_{\mathbf{i}}=(\mathbf{x}-\mathbf{a}{)}^{\top }{\mathbf{u}}_{\mathbf{i}}}{=}\sum _{i=1}^n\frac{y_i^2}{{\lambda }_i}\end{array}$$ 可见 $△=k$ 时，$△={\sum }_{i=1}^n\frac{y_i^2}{{\lambda }_i}=k$ 是 $p$ 维空间中一个超椭圆。注意到 $y_i=(x-a{)}^{\top }{u}_i$，可以看作先把 $x$ 沿 $a$ 方向平移，然后向 $u_i$ 方向上的投影。不妨在 $u_i$ 方向设置坐标轴 $y_i$，二维情况（$n=2$）的示意图如下
  
  可见随着 $k$ 值变化，$△$ 对应到空间中一系列超椭圆，若把 $k$ 看做等高线高度，这一系列椭圆就形成了一个柱状体。进一步考虑整个 n 元高斯分布的概率密度函数 $f(x)=\frac{1}{(2\pi {)}^{n/2}|B{|}^{1/2}}{e}^{-\frac{1}{2}△}$，前面分数部分是个常数，$e^x$ 则是和 $x$ 正相关，所以概率密度函数 $f(x)$ 也是一个相似的柱状体
• 可以如下绘制二维情况下 $f(x)$ 图像，这里设置期望为 $a=\left[\begin{array}{c}0\\ 0\end{array}\right]$，协方差矩阵为 $B=\left[\begin{array}{cc}0.8& 0.2\\ 0.2& 0.2\end{array}\right]$

  %matplotlib notebook
  import numpy as np
  import scipy.stats as st
  import matplotlib.pylab as plt
  from mpl_toolkits.mplot3d import Axes3D
  
  mu = np.array([0,0])
  cov = np.array([[0.8, 0.2], 
                  [0.2, 0.2]])
  
  fig = plt.figure(figsize = (10,5))
  a0 = fig.add_subplot(1,2,1,label='a0',projection='3d') 
  a1 = fig.add_subplot(1,2,2,label='a1',projection='3d') 
  
  x, y = np.mgrid[-2.5:2.5:.1, -2.5:2.5:.1]
  pos = np.empty(x.shape + (2,))
  pos[:, :, 0] = x; pos[:, :, 1] = y
  rv = st.multivariate_normal(mu, cov)   # 生成多元正态分布
  a0.scatter(x, y, rv.pdf(pos),s=1,alpha=0.5,cmap="rainbow") 
  a1.plot_surface(x, y, rv.pdf(pos),alpha=0.5,cmap=plt.cm.cool)
  

2.2 参数数量

• 考察上述 $B$ 矩阵的参数个数，由于 $B$ 是对称矩阵，当尺寸为 $n\times n$ 时，参数有 $1+2+...+n=\frac{n^2+n}{2}$ 个，这时所有参数都非零，意味着 $n$ 元高斯随机变量 $x=\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]$ 中任意两个维度 $x_i,x_j,i\ne j$ 相关。举例来说，期望为 $a=\left[\begin{array}{c}1\\ 2\end{array}\right]$，协方差矩阵为 $B=\left[\begin{array}{cc}0.8& 0.2\\ 0.2& 0.2\end{array}\right]$ 时属于这种情况，此时 $△$ 对应的超椭圆是倾斜的，如下所示

  %matplotlib notebook
  import numpy as np
  import scipy.stats as st
  import matplotlib.pylab as plt
  from mpl_toolkits.mplot3d import Axes3D
  
  mu = np.array([1,2])
  cov = np.array([[0.8, 0.2], 
                  [0.2, 0.2]])
  
  fig = plt.figure(figsize = (10,5))
  a0 = fig.add_subplot(1,2,1,label='a0',projection='3d') 
  a1 = fig.add_subplot(1,2,2,label='a1') 
  
  x, y = np.mgrid[-1.5:3.5:.1, -0.5:4.5:.1]
  pos = np.empty(x.shape + (2,))
  pos[:, :, 0] = x; pos[:, :, 1] = y
  rv = st.multivariate_normal(mu, cov)   # 生成多元正态分布
  a0.scatter(x, y, rv.pdf(pos),s=1,alpha=0.5,cmap="rainbow") 
  a1.contourf(x, y, rv.pdf(pos))         # 等高线
  a1.grid(alpha=0.5)                     # 坐标网格
  

• 希望减少参数数量，可以假设 $n$ 元高斯随机变量 $x=\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]$ 中任意两个维度 $x_i,x_j,i\ne j$ 相互独立，这时 $B$ 除了主对角元素外其他元素都是 $0$，参数减少到 $n$ 个。举例来说，期望为 $a=\left[\begin{array}{c}1\\ 2\end{array}\right]$，协方差矩阵为 $B=\left[\begin{array}{cc}0.8& 0\\ 0& 0.2\end{array}\right]$ 时属于这种情况，此时 $△$ 对应的超椭圆是正的，如下所示

• 进一步减少参数数量，可以假设 $n$ 元高斯随机变量 $x=\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]$ 中任意两个维度 $x_i,x_j,i\ne j$ 相互独立且各向同性，这时 $B$ 除了主对角元素外其他元素都是 $0$，且主对角线元素都为相等正数，参数减少到 $1$ 个。举例来说，期望为 $a=\left[\begin{array}{c}1\\ 2\end{array}\right]$，协方差矩阵为 $B=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$ 时属于这种情况，此时 $△$ 对应的超椭圆变为正圆，如下所示

• 可以使用极大似然估计来估计多维高斯分布的参数，请参考：一文看懂 “极大似然估计” 与 “最大后验估计” —— 极大似然估计篇

3. 特征函数

• 注意到 $n$ 元正态分布函数的概率密度函数很复杂，要计算 $B$ 的逆和行列式，所以我们一般通过特征函数来研究其性质
• 普通一元正态分布 $X\sim N(\mu ,{\sigma }^2)$ 的特征函数为（证明见 4.2 节分量独立性证明）
  $${\phi }_X(t)=e^{it\mu -\frac{1}{2}{t}^2{\sigma }^2}$$
• $n$ 元正态分布 $X\sim N(a,B)$ 的特征函数为
  $${\phi }_{\mathbf{x}}(\mathbf{t})=e^{i{\mathbf{t}}^{\top }\mathbf{a}-\frac{1}{2}{\mathbf{t}}^{\top }\mathbf{Bt}}$$ 详细证明过程如下
  

4. 性质

4.1 边缘分布

1. “多元正态随机向量” 的每个元素是 “一个正态随机变量”：$X_j\sim N(a_j,b_{jj})$
2. “多元正态随机向量” 的部分向量仍为 “多元正态随机向量”：若 $X\sim N(a,B)$，则其第 $k_1,...,k_m$ 分量组成的随机向量满足
   $$X^{\ast }=\left[\begin{array}{c}{X}_{k1}\\ X_{k2}\\ ⋮\\ X_{km}\end{array}\right]\sim N(a^{\ast },B^{\ast })$$ 其中 $B^{\ast }=\left[\begin{array}{cc}{I}_m& 0\end{array}\right]B\left[\begin{array}{c}{I}_m\\ 0\end{array}\right]$ 是保留 $B$ 的第 $k_1,...,k_m$ 行列所得的 $m\times m$ 矩阵，$a^{\ast }=\left[\begin{array}{c}{a}_{k_1}\\ a_{k_2}\\ ⋮\\ a_{k_m}\end{array}\right]$ 是 $a$ 的第 $k_1,...,k_m$ 分量拼成的向量。从特征函数角度证明如下
   

4.2 分量独立性

• 独立性：若 $X=[X_1,X_2,\dots ,X_n{]}^{\top }\sim N(a,B)$，以下陈述等价（注: 随机变量 互不相关 指没有线性关系，即协方差为0；独立 指没有一切关系）

  1. $X_1,X_2,\dots ,X_n$ 相互独立
  2. $X_1,X_2,\dots ,X_n$ 两两独立
  3. $X_1,X_2,\dots ,X_n$ 两两互不相关
  4. $B$ 为对角阵（即不同的两个元协方差为0）

• 证明：$1\Rightarrow 2$ 显然；$2\Rightarrow 3$ 是随机变量性质；$3\Rightarrow 4$ 是互不相关定义；只需证明 $4\Rightarrow 1$ 即得四者等价
  

4.3 线性变换

• 一组正态分布随机变量的线性组合（多元正态随机向量的线性变换）仍然服从正态分布

  1. 设有 $n$ 元随机向量 $X\sim N(a,B)$，则对 $\forall l\ne 0\in {\mathbb{R}}^n$ 有
     $$l^{\top }X\sim N(l^{\top }a,l^{\top }Bl)$$ 这里 $l^{\top }X={\sum }_{i=1}^n{l}_i{X}_i$ 其实就是对 $X$ 中所有正态随机变量 $X_i$ 线性组合得到的一元正态随机变量。上式从特征函数角度可证明如下
     
     注意，这里 ${\phi }_{\mathbf{x}}(tl)$ 是 n 元正态分布的特征函数；${\phi }_{{\mathbf{l}}^{\top }\mathbf{x}}(t)$ 是一元正态分布的特征函数，最后要落到这个上

  2. 设有 $n$ 元随机向量 $X\sim N(a,B)$，$C$ 为 $m\times n$ 矩阵，且行向量线性无关，则
     $$\mathbf{CX}\sim N(\mathbf{Ca},\mathbf{CB}{\mathbf{C}}^{\top })$$ 其实这里 $C$ 中的每一行都对应了一个 1 中的线性组合。上式可以用特征函数证明如下
     
     注意，${\phi }_{\mathbf{x}}(C^{\top }t)$ 表示的是 $n$ 元正态分布的特征函数；${\phi }_{\mathbf{Cx}}(t)$ 表示的是 $m$ 元正态分布的特征函数（$C_{m\times n}$ 将 $n$ 元的 $x$ 变换为 $m$ 维），最后要落到这个上

5. 综合例题