[S] Summary: Wishart r.v. generation

John Castelloe (jcaste@stat.uiowa.edu)
Wed, 08 Apr 1998 09:06:44 -0500


Earlier I asked if there is any code available to generate realizations
from a Wishart distribution. Thanks to Bill Venables, Erin Hodgess and
Coen Bernaards for help. Following is a summary of approaches:

The most straightforward reference I found was pp. 203-204 of Mark E.
Johnson (1987) "Multivariate Statistical Simulation" (which also . Other
good references are:

%AUTHOR = Odell, P. L.
%AUTHOR = Feiveson, A. H.
%TITLE = A numerical procedure to generate a sample covariance matrix
%JOURNAL = Journal of the American Statistical Association
%VOLUME = 61
%PAGES = 199-203
%YEAR = 1966
%KEYWORDS = Monte Carlo; Random number generation; Wishart distribution

%AUTHOR = Smith, W. B.
%AUTHOR = Hocking, R. R.
%TITLE = [Algorithm AS 53] Wishart variate generator
%JOURNAL = Applied Statistics
%VOLUME = 21
%PAGES = 341-345
%YEAR = 1972
%KEYWORDS = Bartlett decomposition; Random normal deviates

Bill Venables provided an elegant S-Plus function, which follows after his
description (it has not been thoroughly tested, but it appears to work fine):

"The following code should generate a single p x p Wishart matrix
realisation with df degrees of freedom. If SqrtSigma is not
specified the parent variance matrix is assumed to be the
identity matrix; if it is specified it must be some matrix
"square root" of the variance matrix, Sigma, e.g. chol(Sigma).
In other words crossprod(SqrtSigma) must be equal to Sigma. (As
this is a big-ticket computation I figured it was better to do it
once outside the function than inside every time.)
Limited testing suggests it is not completely stupid, at any
rate, but that's all. Caveat emptor."

--(cut here)--
rwishart <- function(df, p = nrow(SqrtSigma), SqrtSigma = diag(p)) {
if((Ident <- missing(SqrtSigma)) && missing(p))
stop("either p or SqrtSigma must be specified")
Z <- matrix(0, p, p)
diag(Z) <- sqrt(rchisq(p, df:(df-p+1)))
if(p > 1) {
pseq <- 1:(p-1)
Z[rep(p*pseq, pseq) + unlist(lapply(pseq, seq))] <- rnorm(p*(p-1)/2)
}
if(Ident)
crossprod(Z)
else
crossprod(Z %*% SqrtSigma)
}
--(cut here)--

---------------------------------------------------------------
| John Castelloe phone: (319)335-2009 |
| Computer RA fax: (319)335-3017 |
| Dept. of Statistics web: http://www.stat.uiowa.edu/~jcaste |
| University of Iowa office: 269 SH |
---------------------------------------------------------------
-----------------------------------------------------------------------
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