% A Matlab program applying the Metropolis within Gibbs procedure to 
% sample from a bivariate normal distribution, based on
% Gelman et al., 1st ed. pp. 327-28 or 2nd ed. pp. 288-89 and
% Geweke and Keane (2001) p. 3479.
 
clear
y_1 = 0; y_2 = 0;
rho = .8;
M = 500;
th = zeros(M,2);
th(1,1)=2.5; th(1,2)=-2.5;
for ii=2:M
 th(ii,1)=y_1+rho*(th(ii-1,2)-y_2)+sqrt(1-rho^2)*randn(1,1);
  % Metropolis within Gibbs step
   thstar2 = th(ii-1,2)+2*randn(1,1); % proposal
   thold = [th(ii,1); th(ii-1,2)];
   mu = [y_1; y_2];
   Sigma = [1 rho; rho 1];
   % Joint pdf for (th(ii,1), th(ii-1,2))
   pold  = (2*pi)^(-1)*(det(Sigma))^(-1/2)*exp(-(thold -mu)'*(Sigma^(-1))*(thold -mu));
   thstar = [th(ii,1); thstar2];
   % Joint pdf for (th(ii,1), thstar2)
   pstar = (2*pi)^(-1)*(det(Sigma))^(-1/2)*exp(-(thstar-mu)'*(Sigma^(-1))*(thstar-mu));
   % proposal density for th(ii-1,2)
   qold = normpdf(th(ii-1,2),thstar2,2);
   % proposal density for thstar2
   qstar = normpdf(thstar2,th(ii-1,2),2);
   accept = min(pstar*qold/(pold*qstar),1); % probability of accepting proposal
   urn = rand(1); % draw from uniform distribution on (0,1)
   th(ii,2) = thstar2*(urn<=accept)+th(ii-1,2)*(urn>accept);
end;
plot(th(:,1),th(:,2),'-')
xlabel('\theta_1','fontsize',14)
ylabel('\theta_2','fontsize',14)
pause
scatter(th(251:500,1),th(251:500,2),'filled')
xlabel('\theta_1','fontsize',14)
ylabel('\theta_2','fontsize',14)