[R] simulate data from multivariate normal with pre-specified correlation matrix

[Thomas Harte]

# here's your example correlation matrix:
sigma<- matrix(c(1.00, 0.75, 0,
		 0.75, 1.00, 0,
		 0.00, 0.00, 0), nr=3, byrow=TRUE)
chol(sigma)
# Error in chol(sigma) : the leading minor of order 3 is not positive definite

# DUH!

# let's chop off that dangling row and column:
sigma<- sigma[1:2, 1:2]
N<- 1000
set.seed(1)
samp<- matrix(rnorm(N*nrow(sigma)), N, nrow(sigma))%*% chol(sigma)
# check if the empirical correlation is close to the theoretical sigma:
cor(samp)
#          [,1]      [,2]
#[1,] 1.0000000 0.7502607
#[2,] 0.7502607 1.0000000

# not bad!



>   Message: 50
>   Date: Mon, 13 Aug 2007 17:15:12 -0400
>   From: Jiao Yang <yj316 at gwu.edu>
>   Subject: [R] simulate data from multivariate normal with pre-specified
>   	correlation matrix
>   To: r-help at stat.math.ethz.ch
>   Message-ID: <de6a933c4844.46c091a0 at gwu.edu>
>   Content-Type: text/plain; charset=us-ascii

>   For example, the correlation matrix is 3x3 and looks like 

>   1       0.75  0  0  0
>   0.75    1     0  0  0
>    0       0     0   0  0

>   Can I write the code like  this?

>   p<- 3 >  number of variables per observation
>   N<- 10 >  number of samples

>   >  define population correlation matrix sigma
>   sigma<-matrix(0,p,p) > creates a px p matrix of 0
>   rank<-2
>   for (i in 1:rank){
>   for (j in 1:rank){
>   rho<-0.75
>    sigma[i,j]<-rho^abs(i-j)
>   sigma[i,i]<-1
>    }
>   }

>   >  create a data set of 10 observations from multivariate normalcovariance matrix
sigma  
>   library(MASS)
>   Xprime<-mvrnorm(N,rep(1,p), sigma )

>   Also, if I calculate

>   S<- cov(Xprime) 

>   will S be the sample correlation matrix?

>   Thank you very much for your time!