rmvnorm <- function(n, mean=0, cov=matrix(1,1,1)){
## generoi p*n matrisin, jonka sarakkeet ~N(mean, cov)
B <- chol(cov)
p <- ncol(cov)
Z <- matrix(rnorm(n*p), n, p)
t(Z%*%B) + matrix(mean, p, n)
}

X <- rmvnorm(1000, mean=1:3, cov = matrix(c(1,1,1, 1,2,2, 1,2,3), 3,3))