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;
        x(i+1) = x(i);
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?

  • 2
    $\begingroup$ gamma(x(i)) does not sample from a gamma distributed random variable but computes the Gamma function of the value x(i). Consider using gamrnd() $\endgroup$ – muffin1974 Feb 1 '16 at 17:20
  • 1
    $\begingroup$ Maybe you should consider investing more time on learning about those algorithms than just watching YouTube. If you reflect just a wee bit about your question, how could one suggest simulating from a distribution f to bypass simulating from a distribution f? $\endgroup$ – Xi'an Feb 1 '16 at 20:38
  • 1
    $\begingroup$ @muffin1974: even simulating correctly from a Gamma distribution (with which parameters?) would not work since the support of the Gamma distribution is $\mathbb{R}⁺$ not $\mathbb{R}⁺$. And the Metropolis-Hastings ratio is incorrect since missing the Gamma proposals. $\endgroup$ – Xi'an Feb 1 '16 at 21:02

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.


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

  1. 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).


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.