Sampling from the inverse wishart distribution in R

Question

I am reading the following in the book "Bayesian Data Analysis"

In the book, the matrix S is defined as follows:

For the sake of the question, let's assume that I have S calculated. I would like to be able to sample from the distribution of Sigma|y, i.e., the inverse wishart distribution. I am using the MCMCpack package in R to do so and what I can't quite understand is whether the function expects me to pass it S or the inverse of S, i.e.,

library(MCMCpack)
S <- matrix(c(1,.3,.3,1),2,2)

## Use S?
riwish(4, S)

## or use S^{-1}
riwish(4, solve(S))

It is not clear to me from the documentation which one is expected.

Answer 1

I think the way to figure this out is to think carefully about whether each object/variable is on the variance scale or the precision (inverse-variance) scale. ?riwish says:

The mean of an inverse Wishart random variable with ‘v’ degrees of freedom and scale matrix ‘S’ is (v-p-1)^{-1} S.

We want the final answer (the posterior distribution of the covariance matrix) to be on the variance scale. If the mean of the random variable is is proportional to S, that must mean that S as input to riwish should also be on the covariance scale.

We can also confirm by brute force that using S gives us a sensible value (since the Jeffreys prior is weakly informative we should get a mean approximately equal to S):

S <- matrix(c(1,.3,.3,1),2,2)
set.seed(101)
replicate(10000, MCMCpack::riwish(4, S)) |> apply(c(1,2), mean)

         [,1]      [,2]
[1,] 1.009743 0.2656540
[2,] 0.265654 0.9415342

Sinv <- solve(S)
replicate(10000, MCMCpack::riwish(4, Sinv)) |> apply(c(1,2), mean)

           [,1]       [,2]
[1,]  1.0273954 -0.3092717
[2,] -0.3092717  1.0217894