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))})