Acceptance probability for Metropolis-Hastings MCMC on multinomial-Dirichlet model

Question

As an exercise to learn how to manually code MCMC, I've built a Metropolis-Hastings sampler on top of a multinomial-Dirichlet posterior distribution. Since a closed form solution exists, I can compare results from the MCMC with simulations from the actual posterior distribution.

I'm using a Dirichlet proposal distribution with parameters equal to the latest probabilities in the chain times a scaling constant (~1000), which makes the expected value of the distribution equal to those probabilities, with the scaling constant controlling the variance. Since this distribution certainly isn't symmetric, I tried adding the ratio of the of values from the proposal distribution to the calculation of the acceptance probability.

Doing this, however, seems to bias the results away from the results given by the closed form solution. The only way I've gotten results from the MCMC to match the results from the closed-form solution is to calculate the acceptance probability from the posterior distribution alone, as you would if the proposal distribution were symmetric. R code here: https://github.com/sivadivel/nps_stats/blob/master/manual_mcmc.R

My question is: why is this the case?

Answer 1

You either do a random walk (or not) with symmetric proposal distributions, or use a fixed asymmetric proposal distribution.

According to the answer by David Marx under this question:

I suppose you could perform a "pseudo" random walk with an asymmetric distribution which would cause the proposals to drift in the opposite direction of the skew (a left skewed distribution would favor proposals toward the right).

This is probably why you observed:

The only way I've gotten results from the MCMC to match the results from the closed-form solution is to calculate the acceptance probability from the posterior distribution alone, as you would if the proposal distribution were symmetric.

And to calculate ppropos in your code, you should write:

ppropos = log(ddirichlet(chain[i,1:3], a)) - log(ddirichlet(proposal, a))

with the assumption that the proposal distribution is $Dir(a)$ and fixed.