Contents

% -------------------------------------------------
%    Toy example of Hastings algorithm
%
% This example shows how to sample a 2D normal distribution by MCMC
% sampling method. Specifically, we use standard normal distribution as
% proposal distribution. The target distribution is a correlated 2D normal
% distribution. The proposal distribution here is symmetric. So it reduces
% to Hastings algorithm. If it is not symmetric, it would be the Metropolis
% -Hastings algorithm.
%
% written by Zhiyuan Weng
% -------------------------------------------------

clear
clc

Initialize

N = 70; % number of movement

% target distribution
mu  = [2,2];
sig = [2,0.5
       0.5,2];
p = @(x) 1/(2*pi)/sqrt(det(sig))*exp(-.5*(x-mu)/(sig)*(x-mu)');

% proposal distribution
% it is normal with
% covariance is b^2*eye(2) where b is
b = 1.2;

% pre allocate memory
x_out = zeros(N,2);
rej_flag = zeros(N,1);
x_rej = zeros(size(x_out));

main loop

x_out(1,:) = randn(1,2)*b + mu; % start point
for i = 2:N
	y = randn(1,2)*b + x_out(i-1,:); % sample from q ~ N(x(n-1),b^2)
	alpha = min(p(y)/p(x_out(i-1,:)),1);
	u = rand; % sample u from U[0,1)
	if(u<=alpha) % accept
		x_out(i,:) = y;
		rej_flag(i) = 1;
    else % reject
		x_out(i,:) = x_out(i-1,:);
		x_rej(i,:) = y;
	end
end

plot the data

% plot Gaussian pdf contour
[X,Y]  =meshgrid(-2:0.1:6);
Z = zeros(size(X));
i = 0;
for u = -2:0.1:6
    i = i + 1;
    j = 0;
    for v = -2:0.1:6
        j = j + 1;
        Z(i,j) = p([u,v]);
    end
end
contour(X,Y,Z);
hold on;

% plot moving route
for i = 2:N
    s0 = [x_out(i-1,1),x_out(i,1)];
    t0 = [x_out(i-1,2),x_out(i,2)];
    plot(s0,t0,'-.or','LineWidth',1);
    if(rej_flag(i) == 0)
        s = [x_out(i-1,1),x_rej(i,1)];
        t = [x_out(i-1,2),x_rej(i,2)];
        h = plot(s0,t0,'-.or',x_rej(i,1),x_rej(i,2),'+');
        plot(s,t,'g','LineWidth',1);
    end
end
legend(h,'accepted move','rejected move');
axis([-2,6,-2,6])
axis square