# Metropolis Hastings Algorithm - Prior vs Proposal vs Numerator of Bayes Theorem

I've been using this technique in 'black-box' form for a little while as a physics student.

I have been struggling to understand what's happening under the hood for some time and I think I almost have it - but I'm mixed up about several things. I'm trying to fit together the prior distribution, proposal distribution, and Bayes Theorem

What is the difference between the prior $\mathcal{P}$ and the jumping distribution $\mathcal{Q}$ in the simple MH algorithm? Are these the same thing? I know the jumping distribution is used to draw a candidate point $X_c$ that is 'near' the current point $X_i$. Often this is something like a Gaussian function centered at $X_i$, therefore producing $X_c$ that is usually locally 'close' to $X_i$. That's what I keep finding, using resources such as this youtube video.

Elsewhere, like in this amazing tutorial, I see that the prior is a surface above the parameter space in question shaped in a way that makes sense given the problem, and this surface is used to weight the posterior in conjunction with the appropriate likelihood function. For instance maybe we believe in advance that scores on an particular exam are Gaussian distributed, with known mean. So the posterior is weighted by this assumption and the MCMC will return samples with a higher density in that region for that particular parameter.

That makes sense with Bayes Theorem, as the posterior is proportional to the product of the prior and the likelihood.

$P(A|B) \propto P(B|A)P(A)$

I am not seeing where this happens into the MH algorithm. I know these three things (prior, proposal, and Bayes Theorem) are all wrapped up in the 'acceptance ratio' part of the algorithm. Can someone break it down for me like I'm five years old and show me exactly where and how those pieces fit together?

Finally I'm having a lot of trouble understanding how we define likelihood for an arbitrary statistical model. One explanation from the first video is it's an exponential function that literally evaluates the difference between the proposed fit and the real data:

$\mathcal{L} = e^{-(Data - Model)}$

Okay, so this function takes a higher value in regions of parameter space that produce better fits. But how can we say this is exactly the actual likelihood (i.e. the $P(B|A)$ part of Bayes Theorem)? This doesn't seem like it's guaranteed to be at all points proportional to what we want...

Any help, guidance, or resources greatly appreciated! Thanks so much for reading to this point and bearing with me.

As you say, the three elements used in MH are the proposal (jumping) probability, the prior probability, and the likelihood. Say that we want to estimate the posterior distribution of a parameter $\Theta$ after observing some data $\mathbf x$, that is, $p(\Theta|\mathbf{x})$. Assume that we know the prior distribution $p(\Theta)$, that summarizes our beliefs about the value of $\Theta$ before we observe any data.

