[R] R code for var-cov matrix given variances and correlations
Martin Maechler
maechler at stat.math.ethz.ch
Tue Dec 21 17:21:43 CET 2004
>>>>> "PD" == Peter Dalgaard <p.dalgaard at biostat.ku.dk>
>>>>> on 21 Dec 2004 14:34:59 +0100 writes:
PD> "Haynes, Maurice (NIH/NICHD)" <haynesm at cfr.nichd.nih.gov> writes:
>> If the vector of variances were
>> var.vec <- c(14, 12, 7)
>> and the vector of correlations were
>> cor.vec <- c(.4, .2, .5),
>> then the vector of covariances would be:
>> > covAB(c(.4, .2, .5),c(14, 14, 12), c(12, 7, 7))
>> [1] 5.184593 1.979899 4.582576
>> >
>> and the variance-covariance matrix with covariances rounded to
>> the first decimal place would be:
>> > vmat <- matrix(c(14, 5.2, 2.0, 5.2, 12, 4.6, 2.0, 4.6, 7),
>> + nrow=3)
>> > vmat
>> [,1] [,2] [,3]
>> [1,] 14.0 5.2 2.0
>> [2,] 5.2 12.0 4.6
>> [3,] 2.0 4.6 7.0
>> >
PD> # First fill in the correlation matrix:
PD> V <- matrix(NA,3,3)
PD> diag(V) <- 1
PD> V[lower.tri(V)] <- c(.4, .2, .5)
PD> V[upper.tri(V)] <- t(V)[upper.tri(V)]
PD> # then scale rows and columns
PD> D <- diag(sqrt( c(14, 12, 7)))
PD> D %*% V %*% D
PD> # or, more efficiently
PD> s <- sqrt( c(14, 12, 7))
PD> sweep(sweep(V,1,s,"*"),2,s,"*")
and others give similar answers, not as nice as the sweep one.
While I leave the exact comparisons and even more slick
alternatives to Gabor (:-),
yet another alternative is to look at
cov2cor()
which does Cov |--> Cor, i.e., the inverse of what the original
poster wanted.
The source of that function, e.g. in
https://svn.r-project.org/R/trunk/src/library/stats/R/cor.R
is
cov2cor <- function(V)
{
## Purpose: Covariance matrix |--> Correlation matrix -- efficiently
## ----------------------------------------------------------------------
## Arguments: V: a covariance matrix (i.e. symmetric and positive definite)
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 12 Jun 2003, 11:50
p <- (d <- dim(V))[1]
if(!is.numeric(V) || length(d) != 2 || p != d[2])
stop("`V' is not a square numeric matrix")
Is <- sqrt(1/diag(V)) # diag( 1/sigma_i )
if(any(!is.finite(Is)))
warning("diagonal has non-finite entries")
r <- V # keep dimnames
r[] <- Is * V * rep(Is, each = p)
## == D %*% V %*% D where D = diag(Is)
r[cbind(1:p,1:p)] <- 1 # exact in diagonal
r
}
[ Note that this is the *real source* including comments,
contrary to what you get when you just type 'cov2cor' in R ]
What I've wanted to expose here are the lines
Is <- sqrt(1/diag(V)) # diag( 1/sigma_i )
..
r <- V # keep dimnames
r[] <- Is * V * rep(Is, each = p)
which shows ``how to'' compute Diag %*% Matrix %*% Diag,
efficiently, probably a bit more even than Peter's nice double
use of sweep(.., *) {and still keeping dimnames}.
Martin
More information about the R-help
mailing list