R devel - May 2011 - matrix not positive definite (while it should be)

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

On Thu, May 05, 2011 at 02:31:59PM -0400, Arthur Charpentier
wrote:
> 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
Hi.

If X is a random vector from the 2 dimensional normal distribution
with the covariance matrix

        [,1]  [,2]
  [1,]  0.40 -0.25
  [2,] -0.25  0.30

then the vector X[REGION], which consists of replicated components
of X, has the expanded covariance matrix n times n, which you ask
for. Since the mean and the covariance matrix determine the distribution
uniquely, this is also a description of the required distribution.

The distribution is concentrated in a 2 dimensional subspace, since
the covariance matrix has rank 2.

Hope this helps.

Petr Savicky.