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

[Jarrod Hadfield]

Hi,

The third call to rIW is fix which fixes a (sub) matrix of the matrix  
to be sampled. Replace rIW(Sigma,nu,1) with rIW(Sigma,nu,n=1) in your  
code and you get the desired result.

Jarrod



Quoting Dominick Samperi <djsamperi at gmail.com> on Mon, 9 May 2011  
22:21:57 -0400:

> 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 for the pointer. I'm still confused though, because the
> following numbers computed
> by MCMCglmm, MCMCpack, and bayesm should be the same according to your
> man page. But the MCMCglmm number does not agree with the others.
>
> nu <- 5
> 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(nu,Sigma*nu)
> }
> MCMCpack.mean <- x/1000/(5-2-1)
>
> ## Use rIW from MCMCglmm
> require(MCMCglmm)
> x <- matrix(rep(0,4),2,2)
> for(i in 1:1000) {
>     x <- x + rIW(Sigma,nu,1)
> }
> MCMCglmm.mean <- x/1000/(5-2-1)
>
> ## Use rwishart from bayesm
> require(bayesm)
> x <- matrix(rep(0,4),2,2)
> for(i in 1:1000) {
>     x <- x + rwishart(nu, solve(Sigma*nu))$IW
> }
> bayesm.mean <- x/1000/(5-2-1)
>
> list(MCMCpack.mean=MCMCpack.mean, MCMCglmm.mean=MCMCglmm.mean,
>      bayesm.mean=bayesm.mean)
>
>> 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.
>>
>>
>>
>
>



-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.