在 Matlab 中使用 Metropolis-Hastings 方法计算积分答案

【问题标题】：Using Metropolis-Hastings method to evaluate an integral in Matlab在 Matlab 中使用 Metropolis-Hastings 方法计算积分
【发布时间】：2016-11-02 11:22:39
【问题描述】：

我在使用 Metropolis-Hasting 方法在 Matlab 中评估积分时遇到了一些麻烦。积分是从零到无穷大的 e^(x^-2)。我写的代码没有产生错误，但是， 1）我不太确定它是否在做我希望它做的事情 2）即使它做了我想要的，我也不太确定如何从代码产生的数据中“提取”积分值

clc
clear all

%Parameters for Gaussian proposal distribution N(mu, sigma)
sigma = 1;
mu = 0;

f = @(x) exp(x.^-2); %Target distribution

n = 10000;
step = 1;
x = zeros(1, n); %Storage
x(1) = 1; %Starting point, maximum of function f

for jj = 2:n
xtrial = x(jj-1) + step*normrnd(mu,sigma); %Generates candidate

w = f(xtrial)/f(jj-1);

if w >= 1
    x(jj) = xtrial;
else
    r = rand(1); %Generates uniform for comparison
    if r <= w
        x(jj) = xtrial;
    end
    x(jj) = x(jj-1);
end
end

我觉得这个问题可能很简单，我只是错过了有关此方法的一些基本内容。任何帮助将不胜感激，因为我的编程技能非常基础！

【问题讨论】：

标签： matlab mcmc

【解决方案1】：

您的函数也未在 x = 0 处定义，在您编写的代码中，函数 f 的最大值在 x = 1 处，所以我假设积分是从 1 到 inf。积分的值是 x 的平均值，为此我得到了 2.7183 的结果，代码如下：

    clc
    clear all

    %Parameters for Gaussian proposal distribution N(mu, sigma)
    sigma = 3;
    mu    = 0;
    f     = @(x) double(x>1).*exp(x.^-2); %Target distribution
    n     = 10000;
    step  = 1;
    x     = zeros(1, n); %Storage
    x(1)  = 1;           %Starting point, maximum of function f
    acc   = 0;           % vector to track the acceptance rate

    for jj = 2:n
    xtrial = x(jj-1) + step*normrnd(mu,sigma); %Generates candidate

    w = f(xtrial)/f(x(jj-1));

    if w >= 1
        x(jj) = xtrial;
        acc = acc +1;
    else
        r = rand(1); %Generates uniform for comparison
        if r <= w
            x(jj) = xtrial;
            acc = acc +1;
        end
        x(jj) = x(jj-1);
    end
    end
    plot(x)
    (acc/n)*100
    mean(f(x))

【讨论】：