%Age from stage
% Workshop at Age Determination Workshop
% October 2007, Rostock
% Carol Horvitz and Shripad Tuljapurkar
%
%  this script calls function "Nstochseq_mod2_path" 
%       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 
%initial coding in july05 to implement analytical stuff from Tulja's notes 14june05
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  MODULE 1
%
% in this module, we should 
%  
%  1) load the demographic and habitat trans matrices
%  2) specify the structures
%  3) correct for perfect one-period survival
%  4) choose the subsets of demographic matrices to be analyzed
%  5) choose the habitat dynamics matrix
%  6) get the dimensions of the set you have chosen
%demographic matrices must be in an array
%habitat transition must be in a matrx
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%name the species, set up
names=['Ardisia       '];
codes=['AE'];
   

%Markov chain habitat transition matrix follows 
%import from numeric text file
load ardesc_cmat.txt
cmat=ardesc_cmat
  
sp=1;                             %here there is only 1 sp and 1 hab dynamics
name=names(sp,:);                 %finds the name of the species
code=codes(sp,:);                 %find the two letter code of the species
matfilename=(name);               %uses species name to name output file
nstates=size(cmat,1);             %S in the notes

maxage=300; % pick a max age to which you will estimate survivorship


%import the qs from numeric text files and make into array
for i=1:nstates
eval(['qs(:,:,i)=load (''ardesc_q' num2str(i) '.txt'');' ]);
end
                   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
arraydim=size(qs);
nstates=arraydim(3);     %number of states of the environment
nstages=arraydim(1);     %number of stages

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
stagelegend=['seeds    '
             'seedling '
             'juvenile '
             'pre-rep  '
             'small rep'
             'med rep  '
             'large rep'
             'xlrge rep']; %Ardisia escallonioides


%%% choose habitat transition matrix %%%%%%%%%%%%%%%%

cmat=cmat; %cmat was the input file name in this case

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% get ready to do many sample paths to get the stuff of interest
%%%set up some stuff to make arrays for receiving output
nmatrices=nstates;
runs=20%0;      % number of sample paths
time_max=5000; % for class I reduced the interations 
               % in practice we do 40,000 this parameter used in the function input
%transients=150; %may not be relevant?

%set up arrays to receive output
Ns=zeros(nstages,nstages,runs); %each page is a result from a given sample path
stages_xs=zeros(nstages,time_max,runs); %each page is a result from a given sample path
lifeexps=zeros(runs,nstages);   %each row is a  result from a given sample path
survs=zeros(runs,time_max);     %each row is a  result from a given sample path
agemaxs=zeros(runs,1);           %each row is a result from a given sample path
paths=zeros(runs,time_max);
%agedeaths=zeros(runs,1);         %each row is a result from a given sample path

%output for many runs will be created and stored;
%  many of these are arrays
for r=1:runs
[N,lifeexp,surv,stages_x,agemax,path] = Nstochseq_mod2_path(cmat,qs,time_max); % call the function 
Ns(:,:,r)=N;
lifeexps(r,:)=lifeexp;
survs(r,:)=surv;
agemaxs(r)=agemax;
paths(r,:)=path;
stages_xs(:,:,r)=stages_x;
%agedeaths(r)=agedeath;
end 


%%%%% I have been running just to here and saving the
save oct1107run1_ares %saved the version with 20 runs 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%to add cumulative mortality to the story start here
load oct1107run1_ares


thing1=1:time_max-1; %sets the index for the number of mortality measures
cummort=zeros(runs,time_max-1);    
for i=1:runs;
    growth=survs(i,2:time_max)./survs(i,1:time_max-1); %l(x+1)/l(x)
    loggrowth=log(growth); %logs of the individual time step lambdas
    thing2=-cumsum(loggrowth)./thing1;%the converg
    cummort(i,:)=thing2;
end    


%%%%%Here note ...for figures that have cummort AND stage structure...  

%  use code from makefig_state_stage_mu.m  and the very reduced mat file
%  calculates stoch mort for each path without saving it  over all runs...

firstage=1 %for transients

lastage=time_max;
thing1=1:lastage-1; %sets the index for the number of ages, length of path

for i=[8]; % the 8th run, which is not the same for the two data sets
figure(i)
%colormap gray
%subplot(3,1,1) %for transients
subplot(4,1,1) %for asymptotics
plot(paths(i,firstage:lastage-1),'+k','markersize',8);
%xlabel('Age, yrs','FontSize',20);
ylabel('State','FontSize',20);
axis([0 lastage-firstage 0 nstates])

