How to reduce autocorrelation in Metropolis algorithm?

Question

I've been using a Metropolis/Gibbs sampler combination to generate a joint density for some parameters(it is a hierarchical model, with $y_i\sim Poisson(\lambda_i)$, $\lambda_i\sim Gamma(\alpha,\beta)$). What techniques can I use to lower autocorrelation(it is present in $\alpha$ and $\beta$)? I have been using thinning, but even when I use huge lags(4950, which is reaching the memory limit on my computer to use) there is still significant autocorrelation. Is there something I could do with my step size distribution to help with this? I have been drawing the new values of $\alpha$ and $\beta$ from a normal distribution, with mean the current parameter value, and standard deviation 1. Thanks!

Answer 1

First of all, for your model to be hierarchical, you need hyperpriors for $\alpha$ and $\beta$ (as already explained by Procrastinator). For the sake of simplicity, lets assume uniform priors on the positive part of the real axis. So that have a hierarchical model as follows: $$y_{i}| \lambda_{i}\sim Poisson(\lambda_{i})$$ $$\lambda_{i}|\alpha, \beta \sim Gamma(\alpha,\beta)$$ $$\pi(\alpha,\beta)\propto1_{[0,+\infty]}$$

Now you have a two groups of parameters: {$\lambda_{i}$}$_{i=1,N}$ and {$\alpha, \beta$}. You need to draw random samples of these two groups of parameters from your posterior distribution. Although the model itself is not very complicated, you might stuck with very slowly mixing chains for you parameters, especially if your $N$ is very large (say 1000). You can choose as your proposal normal distribution (just remember to reject properly negative proposed values) for whole set of parameters, in which case you would need ($N+2$)-variate normal distribution as a proposal, or you could use N-variate normal distirbution for lambda's and bivariate normal distribution for $\alpha$ and $\beta$.
I would suggest first of all to go with separately proposing lambda's and gamma parameters - i.e. to use Metropolis within Gibbs sampler. This would allow you to slightly decouple thouse chains. In this investigatory step I would use covariance matrix with no autocorrelations (i.e. diagonal matrix) for $N$-variate normal distribution.
If that does not work, I would introduce for proposal distribution covariance matrix where correlations are not equal to zero - this should improve the mixing. And if this does not produce well-mixing chains, I would turn to the Hamiltonian Monte Carlo. But first of all try to play with different matrices of multivariate normal proposal. I would also suggest to modify your model: instead of Poisson intensity $\lambda_{i}$ to use $exp(\lambda_{i})$ and then to use normal distribution for $\lambda_{i}$ with unknown mean and variance parameters.

Answer 2

Can't seem to figure out how to add a comment, so this is part comment part answer:

First, it sounds from the comments that you are using a Random Walk Metropolis step to update $\alpha$ and $\beta$ jointly (or possibly separately), and then a Gibbs step to update each $\lambda_i$ conditional on the current values for $\alpha$ and $\beta$.

You say you're using a truncated Normal proposal in the (random walk) Metropolis step: it probably won't make much difference in your case, but technically I think you need to propose from a full Gaussian, not truncated - otherwise your proposal distribution isn't symmetric, so $q(y|x) \neq q(x|y) \: \forall x,y$, so the Metropolis acceptance ratio doesn't produce a Markov chain which will converge to the desired density. Alternatively you could keep the truncated $q$ and use the full Metropolis-Hastings acceptance ratio.

Secondly, as Procrastinator says, you presumably therefore have a posterior for $(\alpha, \beta)$, so you must have set a prior for these. What is it? Gamma priors with fixed hyperparameters (something like shape=1, rate=0.001) would do probably do the trick).

Now, to attempt to answer some of your question: assuming you've done everything right and you're getting large autocorrelation because each $(\alpha^{(i)},\beta^{(i)})$ is highly correalated with each $(\alpha^{(i-1)},\beta^{(i-1)})$, then maybe the Random Walk Metropolis isn't the best method to draw $\alpha$ and $\beta$. Since your MH step is in only two dimensions, an independence sampler might work better. Slightly more advanced methods like the Metropolis-adjusted Langevin algorithm and Hamiltonian Monte Carlo are also designed to reduce auto-correlation, but that might be overkill here. So maybe try an independence sampler (don't forget to use the full Metropolis-Hastings ratio $\min\left(1,\frac{\pi(y)q(x|y)}{\pi(x)q(y|x)} \right)$ for acceptance).