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

hold on;

% plot standard-deviation contour
[V,D] = eig(sig);
if((sin(acos(V(1,2)))-V(2,2))<0.0001)
    theta = rad2deg(acos(V(1,2)));
elseif((sin(acos(V(1,2)))+V(2,2))<0.0001)
    theta = rad2deg(-acos(V(1,2)));
end
ecc = sqrt(axes2ecc(D(2,2),D(1,1)));
[elat,elon] = ellipse1(mu(1),mu(2),[D(2,2) ecc], theta );
plot(elat,elon,'-.');

% 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

Published with MATLAB® 7.10