% A Matlab program illustrating Gibbs sampling from a bivariate normal 
% distribution, based on Gelman, Carlin, Stern, and Rubin, Bayesian 
% Data Analysis, 1st ed. pp. 327-28 or 2nd ed. pp. 288-89.  In this 
% example we have a single observation (y_1, y_2) from a bivariate normal 
% distribution with known covariance matrix and unknown mean vector 
% (th(1), th(2)). The covariance matrix has ones on its principal 
% diagonal and rho elsewhere.  
clear
y_1 = 0; y_2 = 0;  % the observed values of y_1 and y_2
rho = .8;          % correlation coefficient 
M = 500;           % number of samples to be drawn
th = zeros(M,2);   % create an Mx2 matrix to store draws
th(1,1)=2.5; th(1,2)=-2.5;  % initial values of parameters
% Now we begin Gibbs sampling, using formulas based on the 
% bivariate normal distribution.
for ii=2:M                   
 th(ii,1)=y_1+rho*(th(ii-1,2)-y_2)+sqrt(1-rho^2)*randn(1,1);
 th(ii,2)=y_2+rho*(th(ii,1)-y_1)+sqrt(1-rho^2)*randn(1,1);
end;
plot(th(:,1),th(:,2),'-')   % plot path through parameter space
xlabel('\theta_1','fontsize',14)
ylabel('\theta_2','fontsize',14)
pause
scatter(th(251:500,1),th(251:500,2),'filled') % plot last 250 draws
xlabel('\theta_1','fontsize',14)
ylabel('\theta_2','fontsize',14)