Metropolis-Hastings 방법을 사용하여 Matlab의 적분 값을 계산하십시오.

Question

Metropolis-Hasting의 방법을 사용하여 Matlab의 적분을 평가하는 데 약간의 문제가 있습니다. 적분은 0에서 무한대까지의 e^(x^-2)입니다. 내가 작성한 코드는 오류가 발생하지 않지만, 1) 내가하고 싶은 일을하고 있는지 잘 모르겠다. 2) 내가 원하는대로해도 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

당신의 코드에서 함수 f의 최대 값이 x = 1이므로 x = 0에서도 함수가 정의되어 있지 않으므로 통합은 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))