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

[Dominick Samperi]

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