Now, it is usually impossible to compute the posterior distribution analytically. Instead, an ingenious method is to create an abstract Markov chain whose states are values of $\Theta$, such that that the stationary distribution of such chain is the desired posterior distribution. Metropolis-Hastings (MH) is a schema (not the only one, e.g. there's Gibbs sampling) to construct such a chain, that requires to carefully select a jumping (or proposal) distribution $q(\Theta|\theta)$. In order to go from one value of $\Theta$, denoted as $\theta$, to the next, say $\theta'$, we apply the following procedure:

1. Sample a candidate (or proposed) $\theta^*$ as the next value, by sampling from $q(\Theta|\theta)$, where $\theta$ is the current value.
2. Accept the candidate value with a probability given by the MH acceptance ratio, given by the formula: $$\alpha(\theta,\theta^*) = \min\left[1,\frac{p(\theta^*|\mathbf{x})\;q(\theta|\theta^*)}{p(\theta|\mathbf{x})\;q(\theta^*|\theta)} \right].$$ By applying Bayes rule the the posterior probability terms in the formula above, we get: $$\alpha(\theta,\theta^*) = \min\left[1,\frac{p(\theta^*)\;p(\mathbf{x}|\theta^*)\;q(\theta|\theta^*)}{p(\theta)\;p(\mathbf{x}|\theta)\;q(\theta^*|\theta)} \right].$$

After iterating this process "enough" times, we are left with a collection of points that approximates the posterior distribution.

A counterintuitive thing about the formula above is that the proposal probability of the candidate value appears at the denominator, while the "reverse" proposal probability (i.e. going from the proposed to the original value) is at the numerator. This is so that the overall transition distribution resulting from this process ensures a necessary property of the Markov chain called detailed balance. I found this paper quite helpful on this topic.

Now, it is perfectly possible to use the prior distribution itself as the proposal distribution: $q(\Theta|\theta)=p(\Theta)$. Note that in this case the proposal distribution is not conditional on the current value of $\Theta$, but that is not a problem in theory. If we substitute this in the formula for $\alpha$ above, and carry out some simplifications, we obtain: $$\alpha(\theta,\theta^*) = \min\left[1,\frac{p(\mathbf{x}|\theta^*)}{p(\mathbf{x}|\theta)} \right].$$ What is left is just the ratio of the likelihoods. This is a very simple approach and usually not very efficient, but may work for simple problems.

Regarding the likelihood, I think it really depends on what your model is. Regarding the formula you write, I don't really understand what is going on. What are $Data$ and $Model$ in there?

• Thank you for your reply. To clarify my question about likelihood: I am unsure how we define the likelihood of a set of parameters given an arbitrary statistical model. For instance, suppose I want to estimate $\mu$ and $\sigma$ for a gaussian function. So, I'm at some point in parameter space $\theta$ in your notation, and that corresponds to some values of $\mu$ and $\sigma$. By exactly what metric do I determine the likelihood of this point versus other points? And what about for some more complicated, arbitrary formula whose parameters I want to estimate? – FaceInvader Jul 29 '15 at 22:51
• I need to compare the gaussian I'm guessing to the one I've observed (the signal) and in some way assign a number that represents 'how good' the fit is. I'm a little confused by this. – FaceInvader Jul 29 '15 at 23:00
• Let's take the case of the Gaussian, where $\theta=(\mu,\sigma)$. If the data $\mathbf{x}=x_1$ is a single point drawn from such Gaussian, then the likelihood is simply the density function evaluated at that point, i.e. $$p(x_1|\mu,\theta)=\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{(x_1-\mu)^2}{2\sigma^2}}.$$ If instead you have a sequence $\mathbf{x}=(x_1,x_2,\ldots,x_N)$ of i.i.d. draws from such Gaussian, the likelihood is just the product of the above. In this simple case, however, you wouldn't really need to use MH, since inference is analytically tractable using a conjugate prior. – breadshe Jul 30 '15 at 0:00
• I actually think I get it now, thanks! At least for the gaussian case. What about for some horrible random formula where we are fitting parameters to some complicated high dimensional model? Is it okay to take the square difference and do something like $Likelihood = e^{-(difference)^2}$ ? How can we be sure this is at all points exactly $P(Data|Parameters)$ ? – FaceInvader Jul 30 '15 at 1:09
• I don't really understand what you mean by $Difference$, or $Data-Model$ in the original question. If $Data$ is a training data sequence, that I called $\mathbf{x}$ above, what does it mean to subtract "$Model$"? A model, even if interpreted as a set of parameters, is not a quantity comparable to the $Data$. Where did you get that formula? In general, the specific formula for the likelihood $p(Data|Parameters)$, or $p(\mathbf{x}|\theta)$ in my notation, should be derived from the Bayesian generative model (e.g. the Gaussian model above.) – breadshe Jul 30 '15 at 4:30

First to explain the MH algorithm, it's used to approximate numerically a target distribution, in this case $p(\theta|D)$. At each stage of the algorithm:

1. A value $\theta_{proposed}$ is proposed using the jumping or proposal distribution.
2. An acceptance ratio is calculated, equal to $\frac{p(\theta_{proposed}|D)}{p(\theta_{current}|D)}$. Because we cannot calculate $p(\theta|D)$ directly, we leverage the proportional expression of Bayes rule and calculate this quotient using the products of likelihood and prior corresponding to $\theta_{proposed}, \theta_{current}$. That is, the acceptance ratio is:

$$\frac{p(D|\theta_{proposed})p(\theta_{proposed})}{p(D|\theta_{current})p(\theta_{current})}$$

1. If this ratio exceeds one—intuitively, if the proposed parameter value is more likely given data and prior—we accept this proposal. If not, we accept it with probability equal to the ratio.

As more and more values are sampled in this way, the trace of $\theta$ values more and more closely approximates the true distribution of $\theta$.