Multivariate normal distribution has the following p.d.f.:

and for a bivariate in specific, it is:

where the covariance $\sigma_{12}^2 = \rho\sigma_1\sigma_2$. See also the note on multivariate Gaussian distribution.

The multivariate normal random variable can be considered as a combination of the mean vector and zero-mean multivariate random variable. Assume $\mathbf{y}$ to be a $n$-vector of iid standard normal random variables, i.e., $\bar{\mathbf{y}}=\mathbf{0}$ and $E[\mathbf{yy}^T]=\mathbf{I}$. Let there be a matrix of constants $\mathbf{S}$ of size $m\times n$, and a $m$-vector of constants (the mean) $\mathbf{m}$, and define

Then $\mathbf{x}$ is a multivariate normal random variable, which the covariance matrix is

Finding $\mathbf{S} = \Sigma_{\mathbf{x}}^{1/2}$ given covariance matrix $\Sigma_{\mathbf{x}}$ is typically using Cholesky decomposition.

Taking again bivariate as an example. If the covariance matrix is

In Julia, we can plot a set of samples of the distribution as follows (assumed zero-mean):

Sigma = [4 9; 9 25]
S = chol(Sigma) # upper triangular part, i.e. S' * S == Sigma
iidnorm = randn(2,1000)
mvnorm = S' * iidnorm
plot(x=mvnorm[1,:], y=mvnorm[2,:], Geom.point)


Alternatively, such sampling can also be done using Monte Carlo method, such as the Gibbs sampling. The algorithm is Gibbs sampling is as follows:

• Assume some sensible value for the multivariate random variable $x^{(0)} = (x_1^{(0)}, x_2^{(0)}, \cdots, x_n^{(0)})$ as initial value
• For $k=0$ up to $k=M+N$
• For $i=1$ up to $i=n$, i.e. each component of the multivariate random variable
• Set $x_i^{(k+1)}$ according to probability $p(x_i\mid x_1^{(k+1)}, \cdots, x_{i-1}^{(k+1)}, x_{i+1}^{(k)}, \cdots, x_n^{(k)})$
• Discard the first $M$ samples of $x^{(k)}$ above and use only the last $N$ as terminal distribution

For the particular case of zero-mean bivariate normal random variables, the loop above means:

• For $k=0$ up to $k=M$
• set $x_1^{(k+1)} = \rho x_2^{(k)}\sigma_1/\sigma_2 + \sqrt{\sigma_1^2(1-\rho^2)} Z_1$ where $Z_1 \sim N(0,1)$
• set $x_2^{(k+1)} = \rho x_1^{(k+1)}\sigma_2/\sigma_1 + \sqrt{\sigma_2^2(1-\rho^2)} Z_2$ where $Z_2 \sim N(0,1)$

The formula $x_1^{(k+1)} = \rho x_2^{(k)} + \sqrt{\sigma_1^2(1-\rho^2)} Z_1$ comes from the following:

And this is the corresponding code:

Sigma = [4 9; 9 25]
sigma1 = sqrt(Sigma[1,1])
sigma2 = sqrt(Sigma[2,2])
rho = Sigma[1,2]/(sigma1*sigma2)
states = []
push!(states, [0,0])
for i = 1:1500
Z1,Z2 = randn(2)
x1, x2 = last(states)
x1 = rho * x2 * sigma1 / sigma2 + sqrt(sigma1^2 * (1-rho^2)) * Z1
x2 = rho * x1 * sigma2 / sigma1 + sqrt(sigma2^2 * (1-rho^2)) * Z2
push!(states, [x1,x2])
end
states = hcat(states...)
plot(x=states[1,501:end], y=states[2,501:end], Geom.point)