Suppose V is a covariance matrix on which you want an inverse Wishart
distribution. Suppose also you want V to have _mean_ R and degrees of
freedom d. Then the _inverse_ of V has a Wishart distribution with the
following parameters: degrees of freedom = d, covariance matrix = inverse
of (d*R). So, just generate V-inverse and invert to obtain samples from an
inverse Wishart distribution. In the code below, you would set "df" to d
and "SqrtSigma" to chol(solve(d*R)), the Cholesky decomposition of the
inverse of d*R, and invert the result.
Now, here is Bill Venables code, along with his commentary:
"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)--
Hope this helps.
John
---------------------------------------------------------------
| 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