[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 Jarrod,

It appears that there is some inconsistency in the notation used by
some authors.
For example, Sorensen-Gianola write T ~ IW_p(Sigma, n) and say
E[T|Sigma,n] = Sigma^-1 / (n - p - 1).

Other authors (B. Walsh, e.g.) write T ~ IW_n(Sigma^-1), where
T = 1/W, and W ~ W_n(Sigma). So the confusion is over what is the
parameter of the Inverse Wishart: Sigma or Sigma^-1?

If the "Sigma" in the Sorensen-Gianola definition is replaced by Sigma^-1,
then you get E[T|Sigma,n] = Sigma/(n - p - 1), the formula that appears
in the MCMCpack man page for riwish.

Dominick

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