% MultiSiteCount. A program to do a multi-site stochastic
% PVA simulation using count data.
% NOTES:
% 1. This model is made to generate random simulations of a
% stochastic  model, generating estimates of the mean and SD of 
% realized growth rates and of extinction probabilities over a 
% specified time horizon.
% 2. The model is specified as the means, variances, and 
% correlations of population growth rates.  The model assumes 
% that log growth rates are entered
% so they can be simulated as normals, then transformed. 
clear all;
%*****************Simulation Parameters***********************
% The parameters for population growth are mu 
% and sig^2 values for three populations, then the
% prob. of leaving, and the prob. of reaching another site.
% Data for the California clapper rail; Harding et al. 2001b
vrmeans= [0.043 -0.002 0 0.2 0.5] ; 	% means
vrvars=[0.051 0.041 0.051 0 0 ];		% variances
% correlations of growth rates:
corrmx1 = ...
 [	1.000	    	0.995   0.896   0   0;
	0.995		1.000   0.938   0   0;
	0.896		0.938   1.000   0   0;
    	0           	0       0       1   0
    	0           	0       0       0   1]; 
% you must also have a separate m-file that takes a vector of
% vital rates, vrs, and creates the matrix elements with them:
makemx='railmxdef';
n0= [70 26 33]';  	% initial population sizes of each pop'n
Nmax = [286 60 58];	% maximum numbers in each population 
Ne = [20 20 20]; 	% quasi-extinction thresholds for each pop'n 
                	% (applied separately to each population) 
tmax = 100; 		% number of years to simulate
np = 5;  		    % number of vital rates
popnum = 3; 		% number of local populations 
runs = 1000;      	% how many trajectories to do
%*************************************************************

rand('state',sum(150*clock)); % M12 is used for correlations
[W,D] = eig(corrmx1);  M12 = W*(sqrt(abs(D)))*W';
vrs = [ [vrmeans(1:3) + 0.5*vrvars(1:3)], vrmeans(4:5)]; % the addition of 0.5*variances 
% is necessary to arrive at the proper means for lognormal variables. see M&D equ.8.8
eval(makemx); 		
[uu,lam1] = eig(mx); 	lam1s = max(lam1);		
[lam0,iilam] = max(lam1s); 		% dominant eigenvalue
uvec = uu(:,iilam)/sum(uu(:,iilam));  % dominant right vector

results = [];
PrExt = zeros(tmax,1); 	% extinction time tracker
logLam = zeros(runs,1); % tracker of stoch-log-lambda values
stochLam = zeros(runs,1); % tracker of stochastic lambda values
for xx = 1:runs;
 if round(xx/5) == xx/5 disp(xx); end;
 nt = n0; % start at initial population size vector
 extinct = zeros(1,popnum); % vector of extinction recorders
 for tt = 1:tmax
m = randn(1,np);
   	rawelems = (M12*(m'))';	% correlated str. normals
	vrs = vrmeans + sqrt(vrvars).*rawelems;	    
  	eval(makemx); % make a matrix with the new vrs values 
  	nt = mx*nt;	 % multiply by the population vector  
    	nt = min([nt';Nmax])'; % applying population cap
  	if sum(extinct) < popnum	% these two ifs check for 
	 for nn = 1:popnum   	% extinction and count it
		  if nt(nn) <= Ne(nn) extinct(nn) = 1;  end; 
        end;     
        if sum(extinct) == popnum PrExt(tt)=PrExt(tt) +1; end; 
end; % if 
 end % tt
 logLam(xx) = (1/tmax)*log(sum(nt)/sum(n0));
 stochLam(xx) = (sum(nt)/sum(n0))^(1/tmax);
end; % xx
CDFExt = cumsum(PrExt./runs);
disp('This is the deterministic lambda value');	disp(lam0);
disp('And this is the mean stochastic lambda'); disp(mean(stochLam));
disp('Below is mean and standard deviation of log lambda');
disp(mean(logLam));		disp(std(logLam));
disp('Next is plotted a histogram of logLams'); Hist(logLam);
disp('And now,the extintion time CDF'); figure; plot(CDFExt);