用Metropolis-Hastings方法在Matlab中求积分

Question

在Matlab中，我在使用Metropolis-Hasting方法计算积分时遇到了一些问题。积分是从零到无穷大的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

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

Answer 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))