9
$\begingroup$

I am trying to figure out how to convert a correlation matrix (R) to a covariance matrix (S) for input into a random number generator that only accepts S (rmvnorm("mvtnorm") in R) 

library("mvtnorm") 

TRUTH= 0.8 # target correlation value between X1 and X2
R <- as.matrix(data.frame(c(1, TRUTH), c(TRUTH, 1)))
V <- diag(c(sqrt(1), sqrt(1))) # diagonal matrix of sqrt(variances)
S <- V %*% R %*% V
cor(rmvnorm(100, sigma=S) )

# repeat this to get an idea of the variance around Pearson's estimator

Instance where variances are not equal to 1

V <- diag(c(sqrt(3), sqrt(2))) 
S <- V %*% R %*% V
cor(rmvnorm(100, sigma=S) )

This seems to be correct, but I would like expert criticism.

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ I can't read your code (I don't know R), but conversion of corr to cov matrix or back is switching an SSCP-type matrix to a new diagonal while preserving the cosine. NewMatrix = Matrix &* (coef*t(coef)) where coef = sqrt(NewDiagonal/Diagonal), * is vector multiplication and &* is usual, elementwise multiplication. $\endgroup$
    – ttnphns
    Commented Jun 28, 2013 at 17:32
  • 2
    $\begingroup$ The code looks ok but its comments don't. V had better be the diagonal matrix of square roots of variances for this to work. $\endgroup$
    – whuber ♦
    Commented Jun 28, 2013 at 17:44
  • 1
    $\begingroup$ Another way to accomplish what you want--and in some cases it might be numerically a little better--is to generate your data using the correlation matrix and post-multiply them by V. It should be obvious that this works because (1) separate linear transformations in the variables do not change their correlation but (2) rescaling a unit variance variable by a constant scales its variance by the square of that constant. Then if you look at what mvtnorm does behind the scenes to factor R, you can see how it effectively carries out the same post-multiplication by V. $\endgroup$
    – whuber ♦
    Commented Jun 28, 2013 at 18:10
  • $\begingroup$ Which can be seen in the 1st example where R = S, because the standard deviations are both 1. So, this has the potential to be different? I am not mathematically inclined, so I will have to "prove" this to myself via coding. Thanks again whuber. $\endgroup$
    – Patrick
    Commented Jun 28, 2013 at 18:40
  • $\begingroup$ @Patrick: Note that your formula V %*% R %*% V is equivalent to what I've suggested above. But, I predict, your formula with two matrix multiplications is slower (which must show on big matrices). Elementwise multiplication of matrices (is there such an operation in R?) is faster, AFAIK. $\endgroup$
    – ttnphns
    Commented Jun 28, 2013 at 19:29

1 Answer 1

Reset to default
10
$\begingroup$

Let $R$ be the correlation matrix and $S$ the vector of standard deviations, so that $S\cdot S$ (where $\cdot$ is the componentwise product) is the vector of variances. Then $$ \text{diag}(S) R \text{diag}(S) $$ is the covariance matrix. This is fully explained here.

This can be implemented in R as

cor2cov_1 <- function(R,S){
    diag(S) %*% R %*% diag(S)
}

but is inefficient. An efficient implementation is

cor2cov <- function(R, S) {
 sweep(sweep(R, 1, S, "*"), 2, S, "*")
 }

and you can test yourself they give the same result.

TRUTH= 0.8 
R <- as.matrix(data.frame(c(1, TRUTH), c(TRUTH, 1)))
S = c(sqrt(1), sqrt(1))

cor2cov_1(R,S)

outer(S,S) * R 

smat = as.matrix(S)
R * smat %*% t(smat)

Here is a microbenchmark showing the efficiency of the functions:

library(microbenchmark)
microbenchmark::microbenchmark(outer(S,S) * R ,cor2cov_1(R,S), cor2cov(R,S), R * smat %*% t(smat), times = 10000)

Unit: microseconds
                 expr     min      lq       mean  median      uq      max neval cld
      outer(S, S) * R   1.968   2.214   2.724639   2.337   2.460 3611.362 10000  a 
      cor2cov_1(R, S)   1.722   1.886   2.778045   1.968   2.091 3743.259 10000  a 
        cor2cov(R, S) 113.037 116.071 125.844711 118.039 120.663 5462.020 10000   b
 R * smat %*% t(smat)   1.066   1.230   1.422712   1.435   1.517   12.177 10000  a
Improve this answer
$\endgroup$
5
  • 4
    $\begingroup$ also outer(S,S) * R since * performs the Hadamard (elementwise) product ... $\endgroup$
    – Ben Bolker
    Commented May 12, 2019 at 14:13
  • $\begingroup$ @Ben Bolker: It would be a good exercise to compare speedwise the efficiencies ... $\endgroup$
    – kjetil b halvorsen ♦
    Commented May 12, 2019 at 15:00
  • 5
    $\begingroup$ On my system (Microsoft's version of R), the outer solution is far faster for large dimensions. It is about eight times faster than the sweep implementation. (If you were to replace the inner sweep by R*S it would be almost twice as fast, but still four times slower than outer.) Replacing outer(S,S) by d <- length(S); (matrix(S,d,1) %*% matrix(S,1,d)) is a little faster. cc @BenBolker $\endgroup$
    – whuber ♦
    Commented May 12, 2019 at 15:48
  • 1
    $\begingroup$ What is the idea behind R * smat %*% t(smat)? $\endgroup$
    – M. Beausoleil
    Commented Jun 8, 2022 at 16:06
  • 1
    $\begingroup$ @M.Beausoleil, it's a fast way of doing the "sweep". Consider that you are multiplying each entry in the correlation matrix by the product of the component standard deviations. For example, the covariance matrix entry at [1, 2] will be the SD[1] * SD[2] * cor[1, 2]. By turning the SD vector into a matrix and multiplying it by itself (tcrossprod(S) would be even faster), you are "pregenerating" the needed 4x4 matrix of multiplications of SD entries each by each other. All that is left to do is to scale them, elementwise, by the correlations, thus the * R. Which could be first or second. $\endgroup$
    – Avraham
    Commented Feb 26 at 9:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.