% A fragment of MATLAB code to compute the extinction time
% cumulative distribution function for Hudsonia montana
% using Tuljapurkar's approximations for the mean and variance
% of the stochastic log growth rate of a structured population 
% and a diffusion approximation for the probability of quasi-
% extinction; 
% see Box 7.6 in Morris and Doak, Quantitative Conservation Biology.

hudmats;				% read matrices from file hudmats.m, and
matrices=[A85(:) A86(:) A87(:) A88(:)];	  % put them in "matrices"
s=sqrt(size(matrices,1));				  % s=number of classes
Abar=reshape(mean(matrices'),s,s); 	 	  % calculate mean matrix Abar
C=cov(matrices');						  % calculate covariances of
										  %    matrix elements
[lambdas,lambda1,W,w,V,v]=eigenall(Abar); % compute lambda1, v, and w for 
										  %    Abar using "eigenall" from 
										  %    Box 7.1
S=v*w'/(v'*w);	% compute matrix of eigenvalue sensitivities for Abar		
S=S(:);			%    and convert it to vector form
sigma2=S'*C*S/lambda1^2; % compute sigma^2 from Equations 7.12 and 7.13
loglams=log(lambda1)-0.5*sigma2;  % compute stochastic log lambda from 
								  %    Equation 7.12
n0=[4264; 3; 30; 16; 25; 5];	% current population vector
v=v/(v'*w);		% rescale v so that both sum(w)=1 and v'*w=1
Nc=v'*n0;		% starting pop.=initial total reproductive value
Nx=500;			% quasi-extinction threshold
tmax=50;		% farthest time horizon
d=log(Nc/Nx);	% log distance to quasi-extinction threshold
cdf=extcdf(loglams,sigma2,d,tmax); 	% calculate CDF using function
									%    "extcdf" defined in Box 3.3
plot(cdf)		% plot the CDF
xlabel('Years into the future')
ylabel('Cumulative probability of quasi-extinction')