% Demo about Metropolis-Hastings sampling with random walk chain;
% Implementation is based on MATLAB function: mhsample;
% Target distribution is a combination of 2 2D Gaussian distribution;
% Proposal distribution is 1 2D Gaussian distribution, which centers on the
% previous accepted sample;
% Author: Nan Li
% Date: 12/02/2010

close all;clear;clc;
% two 2D Gaussian distribution for the target distribution
MU1 = [-1 -1];
SIGMA1 = [1 0; 0 4];
MU2 = [1 1];
SIGMA2 = [4 0; 0 1];
% the covariance matrix of the proposal distribution
SIGMA3 = [4 0; 0 4];

% function handles of target distribution density, proposal distribution
% density and random number generator for the proposal distribution
pdf = @(x) 0.3*mvnpdf(x,MU1,SIGMA1) + 0.7*mvnpdf(x,MU2,SIGMA2);
proppdfRandom = @(x,y) mvnpdf(x,y,SIGMA3);
proprndRandom = @(x) mvnrnd(x,SIGMA3);

startPoint = [0,0];
nsamples = 10000;

xR = mhsample(startPoint,nsamples,'pdf',pdf,'proprnd',proprndRandom,...
    'proppdf',proppdfRandom);

% xRS = mhsample(startPoint,nsamples,'pdf',pdf,'proprnd',proprnd,'symmetric',1);

% evaluations of the target distribution
[X,Y] = meshgrid(-8:.5:8, -8:.5:8);
row = size(X,1);
col = size(X,2);
zColumn = 0.3*mvnpdf([X(:),Y(:)],MU1,SIGMA1) + 0.7*mvnpdf([X(:),Y(:)],MU2,SIGMA2);                                    
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(xR(1:1000,:),{-8:.5:8, -8:.5:8});
xlabel('x1'); ylabel('x2');zlabel('samples');
title('1000 samples');
colormap;
grid on;

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

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