[R] R code for var-cov matrix given variances and correlations

[Martin Maechler]

>>>>> "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