%%
%MH MCMC for bivariate normal example
clear;
mu=[5;5];
ssq=3;
rho=-2/3;

%%
J=2000;
a=5;
x=[0;0];
X=zeros(2,J);

%%
for j=1:J
    xp=x+a*(2*rand(2,1)-1);
    alpha=min(1,p(xp,mu,ssq,rho)/p(x,mu,ssq,rho));
    if (rand<alpha), x=xp; end
    
    X(:,j)=x;
end

mean(X')'
cov(X')

%%
figure(1);
plot(X');
figure(2);
plot(X(1,:),X(2,:),'.'); axis([-2,12,-2,12]); axis equal; axis square; hold on;
m=30; r=linspace(-2,12,m); z=zeros(m,m); w=z; y=z;
for i=1:m, for j=1:m, 
        w(i,j)=r(i); y(i,j)=r(j); z(i,j)=p([r(i);r(j)],mu,ssq,rho); 
    end; end
contour(w,y,z,5,'k');


%%
%Gibbs sampler for bivariate normal example
clear;
mu=[5;5];
ssq=3;
rho=-2/3;

%%
J=1000;
x=[0;0];
X=zeros(2,J);

%%
for j=1:J
    x(1)=mu(1)+rho*(x(2)-mu(2))+randn*sqrt((1-rho^2)*ssq);
    x(2)=mu(2)+rho*(x(1)-mu(1))+randn*sqrt((1-rho^2)*ssq);
    X(:,j)=x;
end

mean(X')'
cov(X')