rmultnorm <- function(n, mu, vmat, tol = 1e-07)  
    # a function to generate random multivariate Gaussians 
    {
       p <- ncol(vmat)  
       if (length(mu)!=p)
           stop("mu vector is the wrong length")  
       if (max(abs(vmat - t(vmat))) > tol) 
           stop("vmat not symmetric") 
       vs <- svd(vmat)  
       vsqrt <- t(vs$v %*% (t(vs$u) * sqrt(vs$d))) 
       ans <- matrix(rnorm(n * p), nrow = n) %*% vsqrt
       ans <- sweep(ans, 2, mu, "+")
       dimnames(ans) <- list(NULL, dimnames(vmat)[[2]])  
       return(ans)
   }
   norms <- rmultnorm(100,c(100,100),matrix(c(14,11,11,14),nrow=2))
   x1 <- norms[,1]
   x2 <- norms[,2]
   plot(x1,x2,xlab="Test 1",ylab="Test 2")