Using R to simulate correlated data which is not standard-normally distributed

gistfile1.r

library(MASS)

## let's assume that
## X ~ N(1.0, 1.5)
## Y ~ N(2.0, 2.1)
## Z ~ N(0.5, 1.0)
## r_XY = 0.2, r_XZ = 0.4, r_YZ = -0.3

## create a vector of standard diviations
std <- c(1.5, 2.1, 1.0)

## the intendet correlation matrix
cormat <- c(1,    0.2,  0.4,
            0.2,  1,   -0.3,
            0.4, -0.3,  1)
cormat <- matrix(cormat, ncol = 3, byrow = TRUE)
cormat


## We know that for r_XY = S_XY / (S_X * S_Y)
## with: 
## S_XY: Covariance between X and Y
## S_X, S_Y: Standard diviation for X and Y, respectively 
## We know r_XY, but we need S_XY, hence: 
## S_XY = r_XY * (S_X * S_Y)

## generate a matrix of all products of S_X, S_Y, S_Z
tmp <- standv %*% t(standv)
tmp[2,1] ## S_X * S_Y = 1.5*2.1

covmat <- tmp * cormat
## the diagonal of the covariance matrix equals the variances of X, Y, Z
diag(covmat) <- standv^2
covmat


## simulate data
df <- mvrnorm(10000, mu = c(1, 2, 0.5), Sigma = covmat)

## estimate pm-corelations
cor(df)

## check means
summary(df)
## check standard deviations
sd(df)