Inverse Wishart: rIW, riwish, and Sorensen-Gianola
On Tue, May 3, 2011 at 2:36 AM, Jarrod Hadfield <j.hadfield at ed.ac.uk> wrote:
Hi, In the note section of ?rIW I have provided the relationships between the inverse wishart samplers in MCMCglmm/MCMCpack/bayesm and the notation used on Wikipedia. MCMCglmm is the odd one out becuase V*nu is the inverse scale matrix whereas in other samplers you pass the scale matrix (or its inverse as in riwish) directly. The reason for choosing to define it this way was because it is more intuitive for users when they want to fix (co)variances at some value because they can just specify a value in V. ? Regarding the conventions of quantitive geneticists I would go with what Sorensen and Gianola use.
Thanks Jarrod, It appears that there is some inconsistency in the notation used by some authors. For example, Sorensen-Gianola write T ~ IW_p(Sigma, n) and say E[T|Sigma,n] = Sigma^-1 / (n - p - 1). Other authors (B. Walsh, e.g.) write T ~ IW_n(Sigma^-1), where T = 1/W, and W ~ W_n(Sigma). So the confusion is over what is the parameter of the Inverse Wishart: Sigma or Sigma^-1? If the "Sigma" in the Sorensen-Gianola definition is replaced by Sigma^-1, then you get E[T|Sigma,n] = Sigma/(n - p - 1), the formula that appears in the MCMCpack man page for riwish. Dominick
riwish.Sigma <- 4*diag(2) riwish.nu <- 5 p<-dim(riwish.Sigma)[1] MCMCglmm.V <- riwish.Sigma/riwish.nu MCMCglmm.nu <-5
colMeans(rIW(MCMCglmm.V, MCMCglmm.nu, n=100000))
[1] ?1.995096455 -0.004143937 -0.004143937 ?2.009750258 riwish.Sigma /(riwish.nu ?- p - 1) ? ? [,1] [,2] [1,] ? ?2 ? ?0 [2,] ? ?0 ? ?2 Jarrod Quoting Dominick Samperi <djsamperi at gmail.com> on Tue, 3 May 2011 00:30:52 -0400:
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
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
-- The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.