% Demo about Gibbs sampling;
% Target distribution is a 2D Gaussian distribution;
% Proposal distribution is the conditional distribution for each variable,  
% which is a 1D Gaussian distribution;
% Author: Nan Li
% Date: 12/02/2010

close all;clear;clc;

% target distribution: 2D Gaussion distribution
mu1 = 0;
mu2 = 0;
sigma1 = 1;
sigma2 = 2;
rho = 0.5;
MU = [0 0];
SIGMA = [1 1; 1 4];

nsamples = 10000;
x = zeros(nsamples,2);
% start point x1: (0,0)
for i = 2:nsamples
    x(i,1) = normrnd(mu1+(x(i-1,2)-mu2)*rho*sigma1/sigma2,...
        (1-rho*rho)*sigma1*sigma1);
    x(i,2) = normrnd(mu2+(x(i,1)-mu1)*rho*sigma2/sigma1,...
        (1-rho*rho)*sigma2*sigma2);
end

% evaluations of the target distribution
[X,Y] = meshgrid(-8:.5:8, -8:.5:8); 
row = size(X,1);
col = size(X,2);
zColumn = mvnpdf([X(:),Y(:)],MU,SIGMA);                                    
Z = reshape(zColumn,row,col);

% show the target distribution
figure();
surf(X,Y,Z);
xlabel('x1'); ylabel('x2');zlabel('pdf');
title('Target Probability Distribution');
colormap;
grid on;

% show bivariate histograms of samples
figure();
hist3(x(1:1000,:),{-8:.5:8, -8:.5:8});
xlabel('x1'); ylabel('x2');zlabel('samples');
title('1000 samples');
colormap;
grid on;

figure();
hist3(x(1:5000,:),{-8:.5:8, -8:.5:8});
xlabel('x1'); ylabel('x2');zlabel('samples');
title('5000 samples');
colormap;
grid on;

figure();
hist3(x(1:10000,:),{-8:.5:8, -8:.5:8});
xlabel('x1'); ylabel('x2');zlabel('samples');
title('10000 samples');
colormap;
grid on;