%CorrRates: a program to generate sets of correlated vital rates
clear all;

%*****************Simulation Parameters***********************
% parameters for 3 vital rates: 
% a beta, a lognormal, and a stretched beta.
vrmeans= [0.77,2,5]; % means
vrvars=	[0.02,0.5,2]; % variances (not standard deviations)

% minimum and maximum values for each vital rate: 
% zeros are place holders for rates that are not stretched betas 
vrmins= [0 0 0];
vrmaxs= [0 0 24];
%then a full correlation matrix.
corrmx = ...
[1 		0.3 	-0.2;
0.3		1		-0.6;
-0.2		-0.6	1];

tmax = 200; 	% number of sets of vital rates to simulate
%*************************************************************

np = length(vrmeans); % finds the number of vital rates 
results = [];
randn('state',sum(100*clock)); % seeding random numbers

% find the eigenvalues (D) and eigenvectors (W)
% of the correlation matrix
[W,D] = eig(corrmx);  

%Calculate C12, the matrix to use to make correlated 
% standard normal variates from uncorrelated ones; see
% the text in Chapter 8 for an explanation of this process   
C12 = W*(sqrt(abs(D)))*W';

for tt=1:tmax % loop to do each year's vital rates
 normvals = randn(np,1); 
 % make a set of random standard normal values. 
 corrnorms = C12*normvals; % make correlated normals		
 % get the beta vital rate with the same Fx as the first normal:
 vrate1 = ...      
   betaval(vrmeans(1),(vrvars(1))^0.5,stnormfx(corrnorms(1)));
 % convert the second normal into a corresponding log normal: 
 vrate2 = lnorms(vrmeans(2),vrvars(2),corrnorms(2));
 %then a stretched beta corresponding to the 3rd normal:
 vrate3 = stretchbetaval(vrmeans(3),(vrvars(3))^0.5,...
				vrmins(3), vrmaxs(3),stnormfx(corrnorms(3)));
 results = [results;vrate1,vrate2,vrate3]; % store the values
end; % tt

meanrates = mean(results) % the means of simulated vital rates
variances = var(results)  % variances of simulated rates
correlations = corrcoef(results) % and the correlation matrix