gibbs_sampler_reml_like.m

%This generates reps samples of mu, s2, s2np according to the likelihood using a gibbs sampler
%updated July 2005
%written by eli.holmes@noaa.gov

function gibbsamps = gibbs_sampler_reml_like(data, varargin)
%this function will return a set number of reps, but will actually generate a default of 10,000
%reps from which it will randomly grab a user defined number
%normally called just with 'reps to return' and reps set at 10,000
%but user can put in 'reps to return',reps to change the reps

if(nargin==3) 
   reps=varargin{2}; 
   reps_to_return = varargin{1};
elseif(nargin==2)
   reps=10000;
   reps_to_return = varargin{1};
elseif(nargin==1) %default values
   reps=10000;
   reps_to_return = 1000;
end

burn = 1000;
reps = burn+reps;

%random draws from parameter space
likemu=zeros(1,reps); likes2=zeros(1,reps);
mu=zeros(1,reps); s2s=zeros(1,reps); s2p=zeros(1,reps);
accept=0; reject=0;

%Get some starting values
[mu_ML tmp sigma2p_ML sigma2s_ML]=remltrend(data);
%don't start at ML value
muR=mu_ML+.01;
sigma2p = sigma2p_ML+.001;
sigma2s = sigma2s_ML+.001;

%Set things up
%  generate U matrix such that U'*X = 0 
%  This annhilates fixed effect (mean) for REML estimation
n = length(data);
w = log(data(2:n)) - log(data(1:n-1));
[wr wc]= size(w);
U = diag(-ones(wr,1)) + diag(ones(wr-1,1),-1);
U = U(:,1:wr-1);
% create time-series (z) with fixed mean annhilated
z = U'*w;
[rz cz]= size(z);

%First calculate the likemu using the ML estimates sigma2s and sigma2p
%using ML ests since mu likelihood should not influence s2 selection

sig_1 = 2*sigma2s_ML+sigma2p_ML;
sig_2 = -sigma2s_ML;
Vhat_obs = diag(sig_1*ones(rz+1,1)) + diag(sig_2*ones(rz,1),1) + diag(sig_2*ones(rz,1),-1);
X = ones(wr,1); 

%Start the Metropolis-Hastings algorithm   
%initial conditions
mu0 = mu_ML+.01;
s2s0 = sigma2s_ML+.001;
s2p0 = sigma2p_ML+.001;
likemu_0 = normpdf( mu_ML, mu0, sqrt(inv(X' / Vhat_obs * X)) );
sig_1 = 2*s2s0+s2p0; sig_2 = -s2s0;
Vhat = diag(sig_1*ones(rz+1,1)) + diag(sig_2*ones(rz,1),1) + diag(sig_2*ones(rz,1),-1);
SIGMA = U'*Vhat*U;
negloglik = .5*rz*log(2*pi) + .5*log(det(SIGMA)) + .5*z'*inv(SIGMA)*z;
likes2_0 = exp(-negloglik);
MLElike = likemu_0*likes2_0;

dmu=0.1; %standard deviation of normal for mu generation
dtau = 0.25; %standard deviation of lognormal for s2 generation

%get a good range for s2p
tmp=[];
therange=0.000001:.001:1;
for(s2axis=therange )
   sig_1 = 2*sigma2s_ML+s2axis; sig_2 = -sigma2s_ML;
   Vhat = diag(sig_1*ones(rz+1,1)) + diag(sig_2*ones(rz,1),1) + diag(sig_2*ones(rz,1),-1);
   SIGMA = U'*Vhat*U;
   negloglik = .5*rz*log(2*pi) + .5*log(det(SIGMA)) + .5*z'*inv(SIGMA)*z;
   tmp=[tmp exp(-1 * negloglik)];
end
where_max_tmp = find(tmp==max(tmp));
where_tmp_small = find(tmp < max(tmp)*1e-4);
maxs2p = therange(min(where_tmp_small(where_tmp_small > where_max_tmp)));
if(isempty(maxs2p)) maxs2p = 1; end

%get a good range for s2s
tmp=[];
therange=0.000001:.001:1;
for(s2axis=therange )
   sig_1 = 2*s2axis+sigma2p_ML; sig_2 = -s2axis;
   Vhat = diag(sig_1*ones(rz+1,1)) + diag(sig_2*ones(rz,1),1) + diag(sig_2*ones(rz,1),-1);
   SIGMA = U'*Vhat*U;
   negloglik = .5*rz*log(2*pi) + .5*log(det(SIGMA)) + .5*z'*inv(SIGMA)*z;
   tmp=[tmp exp(-1 * negloglik)];
end
where_max_tmp = find(tmp==max(tmp));
where_tmp_small = find(tmp < max(tmp)*1e-4);
maxs2s = therange(min(where_tmp_small(where_tmp_small > where_max_tmp)));
if(isempty(maxs2s)) maxs2s = 1; end

iter = 0;
while(iter unifrnd(0,1))
      iter=iter+1;
      mu(iter)=mu1;
      s2s(iter)=s2s1;
      s2p(iter)=s2p1;
      likemu(iter)=likemu_1;
      likes2(iter)=likes2_1;
      accept=accept+1;
      s2s0=s2s1;
      s2p0=s2p1;
      mu0=mu1;
      likemu_0 = likemu_1;
      likes2_0 = likes2_1;
   else 
      reject=reject+1;
      iter=iter+1;
      mu(iter)=mu0;
      s2p(iter)=s2p0;
      s2s(iter)=s2s0;
      likemu(iter)=likemu_0;
      likes2(iter)=likes2_0;
   end
end

%now get a random sample from those
randsamp = burn+unidrnd(reps-burn,1,reps_to_return);

gibbsamps.mu=mu(randsamp); gibbsamps.s2s=s2s(randsamp); gibbsamps.s2p=s2p(randsamp);
likemu=likemu(randsamp); likes2=likes2(randsamp);

gibbsamps.like = likemu.*likes2;
gibbsamps.chisqstat = -2*log(gibbsamps.like/MLElike);