metropolis.m

% 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