% Metropolis sampling demo

C=[1 0.98; 0.98 1];
Ci=inv(C);
m=[0 0];
clf;

x=[-1 1];
px=exp(-(x-m)*Ci*(x-m)'/2);

plot_gaussian(4*C,m,2,60);
axis([-3 3 -3 3]);
set(gcf,'Renderer','zbuffer');
pause;
s=0.1
s=input('width of proposal:');
nstep=input('how many steps:');

for i=1:nstep,
  xold=x;
  pxold=exp(-(xold-m)*Ci*(xold-m)'/2);
  
  x=x+randn(1,2)*sqrt(s); % propose a new data point
  px=exp(-(x-m)*Ci*(x-m)'/2); % compute its probability
  if i<20, pause(1); end;
  if rand<min(1,px/pxold), % accept (plot in blue)
    plot(xold(1),xold(2),'wo');
    plot(x(1),x(2),'b.');
    hold on;
    plot(x(1),x(2),'bo');
    fprintf('a');
  else % reject (plot in red)
    plot(x(1),x(2),'.','Color',[1 0.5 0.5]);
    fprintf('r');
    x=xold;
  end;

  drawnow;
  pause(0.1);
end;
     
