%
% 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 = 1;
b = 2; 

% store the desired distribution we will draw from (proportional to a beta distribution): 
vm = inline( '(1/x)*exp( -0.5*(log(x)^2) )', 'x' );

% generate "n" iterations of the chain: 
n = 2000; 
% this is our value of a "burn-in" parameter: 
b = 100; 

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

% get a starting value: 
X(1) = 0.1; 

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

figure; plot( X(b:end), 'x' ); 
xlabel('iteration value'); ylabel('X(i)'); 
saveas( gcf, '../../WriteUp/Graphics/Chapter11/prob_11_13_chain', 'epsc' ); 
figure; 
%hist( X, 20 ); 
cshistden( X(b:end), 50 ); hold on; 
title( 'univariate density histogram' ); 
x = linspace(0,10,500); y = lognpdf(x,0,1); plot( x, y, '-og' ); axis( [0,10,0,+Inf] ); 
xlabel('x value'); ylabel('pdf value'); 
saveas( gcf, '../../WriteUp/Graphics/Chapter11/prob_11_13_hist', 'epsc' ); 

clear functions;