2016-11-02 5 views
1

Matlabの積分を評価するMetropolis-Hastingの方法を使用する際に問題があります。積分は、ゼロから無限大までのe ^(x^-2)です。私が書いたコードではエラーは発生しませんが、 1)それは私がやっていることをしているのかどうかはわかりません 2)私が望むことをしても、 「抽出」は」データからの積分値がコードは、この問題は、おそらく非常にシンプルで、私はこの方法についての基本的な何かを見逃していることのように私は感じMetropolis-Hastingsメソッドを使用してMatlabの積分を評価する

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 

を生成します。私のプログラミングスキルは非常に基本的なので、どんな助けも大いに評価されます!

答えて

0

あなたの関数では、関数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)) 
関連する問題