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

trueMean = 0.06; 

% generate our "n" standard normals: 
n      = 20; 
z_samp = randn(1,n) + trueMean;
z_bar  = mean(z_samp); 

% store the distribution we desire to draw from this is \pi(\cdot)
% 
vm = inline( '(1/(1+t^2))*exp(-0.5*n*(t-z_bar)^2)', 't', 'n', 'z_bar' ); 
%vm = inline( 'exp(-0.5*n*(t-z_bar)^2)', 't', 'n', 'z_bar' );   % <- a "uniform" prior ... 

% generate this many chains: 
nChains = 5; 

% generate "m" iterations of the chain: 
m = 5000; 
%b = 1; 
%b = round(0.1*m);
b = 2500;      % <- based on the gelman-rubin convergence diagnostic ... 

% set up some memory: 
X = zeros(nChains,m); 
mus = zeros(nChains,m);   % mean
rhat_m = zeros(1,m);    % stores convergence criterion ... 

% get a starting value to draw samples from our posteriori from: 
X(:,1) = z_bar + randn(nChains,1);

for i=2:m,
  for ci=1:nChains,
    % generate a variable from the proposal distribution q take a wide enough uniform RV: 
    y = unifrnd(min(z_samp),max(z_samp));
    %y = unifrnd(0.0,+1.0);
    % calculate alpha: 
    alpha = min([ 1, vm(y,n,z_bar)/vm(X(ci,i-1),n,z_bar) ]); 
    if( rand(1,1) < alpha )
      X(ci,i)=y;             % set the chain to y
    else
      X(ci,i)=X(ci,i-1); 
    end
  end
  mus(:,i)  = mean(X(:,1:i).').';
  rhat_m(i) = csgelrub(mus(:,1:i));
end
t = X(:,b:end); % <- exclude the initial burn-in samples ... 
fprintf('the mean of all our samples from our posteriori distribution is given by %10.6f\n',mean(t(:)));
fprintf('the variance of all our samples from our posteriori distribution is given by %10.6f\n',var(t(:)));
%%%fprintf('the true mean is %10.6f the difference is %10.6f\n',trueMean,trueMean-mean(t(:)));

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 the mean' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter11/prob_11_9_gelman_rubin_conv_m', 'epsc' ); 

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

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

clear functions;