%This generates reps samples of mu, s2, s2np according to the likelihood using a gibbs sampler
%updated Dec 25, 2004
%THIS CODE GENERATE SAMPLES ACCORDING TO THE FREQ YOU WOULD EXPECT BY RANDOM CHANCE
%TO ESTIMATE THIS, IT USES THE ASYMP. CHI-SQ DISTRIBUTION FOR -2*LOG(LAMBDA)
%THIS IS A DECENT APPROX. AS SHOWN IN test_like_ratio_distribution.m

function [mu, s2p, s2s, chisqstat]=gibbs_sampler_kalman_f(data, reps)
numparams = 3;
%set up data for kalman_loglik which needs log transformed data with missing data as -99
tmp = data;
tmp(data==-99) = 1;
y = log(tmp);
y(data==-99) = -99;
initx=log(data(1));

Nratios = data(2:end)./data(1:(end-1));
Nratios = Nratios( Nratios > 0 );
Nratios = log(Nratios);
sigma2init = 0.5 * var(Nratios);
sigma2npinit = 0.5 * var(Nratios);
V1 = .001;

burn=2000;
like=zeros(1,reps);
mu=zeros(1,reps); s2p=zeros(1,reps); s2s=zeros(1,reps);
f = zeros(1,reps); tmp3 = zeros(1,reps);
accept=0; reject=0;

%Start the Metropolis-Hastings algorithm   
%initial conditions
[mu0, s2p0, s2s0, initx, negloglik] = kalman(data);
%redo just to standardize the initx and V1 used
MLElike = exp( -1 * kalman_loglik([mu0 log(s2p0) log(s2s0), initx],y,V1,0) ); %0 is flag to not plot

%get a good range for s2p
tmp=[];
therange=0.000001:.001:1;
for(s2axis=therange )
   tmp=[tmp exp(-1 * kalman_loglik([mu0 log(s2axis) log(s2s0), initx],y,V1,0))];
end
where_max_tmp = find(tmp==max(tmp));
where_tmp_small = find(tmp < max(tmp)*1e-3);
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 )
   tmp=[tmp exp(-1 * kalman_loglik([mu0 log(s2p0) log(s2axis) initx],y,V1,0))];
end
where_max_tmp = find(tmp==max(tmp));
where_tmp_small = find(tmp < max(tmp)*1e-3);
maxs2s = therange(min(where_tmp_small(where_tmp_small > where_max_tmp)));
if(isempty(maxs2s)) maxs2s = 1; end

%don't start right at MLE
mu0 = mu0 + 0.001;
log_like_0 =  -1 * kalman_loglik([mu0 log(s2p0) log(s2s0) initx],y,V1,0);
like_0 = exp( log_like_0 );
f_0 = chi2pdf(-2*(log_like_0 - log(MLElike)),numparams);

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

iter = 0;
while(iter<(burn+reps))
   % use a uniform proposal distribution
   mu1 = normrnd(mu0,dmu);
   %mu1 = unifrnd(-.2,.2);
   %Normal proposition distribution; NOT USED
   %s2s1 = normrnd(s2s0,dtau); %m = s2s0 = exp(log(s2s0)); sigma = dtau
   %pr01 = normpdf(s2s1,s2s0,dtau);
   %pr10 = normpdf(s2s0,s2s1,dtau);
   %if(s2s1 < 1e-10) s2s1 = 1e-10; end
   %qs01 = (pr01/pr10);
   
   %s2p1 = normrnd(s2p0,dtau); %m = s2s0 = exp(log(s2s0)); sigma = dtau
   %pr01 = normpdf(s2p1,s2p0,dtau);
   %pr10 = normpdf(s2p0,s2p1,dtau);
   %if(s2p1 < 1e-10) s2p1 = 1e-10; end
   %qp01 = (pr01/pr10);
   
   %Gamma proposition distribution; NOT USED
   %gamscale = max(s2p0,.001);
   %gamscale = .01;
   %s2p1 = gamrnd(1,gamscale); %
   %prob of going from 0 to 1
   %pr01 = gampdf(s2p1,1,gamscale);
   %prob of going from 1 to 0
   %pr10 = gampdf(s2p0,1,gamscale);
   %should be ratio of going 1 to 0 over going 0 to 1
   %qp01 = pr10/pr01;
   
   %Uniform proposition distribution
   %I've found that this inefficient proposition dist works best
   s2p1 = unifrnd(0,maxs2p);
   qp01 = 1;
   s2s1 = unifrnd(0,maxs2s);
   qs01 = 1;
   
   log_like_1 =  -1 * kalman_loglik([mu1 log(s2p1) log(s2s1) initx],y,V1,0);
   f_1 = chi2pdf(-2*(log_like_1 - log(MLElike)),numparams);
   like_1 = exp( log_like_1 );
   
   logr = log(qs01) + log(qp01) + log(f_1) - log(f_0);
   
   if(exp(logr) > unifrnd(0,1))
      iter=iter+1;
      if(iter>burn) 
         mu(iter-burn)=mu1;
         s2s(iter-burn)=s2s1;
         s2p(iter-burn)=s2p1;
         like(iter-burn)=like_1;
         f(iter-burn)=f_1;
         tmp3(iter-burn) = exp(logr);
      end
      accept=accept+1;
      s2s0=s2s1;
      s2p0=s2p1;
      mu0=mu1;
      like_0 = like_1;
      f_0 = f_1;
   else 
      reject=reject+1;
      iter=iter+1;
      if(iter>burn)
         mu(iter-burn)=mu0;
         s2p(iter-burn)=s2p0;
         s2s(iter-burn)=s2s0;
         like(iter-burn)=like_0;
         f(iter-burn)=f_0;
         tmp3(iter-burn) = exp(logr);
      end
   end
end
disp(reject/(accept+reject));
chisqstat = -2*log(like/MLElike);
   keyboard