13
$\begingroup$

I have been trying to understand the Metropolis-Hastings algorithm in order to write a code for estimating the parameters of a model (i.e. $f(x)=a*x$). According to bibliography the Metropolis-Hastings algorithm has the following steps:

  • Generate $Y_t \sim q(y|x^t)$
  • $X^{t+1}=\begin{cases} Y^t, & \text{with probability} \quad \rho(x^t,Y_t), \\ x^t, & \text{with probability} \quad 1-\rho(x^t,Y_t), \end{cases}$

where $\rho(x,y)=\min \left( \frac{f(y)}{f(x)}*\frac{q(x|y)}{q(y|x)},1 \right)$

How I would like to ask few questions:

  • The bibliography states that if $q$ is a symmetric distribution the ratio $q(x|y)/q(y|x)$ becomes 1 and the algorithm is called Metropolis. Is that correct? The only difference between Metropolis and Metropolis-Hastings is that the first uses a symmetric distribution? What about "Random Walk" Metropolis(-Hastings)? How does it differ from the other two?
  • Most of the example code I find on-line uses a Gaussian proposal distribution $q$ and thus $\rho(x,y)=\min( f(y)/f(x),1)$ where $f(y)/f(x)$ is the likelihood ratio. What if the proposal distribution is a Poisson distribution? I think I understand rationally why that ratio does not become 1 when using an asymmetric distribution but I am not quite sure if understand it mathematically or how to implement it with code. Could someone give me a simple code (C, python, R, pseudo-code or whatever you prefer) example of the Metropolis-Hastings algorithm using a non symmetric proposal distribution?
$\endgroup$
16
$\begingroup$

The bibliography states that if q is a symmetric distribution the ratio q(x|y)/q(y|x) becomes 1 and the algorithm is called Metropolis. Is that correct?

Yes, this is correct. The Metropolis algorithm is a special case of the MH algorithm.

What about "Random Walk" Metropolis(-Hastings)? How does it differ from the other two?

In a random walk, the proposal distribution is re-centered after each step at the value last generated by the chain. Generally, in a random walk the proposal distribution is gaussian, in which case this random walk satisfies the symmetry requirement and the algorithm is metropolis. I suppose you could perform a "pseudo" random walk with an asymmetric distribution which would cause the proposals too drift in the opposite direction of the skew (a left skewed distribution would favor proposals toward the right). I'm not sure why you would do this, but you could and it would be a metropolis hastings algorithm (i.e. require the additional ratio term).

How does it differ from the other two?

In a non-random walk algorithm, the proposal distributions are fixed. In the random walk variant, the center of the proposal distribution changes at each iteration.

What if the proposal distribution is a Poisson distribution?

Then you need to use MH instead of just metropolis. Presumably this would be to sample a discrete distribution, otherwise you wouldn't want to use a discrete function to generate your proposals.

In any event, if the sampling distribution is truncated or you have prior knowledge of its skew, you probably want to use an asymmetric sampling distribution and therefore need to use metropolis-hastings.

Could someone give me a simple code (C, python, R, pseudo-code or whatever you prefer) example?

Here's metropolis:

Metropolis <- function(F_sample # distribution we want to sample
                      , F_prop  # proposal distribution 
                      , I=1e5   # iterations
               ){
  y = rep(NA,T)
  y[1] = 0    # starting location for random walk
  accepted = c(1)

  for(t in 2:I)    {
    #y.prop <- rnorm(1, y[t-1], sqrt(sigma2) ) # random walk proposal
    y.prop <- F_prop(y[t-1]) # implementation assumes a random walk. 
                             # discard this input for a fixed proposal distribution

    # We work with the log-likelihoods for numeric stability.
    logR = sum(log(F_sample(y.prop))) -
           sum(log(F_sample(y[t-1])))    

    R = exp(logR)

    u <- runif(1)        ## uniform variable to determine acceptance
    if(u < R){           ## accept the new value
      y[t] = y.prop
      accepted = c(accepted,1)
    }    
    else{
      y[t] = y[t-1]      ## reject the new value
      accepted = c(accepted,0)
    }    
  }
  return(list(y, accepted))
}

Let's try using this to sample a bimodal distribution. First, let's see what happens if we use a random walk for our propsal:

set.seed(100)

test = function(x){dnorm(x,-5,1)+dnorm(x,7,3)}

# random walk
response1 <- Metropolis(F_sample = test
                       ,F_prop = function(x){rnorm(1, x, sqrt(0.5) )}
                      ,I=1e5
                       )
y_trace1 = response1[[1]]; accpt_1 = response1[[2]]
mean(accpt_1) # acceptance rate without considering burn-in
# 0.85585   not bad

