For Metropolis-Hastings algorithm, should target density and proposal distribution have the same distribution?

Question

I watched some youtube videos about the Metropolis-Hastings algorithm. They used a Gaussian as a proposal function to estimate an unknown Gaussian, or used a Gamma function as the proposal function to estimate an unknown Gamma.

I tested this with some simple MATLAB code and it didn't work.

n = 20000;
x = zeros(n,1);
x(1) = 0.5;

for i = 1:n-1
    % x_c = normrnd(x(i), 0.05); % proposal function (0.5, 0.05)
    x_c = gamma(x(i)); % will it works with a different distribution?
    if rand() < min(1, normpdf(x_c) / normpdf(x(i)))
        x(i+1)= x_c;
    else
        x(i+1) = x(i);
    end
end
hist(x, 100) % to get the standard Gaussian (0,1)

So I wonder when using Metropolis-Hastings algorithm to estimate an unknown distribution, should the proposal function have the same distribution?

Answer 1

The target function's form is often not known and often does not have any particular parametric representation, so the answer is no. I frequently sample from a Gaussian proposal distribution when I don't expect the target distribution to have any particular form.

If you have some knowledge of the domain, such as nonnegative reals or integers or something, that can help you decide which distributions from which to sample.

As others pointed out in the comments, you are not using the correct Matlab function for sampling. See gamrnd.

Answer 2

If you can, then yes it should. When the proposal distribution equals the target distribution, then you are in the optimal situation of Monte Carlo integration. Lets see why.

Formally, what you are describing is that $g(x'|x) = \pi(x')$. In this case:

1. the acceptance probability is always 1 because

$$\frac{g(x|x')}{g(x'|x)}\frac{\pi(x')}{\pi(x)} = 1$$

2. The states are always independent. This is a consequence of $g(x'|x)$ not depending on $x$.

Therefore, the Markov process corresponds to directly draw independent samples $x$ from $\pi(x)$, which is actually the problem you were trying to solve in the first place.

The main issue is that it is often the case that $\pi(x)$ is not an easy function nor it can be directly sampled from. In this case, the best you can do is to correlate the states in order to sample from $\pi(x)$ (through $g(x'|x)$). This will naturally require a different distribution, that often correlates the samples. Therefore, the autocorrelation time of the process will be higher and you will need more samples, making the process suboptimal (and often the best you can do about it).