Dear list,
This little question of mine has to do with simulating multivariate normals.
It's been nagging me for sometime and I'd like to hear the list's comment.
Splus has this build-in rmvnorm function (or the mvrnorm in MASS) which works
just fine, until I realize it should be exactly the same thing as taking the
choleski decomposition of the covariance matrix and multiply it with a sample
of N(0,I). So I did this,
y <- chol(cov) %*% matrix(rnorm(n*10000), nrow=n, ncol=10000)
where my covariance matrix is positive definite.
Now, my y should behave like what the cov says it should. Given a large enough
sample of unit normal variates, I should get var(y) ~ cov, ie. the covariance
matrix of my simulated sample should come close to cov.
However, trial after trial, my var(y) is so unrecognizable that I suspect
something is missing. On the other hand, with rmvnorm,
y1 <- rmvnorm(1000, cov=cov, d=n)
a modest n (~ 1000) already gives a very close approximation of cov.
What is the reason my homemade simulator doesn't do as well as rmvnorm?
Thank you all in advance.
Horace W. Tso
PGE
Portland, Oregon
|