2013-03-28 7 views
2

현재 Metropolis-Hastings 알고리즘 및 일부 수치 예제에 대한 개요를 기반으로하는 수학 학위 과정의 최종 프로젝트를 진행하고 있습니다. 지금까지 Gaussian으로 제안서 배포판을 사용하고 몇 가지 다른 배포판을 샘플링하여 좋은 결과를 얻었지만 다른 제안 배포판을 사용하여 한 걸음 더 나아가려고합니다.Metropolis-Hastings 알고리즘, Matlab의 가우시안 이외의 제안 분포 사용

지금까지 나는 다른 제안을 사용하는 것에 관해 온라인으로 제한된 리소스로이 코드를 사용했습니다 (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 루프가 완전히 잘못되었다는 것에 대한 제 세대에 대해서는 매우 조심 스럽습니다).

감사합니다, 톰 참고로

답변

1

,이, @ripegraph에 의해 다른 사이트에 해결 된 라플라스 분포 난수를 생성하기위한 나의 방법은 올바르지 실제로 사용하여 수행해야합니다 http://en.wikipedia.org/wiki/Laplace_distribution#Generating_random_variables_according_to_the_Laplace_distribution

그는 라플라스 분포는 대칭 적이기 때문에 코드에 전혀 포함될 필요가 없다는 점도 지적했다.

좀 더 연구를 한 후에 X ~ 감마 (v/2, 2)를 사용하면 X ~ ChiSquare (v)가되고 비 가우스 제안을 사용하는 훨씬 좋은 예가된다는 것을 알았습니다. 그러나이 예제를 사용하려면 독립 샘플러 http://www.math.mcmaster.ca/canty/teaching/stat744/Lectures5.pdf (슬라이드 89)을 사용해야합니다.

누군가에게 유용 할 수 있기를 바랍니다.