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

% constants of our distribution: 
a = 5; b = 1;

% store the desired distribution we will draw from (proportional to a gamma distribution): 
vm = inline( '(x^(a-1))*exp(-b*x)', 'x', 'a', 'b' ); 

% generate "n" iterations of the chain: 
n = 10000; 

% set up some memory: 
X = zeros(1,n); 

% get a starting value for our chain: 
X(1) = 1; 

for i=2:n,
  % generate a variable from the proposal distribution: 
  y = unifrnd(0,20);
  % calculate alpha: 
  alpha = min([ 1, vm(y,a,b)/vm(X(i-1),a,b) ]); 
  if( rand(1,1) < alpha )
    X(i)=y;             % set the chain to y
  else
    X(i)=X(i-1); 
  end
end

figure; plot( X, 'x' ); 
xlabel('iteration value'); ylabel('X(i)'); 
saveas( gcf, '../../WriteUp/Graphics/Chapter11/prob_11_3_chain', 'epsc' ); 
figure; 
%hist( X, 20 ); 
cshistden( X, 20 ); hold on; 
title( 'univariate density histogram' ); 
x = linspace(0,20,200); y = gampdf(x,a,b); plot( x, y, '-og' ); 
xlabel('x value'); ylabel('pdf value'); 
saveas( gcf, '../../WriteUp/Graphics/Chapter11/prob_11_3_hist', 'epsc' ); 

clear functions;