Solved – Metropolis-Hastings algorithm, using a proposal distribution other than a Gaussian in Matlab

markov-chain-montecarlomathematical-statisticsMATLABmetropolis-hastingsnormal distribution

I am currently working on my final year project for my mathematics degree which is based on giving an overview of the Metropolis-Hastings algorithm and some numerical examples. So far I have got some great results by using my proposal distribution as a Gaussian, and sampling from a few other distributions, however I am trying to go one step further by using a different proposal distribution.

So far I have got this code (I am using Matlab), however with limited resources online about using different proposals it is hard to tell if I am close at all, as in reality I am not too sure how to attempt this, (especially as this gives no useful data output so far).

It would be fantastic if someone could lend a hand if they know or forward me to some easily accessible information (I understand that I am not just asking coding advice, but Mathematics as well).

So, I want to sample from a Gaussian using a proposal distribution of a Laplace, this my code so far:

n = 1000;       %%%%number of iterations

x(1) = -3;      %%%%Generate a starting point

%%%%Target distribution: Gaussian:

strg = '1/(sqrt(2*pi*(sig)))*exp(-0.5*((x - mu)/sqrt(sig)).^2)';
tnorm = inline(strg, 'x', 'mu', 'sig');

mu = 1;    %%%%Gaussian Parameters (I will be estimating these from my markov chain x)
sig = 3;


%%%%Proposal distribution: Laplace:

strg = '(1/(2*b))*exp((-1)*abs(x - mu)/b)';
laplace = inline(strg, 'x', 'b', 'mu');

b = 2;       %%%%Laplace parameter, I will be using my values for y and x(i-1) for mu


%%%%Generate markov chain by acceptance-rejection

for i = 2:n

    %%%%Generate a candidate from the proposal distribution
    y = laplace(randn(1), b, x(i-1));

    %%%%Generate a uniform for comparison
    u = rand(1);

    alpha = min([1, (tnorm(y, mu, sig)*laplace(x(i-1), b, y))/(tnorm(x(i-1), mu, sig)*laplace(y, b, x(i-1)))]);


    if u <= alpha
        x(i) = y;
    else
        x(i) = x(i-1);
    end 
end

If anyone can tell me if the above is completely wrong/going about it in the wrong way, or there are just a few mistakes (I am very wary about my generation of 'y' in for the for loop being completely wrong) that would be fantastic.

Thanks, Tom

Best Answer

As far as for the code, I think there is only one error which is the generation of random numbers according to Laplace distribution.

Replace 

y = laplace(randn(1), b, x(i-1));

by (See the Laplace distribution wiki to learn how to generate random Laplace distributed numbers before seeing my answer.)

u = rand()-0.5;
uu = b / sqrt(2); 
y = mu - uu * sign(u).* log(1- 2* abs(u));

One more thing, Laplace distribution is symmetric, so laplace(x(i-1), b, y))=laplace(y, b, x(i-1)). Therefore, when calculating your alpha, you can omit laplace(x(i-1), b, y))/laplace(y, b, x(i-1))=1.

There are a lot of books on MCMC approaches. But for a quick start, maybe this paper and this note is good enough.

Related Question