{VERSION 5 0 "SUN SPARC SOLARIS" "5.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0
1 0 0 0 0 1 }{CSTYLE "" -1 256 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 257 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1
258 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 259 "" 0 1 0 0
0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 0 1 0 0 0 0
0 0 0 1 }{CSTYLE "" -1 261 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 262 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 }{CSTYLE "" -1
263 "" 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 }{CSTYLE "" -1 264 "" 0 1 0 0
0 0 2 0 1 0 0 0 0 0 0 1 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 0 0 1 0 0 0
0 0 0 1 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Helvetica" 1 12 0
0 0 1 1 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Headi
ng 1" -1 3 1 {CSTYLE "" -1 -1 "Helvetica" 1 12 0 0 0 1 1 1 2 2 2 2 1
1 1 1 }1 1 0 0 6 6 1 0 1 0 2 2 0 1 }{PSTYLE "R3 Font 0" -1 256 1
{CSTYLE "" -1 -1 "Helvetica" 1 12 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0
0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "R3 Font 2" -1 257 1 {CSTYLE "" -1 -1 "
Courier" 1 14 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0
1 }{PSTYLE "R3 Font 2" -1 258 1 {CSTYLE "" -1 -1 "Courier" 1 12 0 0 0
1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Heading 2
" -1 259 1 {CSTYLE "" -1 -1 "Helvetica" 1 12 0 0 0 1 1 1 2 2 2 2 1 1
1 1 }1 1 0 0 8 2 1 0 1 0 2 2 0 1 }}
{SECT 0 {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1 " " }}}{EXCHG {PARA 0
"" 0 "" {TEXT 256 10 "Leslie.mws" }{TEXT -1 57 " (revised Sept. 8, \+
2001) Maple worksheet to accompany " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }
}{PARA 0 "" 0 "" {TEXT -1 28 "Brault and Caswell, 1993, " }{TEXT
262 40 "Pod-specific demography of killer-whales" }{TEXT -1 10 ", Ecol
ogy " }{TEXT 261 2 "74" }{TEXT -1 12 ": 1444-1454 " }}{PARA 0 "" 0 ""
{TEXT -1 4 "and " }}{PARA 0 "" 0 "" {TEXT -1 21 "Rodrigo Bernal,1998. \+
" }{TEXT 263 40 "Demography of the vegetable ivory palm " }{TEXT 264
19 "Pytelephas seemanii" }{TEXT 265 48 " in Colombia, and the impact o
f seed harvesting." }{TEXT -1 28 " Journal of Applied Ecology " }
{TEXT 260 2 "35" }{TEXT -1 7 ":64-74 " }}{PARA 0 "" 0 "" {TEXT -1 34 "
and other Leslie-Lefkowitch models" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 329 "Stage-based models are a modification of
the Leslie model. In this case, the assumption is that the stage of \+
an organism (size class) is a better indicator of reproductive perform
ance and survival than age is. This assumption works well for many pl
ants, and any organism in which size is correlated with survival and f
ecundity." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1
304 "The general form of the model treats the stages as elements of a \+
vector. The population projection matrix A (sometimes called a Lefkov
ich matrix, or, among mathematicians, transition matrix) is used to ca
lculate the number of individuals in each stage class at time t+1, bas
ed on the numbers at time t. " }}{PARA 0 "" 0 "" {TEXT -1 21 " \+
" }}{PARA 0 "" 0 "" {TEXT -1 10 " " }{TEXT 257
2 "n(" }{TEXT -1 3 "t+1" }{TEXT 258 9 ") = A n(" }{TEXT -1 1 "t" }
{TEXT 259 2 ") " }{TEXT -1 3 " " }}{PARA 0 "" 0 "" {TEXT -1 6 " \+
" }}{PARA 0 "" 0 "" {TEXT -1 44 "The (non-zero) elements of the matri
x A are:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1
68 "Pi the probabilities of surviving and remaining in the same stag
e." }}{PARA 0 "" 0 "" {TEXT -1 63 "Gi the probabilities of surviving a
nd moving to the next stage." }}{PARA 0 "" 0 "" {TEXT -1 90 "Fi the f
ertility, or number of female offspring at time t+1, per adult female
at time t." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "restart: wi
th(linalg): " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 73 "A:= matrix
(4,4, [ [0,F2,F3,F4], [G1,P2,0,0], [0,G2,P3,0], [0,0,G3,P4] ]);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 102 "Here are the respective matrices:
Caswell goes with the orca (whale) model and Bernal with palm trees.
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 111 "Caswell:= matrix(4, 4,
[ [0,0.0043,0.1132,0], [0.9775,0.9111,0,0],\n[0,0.0736,0.9534,0], [0,0
,0.0452,0.9804] ] );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 185 "Be
rnal:= matrix(6,6,[[0.7052,0,0,14.8173,14.9770,14.7732],\n[0.1548,0.63
39,0,0,0,0],\n[0,0.0216,0.9100,0,0,0],\n[0,0,0.0330,0.9614,0,0],\n[0,0
,0,0.0386,0.9220,0],\n[0,0,0,0,0.0188,0.9535]]);" }}}{EXCHG {PARA 0 ">
" 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 30 "Bern
al stable age distribution" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0
22 "A:= evalm( Bernal ); " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0
14 "n:= rowdim(A);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0
"" 0 "" {TEXT -1 257 " The elements on the top row are the reproductiv
e outputs. The subdiagonal elements represent the probability of tran
sition to the next stage. The diagonal elements represent the probabil
ity of remaining in a stage. here are some questions to think about.
" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 111 "*In \+
the palm tree model what fraction of seeds remain in the 'seed bank' e
ach year, that is, do not germinate? " }}{PARA 0 "" 0 "" {TEXT -1 0 "
" }}{PARA 0 "" 0 "" {TEXT -1 182 "**Of 1000 stage 1 individuals, how m
any survive to become stage 2? If there are currently 1000 stage 2 an
d 1000 stage 3, how many stage 1 individuals will be present after 1 y
ear? " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1
143 "***Suppose an individual is in the highest stage. What is the pr
obablility that this individual will still be alive 14 time steps from
now?***" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 ""
{TEXT -1 176 " Let's track what happens to a single burst of reproduct
ion, say 1000 seeds or newborns. Each time step is computed by taking
the matrix A times the current population vector." }}{PARA 0 "" 0 ""
{TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "pop||0:= vector( n \+
, [1000, 0, 0, 0, 0, 0]); " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0
44 "distrib= vector( n , [ 1, 0, 0, 0, 0, 0] ); " }}}{EXCHG {PARA 0 ">
" 0 "" {MPLTEXT 1 0 11 "Digits:=6; " }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 37 "points:= [0, 1000]; lpoints:= [0, 3];" }}}{EXCHG
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 142 "Next we'
re going to track the population as it changes. We'll go for 50 years
, but take a snapshot every 5 years. Feel free to run it longer." }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "for i from 0 by 5 to 50 do"
}}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1
0 38 "pop||(i+5):= evalm ( A^5 &* pop||i ): " }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 54 "total||(i+5):= sum( (pop||(i+5))['j'] , 'j' = 1 .. n)
:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 84 "points:= points, [i+5, total||
(i+5)]: lpoints:= lpoints, [i+5, log10(total||(i+5))]:" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 54 "distrib:= scalarmul( pop||(i+5), 1/ total||(i+
5)): od;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 266 "What app
ears to be happening? Is the population growing? What's happening in
the individual stages? to the overall distribution? Try it over again
with more than 50 years. If you get too much output change od; to od
: and use the following line for the final step." }}}{EXCHG {PARA 0 ">
" 0 "" {MPLTEXT 1 0 150 "pop||i:= evalm(pop||i); total||i:= total||i
; distrib:= evalm(distrib); five_year_rate:= total||i/total||(i-5); a
nnual_rate:= (five_year_rate)^(1/5);" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 66 "plot([points], title = \"Total population as a functi
on of time\" );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "plot([lpoints], title = \"Lo
g10 of total pop. vs. time\");" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT
1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 55 "Is there any pattern h
ere? What about in the long term?" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "
" 0 "" {TEXT -1 66 "Now let's compare this with the eigenvalue - eigen
vector analysis." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "eigsys:
=eigenvects(A): # (output printing suppressed)" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 96 "EVec:= seq( eigsys[i][3], i = 1 .. nops([eigsys])): #
eigenvectors (output printing suppressed)" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 65 "EVal:= seq( eigsys[i][1], i = 1 .. nops([eigsys])); #
eigenvalues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "EM:= seq( eigsys[i]
[2], i = 1 .. nops([eigsys])): # multiplicities" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 64 "magnitudes:=op( map( abs, [EVal] )); # magnitudes of \+
eigenvalues" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 195 "Out of the comple
te list of eigenvalues we find the maximum or \"dominant\" eigenvalue \+
(lambda = population growth rate), and the index to the corresponding \+
eigenvector (stable age distribution). " }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 24 "lambda:=max(magnitudes);" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 94 "member( lambda, [magnitudes], 'position'): indx:= po
sition: # isolate the dominant eigenvalue" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 82 "V:=op( EVec[indx] ); # find an eigenvector that is as
sociated with this eigenvalue" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "
" }}{PARA 0 "" 0 "" {TEXT -1 66 "How does this compare with the annual
rate that you found earlier?" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1
0 12 "annual_rate;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 507 "If the population grows for many \+
generations with constant transition probabilities given by matrix A, \+
the population age or stage vector will approach an eigenvector associ
ated with the dominant eigenvector of A. This is the stable age or s
tage distribution vector. Typically the stable stage vector w is scal
ed so that the sum of the elements is 1.0. To scale the stable stage \+
vector, we have to do some fiddling. Add up the elements in the vecto
r; then divide all elements in the vector by this sum:" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 45 "total:= sum( V['j'] , 'j' = 1 .. rowdim(A) );
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1
56 "We now have the stable age or stage distribution vector:" }}{PARA
0 "> " 0 "" {MPLTEXT 1 0 25 "w:=scalarmul(V, 1/total);" }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT
-1 78 "How does this compare with the distribution vector that you fou
nd iteratively?" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "evalm(di
strib);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 ""
{TEXT -1 34 "Let's check; is A w = (lambda) w? " }}{PARA 0 "" 0 ""
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 337 " Remember this is our de
finition of an eigenvalue and eigenvector - that the eigenvector multi
plied by the matrix is the same as the eigenvector multiplied by the e
igenvalue. If we got the correct eigenvalues and eigenvectors, we sho
uld obtain the same answer by multiplying the stable age distribution \+
by matrix A or by the eigenvalue:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0
41 "evalm( A &* w ); scalarmul( w , lambda);" }}}{EXCHG {PARA 0 "> "
0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 31 "Caswel
l stable age distribution" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35
"restart: with(linalg): with(plots):" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 110 "Caswell:=matrix(4, 4,[ [0,0.0043,0.1132,0], [0.9775,
0.9111,0,0],\n[0,0.0736,0.9534,0], [0,0,0.0452,0.9804] ] );" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "A:= evalm( Caswell ); " }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "n:= rowdim(A);" }}}{EXCHG
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 257 " The ele
ments on the top row are the reproductive outputs. The subdiagonal el
ements represent the probability of transition to the next stage. The \+
diagonal elements represent the probability of remaining in a stage. \+
here are some questions to think about." }}{PARA 0 "" 0 "" {TEXT -1 0
"" }}{PARA 0 "" 0 "" {TEXT -1 183 "*In the whale model what tells you \+
that juveniles do not reproduce? that older females do not reproduce? \+
Does this remind you of any other creature with a similar life stage p
attern?*" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1
182 "**Of 1000 stage 1 individuals, how many survive to become stage 2
? If there are currently 1000 stage 2 and 1000 stage 3, how many stag
e 1 individuals will be present after 1 year? " }}{PARA 0 "" 0 ""
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 143 "***Suppose an individual
is in the highest stage. What is the probablility that this individu
al will still be alive 14 time steps from now?***" }}{PARA 0 "" 0 ""
{TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 176 " Let's track wha
t happens to a single burst of reproduction, say 1000 seeds or newborn
s. Each time step is computed by taking the matrix A times the curren
t population vector." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> "
0 "" {MPLTEXT 1 0 40 "pop||0:= vector( n , [1000, 0, 0, 0 ]); " }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "distrib= vector( n , [ 1, 0,
0, 0] ); " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "Digits:=6; \+
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "points:= [0, 1000]; lpo
ints:= [0, 3];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "for i fro
m 0 by 10 to 100 do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "pop||(i+10):
= evalm ( A^10 &* pop||i ):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "tota
l||(i+10):= sum( pop||(i+10)['j'] , 'j' = 1 .. n):" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 54 "distrib:= scalarmul( pop||(i+10), 1/ total||(i+10) \+
): " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "points:= points, [(i+10), to
tal||(i+10)]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "lpoints:= lpoints,
[(i+10), log10(total||(i+10))]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "
od;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT
-1 383 "What appears to be happening? Is the population growing? Wha
t's happening in the individual stages? to the overall distribution? T
ry it over again with more than 20 steps. Suggestion: change the ; to \+
a : after the od to avoid screenfulls of output, and don't be too gree
dy. Let's just print the last computation and get an idea of the grow
th rate by comparing the last two totals." }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 146 "pop||i:= evalm(pop||i); total||i:= total||i; dis
trib:= evalm(distrib); ten_yr_rate:= total||i/total||(i-10); annual_ra
te:= (ten_yr_rate)^(1/10);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0
17 "plot( [points] );" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }
}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "plot( [lpoints] );" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 ""
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 66 "Now let's compare this wi
th the eigenvalue - eigenvector analysis." }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 54 "eigsys:=eigenvects(A): # (output printing suppress
ed)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 96 "EVec:= seq( eigsys[i][3], i \+
= 1 .. nops([eigsys])): # eigenvectors (output printing suppressed)"
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "EVal:= seq( eigsys[i][1],
i = 1 .. nops([eigsys])); # eigenvalues" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 66 "EM:= seq( eigsys[i][2], i = 1 .. nops([eigsys])): # m
ultiplicities" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "magnitudes:=op( ma
p( abs, [EVal] )); # magnitudes of eigenvalues" }}}{EXCHG {PARA 0 ""
0 "" {TEXT -1 261 "Out of the complete list of eigenvalues we find the
maximum or \"dominant\" eigenvalue (lambda = population growth rate),
and the index to the corresponding eigenvector (stable age distributi
on). How are these eigenvalues different from the ones we got with th
e " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "lambda:=max(magnitude
s);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "member( lambda, [magnitudes]
, 'position'): indx:= position:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20
"V:=op( EVec[indx] );" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 508 "If the population grows for many generat
ions with constant transition probabilities given by matrix A, the pop
ulation age or stage vector will approach the eigenvector associated w
ith the dominant eigenvector of A. This is the stable age or stage d
istribution vector. Typically the stable stage vector w is scaled so \+
that the sum of the elements is 1.0. To scale the stable stage vector
, we have to do some fiddling. Add up the elements in the vector; the
n divide all elements in the vector by this sum:" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 37 "total:= sum( V['j'] , 'j' = 1 .. n );" }}}{EXCHG
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 56 "We now ha
ve the stable age or stage distribution vector:" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 25 "w:=scalarmul(V, 1/total);" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 89 "Compare wit
h the annual growth rate and the long term distribution that we found \+
earlier." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "annual_rate; e
val(distrib);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0
"" {TEXT -1 34 "Let's check; is A w = (lambda) w? " }}{PARA 0 "" 0 ""
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 337 " Remember this is our de
finition of an eigenvalue and eigenvector - that the eigenvector multi
plied by the matrix is the same as the eigenvector multiplied by the e
igenvalue. If we got the correct eigenvalues and eigenvectors, we sho
uld obtain the same answer by multiplying the stable age distribution \+
by matrix A or by the eigenvalue:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0
41 "evalm( A &* w ); scalarmul( w , lambda);" }}}{EXCHG {PARA 0 "> "
0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 259 "" 0 "" {TEXT -1 46 "\"Fi
sh model\" and a look at reproductive value " }}{EXCHG {PARA 0 "> " 0
"" {MPLTEXT 1 0 37 "restart: with(linalg): with(plots):" }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "Fish:= matrix( 3, 3, [ [0, 0.8, 2.2
], [0.5, 0, 0], [0, 0.6, 0] ]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT
1 0 34 "A:= evalm( Fish ); n:= rowdim(A);" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 35 "pop||0:= vector( n , [100, 0, 0]); " }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "distrib= vector( n , [ 1, 0, 0] ); \+
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "total||0:= 100; " }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "points:= [0, 100]; lpoints:=
[0, 2];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "for i from 0 by
1 to 30 do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 35 "pop||(i+1):= evalm ( A &* pop||i ):" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 52 "total||(i+1):= sum( pop||(i+1)['j'] , 'j' = 1 \+
.. n):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "distrib:= scalarmul( pop|
|(i+1), 1/ total||(i+1) ): " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "poin
ts:= points, [(i+1), total||(i+1)]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0
48 "lpoints:= lpoints, [(i+1), log10(total||(i+1))]:" }}{PARA 0 "> "
0 "" {MPLTEXT 1 0 3 "od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 113 "pop||i:= evalm(pop||i)
; total||i:= total||i; distrib:= evalm(distrib); annual_rate:= (tota
l||i)/(total||(i-1));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "plot( [points] );" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 18 "plot( [lpoints] );" }}}{EXCHG {PARA 0 "" 0 ""
{TEXT -1 43 "Here's the stable age distribution anaysis." }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "eigsys:=eigenvects(A): # (output p
rinting suppressed)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 96 "EVec:= seq( \+
eigsys[i][3], i = 1 .. nops([eigsys])): # eigenvectors (output printi
ng suppressed)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "EVal:= seq( eigsy
s[i][1], i = 1 .. nops([eigsys])); # eigenvalues" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 66 "EM:= seq( eigsys[i][2], i = 1 .. nops([eigsys])): # m
ultiplicities" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "magnitudes:=op( ma
p( abs, [EVal] )); # magnitudes of eigenvalues" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 24 "lambda:=max(magnitudes);" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 60 "member( lambda, [magnitudes], 'position'): indx:= po
sition:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "V:=op( EVec[indx] );" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "total:= sum( V['j'] , 'j' = 1 .. n \+
):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT
-1 56 "We now have the stable age or stage distribution vector:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "w:=scalarmul(V, 1/total);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 ""
{TEXT -1 35 "Compare with our 30 year iteration:" }}}{EXCHG {PARA 0 ">
" 0 "" {MPLTEXT 1 0 28 "annual_rate; evalm(distrib);" }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT
-1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 155 "Here we transpose the A matrix
in order to set up calculation of the reproductive value vector, whic
h is the dominant eigenvector of the transposed matrix." }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 30 "AT:=transpose(A); Digits:= 10:" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 57 "eigsysT:=eigenvects(AT): # (output printing \+
suppressed)" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "
" {TEXT -1 120 "Find maximum eigenvalue and index to corresponding eig
envector; this eigenvector is the reproductive value distribution." }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "EVecT:= seq( eigsysT[i][3], i = 1 .
. nops([eigsysT]) ): # eigenvectors" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 69 "EValT:= seq( eigsysT[i][1], i = 1 .. nops([eigsysT]) \+
); # eigenvalues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "EMT:= seq( eigs
ysT[i][2], i = 1 ..nops([eigsysT]) ): # multiplicities" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 137 "magsT:= op( map( abs , [EValT] )); # magnitud
es of eigenvalues (Note: [ ] turns a sequence into a list; op turns a \+
list into a sequence.)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "m
u:= max( magsT); # get dominant eigenvalue" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 53 "member(mu, [magsT], 'position'): indexT:= position:
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1
76 "Convert eigenvector lists to actual vectors for use in sensitivity
analysis." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "VT:= op( EVecT[indexT
] ); " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 ""
{TEXT -1 309 "This is the reproductive value distribution. Note that \+
vector v is scaled so the first element is 1.0. We do a scalar multi
plication of all elements in the vector by the reciprocal of the first
element. (Recall that a scalar multiple of an eigenvector is again a
n eigenvector with the same eigenvalue.) " }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 28 "v:= scalarmul(VT, 1/VT[1] );" }}}{EXCHG {PARA 0 "> "
0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 180 "What in
formation does this vector give us? Which age or stage group is contri
buting the most to reproduction, and why is this plausible in terms of
both fecundity and survivorship?" }}{PARA 0 "" 0 "" {TEXT -1 1 "." }}
}{EXCHG {PARA 0 "" 0 "" {TEXT -1 108 "Let's confirm that we correctly \+
identified the eigenvalue/eigenvector pair by checking that AT v = (mu
) v. " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "evalm( AT &* v ); scala
rmul( v , mu);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}
{SECT 1 {PARA 3 "" 0 "" {TEXT -1 43 "Sensitivity analysis for Bernal \+
(optional)" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 ""
{TEXT -1 831 "Now we will calculate the sensitivity of the system to c
hanges in fecundity. This is the analysis proposed in Caswell. We ar
e analytically determining the partial derivative of lambda (populatio
n growth rate) with respect to each element in the projection matrix. \+
Thus we are asking: given a change in a particular element in the pro
jection matrix, what is the corresponding change in population growth \+
rate? This allows us to determine which components of the life histor
y have the biggest effects on population growth, and which may be most
sensitive to natural selection. The first parameter we need is the i
nner product of the two eigenvectors (stable age distribution and repr
oductive value distribution). This is symbolized as in the Ca
swell handout. This value is used as the denominator of the calculati
on. " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart:" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "with(linalg):" }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 185 "Bernal:= matrix(6,6,[[0.7052,0,0,1
4.8173,14.9770,14.7732],\n[0.1548,0.6339,0,0,0,0],\n[0,0.0216,0.9100,0
,0,0],\n[0,0,0.0330,0.9614,0,0],\n[0,0,0,0.0386,0.9220,0],\n[0,0,0,0,0
.0188,0.9535]]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "A:=eval
m(Bernal):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "eigsys:=eigen
vects(A): # (output printing suppressed)" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 96 "EVec:= seq( eigsys[i][3], i = 1 .. nops([eigsys])): #
eigenvectors (output printing suppressed)" }}}{EXCHG {PARA 0 "> " 0
"" {MPLTEXT 1 0 65 "EVal:= seq( eigsys[i][1], i = 1 .. nops([eigsys]))
; # eigenvalues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "EM:= seq( eigsys
[i][2], i = 1 .. nops([eigsys])): # multiplicities" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 64 "magnitudes:=op( map( abs, [EVal] )); # magnitudes o
f eigenvalues" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1
0 24 "lambda:=max(magnitudes);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "m
ember( lambda, [magnitudes], 'position'): indx:= position:" }}{PARA
0 "> " 0 "" {MPLTEXT 1 0 20 "V:=op( EVec[indx] );" }}}{EXCHG {PARA 0 "
> " 0 "" {MPLTEXT 1 0 45 "total:= sum( V['j'] , 'j' = 1 .. rowdim(A) )
;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 45 "Calculate the scaled stable \+
age distribution:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "w:=scalarmul(V
,1/total);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 145 "Transpose the matr
ix in order to set up calculation of the reproductive value vector, wh
ich is the dominant eigenvector of the transposed matrix." }}{PARA 0 "
> " 0 "" {MPLTEXT 1 0 17 "AT:=transpose(A);" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 57 "eigsysT:=eigenvects(AT): # (output printing suppres
sed)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 ""
{TEXT -1 120 "Find maximum eigenvalue and index to corresponding eigen
vector; this eigenvector is the reproductive value distribution." }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "EVecT:= seq( eigsysT[i][3], i = 1 .
. nops([eigsysT]) ): # eigenvectors" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 69 "EValT:= seq( eigsysT[i][1], i = 1 .. nops([eigsysT]) \+
); # eigenvalues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "EMT:= seq( eigs
ysT[i][2], i = 1 ..nops([eigsysT]) ): # multiplicities" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 137 "magsT:= op( map( abs , [EValT] )); # magnitud
es of eigenvalues (Note: [ ] turns a sequence into a list; op turns a \+
list into a sequence.)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "m
u:= max( magsT); # get dominant eigenvalue" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 54 "member(mu, [magsT], 'position'): indexT:= position
:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1
114 "Convert eigenvector lists to actual vectors for use in sensitivit
y analysis. VT is the reproductive value vector." }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 26 "VT:= op( EVecT[indexT] ); " }}}{EXCHG {PARA 0 "" 0
"" {TEXT -1 78 "v is the scaled reproductive value vector (scaled so t
he first element is 1.0)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "v:=scal
armul(VT,1/VT[1]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 114 "Sensitivit
y to fecundity is related both to the stable age distribution and to t
he reproductive value distribution" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 53 "IP:= innerprod(VT , V); # dot, or scalar, product \+
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 74 "sensitivity_to_fecundity:= [seq
( V[i]*VT[1] / IP , i = 1 .. rowdim(A) )]; " }}}{EXCHG {PARA 0 "> " 0
"" {MPLTEXT 1 0 0 "" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 111 "Observe that
it is immaterial whether we use the normalized vectors or not, as lon
g as we do this consistently." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1
0 59 "[seq( w[i]*v[1] /innerprod( w, v) , i = 1 .. rowdim(A) )]; " }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 215 "In terms of the paper this corres
ponds to the top row of matrix S, equation 7, in Caswell except that t
he first and last elements have been suppressed because the original m
atrix had zero entries in these locations." }}{PARA 0 "" 0 "" {TEXT
-1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 67 "We do a similar calcuation to o
btain the sensitivity to survival: " }}{PARA 0 "> " 0 "" {MPLTEXT 1
0 76 "sensitivity_to_survival:= [seq( V[i]*VT[i+1] / IP , i = 1 .. row
dim(A)-1 )];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 401 "This vector corr
esponds to the subdiagonal of matrix S. The entire sensitivity matri
x can be calcuated by generalizing the above operations. Note there a
re some entries that can be taken to be zero, or otherwise ignored, be
cause the A matrix had zero entries (but this isn't really necessary, \+
since as we shall see in a moment we end up multiplying each entry of \+
S by the corresponding entry of A)." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0
106 "S:=matrix(rowdim(A),rowdim(A), [ seq( seq( V[i] * VT[j] / IP, i =
1.. rowdim(A)), j = 1.. rowdim(A) ) ]);" }}}{EXCHG {PARA 0 "" 0 ""
{TEXT -1 418 "The elasticity values are calculated by scaling the sens
itivities by aij / lambda (Caswell handout). The resulting matrix is \+
reported as Table 4 on p 70. The elasticities are a bit easier to int
erpret than the sensitivities, because they represent the proportional
contributions of matrix elements to lambda. Recall the we can't use \+
the label E in Maple because this means the base of natural logarithms
(2.718...)." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 111 "e:=matrix(rowdim(A
),rowdim(A),[ seq( seq( S[j,i] * A[j,i] / lambda, i = 1 .. rowdim(A)),
j = 1.. rowdim(A)\\)]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 401 "***Try seeing if you believe that the se
nsitivities or elasticities are correct, following the pattern illustr
ated below.*** We will change one of the entries in the original \"A
\" matrix by 10% and see what percent change there is in lambda. Comp
are this change with the corresponding value in the elasticity or sens
itivity matrix. Let's first recall the original matrix, eigenvalues,
eigenvectors:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "A:= evalm(A);" }}
}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "lambda[old]:= lambda; mu[
old]:= mu;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "RV:= evalm(v)
; SAD:= evalm(w);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 428 "Let's chan
ge the P2 value by10% (we'll do this by multiplying P2 by 1.1) and see
what percent change we get in lambda. This is the so-called [2,2] en
try of the matrix A, that is, the number that appears in the first row
(from top to bottom) and second column (from left to right). If the
author's description of elasticity is correct, then we should get a c
hange in lambda that is proportional to the elasticity value for P2."
}}{PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "A[2,2]:=1.1*A[2,2]; # make the \+
change in the affected entry" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1
0 41 "A:=evalm(A); # confirm the matrix is OK" }}}{EXCHG {PARA 0 ""
0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 77 "Now let's calculate \+
the eigenvalues and eigenvectors for this revised matrix." }}{PARA 0 "
> " 0 "" {MPLTEXT 1 0 54 "eigsys:=eigenvects(A): # (output printing s
uppressed)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 67 "EVec:= seq( e
igsys[i][3], i = 1 .. nops([eigsys])): # eigenvectors " }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "EVal:= seq( eigsys[i][1], i = 1 .. \+
nops([eigsys])); # eigenvalues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "m
ags:= op( map( abs, [EVal] )); lambda[new]:=max(mags);" }}}{EXCHG
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 63 "Compare t
his new value for lambda with the old value of lambda:" }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "lambda[old];" }}}{EXCHG {PARA 0 ""
0 "" {TEXT -1 163 "This is the ratio of the new value to the old; the \+
fractional change in lambda is the portion after the decimal. Did you
expect lthe new lambda to be larger? Why?" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 57 "ratio:= lambda[new]/lambda[old]; frac_change:= rati
o - 1;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 77 "Here is the elasticity \+
value. Compare it to the fractional change in lambda." }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "elas:= e[2,2];" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 410 "Note that we made a 10 percent change in one elem
ent of the matrix (a fractional change of 0.1). If you multiply this \+
elasticity value by the fractional change made in the matrix element, \+
you should get a number close to the fractional change in lambda. If
you made a 100% change (a fractional change of 1.0) in the matrix ele
ment, you should see a fractional change in lambda equal to the elasti
city value." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "(elas * 0.1)
- (frac_change); # close to zero?" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 59 "member( lambda[new], [mags], 'position'): indx:= pos
ition:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "evalm(SAD); # here's the \+
original stable age distribution" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1
47 "\nCompare with the new stable age distribution. " }}}{EXCHG {PARA
0 "> " 0 "" {MPLTEXT 1 0 120 "V := op( EVec[indx] ); total:= sum( V['k
'] ,' k' = 1 .. rowdim(A) ): w:= scalarmul(V, 1/total); # compare with
SAD above" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 147 "Transpose the A ma
trix in order to set up calculation of the reproductive value vector, \+
which is the dominant eigenvector of the transposed matrix." }}{PARA
0 "> " 0 "" {MPLTEXT 1 0 17 "AT:=transpose(A);" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 57 "eigsysT:=eigenvects(AT): # (output printing suppres
sed)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "EVecT:= seq( eigsys
T[i][3], i = 1 .. rowdim(A) ): # eigenvectors" }}}{EXCHG {PARA 0 "> "
0 "" {MPLTEXT 1 0 63 "EValT:= seq( eigsysT[i][1], i = 1 .. rowdim(A) )
; # eigenvalues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "magsT:= op( map(
abs , [EValT] )); # magnitudes of eigenvalues " }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 82 "mu[new]:= max( magsT); member(mu[new], [magsT], 'po
sition'): indexT:= position:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1
0 78 "VT:= op( EVecT[indexT] ): v:= scalarmul(VT, 1/VT[1] ); # compare
with RV above" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "evalm(RV)
;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 ""
0 "" {TEXT -1 365 "***Now, how did v and w change, given that 10% chan
ge in a single entry of the matrix A? Considering the change that was
made can you give an explanation for the result? Now change the surv
ival probability G3 by 15% (multiply it by 1.5). (Don't forget to res
et A[1,2] back to its original value first) What effect do you expect
to see? What happens in fact?***" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 252 "The author states on p71 that \"an 86% r
eduction in fecundity of all classes is required to reduce lambda from
1.0589 to 1.0\". Try making reductions in fecundity and see if your \+
observations confirm this: First we will get an original copy of matr
ix A" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "A:=evalm(Bernal);"
}}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT
-1 53 "Then we define our fractional reduction in fecundity." }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "fecundity_reduction:=0.86;"
}}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 271 "We will now reduce the seed pr
oduction entries by this amount. We will not change A[1,1], because t
hat element of the matrix defines the rate of retention of seeds in th
e seed bank, not a fecundity. Therefore we change only elements 2 to \+
6 in the top row of the matrix:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT
1 0 28 "for i from 2 to rowdim(A) do" }}{PARA 0 "> " 0 "" {MPLTEXT 1
0 46 " A[1,i]:=A[1,i]*(1-fecundity_reduction); od:" }}}{EXCHG {PARA
0 "> " 0 "" {MPLTEXT 1 0 9 "evalm(A);" }}}{EXCHG {PARA 0 "" 0 ""
{TEXT -1 235 "Now move back to the top of this section of the workshee
t, and rerun from just below the original definition of matrix A. See
if the values for lambda, stable stage distribution, and reproductive
value distribution change as expected." }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 44 "Sensitivity
analysis for Caswell (optional)" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1
99 "FIrst we need to get the eigenvalues and eigenvectors of both the \+
original matrix and its transpose" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 22 "restart: with(linalg):" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 111 "Caswell:= matrix(4, 4,[ [0,0.0043,0.1132,0], [0.9775
,0.9111,0,0],\n[0,0.0736,0.9534,0], [0,0,0.0452,0.9804] ] );" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "A:= evalm(Caswell):" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "eigsys:=eigenvects(A): # (o
utput printing suppressed)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 96 "EVec:
= seq( eigsys[i][3], i = 1 .. nops([eigsys])): # eigenvectors (output
printing suppressed)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "EV
al:= seq( eigsys[i][1], i = 1 .. nops([eigsys])); # eigenvalues" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "EM:= seq( eigsys[i][2], i = 1 .. no
ps([eigsys])): # multiplicities" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "
magnitudes:=op( map( abs, [EVal] )); # magnitudes of eigenvalues" }
{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "lambda:=max
(magnitudes);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "member( lambda, [m
agnitudes], 'position'): indx:= position:" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 20 "V:=op( EVec[indx] );" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 45 "total:= sum( V['j'] , 'j' = 1 .. rowdim(A) );" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 45 "Calculate the scaled stable age di
stribution:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "w:=scalarmul(V,1/tot
al);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 145 "Transpose the matrix in \+
order to set up calculation of the reproductive value vector, which is
the dominant eigenvector of the transposed matrix." }}{PARA 0 "> " 0
"" {MPLTEXT 1 0 17 "AT:=transpose(A);" }}{PARA 0 "> " 0 "" {MPLTEXT 1
0 57 "eigsysT:=eigenvects(AT): # (output printing suppressed)" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT
-1 120 "Find maximum eigenvalue and index to corresponding eigenvector
; this eigenvector is the reproductive value distribution." }}{PARA 0
"> " 0 "" {MPLTEXT 1 0 70 "EVecT:= seq( eigsysT[i][3], i = 1 .. nops([
eigsysT]) ): # eigenvectors" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0
69 "EValT:= seq( eigsysT[i][1], i = 1 .. nops([eigsysT]) ); # eigenval
ues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "EMT:= seq( eigsysT[i][2], i \+
= 1 ..nops([eigsysT]) ): # multiplicities" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 137 "magsT:= op( map( abs , [EValT] )); # magnitudes of e
igenvalues (Note: [ ] turns a sequence into a list; op turns a list in
to a sequence.)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "mu:= max
( magsT); # get dominant eigenvalue" }}{PARA 0 "> " 0 "" {MPLTEXT 1
0 53 "member(mu, [magsT], 'position'): indexT:= position:" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 114 "C
onvert eigenvector lists to actual vectors for use in sensitivity anal
ysis. VT is the reproductive value vector." }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 26 "VT:= op( EVecT[indexT] ); " }}}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 78 "v is the scaled reproductive value vector (scaled so th
e first element is 1.0)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "v:=scala
rmul(VT,1/VT[1]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 940 "N
ow we will calculate the sensitivity of the system to changes in fecun
dity. This is the analysis proposed on p. 1447 in equation 3. We are
analytically determining the partial derivative of lambda (population
growth rate) with respect to each element in the projection matrix. \+
Thus we are asking: given a change in a particular element in the proj
ection matrix, what is the corresponding change in population growth r
ate? This allows us to determine which components of the life history
have the biggest effects on population growth, and which may be most \+
sensitive to natural selection. The first parameter we need is the in
ner product of the two eigenvectors (stable age distribution and repro
ductive value distribution). This is symbolized as in equatio
n 3 on p. 1447. This value is used as the denominator on the right ha
nd size of equation 3. Observe that it is immaterial whether we use t
he normalized vecotrs or not." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1
0 53 "IP:= innerprod(VT , V); # dot, or scalar, product " }}{PARA
0 "> " 0 "" {MPLTEXT 1 0 66 "sensitivity_to_fecundity:= [seq( V[i]*VT[
1] / IP , i = 1 .. 4 )]; " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 111 "Observe that it is immaterial whether we
use the normalized vectors or not, as long as we do this consistently
." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "[seq( w[i]*v[1] /inner
prod( w, v) , i = 1 .. 4 )]; " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "
" }}{PARA 0 "" 0 "" {TEXT -1 215 "In terms of the paper this correspon
ds to the top row of matrix S on p. 1448, equation 7, except that the \+
first and last elements have been suppressed because the original matr
ix had zero entries in these locations." }}{PARA 0 "" 0 "" {TEXT -1 0
"" }}{PARA 0 "" 0 "" {TEXT -1 67 "We do a similar calcuation to obtain
the sensitivity to survival: " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "
sensitivity_to_survival:= [seq( V[i]*VT[i+1] / IP , i = 1 .. 3 )];" }}
}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 401 "
This vector corresponds to the subdiagonal of matrix S. The entire s
ensitivity matrix can be calcuated by generalizing the above operation
s. Note there are some entries that can be taken to be zero, or other
wise ignored, because the A matrix had zero entries (but this isn't re
ally necessary, since as we shall see in a moment we end up multiplyin
g each entry of S by the corresponding entry of A)." }}{PARA 0 "> " 0
"" {MPLTEXT 1 0 74 "S:=matrix(4,4, [ seq( seq( V[i] * VT[j] / IP, i = \+
1.. 4), j = 1.. 4 ) ]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 420 "The elasticity values are calculated by \+
scaling the sensitivities by aij / lambda (equation 4 on p. 1447). Th
e resulting matrix is reported as equation 8 on p. 1448. The elastici
ties are a bit easier to interpret than the sensitivities, because the
y represent the proportional contributions of matrix elements to lambd
a. Recall the we can't use E in Maple because this means the base of \+
natural logarithms (2.718...)." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 78 "e
:=matrix(4,4,[ seq( seq( S[j,i] * A[j,i] / lambda, i = 1 .. 4), j = 1.
. 4)]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 ""
{TEXT -1 401 "***Try seeing if you believe that the sensitivities or e
lasticities are correct, following the pattern illustrated below.*** \+
We will change one of the entries in the original \"A\" matrix by 10% \+
and see what percent change there is in lambda. Compare this change w
ith the corresponding value in the elasticity or sensitivity matrix. \+
Let's first recall the original matrix, eigenvalues, eigenvectors:" }
}{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "A:= evalm(A);" }}}{EXCHG {PARA 0 "
> " 0 "" {MPLTEXT 1 0 37 "lambda[old]:= lambda; mu[old]:= mu;" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "RV:= evalm(v); SAD:= evalm(
w);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT
-1 416 "Let's change the F2 value from 0.0043 to 0.00473 (a 10% change
) and see what percent change we get in lambda. This is the so-called
[1,2] entry of the matrix A, that is, the number that appears in the \+
first row (from top to bottom) and second column (from left to right).
If Caswell's description of elasticity is correct, then we should g
et a change in lambda that is proportional to the elasticity value for
F2." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "A[1,2]:=.00473; # make th
e change in the affected entry" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT
1 0 41 "A:=evalm(A); # confirm the matrix is OK" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 77 "Now let's calculat
e the eigenvalues and eigenvectors for this revised matrix." }}{PARA
0 "> " 0 "" {MPLTEXT 1 0 54 "eigsys:=eigenvects(A): # (output printin
g suppressed)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "EVec:= seq
( eigsys[i][3], i = 1 .. 4): # eigenvectors " }}}{EXCHG {PARA 0 "> "
0 "" {MPLTEXT 1 0 52 "EVal:= seq( eigsys[i][1], i = 1 .. 4); # eigenva
lues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "mags:= op( map( abs, [EVal]
)); lambda[new]:=max(mags);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "
" }}{PARA 0 "" 0 "" {TEXT -1 118 "Compare this new value for lambda wi
th the \"estimate\" of lambda at bottom of p 1447. This is the old va
lue of lambda:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "lambda[ol
d];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 163 "This is the ratio of the \+
new value to the old; the fractional change in lambda is the portion a
fter the decimal. Did you expect lthe new lambda to be larger? Why?"
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "ratio:= lambda[new]/lambd
a[old]; frac_change:= ratio - 1;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1
77 "Here is the elasticity value. Compare it to the fractional change
in lambda." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "elas:= e[1,2
];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 390 "Note that we made a 10 per
cent change in fecundity (a fractional change of 0.1). If you multipl
y this elasticity value by the fractional change made in the fecundity
, you should get a number close to the fractional change in lambda. \+
If you made a 100% change (a fractional change of 1.0) in the fecundit
y value, you should see a fractional change in lambda equal to the ela
sticity value." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "(elas * 0
.1) - (frac_change); # close to zero?" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 63 "member( lambda[new], [mags], 'position'): location:=
position:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "evalm(SAD); #
here's the original stable age distribution" }}}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 47 "\nCompare with the new stable age distribution. " }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 116 "V := op( EVec[location] ); \+
total:= sum( V['k'] ,' k' = 1 .. 4 ):\nw:= scalarmul(V, 1/total); # co
mpare with SAD above" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 147 "Transpos
e the A matrix in order to set up calculation of the reproductive valu
e vector, which is the dominant eigenvector of the transposed matrix.
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "AT:=transpose(A);" }}{PARA 0 ">
" 0 "" {MPLTEXT 1 0 57 "eigsysT:=eigenvects(AT): # (output printing
suppressed)" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "EVecT:= seq
( eigsysT[i][3], i = 1 .. 4 ): # eigenvectors" }}}{EXCHG {PARA 0 "> "
0 "" {MPLTEXT 1 0 55 "EValT:= seq( eigsysT[i][1], i = 1 .. 4 ); # eige
nvalues" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "magsT:= op( map( abs , [
EValT] )); # magnitudes of eigenvalues " }}{PARA 0 "> " 0 "" {MPLTEXT
1 0 82 "mu[new]:= max( magsT); member(mu[new], [magsT], 'position'):
indexT:= position:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 78 "VT:
= op( EVecT[indexT] ): v:= scalarmul(VT, 1/VT[1] ); # compare with RV \+
above" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "evalm(RV);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 372 "***Now, ho
w did v and w change, given that 10% change in a single entry of the m
atrix A? Considering the change that was made can you give an explana
tion for the result? Now change the survival probability G3 by 15% fr
om 0.0452 to 0.05198. (Don't forget to reset A[1,2] back to its origi
nal value of 0.0043!) What effect do you expect to see? What happens
in fact?***" }}}}}{MARK "0 0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }
{PAGENUMBERS 0 1 2 33 1 1 }