Plotting the marginal prior distributions for the elements of an Inverse-Wishart matrix

ggplot-marginals_Wishart.R

### Plotting the marginal prior distributions for the elements of
### an Inverse-Wishart matrix
library(MCMCpack)   # Wishart distribution
library(actuar)     # inverse gamma distribution
library(ggplot2)

## Hyperparameters of the inverse-Wishart prior
Phi <- matrix(c(1,.3,.3,2),2,2)   # Scale matrix
nu  <- 5                          # Degrees of Freedom

## Marginal priors of the variances
##
## Use the fact that, if S ~ invW(Phi, nu), then
## s_ii ~ invGa(alpha = nu/2, beta = Phi_ii/2)

limits <- sapply(diag(Phi), 
                 function(x) qinvgamma(c(.01, .99), shape = nu/2, scale = x/2))

## We can choose wether to plot the marginal prior variances under
## the same scale or not. If that's the case, we use the maximal interval
limits <- c(min(limits[1, ]), max(limits[2, ]))
x <- seq(limits[1], limits[2], length = 501)
y <- do.call('c',
             lapply(diag(Phi),
                    function(ss) dinvgamma(x,
                                           shape = nu/2,
                                           scale = ss/2)))
priordf.var <- data.frame(comp = factor(c(rep('sigma2_a', length(x)),
                                          rep('sigma2_c', length(x)))),
                          x = c(x, x),
                          y)


## Plot prior curves and true value
ggplot(priordf.var, aes(x, y, col = comp)) + 
  geom_line() +
  theme_bw()


## For the marginal correlation, the analytic derivation of the distribution
## is not very straightforward, so we just sample and plot an empirical density

## Use the fact that A ~ invW(Phi, nu) iff invA ~ W(invPhi, nu)
## and build a function to get a sample of the correlation or covariance element
rinvWish <- function(n, Phi, nu, corr = TRUE) {
  invsample <- lapply(seq.int(n),
                   function(x) rwish(nu, solve(Phi)))
  sample <- lapply(invsample, solve)
  cov <- sapply(sample,
                function(x) x[1,2])
  if( corr ) {
    denom <- apply(sqrt(sapply(sample, diag)), 2, prod)
    stopifnot(abs(cov) < denom)
    cov <- cov/denom
  }
  
  return(cov)
}

N <- 1e4
cor.smpl <- rinvWish(N, Phi, nu, corr = TRUE)
ggplot(as.data.frame(density(cor.smpl, from = -1, to = 1)[1:2]),
       aes(x, y)) + 
  geom_line()