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?