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

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?

• We're not trying to estimate $f$; rather, we're trying to create a sample from the distribution it represents. Can you expand on what you mean by "sample from the distribution directly"? Jan 3 at 18:02
• @JohnMadden Thanks, my confusion is that we seem to be able to sample from other distributions like uniform, Gaussian, etc without having to resort to Metropolis-Hastings. I thought sampling just entails generating a random number between 0 and 1 and plugging that into your density function, no MH needed. So I suppose I should ask, when do we not need a special algorithm like MH to sample from a distribution? Jan 3 at 18:34
• Ahh, here's the confusion: rather than putting the random number into the density function, we actually put it into the quantile function (aka inverse cdf). The quantile function is the functional inverse of the integral of $f$, so it will not be known for general $f$. Jan 3 at 18:39
• Jan 4 at 14:00

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.

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.

• Sorry I don't understand, What does the LLN under the approximately equal sign mean? And how is $cos(x_1)\cdot cos(x_2)$ related to $\frac1c (x^8_1 + x^2_2)$? Jan 4 at 3:48
• And how do we use the samples to calculate the expected value without $c$? Jan 4 at 4:13

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$$

• As it’s currently written, your answer is unclear. Please edit to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers in the help center.
– Community Bot
Jan 4 at 22:57
• Thank you, this helped a lot. My only question is about your statement "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." So why is it that we cannot just use the inverse of the PDF, rather than the inverse of the CDF as in inverse transform sampling? Jan 4 at 23:57
• @xojfqa Glad I could help. I don't know the exact reason and I feel like a good answer would not fit in a comment so perhaps this could be a new question. Jan 5 at 14:26