Lab3.m

%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