Re: Gamma distribution & Markov chain Monte Carlo

On Jul 18, 2:59 am, Maniaoh <n.hoai...@xxxxxxxxx> wrote:
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?

If you are still seeing negative values of y then there is something
wrong in your code. I'm sure truncating y to be positive does improve
the histograms, but that doesn't mean it's correct. :)

I implemented your code myself as follows. It appears to work fine --
no negative y values and a histogram of 10,000 samples is a pretty
good looking gamma. Use this as a reference to check your own code.

function s = gamma_mcmc(n)

s = zeros(1,n);

gampdf = @(x,a,b) x.^(a-1) * exp(-b*x) .* (x>=0);
normpdf = @(x,mu,sig) exp(-0.5*((x-mu)/sig).^2);

a = 1;
b = 1;
sig = 1;

% Initialize chain at the mean of the gamma distribution
x = a * b;

for i = 1:n,

% Generate a proposal from a normal distribution
x_star = x + sig .* randn;

% Calculate the M-H acceptance
q_x_given_x_star = normpdf( x, x_star, sig );
q_x_star_given_x = normpdf( x_star, x, sig );
p_x = gampdf( x, a, b );
p_x_star = gampdf( x_star, a, b );

alpha = (p_x_star * q_x_given_x_star) / (p_x * q_x_star_given_x );

if rand <= alpha
x = x_star;
end

s(i) = x;
end
.

Relevant Pages

  • Re: compare 2 samples
    ... I.e. mean .005 seconds with std dev .02 seconds. ... For some reason I was kind of envisioning a normal distribution just truncated at zero. ...
    (sci.stat.math)
  • Re: Is driving time or average velocity normally distributed?
    ... >>commute is about 10 miles away, and on the average, it takes me 25 ... Not sure what you mean by a histogram of "the average driving times". ... >very well approximated by a normal distribution, ...
    (sci.math)
  • Re: Probability Problem 2
    ... Given three independent draws from a normal distribution, ... You *can't* assume it to be at zero, ... You all think I'm paranoid, ...
    (rec.puzzles)
  • Geometrical analogies to statistical distribution functions?
    ... I just realized this in bed, but the normal distribution corresponds to the ... cross-sectional height of a parallelogram. ... The histogram of the sum of two dice approximates the normal ...
    (sci.math)
  • Re: Normal Distribution
    ... close to zero and the coefficient of kurtosis is close to 3, ... safe assuming the sample is from a normal distribution with the sample mean ... If you use Excel's KURT function, the online help shows an approximate ... indicates approximate normality. ...
    (microsoft.public.excel.worksheet.functions)