% Gibbs 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(sqrt(2)*C,m,2,60);
axis([-3 3 -3 3]);
set(gcf,'Renderer','zbuffer');
pause;
s=0.1

for i=1:200,
xold=x;
pxold=exp(-(xold-m)*Ci*(xold-m)'/2);

x=x+randn(1,2)*sqrt(s);
px=exp(-(x-m)*Ci*(x-m)'/2);

if rand<min(1,px/pxold), % accept
plot(xold(1),xold(2),'wo');
plot(x(1),x(2),'b.');
hold on;
plot(x(1),x(2),'bo');
fprintf('a');
else
plot(x(1),x(2),'r.');
fprintf('r');
x=xold;
end;

drawnow;
pause(0.4);
end;