%
% Example on epage 466
% 
% Problem on epage 472
% 
% 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 = 6000;   % generate a chain of size k
b = 200;    % burn-in will be b 

X1 = zeros(1,k); X1(1) = 0.5; 
X2 = zeros(1,k); X2(1) = 0.5; 
X3 = zeros(1,k); X3(1) = 0.5; 

% Start the Gibbs Sampler.
% 
for i = 2:k,
  x1 = X1(i-1); x2 = X2(i-1); x3 = X3(i-1); 
  
  y1 = (1-x2-x3)*betarnd(5,2);
  y2 = (1-y1-x3)*betarnd(4,2); 
  y3 = (1-y1-y2)*betarnd(3,2); 
  
  X1(i) = y1; X2(i) = y2; X3(i) = y3; 
end

figure; 
subplot( 1,3,1 ); p1=plot( X1(b:end), 'xr' ); ylabel( 'X_1' ); xlabel( 'sample index' ); axis tight; 
subplot( 1,3,2 ); p2=plot( X2(b:end), 'og' ); ylabel( 'X_2' ); xlabel( 'sample index' ); axis tight; 
subplot( 1,3,3 ); p3=plot( X3(b:end), '.k' ); ylabel( 'X_3' ); xlabel( 'sample index' ); axis tight; 
saveas( gcf, '../../WriteUp/Graphics/Chapter11/prob_11_8_xi_samps', 'epsc' );

figure; 
subplot( 1,3,1 ); hist( X1(b:end), 20 ); 
subplot( 1,3,2 ); hist( X2(b:end), 20 ); 
subplot( 1,3,3 ); hist( X3(b:end), 20 ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter11/prob_11_8_hists', 'epsc' );

return;