[R-sig-ME] Inverse Wishart: rIW, riwish, and Sorensen-Gianola
Jarrod Hadfield
j.hadfield at ed.ac.uk
Tue May 3 08:36:02 CEST 2011
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.
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.
More information about the R-sig-mixed-models
mailing list