[R] how to generate a random correlation matrix with restrictions
Ravi Varadhan
rvaradhan at jhmi.edu
Sun Aug 30 02:54:57 CEST 2009
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.
> >
More information about the R-help
mailing list