Inverse Wishart: rIW, riwish, and Sorensen-Gianola

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

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