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

[Dominick Samperi]

Hello,

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
require(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
require(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.

Thanks,
Dominick