Representation of 3d Gaussian as vector of functions of 3iid gaussians.

Question

If I have a multivariate guassian. $X \sim \mathcal{N}(0^3,\Sigma)$,where $\Sigma= \begin{bmatrix} 1 &r_{12}&r_{13}\\ r_{12}&1&r_{23}\\r_{13}&r_{23}&1\end{bmatrix}$

I want to write it as an array of functions of three standard Gaussians. $N_1,N_2,N_3\sim\mathcal{N}(0,1)$

I believe I can use successive marginalization starting from the bottom to write

$X = \left[\begin{array}{c}r_{13}N_3 + (r_{12}-r_{13}r_{23})\sqrt{1-r_{23}^2}N_2+\sqrt{1-r_{13}^2-\frac{1}{1-r_{23}^2}(r_{12}-r_{13}r_{23})^2}N_1 \\r_{23}N_3+\sqrt{1-r_{23}^2}N_2\\N_3\end{array}\right]$

However when I test this by using numpy to create a realization of $X$ and finding the covariance matrix of $X$ it's close but not quite correct. Specifically $r_{12}$ of the realization is off.

I am wondering if anyone could tell me what I am doing wrong. Is the idea of using successive marginalization flawed?

Answer 1

From the properties of multivariate gaussian distribution, if $X \sim \mathcal{N}(\mu, \Sigma)$, then $Y = AX + b \sim \mathcal{N}(\mu + b, A\Sigma A^T)$. In particular, $X \sim \mathcal{N}(\mu, \Sigma)$ can be decomposed as $X = \Sigma^{1/2} N + \mu$, where $N$ has standard gaussian distribution $\mathcal{N}(0, I)$ and $\Sigma^{1/2}$ is the matrix square root of $\Sigma$ (which is defined according to the general theory of function on symmetric matrices) with properties $\Sigma^{1/2} \Sigma^{1/2} = \Sigma$ and $\left(\Sigma^{1/2}\right)^T = \Sigma^{1/2}$.

In your case it's needed to find $\Sigma^{1/2}$ for given $\Sigma$. You found that $$ X = \begin{bmatrix} \sqrt{1 - r_{13} - \frac{1}{1-r_{23}} (r_{12} - r_{13} r_{23})^2} & (r_{12} - r_{13} r_{23}) \sqrt{1 - r_{23}^2} & r_{13} \\ 0 & \sqrt{1- r_{23}^2} & r_{23} \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} N_1 \\ N_2 \\ N_3 \end{bmatrix} = R \times N $$ $R$ is not symmetric and from what I computed, $R R^T \neq \Sigma$.

To find $\Sigma^{1/2}$ explicitelly, you need to diagonalize it, $\Sigma = Q \Lambda Q^T$ where $Q$ is orthogonal and $\Lambda$ is diagonal with eigenvalues on the diagonal, and then calculate $\Sigma^{1/2} = Q \Lambda^{1/2} Q^T$, where $\Lambda^{1/2}$ is obtained from $\Lambda$ by replacing eigenvalues with their positive square roots.

If you need only numerical procedure in Python, you can use Numpy and sqrtm function from Scipy in the following way:

>>> S = np.array([[1, 0.3, 0.5], [0.3, 1, 0.7], [0.5, 0.7, 1]])
>>> S_sqrt = sqrtm(S)
>>> N = np.random.randn(3, 1_000_000) # 1 million samples of 3D standard gaussian
>>> X = S_sqrt @ N
>>> np.cov(X)
array([[0.99990112, 0.2979038 , 0.50002984],
       [0.2979038 , 1.00028636, 0.7001492 ],
       [0.50002984, 0.7001492 , 1.00095704]])