[R] generating random covariance matrices (with a uniform distribution of correlations)

Petr Savicky savicky at praha1.ff.cuni.cz
Sat Jun 4 11:44:52 CEST 2011


On Fri, Jun 03, 2011 at 01:54:33PM -0700, Ned Dochtermann wrote:
> Petr,
> This is the code I used for your suggestion:
> 
> 	k<-6;kk<-(k*(k-1))/2
> 	x<-matrix(0,5000,kk)
> 	for(i in 1:5000){
> 	A.1<-matrix(0,k,k)
> 	rs<-runif(kk,min=-1,max=1)
> 	A.1[lower.tri(A.1)]<-rs
> 	A.1[upper.tri(A.1)]<-t(A.1)[upper.tri(A.1)]
> 	cors.i<-diag(k)
> 	t<-.001-min(Re(eigen(A.1)$values))
> 	new.cor<-cov2cor(A.1+(t*cors.i))
> 	x[i,]<-new.cor[lower.tri(new.cor)]}
> 	hist(c(x)); max(c(x)); median(c(x))
> 
> This, unfortunately, does not maintain the desired distribution of
> correlations.

Hello.

On the contrary to what i thought originally, there are solutions
also for the case of the correlation matrix. The first solution
creates a singular correlation matrix (of rank 3), but the nondiagonal
entries have exactly the uniform distribution on [-1, 1], since
the scalar product of two independent uniformly distributed unit
vectors in R^3 has the uniform distribution on [-1, 1].

  x <- matrix(rnorm(18), nrow=6, ncol=3)
  x <- x/sqrt(rowSums(x^2))
  a <- x %*% t(x)

The next solution produces a correlation matrix of full rank, whose
non-diagonal entries have distribution very close to the uniform on
[-1, 1]. KS test finds a difference only with sample size more
than 50'000.

  w <- c(0.01459422, 0.01830718, 0.04066405, 0.50148488, 0.60330865, 0.61832829)
  x <- matrix(rnorm(36), nrow=6, ncol=6) %*% diag(w)
  x <- x/sqrt(rowSums(x^2))
  a <- x %*% t(x)

Hope this helps.

Petr Savicky.



More information about the R-help mailing list