Simulating random variables with a particular correlation structure

I often want to simulate multivariate random variables with a particular correlation structure (example), so I developed a flexible R package for doing this called rmv. The name rmv is based on R‘s convention for naming random number generators: paste('r', distribution.name, sep = ''), where in this case distribution.name <- 'mv' for generic multivariate.

I think the interface is easy to use:

rmv(n, covmat, rfunc = rnorm,
    method = c("chol", "eigen"), ...)

n is just the number of random vectors. covmat is just a covariance matrix. rfunc is the name of a function that returns a vector of univariate random variables. The only real restrictions on the choice of rfunc is that its first argument must accept the number of random draws and it must produce independent and identically distributed draws. The ... passes arguments to rfunc. And I’ll explain the method argument after.

For most practical purposes…the resulting random variables will have the same correlation structure as random variables with a covariance matrix given by covmat multiplied by the variance of the random numbers produced by rfunc.

Here is why it works:

Finally, the method argument just indicates how to compute the matrix square root of covmat, which is required by the algorithm (see the image above). Different method specifications will create different results but (almost…I think) nothing will change the fact that is in bold above. In any case it has always worked for me!