[Multivariate Normal Sampling]

rmultinorm.R

library(assert)

rmultinorm.fast <- function(k, mu, sigma) {
  p = length(mu)
  t(chol(sigma)) %*% matrix(rnorm(p*k, 0, 1), nrow=p) + mu
}

#' Sample from a multivariate normal distribution. Attaches both `rmultinorm`
#' and `rmultinorm.fast` to the current environment.
#' 
#' @param k Number of samples
#' @param mu Mean vector of the distribution
#' @param sigma Covariance matrix
#' @return A matrix (of dimension length(mu) x k) with columns corresponding to 
#'   independent draws from the multivariate normal distribution.
#'  
#' @examples
#' mu = c(-10, 10)
#' sigma = matrix(c(1,0,0,1), nrow=2)
#' r = rmultinorm(1000, mu, sigma)
#' hist(c(1,1) %*% r)  
rmultinorm <-  function(k, mu, sigma) {
  assert(is.numeric(k),
          length(k) == 1,
          k > 0,
          is.numeric(mu),
          is.matrix(sigma),
          all(length(mu) == dim(sigma)))

  return(rmultinorm.fast(k, mu, sigma))
}