统计模拟(四)——多维正态分布与乔列斯基分解

本文为第一系列的原创文章之四, 目的在于总结“Simulation” -Sheldon M.Ross一书中出现的模拟随机变量的方法.

Multivariate Normal Distribution

在上一篇文章中, Box-Muller 生成了一个二维的正态分布 $X$ 与 $Y$ , 二者相互独立, 但如果想要生成有相关性的多维(二维以及以上)变量, Box-Muller变换就有些捉襟见肘了.

在数学中, 假设向量 $\overrightarrow{Z}=(Z_1,Z_2,Z_3,…,Z_m)^T$ 独立同分布于标准正态分布, 如果对于矩阵 $A_{i\times j},\quad i=1,2,3,…,n,\quad j=1,2,3,…m$与向量 $\overrightarrow{\mu}=(\mu_1,\mu_2,\mu_3,…,\mu_n)^T$ 有:

$$\overrightarrow{X}=A\overrightarrow{Z}+\overrightarrow{\mu}$$

那么我们认为向量 $\overrightarrow{X}=(X_1,X_2,X_3,…,X_n)^T$ 服从多维正态分布.

很明显易得: 向量 $\overrightarrow{X}$ 的均值为 $\overrightarrow{\mu}$ , 其协方差为

$$\begin{align*} Cov(X_i,\ X_j)&=Cov(\sum_{k=1}^{m} a_{ik}Z_k,\ \sum_{r=1}^{m} a_{jr}Z_r)\\\\ &=\sum_{k=1}^{m} \sum_{r=1}^{m} Cov(a_{ik}Z_k,\ a_{jr}Z_r))\\\ &=\sum_{k=1}^{m} \sum_{r=1}^{m} a_{ik} a_{jr} Cov(Z_k,\ Z_r)\\\ &=\sum_{k=1}^{m} a_{ik} a_{jr} \end{align*}$$

上面的推导中:

$$Cov(Z_k,\ Z_r)= \left \{ \begin{matrix} 1,\quad if\ r=k\\ 0,\quad if\ r\neq k \end{matrix}\right.$$

那么可以得到, $\overrightarrow{X}$ 的协方差矩阵可以表示为:

$$C=AA^T$$

该多维正态分布的一个重要的特性是, 其分布由 $\overrightarrow{\mu}$ 与协方差矩阵 $C$ 唯一确定. 可由 $X$ 的联合矩母函数(Joint Moment Generating Function) 证明:

令
$$W=\sum_{t=1}^{n} t_i X_i\\ E[e^{W}]=e^{E(W)+\frac{var(W)}{2}}$$

由于
$$E(W)=E(\sum_{t=1}^{n} t_i X_i)=E(\sum_{t=1}^{n} t_i \mu_i)$$

且

$$\begin{align*} Var(\sum_{t=1}^{n} t_i X_i)&=Cov(\sum_{i=1}^{n} t_i X_i,\ \sum_{j=1}^{n}t_j X_j)\\\ &=\sum_{i=1}^{n} \sum_{j=1}^{n}t_i t_j Cov(X_i,\ X_j)\end{align*}$$

由于矩母函数可以认为是某特定概率函数唯一对应的函数, 所以可以认为该多维正态分布由 $\overrightarrow{\mu}$ 与协方差矩阵 $C$ 唯一确定.

Choleski Decomposition

假设现在已知某多维正态分布的均值 $\overrightarrow{\mu}$ 与协方差 $C$, 并打算模拟该多维正态分布. 由于

$$\overrightarrow{X}=A\overrightarrow{Z}+\overrightarrow{\mu}$$

生成 $n$ 个独立的标准正态随机变量 $Z_i$ , 找到矩阵 $A$ 使得

$$C=AA^T$$

这样就可以对多位正态分布 $\overrightarrow{X}$ 进行模拟.

所以, 方便的做法是, 已知 $C$ 的情况下, 对 $C$ 进行乔列斯基分解 (Choleski Decomposition). 将 $C$ 分解为一个对角阵转置相乘的形式.

乔列斯基分解是对称正定矩阵的LU分解. 其过程与证明不属于 “Simulation” 一书的主要内容, 但是考虑到之前数值分析一课中学习过相关的证明与计算, 所以一并在此附上相关介绍.

定理:
设 $A\in R^{n\times n}$ , 如果 $A$ 所有的顺序主子式均不为 $0$ , 则方阵 $A$ 存在唯一的分解式:
$$A=LDR$$
其中, $L, \ R$ 分别是单位下三角阵与单位上三角阵, $D$ 是非奇异对角阵. (证明略)

那么当 $A$ 为对称正定矩阵, 有

$$\begin{align*} A&=A^T\\\ LDR&=R^TDL^T\\\ \end{align*}$$

可以得到

$$\begin{align*} L&=R^T\\\ A&=LDL^T\\\ \end{align*}$$

如果 $D=diag[d_1, d_2,…,d_n]$ , 由于A正定, 记 $D^{\frac{1}{2}}=diag[\sqrt{d_1},\sqrt{d_2},\sqrt{d_3},…,\sqrt{d_n}]$, 那么有 $A=LD^{\frac{1}{2}}(LD^{\frac{1}{2}})^T$, 记

$$\tilde{L}= LD^{\frac{1}{2}}$$

那么有

$$A=\tilde{L}{\tilde{L}}^T$$

这就是对称正定矩阵的乔列斯基分解. 为了方便 $\tilde{L}$ 常被记为 $L$. 其中

$$\begin{align*} L=\begin{bmatrix}l_{11}& & & \\\ l_{21}& l_{22} & & \\\ \vdots & \vdots & \ddots & \\\ l_{n1}& l_{n2} & \cdots & l_{nn}\end{bmatrix}\end{align*}$$

$$a_{i,j}=l_{i1} \times l_{j1}+ l_{i2}\times l_{j2}+,…,+ l_{ii} \times l_{ji}\quad (i<=j)$$

直观的认为, $a_{i,j}$ 是将矩阵 $L$ 的第 $i$ 行和第 $j$ 行提取出来, 上下相乘并加起来. 解出每一个 $l_{i,j}$ 即可完成分解.

具体的模拟涉及到大量简单但是繁琐的计算, 故不在此列出, 仅给出以上内容作为数学推导说明原理.