% mat5stochseq_08 % May 2008 % uses mat-file "jpdata" % (does not need mat-file "rnd5001m" since the random numbers generated % internally) % (does not need function file "senstuff" since sensitivities etc are calculated % internally) % %Demography lessons % % Stochastic sequence for a Markovian environment and sample paths. % Stochastic population growth rate and what influences it. % Time averages, geometric means, cumulative population growth, lognormal distributions % % MATLAB lessons % arrays: storing data in 3 or more dimensions % using generic code % accumulating sums (or products) within loops % conditioning a calculation on the value of a variable %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This code is based on "stochseq_codeforcolleagues" % which was written at Stanford University 2004-2007 % by Carol C. Horvitz and Shripad Tuljapurkar %authorship rights % the stochastic programs and the stochastic simulation reported here is % described in Tuljapurkar, Horvitz and Pascarella, 2003 AmNat % and in Horvitz, Tuljapurkar and Pascarella, 2005 Ecology %to run this program you need %a mat-file with the demograhy matrices %a mat-file with the habitat transition rules %c(i,j) = Prob{state at time t+1 is i GIVEN that state at time t is j} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %the program is organized in modules %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %(0) before all modules: load the data %(1) module 1: %you choose %tr the number of time steps that are considered transient dynamics %ts the number of time steps that you will actually use %program generates %t2 the number of time steps in the sequence generated by the %simulation %nstage the number of life history stages %nstate the number of different states of the habitat in the system %sets up matrices and vectors to receive %states the sequence of environmental states generated by the %probability rules %growth the sequence of one time-step growth rates generated by the %simulation %uvecs the sequence population structures generated by the %simulation %vvecs the sequence of reproductive value vectors generated by the %simulation %cstar the equilibrium distribution of states %avmat the average matrix, Abar %varmat the variance of each matrix element %sdmat the standard deviation of each matrix element %cvmat the coefficient of variability of each matrix element %cov_ij the covariance of matrix elements %(2) module 2: uses the habitat transition rules to generate %states sequence of environmental states %cstar equilibrium distribution of habitats %avmat the mean matrix, Abar %varmat the variance of each matrix element %sdmat the standard deviation of each matrix element %cvmat the coefficient of variability of each matrix element %cov_ij the covariance of matrix elements %(3) module 3: uses the sequence of environmental states to generate %growth sequence of one-time growth rates %uvecs sequence of population structures %vvecs sequence of reproductive value vectors %uses the sequence of growth rates to calculate %a stochastic growth rate in "r" currency %lam_S stochastic growth rate in "\lambda" currency %(4) module 4: uses the sequences of vectors and growth rates to generate %sensE E^S sensitivity of lam_S to proportional perturbations %of ij^th rate in all states %sensA E^{S,\mu} sensitivity of lam_S to proportional perturbations %of the mean of the ij^th rate %sensV E^{S,\sigma} sensitivity of lam_S to proportional perturbations %of variability of ij^th rate %eijSbs E^S_\beta sensitivity of lam_S to proportional perturbations %of ij^th rate in a given, particular state %ebS summed habitat elasticity %sensitivity of lam_S to proportional perturbations %of all rates in a given, particular state %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% You need to input % demo matrices, one for each environment % one habitat transition matrix, which describes transitions environments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all close all load jpdata %loads the mat-file that contains John Pascarella's dissertation data % reference: % Pascarella, J.B. and C.C. Horvitz. 1998. Hurricane disturbance % and the population dynamics of a tropical understory shrub: % megamatrix elasticity analysis: Ecology 79: 547-563 % some of these parameters were altered in subsequent papers % to make some reducible % matrices irreducible, etc, but the version you have here is % the one published in this reference %for Pascarella's data: make an "array" of the 7 demography matrices. %%%%%%%%%%%%%%%%%%%%%%%%%CREATE AN ARRAY%%%%%%%%%%%%%%%%%%%%%%%%%% % Create an array in which each "page" is a demo matrix (there are 7 pages) % and each demo matrix is nstage X nstage. % (Elements in a 3-WAY array are referenced by 3 subscripts (row, column, % page); for example, the (2,1) element in habitat 7 is "ms(2,1,7)") for i=1:7 eval([ 'ms(:,:,i)=m' num2str(i) ]) %uses "evaluate" and "num2str" end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %MODULE 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% tr=200; %set the transients at 200 (shortened from 2,000 classroom use) ts=5000; %set the number you actually want at 5,000 (shortened from 100,000 for classroom use t2=ts + 2*tr; %simulate 5,000 plus 200 transients at beginning and end, respectively nstage=length(ms(:,:,1)); nstate=length(cmat); states=zeros(1,t2+1); %to receive the sequence of states growth=zeros(1,ts); %to receive the sequence of one-time growth rates uvecs=zeros(nstage,ts); %to receive the sequence of population structures vvecs=zeros(nstage,ts); %to receive the sequence of reproductive value vectors avmat=zeros(nstage); %to receive the average matrix varmat=zeros(nstage); %the variance of each matrix element sdmat=zeros(nstage); %the standard deviation of each matrix element cvmat=zeros(nstage); %the coefficient of variability of each matrix element cov_ij=zeros(nstage^2); %the covariance matrix of matrix elements %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %MODULE 2 %guts of mchain2.m, modified such that %p(i,j) defined as Prob (j->i ); % and calculates Abar, %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %outputs: %states sequence of states %cstar equilibrium distribution %avmat the mean matrix, Abar %varmat the variance of each matrix element %sdmat the standard deviation of each matrix element %cvmat the coefficient of variability of each matrix element %cov_ij the covariance matrix of matrix elements %inputs: %cmat habitat transition probability matrix, %p(i,j) defined as Prob (j->i ); %t2 length of run % determine equilibrium dist of hab states [vec,c] = eig(cmat); [cmax,loc] =max(abs(diag(c))); cstar = vec(:,loc); cstar = cstar/sum(cstar); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % AT THIS POINT you have the long run equilibrium frequency of environments figure(1) bar(cstar) xlabel('Environmental state') ylabel('Frequency') title('Long-term equilibrium frequency of environments') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %calculate Abar %calculate covariances %adapted from morris and doak, p.235, %but uses n rather than n-1 weighting % cov(x,y) = E(x*y) - E(x)*E(y), % you need the means of each parameter % & the mean of the products of parameters Abar=zeros(nstage^2,1); Exy=zeros(nstage^4,1); for i=1:nstate extractmat=ms(:,:,i); A=extractmat(:); %chooses the matrix for state i, represent it as a column vector that is nstages^2 long Exy=Exy+cstar(i)*kron(A,A); %kron(A,A) gives a vector that is nstages^4 long, since A is nstages^2 long % the pairwise products of parameters, % weight by hab freq, accumulate to get mean %try a vector that is [1 2 3 4]' to see how %this works, in this simple example the %first 4 are 1 times [1 2 3 4]' %next 4 are 2 times [1 2 3 4]', %next 4 are 3 times [1 2 3 4]', %next 4 are 4 times [1 2 3 4]' Abar=Abar+cstar(i)*A; %weight environment by its frequency, %accumulate to get mean matrix end %C=(Exy-kron(Abar,Abar))*nstate/(nstate-1); weight this way if you have a sample of states C=(Exy-kron(Abar,Abar)); % weight this way if you have the full universe of states cov_ij=reshape(C,nstage^2,nstage^2); %puts it into the form of the variance-covariance matrix avmat=reshape(Abar,nstage,nstage);%puts it into the form of the population projection matrix varmat=diag(cov_ij); varmat(varmat<0)=0; %sets spurious negatives to zero varmat=reshape(varmat,nstage,nstage); sdmat= sqrt(varmat); cvmat = sqrt(varmat); %intermediate step cvmat(avmat>0) = cvmat(avmat>0)./avmat(avmat>0);%the coefficients of variability for the nonzero Abars %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % AT THIS POINT you have summary stats about the mean, variability of matrix entries figure(2) subplot(2,1,1) bar3(avmat) title('Mean matrix') subplot(2,1,2) bar3(cvmat) title('CV of matrix entries') %GET RANDOM NUMBERS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% rx = rand(1,t2+1); %%%%%%submodule combines "cumvec", "randcomp" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% probvec = cstar; %choose initial state at random from equilibrium vector for i = 1:t2+1 psum = cumsum(probvec); %find the cumulative probability for a given vector j = min(find(psum > rx(i))); %get random number at time i from the rx vector % and decide which state it corresponds to states(i) = j; %put the state into the sequence of states probvec = cmat(:,j); %choose new column from the matrix as the next probability vector end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % AT THIS POINT: you have the FULL sequence of environmental states... you % have a sample path, you can make a figure of 100 time steps after t=3000 % figure(3) subplot(3,1,1) hist(states(1:20)) title('Sample path: first 20 states') axis([0 8 0 20]) subplot(3,1,2) hist(states(1:1000)) title('Sample path: first 1000 states') axis([0 8 0 1000]) ylabel('Frequency') subplot(3,1,3) hist(states(1:5000)) title('Sample path: all 5000 states') axis([0 8 0 5000]) xlabel('Environmental state') figure(4); plot(3001:3100,states(3001:3100),'-o'); xlabel('Time'); ylabel('Environmental state'); title('Section of sample path: sequence of environments') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %MODULE 3 %generate sequences of growth rates, population and reproductive value vectors, % calculate the stochastic growth rate %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %STOCHASTIC GROWTH RATE: CONCEPTUAL AND CALCULATION BY SIMULATION % The spirit of projection... as we move through time % project the population, but each time we change the matrix, % according to the sequence of environments. % Each time step we observe the size of the population and its shape. % Calculate N(t+1)/N(t) for each step...producing a sequence of % one-time-step growth rates...from these we calculcate a % stochastic growth rate for the population in the dynamic environment.. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %NOTE MEANINGS OF TWO RELATED PARAMETERS: % The phrase "stochastic growth rate" % sometimes refers to the parameter % 'a' (which is the same currency as "r") % but other times it refers to the parameter % 'lambda_S' (which is the same currency as "lambda) % These two parameters are related: % lambda_S = e^a (where e is the number e) % log(lambda_S)=a %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Conceptually, we want the time average of % cumulative population growth, N(t)/N(0) for very big t %Computationally, we produce a really long time series, getting % a time step specific N(t+1)/N(t) at each step. % Since population growth is multiplicative, % We want the geometric mean of these, which is obtained by % calculating log (N(t+1)/N(t)) at each step % taking the arithmetic mean of that % then raise e to that number... % This gives the geometric mean of all the entries in the series %TRY a simple series to see how logs and geometric means work %%% SIDE LESSON ABOUT GEOMETRIC MEANS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% YOU CAN COMMENT IT OUT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% series=[1 2 3] logs=log(series) %by computational method, step 1 arithmeanoflogs=mean(logs) %by computational method, step 2 etomean=exp(arithmeanoflogs) %by computational method, step 3 geomean=(prod(series))^(1/3) %by definition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% RETURN TO PROGRAM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% vec1=ones(1,nstage); vec1=vec1/sum(vec1); vec1=vec1'; %the 3 lines create a vector to use as an initial population vector vec2 = vec1; trun=ts+tr; %time to run is desired length of sample path plus transients %set up a new cumulative growth for all trun mgrowth = 0; %mgrowth2 = 0; for i = 1:trun i2=states(i); mat1 = ms(:,:,i2); vec1 = mat1* vec1; growth1 = sum(vec1); vec1= vec1/growth1; mgrowth = mgrowth + log(growth1); %mgrowth2 = mgrowth2 + log(growth1)*log(growth1); if( i > tr) %only store the ones after the transients i1 = i - tr; uvecs(:,i1) = vec1; growth(i1) = growth1; end %in reverse order for the rv's j = (t2 - i+1); %to make v(3000) use X(3000) but this comes from environment 3001 i2 = (states(j+1)); mat1 = ms(:,:,i2); %the next line Carol rewrote to make 'left' vectors vec2 = (vec2'*mat1)'; vec2= vec2/(sum(vec2)); if (i > tr) %only store the ones after the transients vvecs(:, j-tr) = vec2; end end a1=mean(log(growth(1:ts))); lam_S=exp(a1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %AT THIS POINT: % You have a sequence of states, % a sequence of single time-step growth rates, % a sequence of stage distributions, % a sequence of reproductive values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure(5) subplot(3,1,1) plot(3001:3100,growth(3001:3100),'-v'); xlabel('Time step'); ylabel('\lambda_t'); title('Section of sample path: Single time-step growth rates, stages, and rep values') subplot(3,1,2) area(uvecs(:,3001:3100)') ylabel('Stage distribution') subplot(3,1,3) area(vvecs(:,3001:3100)') ylabel('Reproductive values') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %YOU ALSO HAVE THE STOCHASTIC GROWTH RATE OVER THE ENTIRE SIMULATION figure(6) subplot(2,1,1) hist(growth) xlabel('\lambda_t') ylabel('frequency') title([' \lambda_S=' num2str(lam_S) ', log(\lambda_S)=a=' num2str(a1) ]) subplot(2,1,2) hist(log(growth)) xlabel('log(\lambda_t)') ylabel('frequency') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %MODULE 4 % generate elasticities: E^S, E^{S,\mu}, E^{S,\sigma}, E^S_\beta, % and summed habitat elasticity %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %set up matrices to receive the output sensE =zeros(nstage); %E^S, perturb the rate in every state by same percent sensA =zeros(nstage); %E^{S\mu}, perturb the mean only sensV =zeros(nstage); %E^{S\sigma}, pertub the variability only eijabS=zeros(nstate*nstage); %make hab-stage matrices for e's, intermediate step eijSbs=zeros(nstage,nstage,nstate); %E^S\beta %To estimate stochastic elasticity, we calculate an "elasticity" -looking object % at each time step, again this has to be for very big t, % and then average over t to get the desired result %start loop for i = (tr+1):(tr+ts-1) %choose index for states, starting one time step after transients %go as far as the state just before the (tr+ts)^th state %find state -- Keep in mind that: phenology of disturbance determines %whether you pick state(i) or state(i+1) for the matrix % % (note: may still need to be checked which is which again) %In either case that matrix becomes X(t)...Tulja .. u(t)=X(t)*u(t-1) % growth(t) is also named by the end of the interval... % growth(t) is the amount of growth from t-1 to t i2 = (states(i+1)); mat1 = ms(:,:,i2); itime = i+1-tr; % if this is: t i1 = i-tr; % then this is: t-1 mat2S = mat1; mat2A = avmat; mat2V = mat1 - avmat; mat2S = diag(vvecs(:,itime))*mat2S*diag(uvecs(:,i1)); mat2A = diag(vvecs(:,itime))*mat2A*diag(uvecs(:,i1)); mat2V = diag(vvecs(:,itime))*mat2V*diag(uvecs(:,i1)); scale1=( (vvecs(:,itime))' )*(uvecs(:,itime)); %makes scalar product of sd with rv mat2S = (1/(scale1*growth(itime)))*mat2S; %separate it out for hab-stage %use the growth rate at position 3 (t+1) sensE = sensE + mat2S; %accumulate these for E_{ij} mat2A = (1/(scale1*growth(itime)))*mat2A; %separate it out %use the growth rate at position 3 (t+1) sensA = sensA + mat2A; %accumulate these for E_{ij} mat2V = (1/(scale1*growth(itime)))*mat2V; %separate it out %use the growth rate at position 3 (t+1) sensV = sensV + mat2V; %accumulate these for E{ij} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % submodule: hab-stage blocking %lines are from "sensimod2.m" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% blockrow = states(i-1); blockcol = states(i-2); %so we are in block (blockrow, blockcol) rowdex = nstage*(blockrow-1) + (1:nstage); coldex = nstage*(blockcol-1) + (1:nstage); eijabS(rowdex, coldex) = eijabS(rowdex, coldex) + mat2S; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% sensE = sensE/(ts-1); sensA = sensA/(ts-1); sensV = sensV/(ts-1); eijabS = eijabS/(ts-1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%% calculate E^S_{beta}, and summed habitat elasticity eijS = zeros(nstage,nstage); ebS = zeros(nstate); for b = 1:nstate colx = nstage*(b-1); matrix1 = zeros(nstage,nstage); for a = 1:nstate rowx = nstage*(a-1); matrix1 = matrix1 + eijabS( (rowx + (1:nstage)), (colx + (1:nstage)) ); eijS = eijS + eijabS( (rowx + (1:nstage)), (colx + (1:nstage)) ); end eijSbs(:,:,b) = matrix1; %puts the matrix into the array ebS(b) = ebS(b) + sum(sum(matrix1)); end ebS=(ebS(:,1)); %%%%%%%%%%%%%%%%%%%end set for eijS_\beta*************** check1=sum(sum(eijabS)) %should sum to 1 check2=sensE-(sensA+sensV) %should be a matrix of zero's check3=sum(sum(sensE)) %should sum to 1 %%%%%%%%%%%%%ELASTICITY FIGURES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure(7);bar3(sensE);axis ij;title('E_{ij}^S') %perturbs the mean and deviation figure(8);bar3(sensA);axis ij;title('E_{ij}^{S\mu}') %perturbs only mean figure(9);bar3(sensV);axis ij;title('E_{ij}^{S\sigma}') %perturbs only deviation %%%%FIGURE TO COMPARE WITH MEGAMATRIX RESULTS figure(10)%for comparison to megamatrix waterfall(eijabS) axis ij, grid off; AXIS([1 56 1 56 0 0.2]); % manually sets the x, the y and the z axis % the x goes from 1 to 56, because there are 56 habitat-stage classes % the y goes from 1 to 56, " % the z goes from 0 0.2, set this high to be able to compare to another shading faceted; ylabel('Habitat-stages at t+1','FontName','Helvetica','FontSize',10,'FontWeight','bold') xlabel('Habitat-stages at t','FontName','Helvetica','FontSize',10,'FontWeight','bold') zlabel('Elasticity','FontName','Helvetica','FontSize',16,'FontWeight','bold') title('Stochastic sequence','FontName','Helvetica','FontSize',18,'FontWeight','bold') for i=1:nstate my_xvalues(i,:)=num2str(i); my_xticks(i)=nstage*(i-1)+1; % 8*(i-1) + 1 for JP data, try it for yours! end set(gca,'XTick',my_xticks,'FontName','Helvetica','FontSize',10) set(gca,'XTickLabel',{my_xvalues},'FontName','Helvetica','FontSize',10) set(gca,'YTick',my_xticks,'FontName','Helvetica','FontSize',10) set(gca,'YTickLabel',{my_xvalues},'FontName','Helvetica','FontSize',10) %%%FIGURES TO SHOW SHAPES OF HABITAT-STAGE ELASTICITY FOR EACH ENVIRONMENT for i=1:nstate; %make figures of habitat-stage elasticities FOR EACH \BETA figure(10+i); bar3(eijSbs(:,:,i));title(['E_{ij,\beta=' num2str(i) '}^S']); axis ij; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%