# install.packages("mvtnorm")
library(mvtnorm)
set.seed(12345)

## Assign mu, Sigma   
mu <- c(60, 100, 80, 75)
Sigma <- matrix(c(12,2,5,-4, 2,16,3,-2, 5,3,20,-2, -4,-2,-2,18),ncol=4)

## Draw Random sample of Multivariate Normal Distribution
X <- rmvnorm(5000, mu, Sigma)

## Compute sample means, variance-covariance matrix
colMeans(X)
cov(X)

######################################################################
## Brute-force
## Draw a random sample of 4*n of N(0,1) and put in 4 columns of n rows
X1 <- matrix(rnorm(20000,0,1),ncol=4)

## Obtain the Cholesky "square-root" matrix (upper triangular)
(Sigma12.c <- chol(Sigma))

## Obtain the "square-root" matrix from eigenvalues/vectors
P <- eigen(Sigma)$vector
L12 <- diag(sqrt(eigen(Sigma)$value))
(Sigma12.e <- P %*% L12 %*% t(P))


## Confirm Sigma12'Sigma12 = Sigma
t(Sigma12.c) %*% Sigma12.c
t(Sigma12.e) %*% Sigma12.e

## Create Mu matrix  (n x 4)
mu1 <- matrix(rep(mu,5000),byrow=T,ncol=4)

## Transform N(0,1) matrix to MVN(mu, Sigma)
X2.c <- mu1 + X1 %*% Sigma12.c
X2.e <- mu1 + X1 %*% Sigma12.e

## Compute sample means, variance-covariance matrix
colMeans(X2.c)
cov(X2.c)
colMeans(X2.e)
cov(X2.e)