[R-sig-ME] Inverse Wishart: rIW, riwish, and Sorensen-Gianola

Dominick Samperi djsamperi at gmail.com
Tue May 3 06:30:52 CEST 2011


The conventions used for the Inverse-Wishart generator included
with MCMCglmm (rIW) and MCMCpack (riwish) seem to be different,
and both do not seem to agree with the notation used in the
Sorensen-Gianola book (Likelihood, Bayesian, and MCMC methods
in Quantitative Genetics). Can somebody say a few words about
the conventions used in Quantitative Genetics?

Here is a simple test:

Sigma <- 4*diag(2)

## Use riwsh from MCMCpack
x <- matrix(rep(0,4),2,2)
for(i in 1:1000) {
    x <- x + riwish(5,Sigma)
MCMCpack.mean <- x/1000

## Use rIW from MCMCglmm
x <- matrix(rep(0,4),2,2)
for(i in 1:1000) {
    x <- x + rIW(Sigma,5,1)
MCMCglmm.mean <- x/1000

list(MCMCpack.mean=MCMCpack.mean, MCMCglmm.mean=MCMCglmm.mean)

MCMCpack man page says the mean should be Sigma/(n - p - 1) = Sigma/(5
- 2 - 1), and
this is what we get when using riwish, but not when we use rIW.

Also, Sorensen-Gianola Eqn. (1.107) says the mean should be
Sigma^(-1)/(n-p-1), but
this is not what we get.


More information about the R-sig-mixed-models mailing list