[Math] How to sample a bivariate Gaussian distribution using Gibbs sampling

samplingstatistics

I'm trying to sample a bivariate Gaussian distribution using Gibbs sampling, but I think I don't have the correct conditional probabilities. According to this lecture slides, the conditional expectation and variance of the bivariate Gaussian distribution are:

$E[X|Y=y]=\mu_X+\sigma_X\rho(\frac{\displaystyle y-\mu_Y}{\displaystyle \sigma_Y})$

and

$Var[X|Y=y]={{\sigma}_X}^2(1-{\rho}^2)$

So I'm sampling the bivariate Gaussian distribution using the standard normal distribution as follows:

$f(X|Y=y)=\mu_X+\sigma_X\rho(\frac{\displaystyle y-\mu_Y}{\displaystyle \sigma_Y})+{{\sigma}_X}^2(1-{\rho}^2)\mathcal{N}(0,1)$

$f(Y|X=x)=\mu_Y+\sigma_Y\rho(\frac{\displaystyle x-\mu_X}{\displaystyle \sigma_X})+{{\sigma}_Y}^2(1-{\rho}^2)\mathcal{N}(0,1)$

Are these equations ok to sample the bivariate Gaussian?

Best Answer

Conditionally on $[Y=y]$, $X$ is distributed like $$ \mu_X+\sigma_X\rho\sigma_Y^{-1}(y-\mu_Y)+\sigma_X\sqrt{1-\rho^2}Z, $$ where $Z$ is standard normal (note the factor $\sigma_X\sqrt{1-\rho^2}$ instead of its square). Likewise for $Y$ conditionally on $[X=x]$. But to sample $(X,Y)$, one should not use these two equations, rather, one can generate $Y$ by $$ Y=\mu_Y+\sigma_YT, $$ where $T$ is standard normal, then generate $X$ by $$ X=\mu_X+\sigma_X\rho\sigma_Y^{-1}(Y-\mu_Y)+\sigma_X\sqrt{1-\rho^2}Z, $$ where $Z$ is standard normal independent of $T$, or, equivalently, $$ X=\mu_X+\sigma_X\rho T+\sigma_X\sqrt{1-\rho^2}Z. $$