## Goal: sample from a bivariate normal 
## distribution

set.seed(99355)
# number of samples 
N = 5000

# simulate 2*N random variables
# from a standard normal distribution
x = rnorm(2*N)
x = matrix(x, ncol=N)

# we are interested in getting samples from 
# N(mu, Sigma), where Sigma is the covariance 
# matrix.
Sigma <- matrix(c(10,3,3,2),2,2)
mu <- c(1,2)

# compute the Cholesky decomposition returns A, 
# so that t(A)*A = Sigma
A = chol(Sigma)
# check
(t(A) %*% A)

y = mu + t(A) %*% x
y = t(y)

colMeans(y)
var(y)

## compare with integrated multivariate normal 
## distribution in R
library(MASS)
xR = mvrnorm(n=N, mu, Sigma)

var(xR)
colMeans(xR)

## use a kernel density estimator to plot 
## the distribution
library(gplots)
y.kde = kde2d(y[,1], y[,2], n=100)
xR.kde = kde2d(xR[,1], xR[,2], n=100)

par(mfrow=c(1,2))
image(y.kde, xlim=c(-10,10), ylim=c(-2,6))
contour(y.kde, add = T)
image(xR.kde, xlim=c(-10,10), ylim=c(-2,6))
contour(xR.kde, add = T)