1
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ A scaling constant of 1000 seems enormous and may have the consequence of turning the proposal into a Dirac mass. It has also the side effect of producing a normalising constant in the Dirichlet that may create an overflow: in R gamma(1000) returns Inf $\endgroup$ – Xi'an Jul 12 '17 at 7:07
  • $\begingroup$ You're right Xi'an, using a scaling constant of 1000 produces an acceptance rate of 91% with a lot of autocorrelation in the chain, and lowering it to around 35 produces an acceptance rate of ~45% with much less autocorrelation. $\endgroup$ – Levi Davis Jul 12 '17 at 14:47
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.