The Gamma distribution is log concave if the shape parameter is
at least one; the log Gamma distribution for any value of
the parameter. If the parameter is at least 1/3, the cube
root of a Gamma random variable is "dominated" by the normal
density fitting it at the mode; this is the Wilson-Hilferty
transformation.  Acceptance rejection methods are quick.

MKMC does not give independent Gamma random variables; this
does with no "burn in".

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

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

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

% draw u from normal distribution
u = randn(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)));

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

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

Hi guys,

Thanks in advance. Any comment is highly appreciated.
Iaoh.

--------------------

clear all

% number of iterations
n = 5000;
% burn-in rounds
burnin = 2000;
% preallocate the array
vari = zeros(1, n);

str = 'x^(a - 1)*exp(-b*x)'; % gamma distribution
gamDist = inline(str, 'x', 'a', 'b');
a = 1;
b = 1;

str = 'exp(-0.5*((x-mu)/sig)^2)';
norm = inline(str, 'x', 'mu', 'sig');
sig = 1;

vari(1) = randn(1); % initial value
for index = 2:n
% normal distribution is used as proposal distribution
y = vari(index - 1) + randn(1)*sig;
% draw u from normal distribution
u = randn(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)));
% by Metropolis-only method
% alpha = min(1, (gamDist(y, a, b))/(gamDist(vari(index - 1), a,
b)));
if u <= alpha
vari(index) = y;
else
vari(index) = vari(index - 1);
end
end

% remove burn-in elements
vari = vari(burnin:n);
hist(vari);

On Jul 17, 1:23 am, Maniaoh <n.hoai... at (no spam) gmail.com> 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

Hi Lucas,

Thank you very much. I was able to code it right. Use of normal
distribution sometimes create value < 0 that causes terrible pain. I
also did add a guard to check it and change the way I draw u.

On Jul 18, 12:47 am, Maniaoh <n.hoai... at (no spam) gmail.com> wrote:

Hi Lucas,

Thank you very much. I was able to code it right. Use of normal
distribution sometimes create value < 0 that causes terrible pain. I
also did add a guard to check it and change the way I draw u.

Just to be clear, there's nothing wrong with your Normal proposal
generating values, y, less than zero. The gamma density, gamDist(y,
a, b), should return zero for negative arguments. Since this factor
is in the numerator of the acceptance ratio, alpha becomes zero and
the the negative value will never be accepted.

It sounds like you may have "fixed" your code by preventing the
proposal distribution from proposing negative numbers. If you did this
by adding a loop like

% normal distribution is used as proposal distribution
y = vari(index - 1) + randn(1)*sig;
while y < 0,
y = vari(index - 1) + randn(1)*sig;
end

then you have *changed your proposal distribution* from a Normal to a
truncated Normal distribution and your acceptance ratio calculation is
no longer correct. In fact, you would have to use only Metropolis-
Hastings since this proposal is not symmetric and the Metropolis
algorithm no longer applies.

It's very easy to invalidate a sampling algorithm by making "obvious"
changes.

-Lucas

Hi Lucas,

Just as you say, a proposal of normal distribution does nothing
harmful to the target distribution if it generates value less than
zero. It can be confirmed through explicit expression of alpha.
Nevertheless, when I allow normal distribution to generate less-than-
zero values, the target distribution does not seem right. When I draw
histogram, on the negative side of axis, some values of y (target) are
generated.

You exactly state what I did, even the code , I took a guard and
therefore the proposal normal distribution became a truncated one. In
this case, the result is better, i.e., histogram. Do you have any
comment on this?