Try:
y <- matrix(rnorm(n*10000), ncol=n) %*% chol(cov)
--
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.snow@intermountainmail.org
(801) 408-8111
> -----Original Message-----
> From: s-news-owner@lists.biostat.wustl.edu
> [mailto:s-news-owner@lists.biostat.wustl.edu] On Behalf Of Horace Tso
> Sent: Friday, February 03, 2006 12:19 PM
> To: s-news@lists.biostat.wustl.edu
> Subject: [S] Simulating multivariate normals
>
> 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
>
>
>
>
> --------------------------------------------------------------------
> This message was distributed by
> s-news@lists.biostat.wustl.edu. To unsubscribe send e-mail
> to s-news-request@lists.biostat.wustl.edu with the BODY of
> the message: unsubscribe s-news
>
|