In R, I would use the MASS package: library(MASS) mu <- c(0,0) # Mean Sigma <- matrix(c(1, .5, .5, 1), 2) # Covariance matrix bivn <- mvrnorm(5000, mu = mu, Sigma = Sigma ) # from Mass package