4
$\begingroup$

If I understood it correctly, the Metropolis-Hastings algorithm allows one to sample from a distribution without an analytical representation, which comes in handy, for instance in the Bayesian statistics literature.

Quite generally, the algorithm seems to work as follows (this is similar to Koop (2004, p. 93): Bayesian Econometrics):

  1. Choose a starting value, $\theta^{(0)}$.
  2. Take a candidate draw, $\theta^{*}$, from the candidate generating density, $q(\theta^{(s-1)}; \theta)$.
  3. Calculate an acceptance probability based on the candidate draw $\theta^{*}$ and the previous draw $\theta^{(s-1)}$: $\alpha(\theta^{(s-1)}, \theta^*)$.
  4. Set $\theta^{(s)}=\theta^*$ with probability $\alpha$ and set $\theta^{(s)}=\theta^{(s-1)}$ with probability $1-\alpha$.
  5. Repeat steps 2 to 4 often enough to reach convergence.

Later on, they argue that the acceptance probability is a function of the posterior density evaluated at point $\theta^*$ or point $\theta^{(s-1)}$. This, I don't understand. It seems that you need to know the very same posterior density function to compute the acceptance probability that the algorithm is meant to sample from because it is unknown. This seems circular to me.

Also, I looked at the Metropolis-Hastings chapter from Robert/Casella (2010): Introducing Monte Carlo Methods with R. Here, I found for instance this example to draw from a beta distribution:

a=2.7; b=6.3; c=2.669 # initial values
Nsim=500000
X=rep(runif(1),Nsim) # initialize the chain
for (i in 2:Nsim){
  Y=runif(1)
  rho=dbeta(Y,a,b)/dbeta(X[i-1],a,b)
  X[i]=X[i-1] + (Y-X[i-1])*(runif(1)<rho)
}

#Here, the algorithm is just an alternative to sampling directly from rbeta
summary(X)
  Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00255 0.18970 0.28480 0.30010 0.39520 0.92870 
summary(rbeta(Nsim, a, b))
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.001079 0.189500 0.284500 0.299900 0.395000 0.916600

Note that they have to call dbeta to compute the acceptance probability. They even mention in their book that this example is just for educational purposes because you can directly sample from rbeta in this case. But what are real cases in which I would know a density function (here dbeta) but couldn't sample from it (here with the function rbeta)?

$\endgroup$
  • $\begingroup$ I posted an answer below, but also think about how you would be able to answer this question if the functions dbeta and rbeta did not exist. Those functions are for convenience but you do not actually need them to solve the problem. $\endgroup$ – user25658 Sep 10 '13 at 16:17
  • $\begingroup$ @user25658: you need a version of dbeta to compute the acceptance probability. If the target density is unavailable, Metropolis-Hastings cannot operate. $\endgroup$ – Xi'an Apr 5 '16 at 12:45
  • $\begingroup$ you put a link to a chapter of our book, which seems like a copyright infringement, unless there is some provision unknown to me for this kind of posting... $\endgroup$ – Xi'an Apr 5 '16 at 14:11
6
$\begingroup$

I think the point you may be missing is that you may know the density but you still don't know how to sample from it and that is where the Metropolis-Hastings (M-H) algorithm comes in. So you can't sample from the density but you can evaluate it. To try and make these points more concrete consider the following density

$$f_X(x)=\frac{100}{x^2}\text{ for }x>100$$

You don't recognize this distribution as something you know how to sample from so you can use the M-H algorithm to obtain samples from it. And clearly you will be using this density to evaluate candidate points for acceptance/rejection even though you are not yet sampling from it.

Also, something that is pretty neat about the M-H algorithm is that all you need to know if the density you want to sample from up to a normalizing constant. For example, if you wanted to sample from the Beta distribution $\text{Beta}(\alpha,\beta)$ you do not need to have the normalizing constant, you just need the kernel of the density, i.e.,

$$f_X(x)\propto x^{\alpha-1}(1-x)^{\beta-1}$$ So just a proportional density is all that you need.

$\endgroup$
  • 2
    $\begingroup$ Thanks a lot. It's true, I never gave it much thought how one actually samples from a density function. I just implicitely assumed that this is the simple part. $\endgroup$ – Christoph_J Sep 10 '13 at 16:24
  • $\begingroup$ No problem. It is a very reasonable assumption that you made. $\endgroup$ – user25658 Sep 10 '13 at 16:29

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.