%
% Gibbs sampling for the model X ~ N(mu,1/phi) where phi and mu have
% independent normal N(m,1/delta) and gamma G(a/2,b/2) priors respectively
% and a sample of size n is observed.
%
% -------------------------------------------------------------------------
% -------------------------------------------------------------------------
% PROGRAM
% -------------------------------------------------------------------------
% -------------------------------------------------------------------------
%
% Prior parameters.
%
m=100;
delta=0.01;
a=1;
b=1;
%
% Data and summary statistics.
%
x=[140,103,73,171,137,91,81,157,146,69,121,134];
n=size(x,2);
xbar=mean(x);
s2=var(x);
s=std(x);
%
% Gibbs sampler iterations.
%
iters=10000;
%
% Initial value for Gibbs sampler.
%
mu(1)=xbar;
%
% Some fixed parameters for speeding up algorithm.
%
par1=0.5*(a+n);
par2=b+(n-1)*s2;
%
% Grid for drawing the posterior density fmu for mu.
%
mumin=xbar-2*s;
mumax=xbar+2*s;
mupts=100;
mugrid=mumin+[0:mupts]*(mumax-mumin)/mupts;
fmu=zeros(1,(mupts+1));
%
% pmuc is P(mu <= c|data).
%
c=100;
%
% Gibbs algorithm.
%
for t=1:iters
    phi(t+1)=gamrnd(par1,2/(par2+n*((mu(t)-xbar)^2)));
    mu(t+1)=normrnd((delta*m+n*xbar*phi(t+1))/(delta+n*phi(t+1)),1/sqrt(delta+n*phi(t+1)));
    fmu=fmu+normpdf(mugrid,(delta*m+n*xbar*phi(t+1))/(delta+n*phi(t+1)),1/sqrt(delta+n*phi(t+1)));
end
%
% Diagram of mu against iterations.
%
figure
plot([1:(iters+1)],mu)
xlabel('t')
ylabel('\mu{(t)}')
xlim([1,(iters+1)])
%
% Diagram of convergence of the posterior mean of mu with the data.
%
emu=cumsum(mu)./[1:(iters+1)];
figure
plot([1:(iters+1)],emu)
xlabel('t')
ylabel('E[\mu|data]')
xlim([1,(iters+1)])
%
% Combined diagram of mu and mean vs iterations.
%
figure
plot([1:(iters+1)],mu)
xlabel('t')
ylabel('\mu{(t)}')
xlim([1,(iters+1)])
hold on
plot([1:(iters+1)],emu,'r')
hold off
%
% Histogram of mu's generated and normalized posterior density. Uses the
% package histadjusted to give a relative frequency histogram.
% 
fmu=fmu/iters;
figure
histadjusted(mu,20)
xlabel('\mu')
ylabel('f(\mu|data)')
hold on
plot(mugrid,fmu,'r')
hold off
%
% Final estimated mean of mu. 
% 
emu(iters+1)
%
% and probability that mu < 100.
%
pmuc=sum(normcdf(c,(delta*m+n*xbar*phi([2:(iters+1)]))./(delta+n*phi([2:(iters+1)])),...
    1./sqrt(delta+n*phi([2:(iters+1)]))))/iters
%
% A different estimate.  The proportion of sampled mu's < c.
%
pmuc2=sum(mu<100)/(iters+1)