[R] how to generate a random correlation matrix with restrictions

[Ravi Varadhan]

Yes, Kingsford.  This problem does not appear to be trival.  I am not sure if there is any literature on this.  I have seen a paper by Davies and Higham (BIT 2000) on "Stable generation of correlation matrices.  There is also a paper by Harry Joe on generating correlation matrices with given partial correlation structure.  I am not sure if these papers would be helpful for this problem.  Of course, one can readily cook up an ad-hoc, trial-error procedure to generate the type of matrices that Ning wants.

By the way, the easiest way to generate a PD matrix is:

N <- 10
amat <- matrix(rnorm(N*N), N, N)
A <- amat %*% t(amat)  # A is PD 

Best,
Ravi

____________________________________________________________________

Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvaradhan at jhmi.edu


----- Original Message -----
From: Kingsford Jones <kingsfordjones at gmail.com>
Date: Saturday, August 29, 2009 7:43 pm
Subject: Re: [R] how to generate a random correlation matrix with restrictions
To: Ravi Varadhan <rvaradhan at jhmi.edu>
Cc: Ning Ma <pningma at gmail.com>, r-help at r-project.org


> Thanks Ravi -- with my limited linear algebra skills I was assuming
>  invertible symmetric was sufficient (rather than just necessary) for
>  positive definiteness.  So, the open challenge is to generate a pd
>  matrix of dimension 100 with r_ij = 1 if i=j else -1< r_ij< 1, with
>  about 10% of the elements >.9.
>  
>  I tried for a bit without luck, but I can offer another way to
>  generate a pd matrix:
>  
>  ev <- runif(100) # choose some positive eigenvalues
>  U <- svd(matrix(runif(10000), nc=100))$u  # an orthogonal matrix
>  pdM <- U %*% diag(ev) %*% t(U)
>  
>  
>  Kingsford
>  
>  
>  
>  On Fri, Aug 28, 2009 at 9:41 PM, Ravi Varadhan<rvaradhan at jhmi.edu> wrote:
>  > Hi Kingsford,
>  >
>  > There is more structure to a correlation matrix than that meets the 
> eye!  Your method will produce a matrix R that looks "like" a 
> correlation matrix, but beware - it is an impostor!
>  >
>  > You can obtain a valid correlation matrix, Q, from the impostor R 
> by using the `nearPD' function in the "Matrix" package, which finds 
> the positive definite matrix Q that is "nearest" to R.  However, note 
> that when R is far from a positive-definite matrix, this step may give 
> a Q that does not have the desired property.
>  >
>  > require(Matrix)
>  >
>  > R <- matrix(runif(16), ncol=4)
>  >
>  > R <- (R * lower.tri(R)) + t(R * lower.tri(R))
>  >
>  > diag(R) <- 1
>  >
>  > eigen(R)$val
>  >
>  > Q <- nearPD(R, posd.tol=1.e-04)$mat
>  >
>  > eigen(Q)$val
>  >
>  > max(abs(Q - R))  # maximum discrepancy between R and Q
>  >
>  > Another easy way to produce a valid correlation matrix is:
>  >
>  > R <- matrix(runif(36), ncol=6)
>  >
>  > RtR <- R %*% t(R)
>  >
>  > Q <- cov2cor(RtR)
>  >
>  > But this does not have the property that the correlations are 
> uniformly distributed.
>  >
>  > Hope this helps,
>  > Ravi.
>  >
>  > ____________________________________________________________________
>  >
>  > Ravi Varadhan, Ph.D.
>  > Assistant Professor,
>  > Division of Geriatric Medicine and Gerontology
>  > School of Medicine
>  > Johns Hopkins University
>  >
>  > Ph. (410) 502-2619
>  > email: rvaradhan at jhmi.edu
>  >
>  >
>  > ----- Original Message -----
>  > From: Kingsford Jones <kingsfordjones at gmail.com>
>  > Date: Friday, August 28, 2009 10:12 pm
>  > Subject: Re: [R] how to generate a random correlation matrix with   
>     restrictions
>  > To: Ning Ma <pningma at gmail.com>
>  > Cc: r-help at r-project.org
>  >
>  >
>  >> Ahh -- Mark Leeds just pointed out off-list that it was supposed 
> to be
>  >> a correlation matrix.
>  >>
>  >> Perhaps
>  >>
>  >> R <- matrix(runif(10000), ncol=100)
>  >> R <- (R * lower.tri(R)) + t(R * lower.tri(R))
>  >> diag(R) <- 1
>  >>
>  >> These are of course uniformly distributed positive correlations.
>  >>
>  >> Kingsford
>  >>
>  >>
>  >> On Fri, Aug 28, 2009 at 7:36 PM, Kingsford
>  >> Jones<kingsfordjones at gmail.com> wrote:
>  >> >  R <- matrix(runif(10000), ncol=100)
>  >> >
>  >> > hth,
>  >> >
>  >> > Kingsford Jones
>  >> >
>  >> > On Fri, Aug 28, 2009 at 7:26 PM, Ning Ma<pningma at gmail.com> wrote:
>  >> >> Hi,
>  >> >>
>  >> >> How can I generate a random 100x100 correlation matrix, R={r_ij},
>  >> >> where about 10% of r_ij are greater than 0.9
>  >> >>
>  >> >> Thanks in advance.
>  >> >>
>  >> >> ______________________________________________
>  >> >> R-help at r-project.org mailing list
>  >> >>
>  >> >> PLEASE do read the posting guide
>  >> >> and provide commented, minimal, self-contained, reproducible code.
>  >> >>
>  >> >
>  >>
>  >> ______________________________________________
>  >> R-help at r-project.org mailing list
>  >>
>  >> PLEASE do read the posting guide
>  >> and provide commented, minimal, self-contained, reproducible code.
>  >