%subplot(3,1,2) %for transients
subplot(4,1,2) %for asymptotics
stagestuff=stages_xs(:,:,i);
%bar(stagestuff(:,firstage:lastage-1)','stacked')
area(stagestuff(:,firstage:lastage-1)')
legend(stagelegend)
ylabel('Stages','FontSize',20)
axis([0 lastage-firstage 0 1])

%subplot(3,1,3) %for transients
subplot(4,1,3) %for asymptotics
mu=log(survs(i,1:lastage-1))-log(survs(i,2:lastage));
plot(mu(firstage:lastage-1),'k','linewidth',2);
ylabel('Mortality','FontSize',20)
%xlabel('Age, yrs','FontSize',20); %for transients


%end; % for transients

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% go further for asymptotics, cumulative
    growth=survs(i,2:time_max)./survs(i,1:lastage-1); %l(x+1)/l(x)
    loggrowth=log(growth); %logs of the individual time step lambdas
    thing2=-cumsum(loggrowth)./thing1;%the convergence of the long run mean of log(growth)
    
    subplot(4,1,4) %for cumulative
    plot(thing2(firstage:lastage-1),'k','linewidth',2);
    ylabel('\it{\mu_{omega}}', 'FontSize',20)
    text(600,0.0155,['Stochastic mortality_{' num2str(lastage-1) '} = ' num2str(thing2(lastage-1))]);
    xlabel('Time step','FontSize',20); %for asymptotics
    axis([0 lastage-firstage 0.014 0.020])


end; % for asymptotics




%%%%%otherwise...go on


%%%%%to go sort by ini, we need a reduced file focused on mortaltity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Save a different file with only cummort and paths, manually
%      add other needed variables
%save  oct10pmrun_ares_simmort cummort paths
%runs=; nstates=; there may be others....forgot to record it
%load oct10pmrun_ares_simmort

ini=paths(:,1);   %get the initial state for each path
for i=1:nstates;
eval(['ini_' num2str(i) '=find(ini==' num2str(i) ');']) %find out which paths start with initial state 1, get their indices
end;

for j=1:nstates;
    voi=(['ini_' num2str(j)]);
    pathnos=eval(voi); %list of the path numbers (indices) for each initial state 
    noi(j)=length(pathnos); %sample size for each initial state
    
for i=1:noi(j); %for each of the corresponding paths with the same initial state
pathno = pathnos(i); % look up a path by its index 
%eval(['Ns_' num2str(j) '(:,:,i)=Ns(:,:,pathno);'])  %the number of N's will be the number of paths that start with 
%eval(['stages_xs' num2str(j) '(:,:,i)=stages_xs(:,:,pathno);'])
%eval(['survs_' num2str(j) '(i,:)=survs(pathno,:);'])
eval(['cummort_' num2str(j) '(i,:)=cummort(pathno,:);'])
end;                      

end;                               % 

%save here again    
%save oct10pmrun_ares_simmort
%load oct10pmrun_ares_simmort

for i=1:nstates
eval(['mu_s=cummort_' num2str(i) '(:,end);']);
mu_smean(i)=mean(mu_s);
mu_sstderr(i)=std(mu_s)/(sqrt(noi(i))); % noi(i)=no. of observations for each initial state
end
 
cummort_mean=zeros(nstates,length(cummort));
cummort_stderr=zeros(nstates,length(cummort));
for i=1:nstates
    eval(['cummort_mean(i,:)=mean(cummort_' num2str(i) ');']);
    eval(['cummort_stderr(i,:)=std(cummort_' num2str(i) ')/(sqrt(noi(i)));']);
end

for i=1:nstates
    if i<10
stuff1(i,:)=['initial state  ' num2str(i)]; %pads the spaces so we can then have two digits
    else 
stuff1(i,:)=['initial state ' num2str(i)];
    end
    end;

%save here again as oct10pmrun_ares_simmort
%save oct10pmrun_ares_simmort  


%%%oct 12%%%%% you can start here
load oct10pmrun_ares_simmort

%2006 october figures of stochastic mortality, 
%    figures saved oct12 as 'stochmort*ardesc.fig'    
figure(1)% the rate and error bars at the last step for all 7 initial states
errorbar(mu_smean,mu_sstderr,'o','linewidth',2);hold on;
plot(mean(cummort(:,end)*ones(1,7)),'--','linewidth',2);
ylabel('Long run stochastic mortality, \it{\mu_{s}} \pm SE','FontSize',16)
xlabel('Initial state','FontSize',16)
axis([0.5 7.5 0.0146 0.0152])
hold off;

figure(2) %averaged by initial states, 7 lines very close
semilogy(cummort_mean','linewidth',2)
ylabel('Mean cumulative mortality','FontSize',16)
xlabel('Age, yrs','FontSize',16)
legend(stuff1)

figure(3) %all 200 sample paths are shown
subplot(2,1,2)
first250=cummort(:,1:250);
%semilogy(cummort')
semilogy(first250')
ylabel('Mean cumulative mortality','FontSize',14)
xlabel('Age, yrs','FontSize',14)
%text(.30*time_max,1.3*cummort(1,end), ...
%    ['\it{\mu_s} \pm SE = ' num2str(mean(cummort(:,end))) '\pm' ...
%    num2str( std(cummort(:,end))/(sqrt(runs) ) ) ], 'FontSize',14);  
%%cumulative stochastic mortality at 40,000 time steps, +/- the standard error
text(.30*250,2.8*first250(1,end), ...
    ['\it{\mu_{s_{250}}} \pm SE = ' num2str(mean(first250(:,end))) '\pm' ...
    num2str( std(first250(:,end))/(sqrt(runs) ) ) ], 'FontSize',14);
subplot(2,1,1)
first50=cummort(:,1:50);
%semilogy(cummort')
semilogy(first50')
%ylabel('Mean cumulative mortality','FontSize',14)
%xlabel('Age, yrs','FontSize',14)
%text(.30*time_max,1.3*cummort(1,end), ...
%    ['\it{\mu_s} \pm SE = ' num2str(mean(cummort(:,end))) '\pm' ...
%    num2str( std(cummort(:,end))/(sqrt(runs) ) ) ], 'FontSize',14);  
%%cumulative stochastic mortality at 40,000 time steps, +/- the standard error
text(.30*50,2.8*first50(1,end), ...
    ['\it{\mu_{s_{50}}} \pm SE = ' num2str(mean(first50(:,end))) '\pm' ...
    num2str( std(first50(:,end))/(sqrt(runs) ) ) ], 'FontSize',14);



figure(4) %histogram 
hist(cummort(:,end))
ylabel('Frequency','FontSize',16)
xlabel('Long run stochastic mortality, \it{\mu_{s}}','FontSize',16)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%STOP STOP STOP STOP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%you would need to reload the file with the N's TO USE what follows,
% BUT THE CODE IS NOT YET IN GENERIC FORM

%clear all
%close all
%load oct10pmrun_ares_sim
%%%%%%%%%%%%%%%%%%%%%%%%%older stuff still needs to be generalized 
lifeexp_Ns_1=squeeze(sum(Ns_1));
lifeexp_Ns_2=squeeze(sum(Ns_2));
lifeexp_Ns_3=squeeze(sum(Ns_3));
lifeexp_Ns_4=squeeze(sum(Ns_4));
lifeexp_Ns_5=squeeze(sum(Ns_5));
lifeexp_Ns_6=squeeze(sum(Ns_6));
lifeexp_Ns_7=squeeze(sum(Ns_7));
lifeexp_Ns_8=squeeze(sum(Ns_8));
lifeexp_Ns_9=squeeze(sum(Ns_9));
lifeexp_Ns_10=squeeze(sum(Ns_10));

lifeexp_byini= [sum(mean(Ns_1,3))  %average the array over the third dimension and sum the result
               sum(mean(Ns_2,3))
               sum(mean(Ns_3,3))
               sum(mean(Ns_4,3))
               sum(mean(Ns_5,3))
               sum(mean(Ns_6,3))  
               sum(mean(Ns_7,3))
               sum(mean(Ns_8,3))
               sum(mean(Ns_9,3))
               sum(mean(Ns_10,3))];

surv_byini=[mean(survs_1,1)
            mean(survs_2,1)
            mean(survs_3,1)
            mean(survs_4,1)
            mean(survs_5,1)
            mean(survs_6,1)
            mean(survs_7,1)
            mean(survs_8,1)
            mean(survs_9,1)
            mean(survs_10,1)];
        
       

  %sept 18, 2006...mean of the difference between the logs?
  %  is this what we want, start with the mus and take the mean.OR..see below
  mu1_byini=[mean(mus_1,1)
            mean(mus_2,1)
            mean(mus_3,1)
            mean(mus_4,1)
            mean(mus_5,1)
            mean(mus_6,1)
            mean(mus_7,1)
            mean(mus_8,1)
            mean(mus_9,1)
            mean(mus_10,1)];
        
 % or do we want the median? 
 mumed_byini= [median(mus_1,1)
            median(mus_2,1)
            median(mus_3,1)
            median(mus_4,1)
            median(mus_5,1)
            median(mus_6,1)
            median(mus_7,1)
            median(mus_8,1)
            median(mus_9,1)
            median(mus_10,1)];
        
            
 % ?OR is this what we want: start with the mean survivorship curves...
 % mean of the logs .NE. to log of the means
 
 lx_byini=surv_byini;
 lxplus1_byini=surv_byini(:,2:size(surv_byini,2));

 mu2_byini=log(lx_byini(:,1:size(surv_byini,2)-1))-log(lxplus1_byini);

 
 
 mumed   = [median(mus_1,1)
            median(mus_2,1)
            median(mus_3,1)
            median(mus_4,1)
            median(mus_5,1)
            median(mus_6,1)
            median(mus_7,1)
            median(mus_8,1)
            median(mus_9,1)
            median(mus_10,1)];
        
 
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
        
%stuff from 05        
%get some extremes and percentiles
survs_mean=mean(survs);
survs_max=max(survs);


survs_min=min(survs);
survs_prctl=prctile(survs,[5 25 50 75 95]);
survs_var=var(survs);
survs_cv=(survs_var.^0.5)./survs_mean;        

%next calculate means...
lifeexp_mean=mean(lifeexps);
lifeexp_median=median(lifeexps);
agemax_max=max(agemaxs);
survs_agemax=survs(:,1:agemax_max); %choose only the ages which have at least once been maxage
survs_agemax_mean = mean(survs_agemax);
survs_agemax_median=median(survs_agemax);
agemax_mean = mean(agemax);
agemax_median = median(agemax);


mus_prctl=prctile(mus,[5 25 50 75 95]);
mus_mean=mean(mus);
mus_max=max(mus);
mus_min=min(mus);
mus_var=var(mus);
mus_cv=(mus_var.^0.5)./mus_mean;        




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%STOP STOP STOP STOP HERE IN 2006%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        

%%%%%%%%%%%%%%%%%%%%%%%%%% figures
figure(1) %look at the paths grouped by initial state
plot(lifeexp_Ns_1,'b');
hold on;
plot(lifeexp_Ns_2,'c');
plot(lifeexp_Ns_3,'y');
plot(lifeexp_Ns_4,'r');
plot(lifeexp_median,'k:','linewidth',4)
ylabel('E[\eta_{j}]');  %this IS the _expected_ life span for a given sample path of environments
                        % It is not "eta" itself
                        % eta is the number of years an individual is alive before being
                        % absorbed...
                        % this is the expectation of that number, averaged over individuals 
title(['Markov case ', cases(c,:), ': conditional life expectancy, E[\eta_{j}], numerical simulation'])
%title(['Iid case, equal probabilities: conditional life expectancy, E[\eta_{j}], numerical simulation']);
xlabel('initial stage, j');

    
%%% figure that shows the distribution of each E[\eta]'s...
%  the N for each is the  number of sample paths %%%%%%%%%%%%%%%%
figure(2);
colormap white;
for i=1:nstages;
subplot(nstages/2,2,i);
hist(lifeexps(:,i)); %these are arranged : rows = runs, columns = stages
ylabel('Frequency');
%xlabel(['e_' num2str(i)])
title(['all states, stage=' num2str(i)]);
end
figure(2);
for i=7:8;
 subplot(nstages/2,2,i);
 xlabel('E[\eta_{j}]')
end


%%%%  the same data subdivided by initial state, it takes four figures

figure(3);
colormap white;
for i=1:nstages;
subplot(nstages/2,2,i);
hist(lifeexp_Ns_1(i,:)); %here the rows are stages and the columns are the runs
ylabel('Frequency');
%xlabel(['e_' num2str(i)])
title(['initial state=1, stage=' num2str(i)]);
end
figure(3);
for i=7:8
 subplot(nstages/2,2,i);
 xlabel('E[\eta_{j}]');
end

figure(4);
colormap white;
for i=1:nstages
subplot(nstages/2,2,i);
hist(lifeexp_Ns_2(i,:));
ylabel('Frequency');
%xlabel(['e_' num2str(i)])
title(['initial state=2, stage=' num2str(i)]);
end
figure(4);
for i=7:8
 subplot(nstages/2,2,i);
 xlabel('E[\eta_{j}]');
end

figure(5)
colormap white;
for i=1:nstages
subplot(nstages/2,2,i);
hist(lifeexp_Ns_3(i,:));
ylabel('Frequency');
%xlabel(['e_' num2str(i)])
title(['initial state=3, stage=' num2str(i)]);
end
figure(5)
for i=7:8
 subplot(nstages/2,2,i);
 xlabel('E[\eta_{j}]');
end

figure(6)
colormap white;
for i=1:nstages
subplot(nstages/2,2,i);
hist(lifeexp_Ns_4(i,:));
ylabel('Frequency');
%xlabel(['e_' num2str(i)])
title(['initial state=4, stage=' num2str(i)]);
end
figure(6)
for i=7:8
 subplot(nstages/2,2,i);
 xlabel('E[\eta_{j}]');
end

%%%%%   now take the mean over all each sample paths that begins with the
%%%%%   same environmental state

figure(7)
colormap gray;
bar(lifeexp_byini');
ylabel('Mean of E[\eta_{j}]');
title(['Markov case ', cases(c,:), ': conditional life expectancy,E[\eta_{j}], numerical simulation']);
%E[\eta_{j}], numerical simulation']);
%title(['Iid case, equal probabilities: conditional life expectancy, E[\eta_{j}], numerical simulation']);
legend('initial state=1','initial state=2','initial state=3','initial state=4');
xlabel('initial stage, j');

%%%%%%%%%%%%%%%%%%%%%%%%looks identical to the analytical figure



figure(8);
%subplot(2,1,1)
plot(lifeexps');
hold on;
plot(lifeexp_median,'k:','linewidth',4);
ylabel('Life expectancy, yrs','FontName','Helvetica','FontSize',16);
xlabel('Stage','FontName','Helvetica','FontSize',16);
%title(['\it {Calathea ovandensis}, Site 3, median life expectancy = ' num2str(lifeexp_median(1))])
title(['Markov case ', cases(c,:), ': median life expectancy =' num2str(lifeexp_median(1))])';
%title(['Iid case, equal probabilities: median life expectancy = ' num2str(lifeexp_median(1))]);
%axis([1 8 0 20])
%legend(['Markov case ', cases(c,:)])
legend(['Iid case']);
%text(1.5,8,['\lambda_S =' num2str(lambda_s)]) %got the lambda separately 

figure(9);
semilogy(survs(:,1:agemax_max)');
hold on;
plot(survs_agemax_median(1:agemax_max),'k:','linewidth',4);
ylabel('Survivorship (\it l_x)','FontName','Helvetica','FontSize',16);
xlabel('Age, yrs','FontName','Helvetica','FontSize',16);
title(['\it {Calathea ovandensis}, Site 3']);
axis([0 90 10^-7 10^0]);
legend(['Markov case ', cases(c,:)]);
%legend(['Iid case, equal probabilities']);
%text(10,.1,['\lambda_S = ' num2str(lambda_s)]) %got the lambda separately 
text(10,.05,['median of maximum age of interest = ' num2str(agemax_median) ]); %last age for which survivorship > 0.00001
%text(10,.01,['median of mean age of death = ' num2str(agedeath_median)])

figure(10); 
semilogy(surv_byini(:,1:agemax_max)');
ylabel('Survivorship (\it l_x)','FontName','Helvetica','FontSize',16);
xlabel('Age, yrs','FontName','Helvetica','FontSize',16);
title(['\it {Calathea ovandensis}, Site 3']);
axis([0 90 10^-7 10^0]);
legend('initial state = 1','initial state = 2', 'initial state = 3', 'initial state = 4');
text(10, .10,['Markov case ', cases(c,:), ': average by initial state']);
%text(10, .10,['Iid case: average by initial state']);
%text(10,.1,['\lambda_S = ' num2str(lambda_s)]) %got the lambda separately 
text(10,.05,['median of maximum age of interest = ' num2str(agemax_median) ]); %last age for which survivorship > 0.00001
%text(10,.01,['median of mean age of death = ' num2str(agedeath_median)])


figure(11) %median only
%subplot(3,1,2)
colormap gray;
plot(mu,'k:','linewidth',4);
title('\it {Calathea ovandensis}, Site 3, median mortality from 1000 sample paths');
%legend(['Markov case ', cases(c,:)]);
legend(['Iid case ']);
ylabel('\mu, deaths capita^-1 yr^-1','FontName','Helvetica','FontSize',16);
xlabel('Age','FontName','Helvetica','FontSize',16);


figure(12) %median and 1000 sample paths
plot(mus');
hold on;
plot(mu,'k:','linewidth',4);
title('\it {Calathea ovandensis}, Site 3, 1000 sample paths');
legend(['Markov case ', cases(c,:)]);
%legend(['Iid case']);
ylabel('\mu, deaths capita^-1 yr^-1','FontName','Helvetica','FontSize',16);
xlabel('Age','FontName','Helvetica','FontSize',16);

figure(13)
semilogy(surv_byini_max');hold on; semilogy(surv_byini_min')
ylabel('Survivorship (\it l_x)','FontName','Helvetica','FontSize',16);
xlabel('Age, yrs','FontName','Helvetica','FontSize',16);
title(['{\it Calathea ovandensis}, Site 3, Numerical simulation, no. of runs =' num2str(runs) ]);
axis([0 90 10^-7 10^0]);
legend('initial state = 1','initial state = 2', 'initial state = 3', 'initial state = 4');
%text(10, .10,['Markov case ', cases(c,:), ': average by initial state']);
text(10, .10,['Iid case: minimum and maximum for each initial state']);
%text(10,.1,['\lambda_S = ' num2str(lambda_s)]) %got the lambda separately 
%text(10,.05,['median of maximum age of interest = ' num2str(agemax_median) ]); %last age for which survivorship > 0.00001
%text(10,.01,['median of mean age of death = ' num2str(agedeath_median)])

figure(14)
plot(surv_byini_cv')
ylabel('CV of Survivorship (\it l_x)','FontName','Helvetica','FontSize',16);
xlabel('Age, yrs','FontName','Helvetica','FontSize',16);
title(['{\it Calathea ovandensis}, Site 3, Numerical simulation, no. of runs =' num2str(runs) ]);
%axis([0 90 10^-7 10^0]);
legend('initial state = 1','initial state = 2', 'initial state = 3', 'initial state = 4');
%text(10, .10,['Markov case ', cases(c,:), ': average by initial state']);
text(10, 3,['Iid case']);

figure(15)
plot(mu_byini_max');hold on; plot(mu_byini_min')
ylabel('Mortality, {\it \mu_x}','FontName','Helvetica','FontSize',16);
xlabel('Age, yrs','FontName','Helvetica','FontSize',16);
title(['{\it Calathea ovandensis}, Site 3, Numerical simulation, no. of runs =' num2str(runs) ]);
axis([0 90 10^-7 10^0]);
legend('initial state = 1','initial state = 2', 'initial state = 3', 'initial state = 4');
%text(10, .10,['Markov case ', cases(c,:), ': average by initial state']);
text(10, .10,['Iid case: minimum and maximum for each initial state']);
%text(10,.1,['\lambda_S = ' num2str(lambda_s)]) %got the lambda separately 
%text(10,.05,['median of maximum age of interest = ' num2str(agemax_median) ]); %last age for which survivorship > 0.00001
%text(10,.01,['median of mean age of death = ' num2str(agedeath_median)])

figure(14)
plot(mu_byini_cv')
ylabel('CV of Mortality {\it \mu_x}','FontName','Helvetica','FontSize',16);
xlabel('Age, yrs','FontName','Helvetica','FontSize',16);
title(['{\it Calathea ovandensis}, Site 3, Numerical simulation, no. of runs =' num2str(runs) ]);
%axis([0 90 10^-7 10^0]);
legend('initial state = 1','initial state = 2', 'initial state = 3', 'initial state = 4');
%text(10, .10,['Markov case ', cases(c,:), ': average by initial state']);
text(10, 3,['Iid case']);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%