I do have some trouble with matrices. I want to build up a covariance matrix
with a hierarchical structure). For instance, in dimension n=10, I have two
subgroups (called REGION).
NR=2; n=10
CORRELATION=matrix(c(0.4,-0.25,
-0.25,0.3),NR,NR)
REGION=sample(1:NR,size=n,replace=TRUE)
R1=REGION%*%t(rep(1,n))
R2=rep(1,n)%*%t(REGION)
SIGMA=matrix(NA,n,n)
for(i in 1:NR){
for(j in 1:NR){
SIGMA[(R1==i)&(R2==j)]=CORRELATION[i,j]
}}
If I run quickly some simulations, I build up the following matrix
> CORRELATION
[,1] [,2]
[1,] 0.40 -0.25
[2,] -0.25 0.30
> REGION
[1] 2 2 1 1 2 1 2 1 1 2
> SIGMA
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0.30 0.30 -0.25 -0.25 0.30 -0.25 0.30 -0.25 -0.25 0.30
[2,] 0.30 0.30 -0.25 -0.25 0.30 -0.25 0.30 -0.25 -0.25 0.30
[3,] -0.25 -0.25 0.40 0.40 -0.25 0.40 -0.25 0.40 0.40 -0.25
[4,] -0.25 -0.25 0.40 0.40 -0.25 0.40 -0.25 0.40 0.40 -0.25
[5,] 0.30 0.30 -0.25 -0.25 0.30 -0.25 0.30 -0.25 -0.25 0.30
[6,] -0.25 -0.25 0.40 0.40 -0.25 0.40 -0.25 0.40 0.40 -0.25
[7,] 0.30 0.30 -0.25 -0.25 0.30 -0.25 0.30 -0.25 -0.25 0.30
[8,] -0.25 -0.25 0.40 0.40 -0.25 0.40 -0.25 0.40 0.40 -0.25
[9,] -0.25 -0.25 0.40 0.40 -0.25 0.40 -0.25 0.40 0.40 -0.25
[10,] 0.30 0.30 -0.25 -0.25 0.30 -0.25 0.30 -0.25 -0.25 0.30
Here are eigenvalues (or at least the lower ones)
> min(eigen(SIGMA)$values)
[1] -2.109071e-17
> min(eigen(CORRELATION)$values)
[1] 0.09504902
theoretically, it should be 0 for SIGMA, but here it is not the case.
The problem here is that here, I cannot generate a Gaussian random vector
with that specific covariance matrix, since Cholesky does not work...
> library(mnormt)
> RES=rmnorm(n,mean=MU,varcov=SIGMA)
Error in chol.default(varcov) :
the leading minor of order 4 is not positive definite
any suggestions ? perhaps a better way of building my covariance matrix ? or
another way of generating a random vector with that specific covariance
matrix ?
Thanks for your help....
[[alternative HTML version deleted]]
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
