% TEXT page 74, 75, 76 Simulation from the NIG Posterior
%
rand('seed',1);
randn('seed',1);
%
sample = sqrt(exp(1)) * (randn(1,20)) + pi;
%----------------------
% "Observations" are 20 normals with mean
% mu = pi, and variance sigma^2 = e.
% Mean and variance are unknown and will
% be estimated by simulating from the bivariate 
% posterior...
%----------------------
close all;
n=length(sample);
xbar = mean(sample);      % sufficient statistics for mu
s2 = 1/(n-1) * sum( (sample - xbar).^2 ); %s2=std(sample)
% sufficient statistics for sigma^2 if mu is not known.
%
NN=10000;     %NN simulated realizations
theta1s=[];
theta2s=[];
for i = 1:NN
   theta2 = 1/(rand_gamma((n-1)/2, (n-1)*s2/2));
   theta1 = theta2/n * randn(1,1) + mean(sample);
   theta1s=[theta1s theta1];
   theta2s=[theta2s theta2];
end
figure(1) %scatterplot of (theta1, theta2) simulated 
% from the bivariate posterior
  plot(  theta1s, theta2s, '.')
  %----------------------------
figure(2)
subplot(2,1,1)
hist(theta1s, 70)
subplot(2,1,2)
hist(theta2s, 70)
% 
% Empirical and Theoretical Posterior means
sprintf(' Posterior sample mean for mu = %0.5g',mean(theta1s)) %posterior averages
sprintf(' Posterior theoretical mean for mu = %0.5g', xbar) %posterior averages
sprintf(' Posterior sample mean for sigma^2 = %0.5g',mean(theta2s)) %posterior averages
sprintf(' Posterior theoretical mean for sigma^2 = %0.5g',1/(n-3) * (n-1)*s2) %posterior averages
