matlab中的Metropolis-Hastings答案

【问题标题】：Metropolis-Hastings in matlabmatlab中的Metropolis-Hastings
【发布时间】：2020-01-03 00:16:08
【问题描述】：

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

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+');

【问题讨论】：

标签： matlab random mcmc random-walk

【解决方案1】：

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

关于您的代码，缺少一些内容： - 函数alpha, - 循环变量i（它可能只是n，但它不适合索引，因为它从零开始）。如果你想动态填充它，你应该总是在matlab中分配内存，即X在你的情况下。

【讨论】：

我更新了我的代码，但我现在看不出它有什么问题。我需要一个 for 循环吗？

【解决方案2】：

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

n = 0;

与

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

最好将X(:,1) 分配给您接受区域内的随机值（使用您稍后使用的相同代码），和/或包括一个老化期。

根据您要对此执行的操作，在 f 函数中评估 sin 的参数以将其保持在 0 到 2 pi 范围内（可能通过将值移动 2 pi 如果它超出了这些界限）

【讨论】：