% Gibbs sampling demo

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

C2=[0.2 0; 0 0.2];
Ci2=inv(C);
m2=[1.5 -1.5];

clf;

x=[-1 1];

px1=exp(-(x-m)*Ci*(x-m)'/2);
px2=exp(-(x-m2)*Ci2*(x-m2)'/2);
px=0.5*px1+0.5*px2;

plot_gaussian(sqrt(2)*C,m,2,60);
plot_gaussian(sqrt(2)*C2,m2,2,60);
axis([-3 3 -3 3]);
axis equal;
set(gcf,'Renderer','zbuffer');
pause;
s=0.1

for i=1:200,
xold=x;

px1=exp(-(xold-m)*Ci*(xold-m)'/2);
px2=exp(-(xold-m2)*Ci2*(xold-m2)'/2);
pxold=0.5*px1+0.5*px2;

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


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;