2

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

4

1 回答 1

1

作为参考，@ripegraph 在另一个站点上解决了这个问题，我从拉普拉斯分布生成随机数的方法不正确，实际上应该使用：http ://en.wikipedia.org/wiki/Laplace_distribution#Generating_random_variables_according_to_the_Laplace_distribution

他还注意到拉普拉斯分布是对称的，因此根本不需要包含在代码中。

在做了更多的研究之后，我发现如果你有 X~Gamma(v/2, 2) 它变成了 X~ChiSquare(v) 并且是使用非高斯提议的一个更好的例子。但是，要使用此示例，您需要使用独立采样器http://www.math.mcmaster.ca/canty/teaching/stat744/Lectures5.pdf（幻灯片 89）。

希望这可能对某人有用。

于 2013-03-29T12:41:55.283 回答