[R] Warning - Naive Question Alert
Prof Brian D Ripley
ripley at stats.ox.ac.uk
Sun Aug 27 08:00:12 CEST 2000
On Sun, 27 Aug 2000, Bill Venables wrote:
> Extending this to the multivariate case should not be all that difficult.
> One way you could do it is generate any old multivariate normal sample of
> the right dimension, "standardize" it by subtracting the mean and
> multiplying by the inverse of a symmetric square root of the variance
> matrix, then multiply by the same sort of square root of the given summary
> statistic variance matrix and add back on the given sample mean. [Warning:
> I am doing this at the keyboard, ideally I would like a bit more time to
> think about it...but I don't have that sort of time right now. Caveat
> emptor!]
Here's the code I promised, simplified:
mvrnorm <- function (n = 1, mu, Sigma, tol = 1e-06, empirical=FALSE)
{
p <- length(mu)
if (!all(dim(Sigma) == c(p, p))) stop("incompatible arguments")
eS <- eigen(Sigma, sym = TRUE)
ev <- eS$values
if (!all(ev >= -tol * abs(ev[1]))) stop("Sigma is not positive definite")
X <- matrix(rnorm(p * n), n)
if(empirical) {
X <- scale(X, T, F) # remove means
X <- X %*% svd(X, nu = 0)$v # rotate to PCs
X <- scale(X, F, T) # rescale PCs to unit variance
}
X <- mu + eS$vectors %*% diag(sqrt(pmax(ev, 0))) %*% t(X)
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)
}
With empirical=TRUE you get the exact mu and Sigma in the sample.
--
Brian D. Ripley, ripley at 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
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list