% PROGRAM extprob
%   Calculates the probability of extinction with bootstrap confidence 
%   intervals for a density-independent model using a diffusion approximation; 
%   See Dennis et al. 1991, Ecological Monographs 61: 115-143;
%   REQUIRES THE FILES extcdf.m, chi2rv.m, stdnormcdf.m, gammarv.m, and betarv.m

%*************************************************************************
% Change the following user-defined parameters:
%*************************************************************************
mu=0.02134;     % Enter the estimated value of mu
sig2=0.01305;   % Enter the estimated value of sigma^2
CI_mu=[-.01621,.05889];  % Enter confidence interval for mu
CI_sig2=[.00867,.02184]; % Enter confidence interval for sigma^2
q=38;           % Enter the number of transitions in the data set
tq=38;          % Enter the length of the census (in years)
Nc=99;          % Enter the current population size
Nx=20;          % Enter the quasi-extinction threshold
tmax=50;        % Enter latest time to calculate extinction probability
Nboot=500;      % # of bootstrap samples for calculating confidence intervals
                %       for extinction probabilities
%**************************************************************************

d=log(Nc/Nx);	        % calculate log distance to quasi-extinction threshold
SEmu=sqrt(sig2/tq);     % calculate standard error of mu
Glo=ones(1,tmax);       % initialize array to store lower bootstrap confidence 
                        %    limits for extinction probabilities
Gup=zeros(1,tmax);      % initialize array to store upper bootstrap confidence 
                        %    limits for extinction probabilities

rand('state',sum(100*clock));  % seed the random number generator
                        
% Calculate the extinction time CDF for the best parameter estimates
Gbest=extcdf(mu,sig2,d,tmax);

% Calculate bootstrap confidence limits for extinction probabilities
for i=1:Nboot,
    
    % generate a random mu within its confidence interval
    murnd=inf;  
    while murnd<CI_mu(1)|murnd>CI_mu(2)
        murnd=mu+SEmu*randn;
    end;
    
    % generate a random sigma^2 within its confidence interval
    sig2rnd=inf;  
    while sig2rnd<CI_sig2(1)|sig2rnd>CI_sig2(2)
        sig2rnd=sig2*chi2rv(q-1)/(q-1);
    end;

    % Calculate extintion probabilities given murnd and sig2rnd
    G=extcdf(murnd,sig2rnd,d,tmax);

    % Store extreme values
    for t=1:tmax
        if G(t)<Glo(t)
            Glo(t)=G(t);
        end;
        if G(t)>Gup(t)
            Gup(t)=G(t);
        end;
    end;
    
end;    %for i=1:Nboot

% Plot G (log-scale) vs t, using a solid line for the best estimate
% of G and dotted lines for the confidence limits
t=[1:tmax];
semilogy(t,Gup,'k:',t,Gbest,'k-',t,Glo,'k:')
axis([0,tmax,0,1])
xlabel('Years into the future')
ylabel('Cumulative probability of quasi-extinction')      

disp(' ');
disp('The four columns of the following table correspond to:');
disp('  time into the future;');
disp('  best estimate of the quasi-extinction probability;');
disp('  lower confidence limit for quasi-extinction probability;');
disp('  upper confidence limit for quasi-extinction probability');
format short e;
[t' Gbest' Glo' Gup']