2

私は現在、メトロポリス-ヘイスティングス アルゴリズムの概要といくつかの数値例を示すことに基づいて、数学の学位を取得するための最終年度のプロジェクトに取り組んでいます。これまでのところ、提案分布をガウスとして使用し、他のいくつかの分布からサンプリングすることでいくつかの素晴らしい結果を得ていますが、別の提案分布を使用してさらに一歩進んでいます。

これまでのところ、私はこのコードを取得しています (私は Matlab を使用しています) が、さまざまな提案の使用に関するオンラインのリソースが限られているため、実際にはこれを試みる方法がよくわからないため、私がまったく近いかどうかを判断するのは困難です (特にこれまでのところ、有用なデータ出力が得られないためです)。

誰かが私を知っているか、簡単にアクセスできる情報に転送してくれたら、手を貸してくれると素晴らしいでしょう (私はコーディングのアドバイスだけでなく、数学についても尋ねていることを理解しています)。

したがって、ラプラスの提案分布を使用してガウスからサンプリングしたいのですが、これまでのコードは次のとおりです。

n = 1000;       %%%%number of iterations

x(1) = -3;      %%%%Generate a starting point

%%%%Target distribution: Gaussian:

strg = '1/(sqrt(2*pi*(sig)))*exp(-0.5*((x - mu)/sqrt(sig)).^2)';
tnorm = inline(strg, 'x', 'mu', 'sig');

mu = 1;    %%%%Gaussian Parameters (I will be estimating these from my markov chain x)
sig = 3;


%%%%Proposal distribution: Laplace:

strg = '(1/(2*b))*exp((-1)*abs(x - mu)/b)';
laplace = inline(strg, 'x', 'b', 'mu');

b = 2;       %%%%Laplace parameter, I will be using my values for y and x(i-1) for mu


%%%%Generate markov chain by acceptance-rejection

for i = 2:n

    %%%%Generate a candidate from the proposal distribution
    y = laplace(randn(1), b, x(i-1));

    %%%%Generate a uniform for comparison
    u = rand(1);

    alpha = min([1, (tnorm(y, mu, sig)*laplace(x(i-1), b, y))/(tnorm(x(i-1), mu, sig)*laplace(y, b, x(i-1)))]);


    if u <= alpha
        x(i) = y;
    else
        x(i) = x(i-1);
    end 
end

上記が完全に間違っているか、間違った方法で行っているか、またはいくつかの間違いがあるかどうかを誰かが教えてくれれば(forループの「y」の生成が完全に間違っていることに非常に警戒しています)。素晴らしいこと。

ありがとう、トム

4

1 に答える 1

1

参考までに、これは@ripegraphによって別のサイトで解決されました。ラプラス分布から乱数を生成する私の方法は正しくなく、実際には以下を使用して行う必要があります。

また、ラプラス分布は対称であるため、コードに含める必要がまったくないことも指摘されました。

もう少し調査を行った後、X~Gamma(v/2, 2) がある場合、それは X~ChiSquare(v) になり、非ガウス提案を使用するより良い例であることがわかりました。ただし、この例を使用するには、独立サンプラーhttp://www.math.mcmaster.ca/canty/teaching/stat744/Lectures5.pdf (スライド 89)を使用する必要があります。

これが誰かに役立つことを願っています。

于 2013-03-29T12:41:55.283 に答える