%
% 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); 

mu      = [1,2].'; sigma      = [1,0.9;0.9,1];  % <- for the target distribution ...
mu_prop = [0,0].'; sigma_prop = [0.6,0;0,0.4];  % <- for the proposal distribution ... 

% generate "n" iterations of the chain: 
n = 10000;
% 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 = X(i-1,:) + mvnrnd(mu_prop,sigma_prop);
  % 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_12_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_12_xi_pdfs', 'epsc' ); 

clear functions;