本文为第一系列的原创文章之四, 目的在于总结“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}\) 即可完成分解.
具体的模拟涉及到大量简单但是繁琐的计算, 故不在此列出, 仅给出以上内容作为数学推导说明原理.