[R] faster mvrnorm alternative

[Ben Bolker]

fantomas <tomas.iesmantas <at> gmail.com> writes:

> 
> 
> Hello,
> does anybody know another faster function for random multivariate normal
> variable simulation? I'm using mvrnorm, but as profiling shows, my algorithm
> spends approximately 50 % in executing mvrnorm function. 
> Maybe some of you knows much faster function for multivariate normal
> simulation?
> I would be very gratefull for advices.

  If you have to use a different variance-covariance matrix for every
iteration of the simulation then I think you may be stuck.  If, however,
you are simulating with the same Sigma many times, you may be better off
looking at the guts of the mvrnorm() function and seeing what can
be abstracted.

  However, upon experimenting I don't actually find that you
can save very much time this way (see below).  My next suggestion
would be that you see about getting some optimized linear algebra
libraries for your system (see e.g. the gcbd package on CRAN).
=================


m0 <- function (mu, Sigma, tol = 1e-06) {
  p <- length(mu)
  if (!all(dim(Sigma) == c(p, p))) 
    stop("incompatible arguments")
  eS <- eigen(Sigma, symmetric = TRUE, EISPACK = TRUE)
  ev <- eS$values
  if (!all(ev >= -tol * abs(ev[1L]))) 
    stop("'Sigma' is not positive definite")
  list(eSv=eS$vector,ev=ev, dd=diag(sqrt(pmax(ev, 0)), p))
}

m1 <- function(n,mu,eL,fancy=FALSE) {
  p <- length(mu)
  X <- matrix(rnorm(p * n), n)
  X <- drop(mu) + eL$eSv %*% eL$dd %*% t(X)
  if (fancy) {
    nm <- names(mu)
    if (is.null(nm) && !is.null(dn <- dimnames(Sigma))) 
      nm <- dn[[1L]]
    dimnames(X) <- list(nm, NULL)
    if (n == 1) 
      drop(X)
    else t(X)
  } else t(X)
}

mu <- rep(2,20)
Sigma <- matrix(0.1,nrow=20,ncol=20)
diag(Sigma) <- 4

set.seed(1001)
x1 <- mvrnorm(1000,mu=mu,Sigma=Sigma)
set.seed(1001)
e <- m0(mu,Sigma); x2 <- m1(1000,mu=mu,eL=e,fancy=TRUE)
identical(x1,x2)

system.time(replicate(1000,mvrnorm(1000,mu=mu,Sigma=Sigma)))
system.time({e <- m0(mu,Sigma); replicate(1000,m1(1000,mu=mu,eL=e))})
system.time({e <- m0(mu,Sigma); replicate(1000,m1(1000,mu=mu,eL=e,fancy=FALSE))})