Covariance matrix proposal distribution

Question

In a MCMC implementation of hierarchical models, with normal random effects and a Wishart prior for their covariance matrix, Gibbs sampling is typically used.

However, if we change the distribution of the random effects (e.g., to Student's-t or another one), the conjugacy is lost. In this case, what would be a suitable (i.e., easily tunable) proposal distribution for the covariance matrix of the random effects in a Metropolis-Hastings algorithm, and what should be the target acceptance rate, again 0.234?

Thanks in advance for any pointers.

Answer 1

Well, if you are looking "for any pointers"...

The (scaled)(inverse)Wishart distribution is often used because it is conjugate to the multivariate likelihood function and thus simplifies Gibbs sampling.

In Stan, which uses Hamiltonian Monte Carlo sampling, there is no restriction for multivariate priors. The recommended approach is the separation strategy suggested by Barnard, McCulloch and Meng: $$\Sigma=\text{diag_matrix}(\sigma)\;\Omega\;\text{diag_matrix}(\sigma)$$ where $\sigma$ is a vector of std devs and $\Omega$ is a correlation matrix.

The components of $\sigma$ can be given any reasonable prior. As to $\Omega$, the recommended prior is $$\Omega\sim\text{LKJcorr}(\nu)$$ where "LKJ" means Lewandowski, Kurowicka and Joe. As $\nu$ increases, the prior increasingly concentrates around the unit correlation matrix, at $\nu=1$ the LKJ correlation distribution reduces to the identity distribution over correlation matrices. The LKJ prior may thus be used to control the expected amount of correlation among the parameters.

However, I've not (yet) tried non-normal distributions of random effects, so I hope I've not missed the point ;-)

Answer 2

I personally use Wishart proposals. For instance, if I want a proposal $\Sigma^*$ around $\Sigma$, I use: $$ \Sigma^* \sim \mathcal{W}(\Sigma/a,a), $$ where $a$ is a large number, like 1000. With that trick you will get $E[\Sigma^*]=\Sigma$ and you can adjust the variance with $a$. If I am not mistaken, the ratio of proposals for $(p\times p)$ matrices has a closed form: $$ \frac{q(\Sigma\to\Sigma^*)}{q(\Sigma^*\to\Sigma)} = \left(\frac{|\Sigma^*|}{|\Sigma|}\right)^{a-(p-1)/2} \cdot e^{[tr({\Sigma^*}^{-1}\Sigma)-tr({\Sigma}^{-1}\Sigma^*)] \cdot a/2}$$

Answer 3

It is well known that if you use non-Gaussian distributions, the conjugacy of the model is lost, see:

http://www.utstat.toronto.edu/wordpress/WSFiles/technicalreports/0610.pdf

Then, you need to use other MCMC methods, such as Metropolis within Gibbs sampling or some adaptive version of it. Fortunately, there is an R package for doing so:

http://cran.r-project.org/web/packages/spBayes/index.html

The recommended acceptance rate is 0.44 but, of course, there are some assumptions behind this number, similarly as in the case of the 0.234.

Are you THE Dimitris Rizopoulos?

Answer 4

Any proposal can be used if you define your log-posterior properly. You just need to use some tricks to implement it and properly define the support of your posterior, see:

How to find the support of the posterior distribution to apply Metropolis-Hastings MCMC algorithm?

There are tons of examples where a Gaussian proposal can be used for truncated posteriors. This is just an implementation trick. Again, you are asking a question with no general solution. Some proposals even have different performance for the same model and different data sets.

Good luck.