%
% Example on epage 445 
% 
% Problem on epage 471
% 
% 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); 

mu = [1,2].'; sigma = [1,0.9;0.9,1]; 

% generate "n" iterations of the chain: 
n = 6000;
% remove the initial "b" variates from our chain: 
b = 4000; 
% Set up the vectors to store the samples.
X = zeros(n,2);
% Get the starting values for the chains.
X(:,1) = 1.0;

% the "sigma" of the uniform random variate:
sig = 2; 

% Run the first chain.
for i = 2:n
  % Generate variate from proposal distribution.
  y(1) = unifrnd(X(i-1,1)-0.5*sig,X(i-1,1)+0.5*sig);
  y(2) = unifrnd(X(i-1,2)-0.5*sig,X(i-1,2)+0.5*sig);
  % Generate variate from uniform.
  u = rand(1);
  % Calculate alpha.
  alpha = min( [ 1, mvnpdf(y.',mu,sigma)/mvnpdf(X(i-1,:).',mu,sigma) ] ); 
  if u <= alpha
    % Then set the chain to the y.
    X(i,:) = y;
  else
    X(i,:) = X(i-1,:);
  end
end

figure; plot( X(b:end,1), X(b:end,2), 'x' ); 
xlabel('X(1)'); ylabel('X(2)'); axis( [-2 4 -1 5] ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter11/prob_11_11_scatter', 'epsc' ); 

figure; 
subplot( 1,2,1 ); cshistden( X(b:end,1), 20 ); 
subplot( 1,2,2 ); cshistden( X(b:end,2), 20 ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter11/prob_11_11_xi_pdfs', 'epsc' ); 

clear functions;