%% A simple Markov Chain
% Uses Metropolis-Hastings to draw from a mixture of bivariate normal
% densities.
% Uses a N(0,tau I) candidate density.  Runs the chain for various
% values of tau to compare results.
% From 14.384 2008 ps5
clear;
% set parameters
mu = {[1.5  1.5], [-1.5 -1.5]};
K = length(mu{1});
sig = {[1 0.5; 0.5 1], [1 0.5; 0.5 1]};
p = [0.5,0.5];
f = @(x) mvnpdf(x,mu{1},sig{1})*p(1) + mvnpdf(x,mu{2},sig{2})*p(2);
n = 2000;
tau = [0.1^2, 1, 8^2];
theta = zeros(length(tau),n,K);
for t=1:length(tau)
  theta(t,1,:) = [10,-10];
  fold = f(squeeze(theta(t,1,:))');
  muTheta(1,:) = mean(theta(t,1:1,:),2);
  varTheta(1,:) = var(theta(t,1:1,:),0,2);
  accept = 0;
  for i=2:n
    % draw candidate
    x = mvnrnd(squeeze(theta(t,i-1,:))',eye(K)*tau(t));
    u = rand();
    fnew = f(x);
    if (u<fnew/fold) 
      theta(t,i,:) = x;
      fold = fnew;
      accept = accept+1;
    else
      theta(t,i,:) = theta(t,i-1,:);
    end
    muTheta(i,:) = mean(theta(t,1:i,:),2);
    varTheta(i,:) = var(theta(t,1:i,:),0,2);
  end
  figure;
  ezcontour(@(t0,t1) f([t0 t1]),[-4,4])
  h1 = gca;
  h2 = axes('Position',get(h1,'Position'));
  plot(squeeze(theta(t,:,1)),squeeze(theta(t,:,2)),'.')
  set(h2,'Ylim',get(h1,'YLim'),'Color','none')
  set(h2,'XLim',get(h1,'XLim'),'Layer','top','Color','none')  
  fprintf('acceptance rate=%g\n',accept/n);
  figure;
  plot(muTheta)
  figure;
  plot(varTheta)
end
