% -----------------------------
% metrohast.m
% Metropolis-Hastings algorithm
% -----------------------------
clear all; close all;
x = rand(100,2).*0.1;
psi = @(x) exp(-(x(1)-2).^2-(x(2)-4).^2);
t = 0;
s = 1;
figure (1);
plot(x(:,1),x(:,2),'bo');  
old = psi(x);
ok = 0;
for tries = 1: 100000
    r = randn(100,2);
    xp = x + r.*s;
    new = psi(xp);
    ratio = new/old;
    if(ratio>1) 
           ok = ok+1;
           x = xp;
           old = new;
    else
        tmp =rand(1);
        if(ratio>tmp)
            ok = ok + 1;
            x = xp;
            old = new;
        end
    end
    accratio = ok/tries*100;
    if accratio<10
        s = 0.9*s;
    end
    if accratio >90
        s = 1.1*s;
    end
    printf("%f %f %f %f\n",old, new, s,accratio)
    
    if(mod(tries,10)==0) 
	plot(x(:,1),x(:,2),'bo'); 
	drawnow; 
    end	
    %input 'cont'
end