# looks about how we'd expect
plot(density(y_trace1))
abline(v=-5);abline(v=7) # Highlight the approximate modes of the true distribution

enter image description here

Now let's try sampling using a fixed proposal distribution and see what happens:

response2 <- Metropolis(F_sample = test
                            ,F_prop = function(x){rnorm(1, -5, sqrt(0.5) )}
                            ,I=1e5
                       )

y_trace2 = response2[[1]]; accpt_2 = response2[[2]]
mean(accpt_2) # .871, not bad

This looks ok at first, but if we take a look at the posterior density...

plot(density(y_trace2))

enter image description here

we'll see that it's completely stuck at a local maximum. This isn't entirely surprising since we actually centered our proposal distribution there. Same thing happens if we center this on the other mode:

response2b <- Metropolis(F_sample = test
                        ,F_prop = function(x){rnorm(1, 7, sqrt(10) )}
                        ,I=1e5
)

plot(density(response2b[[1]]))

We can try dropping our proposal between the two modes, but we'll need to set the variance really high to have a chance at exploring either of them

response3 <- Metropolis(F_sample = test
                        ,F_prop = function(x){rnorm(1, -2, sqrt(10) )}
                        ,I=1e5
)
y_trace3 = response3[[1]]; accpt_3 = response3[[2]]
mean(accpt_3) # .3958! 

Notice how the choice of the center of our proposal distribution has a significant impact on the acceptance rate of our sampler.

plot(density(y_trace3))

enter image description here

plot(y_trace3) # we really need to set the variance pretty high to catch 
               # the mode at +7. We're still just barely exploring it

We still get stuck in the closer of the two modes. Let's try dropping this directly between the two modes.

response4 <- Metropolis(F_sample = test
                        ,F_prop = function(x){rnorm(1, 1, sqrt(10) )}
                        ,I=1e5
)
y_trace4 = response4[[1]]; accpt_4 = response4[[2]]

plot(density(y_trace1))
lines(density(y_trace4), col='red')

enter image description here

Finally, we're getting closer to what we were looking for. Theoretically, if we let the sampler run long enough we can get a representative sample out of any of these proposal distributions, but the random walk produced a usable sample very quickly, and we had to take advantage of our knowledge of how the posterior was supposed to look to tune the fixed sampling distributions to produce a usable result (which, truth be told, we don't quite have yet in y_trace4).

I'll try to update with an example of metropolis hastings later. You should be able to see fairly easily how to modify the above code to produce a metropolis hastings algorithm (hint: you just need to add the supplemental ratio into the logR calculation).

$\endgroup$
  • $\begingroup$ Awesome answer! Thank you so much! In my case I have a 6-7 parameter model and I have no idea how the posterior distribution might look like (but it might be bimodal) because my datasets are sometimes completely different. Based on what you said I can either use Metropolis(-Hastings) using a huge variance in the proposal distribution or use Random Walk Metropolis(-Hastings) with a smaller variance in the proposal distribution. Under no special circumstances the second solution should converge faster to the target distribution. Right? $\endgroup$ – AstrOne Jul 16 '13 at 5:41
  • $\begingroup$ Now related to the Metropolis-Hastings code I was thinking to replace R=exp(logR) with this: R=exp(logR)*(dnorm(y[t-1],y.prop,my_sigma)/dnorm(y.prop,y[t-1],my_sigma)) for both random and non-random walk MH. Is that correct? $\endgroup$ – AstrOne Jul 16 '13 at 5:43
  • 1
    $\begingroup$ Basically, but as I mention in the metropolis code: you want to perform your calculations in log space. Likelihood calculations tend to operate on very small values, so you generally get much better results adding logarithms and exponentiating your results than multiplying raw values together. $\endgroup$ – David Marx Jul 16 '13 at 9:22
  • 1
    $\begingroup$ You don't need to worry about the current state of the chain (i.e. $y_{t-1}$) when you're using a fixed proposal distribution because: $q(y_t|y_{t-1})=q(y_t)$. A fixed proposal distribution generates independent propositions. We take $y_{t-1}$ into account in the metropolis ratio. $\endgroup$ – David Marx Jul 17 '13 at 4:49
  • 1
    $\begingroup$ You state "In a non-random walk algorithm, the proposal distributions are fixed. In the random walk variant, the center of the proposal distribution changes at each iteration" which is not correct. MH versions that are not random walks mostly have proposals that depend on the current state of the Markov chain, sometimes even centred at this state. A major example is the Langevin MCMC algorithm. When the proposal is fixed, this is an independent MH algorithm. $\endgroup$ – Xi'an Feb 13 '16 at 15:49

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.