[R-sig-ME] Inverse Wishart: rIW, riwish, and Sorensen-Gianola
Dominick Samperi
djsamperi at gmail.com
Tue May 3 06:30:52 CEST 2011
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
More information about the R-sig-mixed-models
mailing list