% PROGRAM RandDraw.m
% Does multiple simulations of discrete exponential growth  
% trajectories with a set of observed lambda values.
clear all
%**************** SIMULATION PARAMETERS **********************
lams = [1.00, 1.98, 1.02, 0.92, 0.53];   	
		% The different lambda values to use. There can be as
		% many as you want. In the simulations each will be
		% drawn with equal probability.                  
numzero = 29;	% starting population size
tmax = 100;		% length of simulations
numreps = 100;   % number of replicate trajectories to simulate
outname = 'temp24';	% name of file in which to store results 
%*************************************************************

rand('state',sum(100*clock)); %this seeds the random number
  % generator, which most be done in all Matlab programs
numlams = length(lams);	% how many lambda values are there? 
intvalues = floor(numlams*rand(tmax,numreps)) + 1; 	
			% generates a matrix of random index numbers that 
			% will be used to decide which lams value to use 
			% for each year and each simulation
lambdas = lams(intvalues);										
           % uses the indices to create a matrix of lambdas 
			% to use for each year and simulation
initpopsize = numzero*ones(1,numreps); 
			% sets initial population sizes
popsizes = [initpopsize;numzero*cumprod(lambdas)];						
				% generates the population size every year 
				% for all simulations, using the cumulative
				% products of the lambda values
% the next three lines find extinctions and set subsequent
% population sizes to zero
notextinct= popsizes > ones(tmax+1,numreps); 
notextinct = cumprod(notextinct);
popsizes = popsizes.*notextinct;
      
maxN = 0.5*max(max(popsizes));	
	% the biggest population size to graph. Max pop'n size
	% will often be too high to see most other dynamics
stochL = prod(lams)^(1/numlams)		
		% geometric mean of the lambdas or stochastic lambda 
expectedNt = numzero*(stochL^tmax) 
				% this is the median population size at tmax
wk1write(outname,popsizes)	 % writes results to Lotus spreadsheet
% the next lines make a figure of the population sizes vs. time
plot([0:tmax], popsizes)	 
xlabel('Time');
ylabel('Population size');
axis([0 tmax 0 maxN]);

% the next lines make a histogram of final population sizes
figure % make a new figure window
Y=[0 1E-10 10:10:200 inf]; % set the category boundaries  
N=histc(popsizes(tmax,:),Y); %count up the numbers of cases 
N=100*N/numreps;    % convert numbers to percentages
Nmax=max(N);        
X=(-5:10:215);      % values for x-axis of histogram
bar(X,N)            % plot a bar graph of the histogram
axis([-10 210 0 Nmax+5])
xlabel('Population size (<0=extinct, last bar >200)');
ylabel('Percent of populations');