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

and for a bivariate in specific, it is:

where the covariance . 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 to be a -vector of iid standard normal random variables, i.e., and . Let there be a matrix of constants of size , and a -vector of constants (the mean) , and define

Then is a multivariate normal random variable, which the covariance matrix is

Finding given covariance matrix 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
using Gadfly
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 as initial value
- For up to
- For up to , i.e. each component of the multivariate random variable
- Set according to probability

- For up to , i.e. each component of the multivariate random variable
- Discard the first samples of above and use only the last as terminal distribution

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

- For up to
- set where
- set where

The formula 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)
```