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

[Dominick Samperi]

On Tue, May 10, 2011 at 1:04 AM, Jarrod Hadfield <j.hadfield at ed.ac.uk> wrote:
> 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.

Thanks, I should have seen that. There is another small mistake in my code: I
need to multiply by (nu-p-1) = (5-2-1), not divide, to calculate the
mean correctly.

On my earlier comments about an inconsistency between Sorensen-Gianola
and B. Walsh, I realize now that the inverse Wishart is naturally
parametrized in
terms of the precision matrix (inverse covariance matrix), and the expressions
in Sorensen-Gianola are consistent with this convention.

Dominick

> 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.
>
>
>