%
% Example on epage 445 
% 
% Problem on epage 473
% 
% Written by:
% -- 
% John L. Weatherwax                2008-02-20
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

clear all; 
close all; 
clc; 

addpath('../../Code/CSTool');

randn('seed',0); rand('seed',0); 

% Set up some constants and arrays to store things.
n  = 5000;
bi = 500; 
xar = zeros(n,3);                        % to store samples 

% Set up the function to evaluate alpha for this problem. Note that the constant has been canceled.
strg = 'exp(-0.5*(x-mu)''*inv(covm)*(x-mu))';
norm = inline(strg,'x','mu','covm');

% Generate starting point.
xar(1,:) = [10,10,10].';
for i = 2:n,
  % Get the next proposd variate for this chain.
  y = mvnrnd( xar(i-1,:), 9*eye(3) ).'; 
  num = norm(y,zeros(3,1),eye(3))*mvnpdf(xar(i-1,:).',y,9*eye(3));
  den = norm(xar(i-1,:).',zeros(3,1),eye(3))*mvnpdf(y,xar(i-1,:).',9*eye(3)); 
  alpha = min([1, num/den]);
  if rand(1) <= alpha
    xar(i,:) = y;
  else
    xar(i,:) = xar(i-1,:);
  end
end
%

figure; plot( xar, 'x' ); 
title( 'samples from all chains' ); 
xlabel('iteration value'); ylabel('X(i)'); 
saveas( gcf, '../../WriteUp/Graphics/Chapter11/prob_11_6_chains', 'epsc' ); 

figure; 
%hist( X, 20 ); 
cshistden( xar(bi:n,1), 50 ); hold on; 
title( 'univariate density histogram (all samples)' ); 
xlabel('x value'); ylabel('pdf value'); 
saveas( gcf, '../../WriteUp/Graphics/Chapter11/prob_11_6_x_1_hist', 'epsc' ); 

clear functions;