%
% 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 = 1000; 
nChains = 10; 
xar = zeros(nChains,n,3);                        % to store samples 

USE_EXP = 1; 

% Set up the function to evaluate alpha for this problem. Note that the constant has been canceled.
strg = 'exp(-( x1 + x2 + x3 + x1*x2 + x1*x3 + x2*x3 ));';
norm = inline(strg,'x1','x2','x3');

% Generate a starting point (all chains have the same starting point).
xar(:,1,:) = 1;
for i = 2:n,
  for ci=1:nChains,
    x   = squeeze(xar(ci,i-1,:)); 
    
    % Get the next proposed variate for this chain.
    % 
    if( USE_EXP ) 
      y   = exprnd(1,3,1);
      num = norm(y(1),y(2),y(3))*exppdf(x(1),1)*exppdf(x(2),1)*exppdf(x(3),1);     % <- the numerator and denominator in computing alpha ... 
      den = norm(x(1),x(2),x(3))*exppdf(y(1),1)*exppdf(y(2),1)*exppdf(y(3),1); 
    else
      y   = unifrnd(0,20,3,1);        % <- this proposal distribution is VERY BAD
      %y   = unifrnd(0,5,3,1);         % <- this proposal distribution is DECENT ... 
      num = norm(y(1),y(2),y(3));     % <- the numerator and denominator in computing alpha ... q cancils ...
      den = norm(x(1),x(2),x(3)); 
    end
  
    alpha = min([1, num/den]);
    if rand(1) <= alpha
      xar(ci,i,:) = y;
    else
      xar(ci,i,:) = x;
    end
    t1 = squeeze(xar(ci,1:i,:)).'; % <- [3,i] all samples from 1 to i
    t2 = prod(t1);                 % <- [1,i] compute the product of x_1 x_2 x_3
    m(ci,i) = mean(t2);            % <- compute the mean (expectation) to get E[X_1 X_2 X_3] of samples from 1:i 
  end %end chain loop 
  % lets compare how well we are approximating this expectation ... 
  rhat_m(i) = csgelrub(m(:,1:i));
end

% estimate the desired expectation using all chains: 
ex1x2x3 = zeros(1,nChains);
for ci=1:nChains,
  t1 = squeeze(xar(ci,bi:n,:)).'; % <- [3,i] ... and estimate this with just one chain ... 
  t2 = prod(t1);                  % <- [1,i]
  ex1x2x3(ci) = mean(t2); 
end
fprintf('the calculated value for E[x_1 x_2 x_3] is %10.6f\n',mean(ex1x2x3));

if( 1 ) 

  figure; plot( rhat_m, '-kx' ); grid on; hold on; 
  plot( 1.1*ones(size(rhat_m)), '-g' );  
  plot( 1.2*ones(size(rhat_m)), '-.g' ); 
  title( 'gelman-rubin convergence diagnostic for E[X_1 X_2 X_3]' ); 
  saveas( gcf, '../../WriteUp/Graphics/Chapter11/prob_11_7_gelman_rubin_conv_m', 'epsc' ); 

  figure; plot( squeeze(xar(:,bi:n,1)).', 'x' ); axis tight; 
  title( 'x_1 samples from all chains after burnin' ); 
  xlabel('iteration value'); ylabel('X(i)'); 
  saveas( gcf, '../../WriteUp/Graphics/Chapter11/prob_11_7_chains', 'epsc' ); 
  %%%saveas( gcf, '../../WriteUp/Graphics/Chapter11/prob_11_7_poor_uniform_xi', 'epsc' ); % <- with a very poor proposal distribution ... 

end

clear functions; 

return; 

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_7_x_1_hist', 'epsc' );