%
% Example on epage 466
% 
% 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); 

% Set up preliminaries.
%
% Here we use k for the chain length and b as the burn-in parameter:
%
k = 3000;   % generate a chain of size k
b = 100;    % burn-in will be b 
B = 10;     % the upper bound of the values of $x_{,1}$ and $x_{,2}$ in the marginal distributions 

nChains = 3;             % run this many chains ...  
X1 = zeros(nChains,k);   % this will represent x_1
X2 = zeros(nChains,k);   % this will represent x_2 
m  = zeros(nChains,k);   % mean
v  = zeros(nChains,k);   % variance
s  = zeros(nChains,k);   % skewness

rhat_m = zeros(1,k); rhat_v = zeros(1,k); rhat_s = zeros(1,k);

% specify some relativly "disparate" starting points: 
% ... while initializing each chain: 
X1(1,1) = 0.1*B; X2(1,1) = min(exprnd(1/X1(1,1),1,1),B);
X1(2,1) = 0.5*B; X2(2,1) = min(exprnd(1/X1(2,1),1,1),B);
X1(3,1) = 0.7*B; X2(3,1) = min(exprnd(1/X1(3,1),1,1),B);

% run each chain ... 
for i = 2:k
  for nc=1:nChains,
    X1(nc,i) = min(exprnd(1/X2(nc,i-1),1,1),B);
    X2(nc,i) = min(exprnd(1/X1(nc,i-1),1,1),B);
  end
  m(:,i) = mean(X1(:,1:i).').';
  v(:,i) = var(X1(:,1:i).').'; 
  s(:,i) = skewness(X1(:,1:i).').'; 
  rhat_m(i) = csgelrub(m(:,1:i));
  rhat_v(i) = csgelrub(v(:,1:i));
  rhat_s(i) = csgelrub(s(:,1:i));
end

% plot the Gelman-Rubin convergence criterion for {\em each} of the statistics
figure; 
pm=plot( rhat_m, 'x-r' ); grid on; hold on; 
pv=plot( rhat_v, 'o-g' ); 
ps=plot( rhat_s, 'd-k' ); axis( [0,k,0.6,1.6] ); 
legend( [pm,pv,ps], {'GR mean', 'GR variance', 'GR skewness'}, 'location', 'northeast' ); 
title( 'gelman-rubin convergence diagnostic for the mean' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter11/prob_11_5_gelman_rubin_conv_m', 'epsc' ); 

% estimate the marginal distribution using each chain: 
% the marginal f(x1) is obtained by averaging the conditional f(x1|x2)
%
n_x1_samps = 500; fhat = zeros(nChains,n_x1_samps); 

x1_sample = linspace(0,B,n_x1_samps);
for nc=1:nChains
  for i = 1:length(x1_sample),
    fhat(nc,i) = mean( exppdf(x1_sample(i),1./X2(nc,b:k)) );
  end
end

figure; hold on; 
p1=plot( x1_sample, fhat(1,:), '-xr' ); 
p2=plot( x1_sample, fhat(2,:), '-og' ); 
p3=plot( x1_sample, fhat(3,:), '-.k' ); axis( [0,1,0,Inf] ); 
xlabel( 'X_1' ); ylabel( 'marginal PDF' ); 
title( 'estiamtes of the marginal distribution X_1' ); 
legend( [p1,p2,p3], {'X1(0)=0.1*B', 'X1(0)=0.5*B', 'X1(0)=0.7*B'},  'location', 'northeast' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter11/prob_11_5_x1_marg', 'epsc' );