%Written by Eli Holmes July 2005
%This lab estimates the A2 IUCN Red List criteria: 80,50 or 20% decline in 10 yrs or 3 generations
% or probability of hitting CR, EN, VU criteria in 10 yrs or 3 generations'
%This lab estimates risk metrics assuming that s2s is extraneous variability and that we are 
%interested in modeling the risks of the unobserved x.  
%This lab also focuses on estimating % declines.  That way we don't have to factor in estimation of 
%x(end).  This is straight-forward but time-consuming and thus not ideal in a teaching situation.

%NOTE THIS IS A TEACHING TOOL; NOT CODED FOR REAL RISK ASSESSMENT
%Technical comments regarding the algorithms used are at the bottom of the script

function Lab3()

%Call getdata function to pick a data set
%Set up variables from the data
[data, dataname, datatype] = getdata(-99,20);

%plot diagnostics
diagnosticfig(data,dataname);

%Get the gibbs reps; use reml
load saved_gibbs_reml;
saved_gibbs = saved_gibbs_reml{datatype};
mu = saved_gibbs.mu;
s2p = saved_gibbs.s2p;
s2s = saved_gibbs.s2s;
chisqstat = saved_gibbs.chisqstat;

ansdone = 'N';
fignum = 2;
while(~(strcmp(ansdone,'Y') | strcmp(ansdone,'y')))
   %Get user inputs about what analyses to do
   nelevel = input('Choose 80%, 50% or 20% decline (enter 80, 50 or 20): ');
   nelevel = (100-nelevel)*.01; %change to proportion
   years = input('Choose extinction time horizon: ');
   
   %get the x estimate (x part of cda); remember this is logged
   %need the estimate of the x part of the CDA
   kalman = get_kalman_ests(data);
   Nend = exp(kalman.x(end)); %this is the ML estimate of x(end) from kalman filter
   
   ne=nelevel*Nend;
   nq=Nend;
   
   %Calculate extprob and decprob for each gibbs sample
   %the matlab functions work on vectors so no need to go through each sample one-by-one
   
   %DECPROB extinction prob at end of years      
   %risk.decprob.da is the probability that DA is below below ne at year proj
   risk.decprob.da=(1-normcdf( (log(nq/ne) + mu*years).*((s2p*years).^-0.5) ));
   
   %EXTPROB extinction prob during years
   %risk.extprob.da is based on the diffusion approximation for the DA (i.e. leaving out s2s).
   %it makes sense if the s2s is viewed as extraneous and not reflecting variability in the
   %total population
   tmp=2*log(nq/ne)*abs(mu)./s2p;
   tmp=tmp/700; %this is a trick to get rid of NAs when tmp is big
   tmp(find(tmp<1))=1;
   G = normcdf( (-log(nq/ne) + abs(mu)*years)./sqrt(s2p*years)) +...
      (exp((1./tmp).*(2*log(nq/ne)*abs(mu)./s2p)).*...
      normcdf( (-log(nq/ne) - abs(mu)*years)./sqrt(s2p*years)).^(1./tmp)).^tmp;
   probext = (ne/nq).^(2*mu./s2p);
   probext(mu<0)=1;
   risk.extprob.da = G.*probext;
   risk.everextinct.da = probext;
   
   disp('Getting median time to extinct.  Takes time.');
   %MEDIAN TIME TO EXTINCTION
   yearsvec = 1:100;
   tmp=2*log(nq/ne)*abs(mu)./s2p;
   tmp=tmp/700; %this is a trick to get rid of NAs when tmp is big
   tmp(find(tmp<1))=1;
   tmp=ones(length(yearsvec),1)*tmp;
   G = normcdf( (-log(nq/ne) + yearsvec'*abs(mu))./sqrt(yearsvec'*s2p)) +...
      (exp((1./tmp).*(ones(length(yearsvec),1)*(2*log(nq/ne)*abs(mu)./s2p))).*...
      normcdf( (-log(nq/ne) - yearsvec'*abs(mu))./sqrt(yearsvec'*s2p)).^(1./tmp)).^tmp;
   yearmat = yearsvec'*ones(1,length(mu));
   tmp = yearmat.*(G<.5);
   risk.mediantime.da=max(tmp);
   
   alpha = chi2inv(.95,3); %3 parameters; this is a bit harsh since there are not 3 orthogonal parameters
   [sorted_chi2, sorti] = sort(chisqstat);
   sorted_prob = risk.extprob.da(sorti);
   like_ext = [sorted_prob(1)   max(sorted_prob(find(sorted_chi2<alpha))) min(sorted_prob(find(sorted_chi2<alpha)))];
   
   sorted_prob = risk.decprob.da(sorti);
   like_dec = [sorted_prob(1)   max(sorted_prob(find(sorted_chi2<alpha))) min(sorted_prob(find(sorted_chi2<alpha)))];
   
   h=figure(fignum);
   set(h,'Position',[682 218 574 626]);
   subplot(2,2,1)
   %plot(sort(mu),(1:length(mu))/length(mu))
   hist(mu,-.15:.01:.15);
   title('Long-term trend');
   subplot(2,2,2)
   %plot(sort(risk.extprob.da),(1:length(risk.extprob.da))/length(risk.extprob.da))
   hist(risk.extprob.da,0.125:.25:.875);
   set(gca,'XLim',[0 1]);
   title(['Probability of ' num2str((1-nelevel)*100) '% decline during next ' num2str(years) ' years']);
   subplot(2,2,3)
   hist(risk.everextinct.da,.05:.1:.95);
   title(['Probability of eventual ' num2str((1-nelevel)*100) '% decline']);
   set(gca,'XLim',[0 1]);   
   subplot(2,2,4)
   hist(risk.mediantime.da);
   title(['Median time to decline given it occurs']);
   set(gca,'XLim',[0 100]);
   
   %Put title at the top of the whole graph
   h=axes('Position',[0 0 1 1],'Visible','off');
   text(.5,.98,dataname,'HorizontalAlignment','center','FontSize',16)
   
   set(gcf,'PaperPosition',[.5 1 7.5 9]);
   ansdone = input('Done? (Y or N) ','s');
   fignum = fignum + 1;
end

%Technical stuff:
%This lab uses mu,s2p,s2s samples from the joint-likelihood function using a gibbs sampler
%the risk metrics are calculated from the parameter sets; the histogram of those metrics is then
%proportional to the posterior (using uniform priors)
%Inference is based on a bayesian approach using essentially a uniform (improper) prior
% and intergrating over all parameters for which risk = theta (July 2005)
%In my actual applications, I tend to prefer inference based on likelihood ratios and the 
% asymptotic distribution of those ratios.  But that is harder to explain than the bayesian approach,
% so I used the bayesian approach in the workshop.