Re: Gamma distribution & Markov chain Monte Carlo

On Jul 17, 1:23 am, Maniaoh <n.hoai...@xxxxxxxxx> wrote:

Matlab.) None of Metropolis - Hastings or Metropolis method works.
Please kindly provide me some hints about this problem, I'm totally
lost.

Some comments. Take them as you will...

vari(1) = randn(1); % initial value

If you're trying to sample gamma variates, make sure you start your
chain with a value > 0. Sampling from the standard normal could give
you a negative value which causes the MCMC to blow up since the
density of a gamma < 0 is 0

for index = 2:n
        % normal distribution is used as proposal distribution
        y = vari(index - 1) + randn(1)*sig;

Ok, your proposal distribution is

p(x* | x) = N( x, sig^2 )

% draw u from normal distribution
u = randn(1);

This is wrong. Draw u from a *Uniform* distribution on [0, 1]

        % calculate accepting factor alpha
        % by Metropolis - Hastings method
        alpha = min(1, (gamDist(y, a, b))*norm(vari(index - 1), y,
sig)/...
                ((gamDist(vari(index - 1), a, b))*norm(y, vari(index -
1), sig)));

This looks ok, but you might want to compute each term for clarity
just is case there's something subtle, e.g.

x = vari(index - 1);
p_y = gamDist(y, a, b);
p_x = gamDist(x, a, b);
q_y_x = norm(y, x, sig);
q_x_y = norm(x, y, sig);

alpha = (p_y * q_x_y) / (p_x * q_y_x);

        % by Metropolis-only method
%       alpha = min(1, (gamDist(y, a, b))/(gamDist(vari(index - 1), a,
b)));

Again, looks ok....

        if u <= alpha
                vari(index) = y;
        else
                vari(index) = vari(index - 1);
        end

This is OK as long as you fix how you draw u. Add in some warning()
or error() statements to catch any case where you get bad values, e.g.
Infinite or NaN likelihood and likelihood ratios or if alpha < 0 or if
y every becomes negative.

-Lucas
.