Re: [S] generate multivariate normal data

Prof Brian Ripley (ripley@stats.ox.ac.uk)
Fri, 13 Feb 1998 07:55:26 +0000 (GMT)


Georges Monette wrote:
>
>
> Here's how to create a matrix whose rows are independent
> multivariate normal vectors with mean 0 and variance matrix SIG:
>
> > mat <- matrix(rnorm(n*p),ncol=p) %*% chol(SIG)
>
> where n is the number of observations and SIG is p x p.
>
> If you want a mean of MU:
>
> > mat <- t(t(mat) + MU)
>
This will come to grief if SIG is singular or near-singular.
Adding pivoting to the chol() call would help, but also
complicates the code considerably. My preference is to use the
eigendecomposition, which allows a numerically stable algorithm.

> >
> > Do anybody know how to generate multivariate normal data from the
> > correlation matrix?
> > Is there any standard Splus program can do it?
> > Thank you for any comment.
> >
> > David Wang

In S-PLUS 4.0 there is the function rmvnorm to do this (although for the
taks to be well-defined you need a covariance matrix, not just a
correlation matrix). There is also a function mvrnorm in the V&R2 MASS
library, which is

# file MASS/mvrnorm.q
# copyright (C) 1994-8 W. N. Venables and B. D. Ripley
#
mvrnorm <- function(n=1, mu, Sigma, tol=1e-6)
{
p <- length(mu)
if(!all(dim(Sigma) == c(p,p))) stop("incompatible arguments")
eS <- eigen(Sigma, sym = T)
ev <- eS$values
if(!all(ev >= -tol*abs(ev[1]))) stop("Sigma is not positive definite")
X <- mu + eS$vectors %*% diag(sqrt(pmax(ev, 0))) %*% matrix(rnorm(p*n),p)
nm <- names(mu)
if(is.null(nm) && !is.null(dn <- dimnames(Sigma))) nm <- dn[[1]]
dimnames(X) <- list(nm, NULL)
if(n == 1) drop(X) else t(X)
}

As you will see, a lot of the code is devoted to niceties such as adding
variable names, treating the case n=1 and dealing with singularities.

(By the way, the copyright notice is there to help a particular group
of users: the corresponding licence (in the distribution) is not at all
restrictive.)

-- 
Brian D. Ripley,                  ripley@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu.  To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message:  unsubscribe s-news