% mcmc2.m
%-----------------------
% Y_ij = beta + U_i + epsilon_ij, i=1,...,I;
%                                 j=1,...,J.
% U_i ~ N(0, sigma^2)
% epsilon_ij ~ N(0, tau^2)
% pi(beta, sigma^2, tau^2) = 1/(sigma^2 tau^2)
% Jeffrey's Prior on beta, sigma, tau
% ------
clear all
close all
rand('state',2); randn('state',2);
%-----------------figure defaults
lw = 2;  
set(0, 'DefaultAxesFontSize', 17);
fs = 14;
msize = 5;
%------------simulate data----------
 I=5; J=10;
 eps = randn(I,J);
 u_i = 3*randn(I,1);
 beta = 7;
%------------------------------------
 for i=1:I
   for j=1:J
     data(i,j)=beta+u_i(i) + eps(i,j);
   end
 end
%-------------------------------------

%------------------------------------------
M = 20000; bin=100;
betas = [];
sigma2s=[];
tau2s = [];
gmeans = mean(data');
gmean = mean(gmeans);

tau2 = 17;      % = 3^2 in the simul
sigma2  = 5;    % = 1 in the simul
beta = 3;       % = 7 in the simul

%-------------------------------------------
    h=waitbar(0,'Simulation in progress');
for r = 1: M
 su=0;
 su2=0;
  for i = 1:I
           u_i(i) = 1/sqrt(J/sigma2 + 1/tau2) * randn + ...
           J*(gmeans(i)-beta)/(J + sigma2/tau2 );
           for j=1:J
                su =su + (data(i,j) - u_i(i) - beta)^2;
           end
           su2=su2+u_i(i)^2;
  end
  ubar = mean(u_i);
 beta = sqrt(sigma2/(I*J))  * randn + (gmean - ubar);
 tau2 = ( rand_gamma( I/2, su2/2 ) )^(-1);       
 sigma2   = ( rand_gamma( I*J/2, su/2 ) )^(-1);
 sigma2s = [sigma2s sigma2];
 tau2s = [tau2s tau2];
 betas = [betas beta];
  waitbar(r/M)
end
close(h)

figure(1)
subplot(3,1,1)
plot((M-500:M), betas(M-500:M))
subplot(3,1,2)
plot((M-500:M), tau2s(M-500:M))
subplot(3,1,3)
plot((M-500:M), sigma2s(M-500:M))
print -depsc 'C:\Brani\Courses\Bayes\Handouts\Working11\Figs\mcmc2_1.eps'


figure(2)
burnin = 2000;
subplot(3,1,1)
hist(betas(burnin:M), bin)
subplot(3,1,2)
tauss = tau2s(burnin:M);
tausss = tauss(tauss < 30); %plenty of outliers!!!
hist(tausss, bin)
subplot(3,1,3)
hist(sigma2s(burnin:M), bin)
print -depsc 'C:\Brani\Courses\Bayes\Handouts\Working11\Figs\mcmc2_2.eps'
median(betas(burnin:M))
median(tau2s(burnin:M))
median(sigma2s(burnin:M))