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

[Jarrod Hadfield]

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.