### ### Random number generation from a bivariate normal distribution ### ### AUTHOR: Arnost Komarek ### ### LOG: 20130317 created ### ### =============================================================================== rMVN2 <- function(n, mean = c(0, 0), sd = c(1, 1), rho = 0) { if (any(sd <= 0)) stop("all sd's must be positive") if (rho <= -1 | rho >= 1) stop("rho must be strictly between -1 and 1") ## Covariance matrix Sigma <- matrix(c(sd[1]^2, rho * sd[1] * sd[2], rho * sd[1] * sd[2], sd[2]^2), nrow = 2, ncol = 2) ## Cholesky decomposition of a covariance matrix ## Sigma = t(SigmaHalf) %*% SigmaHalf SigmaHalf <- chol(Sigma) ## Random sample from N(0, I_2) z <- matrix(rnorm(2 * n), nrow = n, ncol = 2) ## t(y[i,]) = mean + t(SigmaHalf) %*% t(z[i,]) ## y[i,] = t(mean) + z[i,] %*% SigmaHalf y <- matrix(mean, nrow = n, ncol = 2, byrow = TRUE) + z %*% SigmaHalf return(y) }