% RECITATION WK2:
% Continuous conjugate models.
%
% Jason L. Pacheco
% 2/9/11
%

randn('seed',sum(100*clock));

%
% NORMAL DENSITY
%
dx = 0.1;
x = -10:dx:10;
p = normpdf( x, 0, 1 );

% remember to normalize
p = p ./ sum( p ) ./ dx;
plot( x, p, '-b', 'LineWidth', 5 );

% we can do this analytically too
p = 1/sqrt(2*pi*1) * exp( -0.5 .* x.^2 ./ 1 );
p = p ./ sum( p ) ./ dx;
plot( x, p, '-k', 'LineWidth', 5 );

% now the CDF
phi = normcdf( x );     % Defaults to mu=0, sigma=1 
plot( x, phi, '-b', ... % No normalizing!  Not a density.
  'LineWidth', 5 );   

% numerical CDF
plot( x, cumsum( p ) * dx, '-k', ...  % Make sure to multiply by dx
  'LineWidth', 5);
close all;

%
% NORMAL-GAMMA (w/ known mean, unknown precision)
%
close all;
figure;
hold on;

% Prior
dx = 0.1;
t = 1:dx:10;
a = 2; b = 2;
prior = gampdf( t, a, b );
prior = prior / sum( prior ) / dx;
log_prior = log( prior );

% plot
plot( t, prior, '-b', 'LineWidth', 5 );
xlabel('Predcition (\tau)','FontSize',14);

% Sample i.i.d ~ N(0,t)
n = 10;
mu = 1;
t_true = 1/3;
X = mu + sqrt( 1/t_true ) * randn( [1,n] );
% load('data.mat','X');

% Likelihood
nu = 1/n * sum( X(:) );
like = t.^(n/2) .* exp( -1/2 .* t .* n .* nu );
log_like = log( like );
plot( t, like, '--b', 'LineWidth', 5 );

% Posterior
log_post = log_prior + log_like;
max_val = max( log_post );
log_post = log_post - max_val;
post = exp( log_post ) ./ sum( exp( log_post ) ) ./ dx;
plot( t, post, '-r', 'LineWidth', 5);
legend('Prior p(\tau)', 'Likelihood p(X|\tau)', 'Posterior p(\tau|X)' );

% Mode/Mean
t_mean = t * post' * dx;
[mode,I_mode] = max( post );
t_map = t( I_mode );
fprintf('Posterior Mean: %f\n', t_mean );
fprintf('MAP: %f\n', t_map );

% Predictions
x = -10:dx:10;
x_map = normpdf( x, mu, t_map );
x_mean = normpdf( x, mu, t_mean );
sig = 1; nu = 2*a;
x_pred = ( 1 + 1/nu * ( (x-mu)/sig ).^2 ).^(- (nu+1)/2 ) ;
x_pred = x_pred / sum( x_pred ) / dx;

% plot prediction distributions
figure;
plot( x, x_map, '-b', x, x_mean, '--b', x, x_pred, '-k', ...
  'LineWidth', 5 );
title('Prediction Distributions');
xlabel('X');
legend('MAP Plug-in', 'Posterior Mean plug-in', 'Posterior Predictive' );








