[R] generating correlated random variables from different distributions

Richard and Barbara Males rbmales at gmail.com
Thu Apr 29 17:17:50 CEST 2010


I need to generate a set of correlated random variables for a Monte
Carlo simulation.   The solutions I have found
(http://www.stat.uiuc.edu/stat428/cndata.html,
http://www.sitmo.com/doc/Generating_Correlated_Random_Numbers), using
Cholesky Decomposition, seem to work only if the variables come from
the same distribution with the same parameters.  My situation is that
each variable may be described by a different distribution (or
different parameters of the same distribution).  This approach does
not seem to work, see code and results below.  Am I missing something
here?  My math/statistics is not very good, will I need to generate
correlated uniform random variables on (0,1) and then use the inverse
distributions to get the desired results I am looking for?  That is
acceptable, but I would prefer to just generate the individual
distributions and then correlate them.  Any advice much appreciated.
Thanks in advance

R. Males
Cincinnati, Ohio, USA

Sample Code:
# Testing Correlated Random Variables

# reference http://www.sitmo.com/doc/Generating_Correlated_Random_Numbers
# reference http://www.stat.uiuc.edu/stat428/cndata.html
# create the correlation matrix
corMat=matrix(c(1,0.6,0.3,0.6,1,0.5,0.3,0.5,1),3,3)
cholMat=chol(corMat)
# create the matrix of random variables
set.seed(1000)
nValues=10000

# generate some random values

matNormalAllSame=cbind(rnorm(nValues),rnorm(nValues),rnorm(nValues))
matNormalDifferent=cbind(rnorm(nValues,1,1.5),rnorm(nValues,2,0.5),rnorm(nValues,6,1.8))
matUniformAllSame=cbind(runif(nValues),runif(nValues),runif(nValues))
matUniformDifferent=cbind(runif(nValues,1,1.5),runif(nValues,2,3.5),runif(nValues,6,10.8))

# bind to a matrix
print("correlation Matrix")
print(corMat)
print("Cholesky Decomposition")
print (cholMat)

# test same normal

resultMatNormalAllSame=matNormalAllSame%*%cholMat
print("correlation matNormalAllSame")
print(cor(resultMatNormalAllSame))

# test different normal

resultMatNormalDifferent=matNormalDifferent%*%cholMat
print("correlation matNormalDifferent")
print(cor(resultMatNormalDifferent))

# test same uniform
resultMatUniformAllSame=matUniformAllSame%*%cholMat
print("correlation matUniformAllSame")
print(cor(resultMatUniformAllSame))

# test different uniform
resultMatUniformDifferent=matUniformDifferent%*%cholMat
print("correlation matUniformDifferent")
print(cor(resultMatUniformDifferent))

and results

[1] "correlation Matrix"
     [,1] [,2] [,3]
[1,]  1.0  0.6  0.3
[2,]  0.6  1.0  0.5
[3,]  0.3  0.5  1.0
[1] "Cholesky Decomposition"
     [,1] [,2]      [,3]
[1,]    1  0.6 0.3000000
[2,]    0  0.8 0.4000000
[3,]    0  0.0 0.8660254
[1] "correlation matNormalAllSame" <== ok
          [,1]      [,2]      [,3]
[1,] 1.0000000 0.6036468 0.3013823
[2,] 0.6036468 1.0000000 0.5005440
[3,] 0.3013823 0.5005440 1.0000000
[1] "correlation matNormalDifferent" <== no good
          [,1]      [,2]      [,3]
[1,] 1.0000000 0.9141472 0.2676162
[2,] 0.9141472 1.0000000 0.2959178
[3,] 0.2676162 0.2959178 1.0000000
[1] "correlation matUniformAllSame" <== ok
          [,1]      [,2]      [,3]
[1,] 1.0000000 0.5971519 0.2959195
[2,] 0.5971519 1.0000000 0.5011267
[3,] 0.2959195 0.5011267 1.0000000
[1] "correlation matUniformDifferent" <== no good
          [,1]      [,2]      [,3]
[1,] 1.0000000 0.2312000 0.0351460
[2,] 0.2312000 1.0000000 0.1526293
[3,] 0.0351460 0.1526293 1.0000000
>



More information about the R-help mailing list