Metropolis-Hastings 算法，在 Matlab 中使用除高斯分布之外的提议分布答案

【问题标题】：Metropolis-Hastings algorithm, using a proposal distribution other than a Gaussian in MatlabMetropolis-Hastings 算法，在 Matlab 中使用除高斯分布之外的提议分布
【发布时间】：2013-03-28 18:30:16
【问题描述】：

我目前正在完成我的数学学位的最后一年项目，该项目基于对 Metropolis-Hastings 算法的概述和一些数值示例。到目前为止，通过将我的建议分布用作高斯分布并从其他几个分布中采样，我已经获得了一些很好的结果，但是我正试图通过使用不同的建议分布更进一步。

到目前为止，我已经得到了这段代码（我正在使用 Matlab），但是由于在线资源有限，关于使用不同的建议很难判断我是否接近，因为实际上我不太确定如何尝试这个, （特别是因为到目前为止这没有提供有用的数据输出）。

如果有人知道或转发给我一些易于访问的信息，如果有人能伸出援手，那就太好了（我知道我不仅在询问编码建议，还询问数学）。

所以，我想使用拉普拉斯的提议分布从高斯样本中采样，这是我目前的代码：

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

如果有人能告诉我上述是否完全错误/以错误的方式进行，或者只是有一些错误（我非常警惕我在 for 循环中生成的“y”是完全错误的) 那太棒了。

谢谢，汤姆

【问题讨论】：

标签： matlab math gaussian mcmc

【解决方案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）。

希望这可能对某人有用。

【讨论】：