%
%Example MATLAB Code for running a simple Gibbs sampler.
%

%Set mu and phi(precision = 1/variance)
mu = 5; sigma2 = 2; phi = 1/sigma2;

%Draw a Sample of size 53 from the N(mu,sigma2) distribution
x = normrnd(mu,sqrt(sigma2),53,1);

%Write the data x to a file
dlmwrite('./sampledata.dat',x,',');

%Read the data file and save it as an object
mydata = load('./sampledata.dat');
xtot=sum(mydata);

%Prior Values
m = 0;  s2 = 10;  ph=1/s2; mdivs2=m/s2;
a = 2;  b = 3;
n = size(mydata,1);

%Allocate Space to Save Gibbs Sampling Draws
totdraws = 5000;
mudraws = zeros(totdraws,1); %Initial Value of mu is 0
sigma2draws = ones(totdraws,1);

%Begin the Gibbs Sampler
for i=2:totdraws

   %Generate a Draw for mu
   s2star = 1/(ph+n/sigma2draws(i-1,:));
   mstar = s2star*(mdivs2 + xtot/sigma2draws(i-1,:));
   mudraws(i,:) = normrnd(mstar,sqrt(s2star),1,1);

   %Generate a Draw for sigma2
   astar = a+n/2;
   bstar = b+sum((mydata-mudraws(i,:)).^2)/2;
   phitemp = gamrnd(astar,1/bstar,1,1);
   sigma2draws(i,:) = 1/phitemp;

end %End Gibbs Sampler For Loop

%Discard the first 500 draws as Burn-in
mudraws = mudraws(501:size(mudraws,1),:);
sigma2draws = sigma2draws(501:size(sigma2draws,1),:);

%Plot Trace Plots
subplot(2,2,1)
plot(mudraws,'k-'); line([0 totdraws-499],[mu mu])
subplot(2,2,2)
plot(sigma2draws,'k-'); line([0 totdraws-499],[sigma2 sigma2])

%Plot Marginal Posterior Histograms
subplot(2,2,3); hist(mudraws)
title('Marginal Posterior Histogram for \mu'); xlabel('\mu')
subplot(2,2,4); hist(sigma2draws)
title('Marginal Posterior Histogram for \sigma^2'); xlabel('\sigma^2')