Metropolis-Hastings在matlab中的应用

Question

我正在尝试使用Metropolis Hastings算法和随机游走采样器来模拟matlab中$$函数的样本，但我的代码出现了一些问题。建议密度是椭圆2s^2 + 3t^2≤1/4上的均匀PDF。我可以使用拒绝验收的方法从建议密度中抽样吗？

N=5000;
alpha = @(x1,x2,y1,y2) (min(1,f(y1,y2)/f(x1,x2)));
X = zeros(2,N);
accept = false;
n = 0;
while n < 5000
    accept = false;
    while ~accept
        s = 1-rand*(2);
        t = 1-rand*(2);
        val = 2*s^2 + 3*t^2;
        % check acceptance
        accept = val <= 1/4;
    end
    % and then draw uniformly distributed points checking that u< alpha?
    u = rand();
    c = u < alpha(X(1,i-1),X(2,i-1),X(1,i-1)+s,X(2,i-1)+t);
    X(1,i) = c*s + X(1,i-1);
    X(2,i) = c*t + X(2,i-1);
    n = n+1;
end
figure;
plot(X(1,:), X(2,:), 'r+');

Answer 1

您可能只想使用matlab mhsample的本机实现。

关于你的代码，有一些东西遗漏了：-函数alpha，-循环变量i (它可能只是n，但它不适合索引，因为它从零开始)。如果你想动态填充内存，你应该总是在matlab中分配内存，也就是你的例子中的X。

Answer 2

为了扩展@max的建议，如果您将i索引更改为n并替换，则代码似乎可以工作

n = 0;

使用

n = 2; X(:,1) = [.1,.1];

将X(:,1)分配给接受区域内的随机值(使用与稍后使用的代码相同的代码)和/或包含一个老化时间段可能会更好。

根据您要做的事情，在f函数中对sin的参数求值以使其保持在0到2 pi内(如果超出这些界限，则可能将值移位2 pi )，这可能会使事情变得更清晰