How can Metropolis-Hastings use the function it is trying to approximate?

Question

The MH algoithm is used to obtain samples from a probability distribution $f$ that is difficult to sample from directly. The process as described in this answer is:

1. Pick a initial random state $x_0$.
2. Generate a candidate $x^{\star}$ from $Q(\cdot\vert x_0)$.
3. Calculate the ratio $\alpha=f(x^{\star})/f(x_0)$.
4. Accept $x^{\star}$ as a realisation of $f$ with probability $\alpha$.
5. Take $x^{\star}$ as the new initial state and continue sampling until you get the desired sample size.

But how can one compute $\alpha=f(x^{\star})/f(x_0)$ if we are trying to estimate $f$ in the first place? If we are able to calculate this ratio, why do we need to use Metropolis-Hastings?

A similar question was asked here and it linked to this tutorial where the author shows an example of sampling from an exponential distribution, but they use that exact exponential distribution in the sampling process, hence my confusion. Why can't we sample from the distribution directly if it is exactly known?

Answer 1

You misunderstood it, $f$ is known, just we don't know how to sample from it. For some distributions, we have dedicated algorithms to sample from them. There are also generic algorithms like the inverse transform sampling, but to use it you need to know the inverse cumulative probability function (or quantile function), but we don't always know it. Algorithms like Metropolis-Hastings are even broader, you can use them to sample from any distribution. More than this, you only need to know the probability density function $f$ up to a constant.

The caveat is that algorithms like Metropolis-Hastings are usually more costly computationally than the dedicated algorithms, so there's no free lunch.

Answer 2

The object of Metropolis-Hastings is to generate samples from $f$. The function $f$ could be complicated and need numerical evaluation. Metropolis-Hastings addresses the question:

How can I sample from $f$?

But, why would we need samples from $f$ if it is possible to evaluate the function? One example is Monte Carlo integration: $$ E[h(x)] = \int h(x) f(x)\,dx \underset{LLN}{\approx} \frac{1}{N}\sum h(X_i), \quad\quad X_i \sim f $$

As an example, consider $f = \frac{1}{c}(x_1^8 + x_2^2)$ and an expected value $E[\cos(x_1) \cdot \cos(x_2)]$. Given a sample, it is very easy to evaluate the expected value. Also, there is no need to calculate $c$ which can be intractable.

Answer 3

I will show you more visually the difference between evaluating the PDF (probability density function, or $f$) and sampling the PDF. Sampling a PDF is usually much harder and this problem gets worse the more dimensions your PDF has.

Using the following Mathematica code we can draw ('evaluate') the PDF. I chose a kind of arbitrary distribution whose PDF is just the sum of two normal distributions.

n = 1000;
distr = MixtureDistribution[{2, 1}, {NormalDistribution[], 
    NormalDistribution[3, 1/2]}];
Plot[PDF[distr][x], {x, -3, 5}, PlotRange -> Full]

Now we can sample the distribution. This means generating a random real number whose probability is proportional to the PDF at that real number. Or in simpler words: more points where the PDF is high. Note that the height y in the code below is just for visual purposes.

x = RandomVariate[distr, n];
y = RandomReal[{1, 2}, n];
ListPlot[{x, y}\[Transpose], PlotRange -> {{-3, 5}, {.5, 2.5}}, 
 AspectRatio -> Automatic, Axes -> {True, False}] 

If you draw a histogram of x it will look just like the PDF if you sample enough points.

There is also a 'direct' way to sample a distribution, but this is only possible when the CDF increases monotonically. See inverse transform sampling.

In physics there are two common use cases for Metropolis-Hastings. If we can write down an expression for the energy of a certain configuration of the system, we immediately know the probability distribution for that system. The probability is $p(\text{configuration})\propto e^{-E/kT}$, where $E$ is the energy of the configuration, $k$ is a constant and $T$ is the temperature. If we now perform Metropolis-Hastings on this system it would produce different configurations as if you were looking at a simulation. So MH can be used to simulate thermal systems, which is the first use case I was telling about. The second use case is to calculate expectation values, which was also told by Hunaphu. Naively one would calculate an expectation value as follows: $$E[g(x)]\approx\frac{\sum_i g(x_i)f(x_i)}{\sum_i f(x_i)}$$ where $x_i$ are evenly spaced on the number line. A lot of time is now spent calculating values that are basically zero. For a normal distribution the PDF quickly goes to zero when $x$ gets large. Imagine we can now efficiently sample from the PDF to produce samples $X_i$. In that case we can also calculate the expectation value as follows $$E[g(x)]\approx\frac{1}{N}\sum_i g(X_i)$$ wherever $f$ is high, $g$ gets sampled many times and contributes more to the expectation value.

Perhaps it is easier to grasp this expectation using a discrete distribution. Imagine a large number of students have taken test which has integer outcomes in the range $[1,10]$. If we know the probability distribution $f$ of the grades we could calculate the average grade using $$\frac{\sum_{x=1}^{10} x \cdot f(x)}{\sum_{x=1}^{10}f(x)}.$$ But, if a large number of students has taken the test we could equivalently calculate the average of these sample grades $$\frac 1 N\sum_i X_i$$