matpopmod.model

This module implements the class MPM used to represent matrix population models. It provides a flexible way to create matrix population models and a high-level interface to study their mathematical properties and biological descriptors (growth rate, stable distribution, reproductive values, etc). The descriptors are computed when they are first needed, and then stored as immutable attributes.

>>> m = MPM(A = "0.1 2; 0.4 0.8")
>>> m.primitive # is the projection matrix primitive?
True
>>> m.lmbd # asymptotic growth rate
1.4104686356149274
>>> m.w # stable distribution
array([0.60414407, 0.39585593])
>>> m.lmbd = 0 # trying to change a descriptor throws an error.
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: can't set attribute

In order to be computed, many descriptors require knowing which entries of the projection matrix correspond to survival probabilities and which correspond to fertilities, i.e. its decomposition into the survival matrix S and the fertility matrix F. Trying to access such descriptors without having specified the matrices S and F when creating the model will raise NotAvailable:

>>> m1 = MPM(A = "0.1 2; 0.4 0.8")
>>> m1.lmbd
1.4104686356149274
>>> m1.R0
Traceback (most recent call last):
  ...
matpopmod.utils.NotAvailable: 'A = S + F decomposition required.'
>>> m2 = MPM(S = "0.1 0; 0.4 0.8", F = "0 2; 0 0")
>>> m2.lmbd
1.4104686356149274
>>> m2.R0
4.4444444444444455

In what follows, we detail how to create MPM objects and we give the complete list of descriptors that can be computed with them.

class matpopmod.model.MPM(A=None, S=None, F=None, fertilities=None, delimiter=None, metadata=None)

A matrix population model is defined by its projection matrix A. Some descriptors (such the growth rate, the stable distribution or the reproductive values) can be computed directly from this projection matrix. However, other ones (such as the net reproductive rate or the generation time) can be computed only if the A = S + F decomposition of the projection matrix into its survival and fertility components is known. Thus, there are two ways to create a matrix population model:

• Using the projection matrix A only:

  MPM(A = "0.1 2; 0.4 0.8")
  

  If the model is created that way, trying to access descriptor that require knowing S and F will raise NotAvailable.

• By specifying any two of the three matrices A, S and F, for instance:

  MPM(S = "0.1 0; 0.4 0.8", F = "0 2; 0 0")
  

  If the three matrices are given, A will be silently ignored (that is, only S and F will be used, without checking whether A = S + F).

  For convenience, instead of specifying two matrices it is also possible to provide A and the list of entries that correspond to reproduction. This is done with the fertility argument:

  MPM(A = "0.1 2; 0.4 0.8", fertilities = [(0, 1)])
  

  However this is less flexible, since this does not allow for overlapping entries between S and F. Note that the indexing of the entries starts at 0.

In addition to this, there are several ways to specify any of the matrices in the arguments:

• Using a string, as we saw in the examples above. Rows have to be separated by a semicolon and columns by a space and/or a comma.

• Using Python iterables, such as a list of lists:

  MPM(A = [[0.1, 2], [0.4, 0.8]])
  

  Valid arguments are the same as for NumPy’s array() function. In particular, NumPy arrays can be used.

• Reading the matrix from a text file. In that case, the name of the file containing the matrix must be passed as a Path from Python’s pathlib module.

  file_A = pathlib.Path("~/myfolder/myfile.txt")
  MPM(A = file_A)
  

  Valid files contain one row of the matrix per line, with columns separated by whitespace (the character separating columns can be specified using the delimiter argument – use "," to load CSV files). Lines starting with #, % or // are ignored.

Finally, each model has a metadata attribute. It is simply a dictionary that can be used to store relevant information about the model.

Complete list of available descriptors

Unless a specific reference is given, the formulas used are classic and can be found in textbooks such as [Casw00] or [KeCa05].

1. Descriptors of the projection matrix
   1. Structural properties
      1. The matrix A
      2. Irreducibility and related quantities
      3. Aperiodicity and related quantities
      4. Primitivity and quasi-primitivity
   2. Eigen-elements and related quantities
      1. Asymptotic dynamics
      2. Second order dynamics
      3. Other eigen-elements
      4. Derivatives (sensitivities and elasticities)
   3. The genealogical Markov chain
      1. Transition matrix and stationary distribution
      2. Entropy rate and other statistics
2. Descriptors of the survival and fertility matrices
   1. The A = S + F decomposition
      1. Matrices
      2. Leslie and Usher models
      3. Newborn and reproductive classes
      4. Steady state distributions
      5. Survival and fertility excesses
   2. Biological descriptors
      1. Reproductive output
      2. Net reproductive rate
      3. Generation time
      4. Life expectancy
      5. Age structure
      6. Entropies
3. Functions of the initial population
   1. Distance to the stable population
   2. Trajectories

I.  Descriptors of the projection matrix

I.1.  Structural properties

I.1.a.  The matrix A

A

The projection matrix $\mathbf{A}=(a_{ij})$, where $a_{ij}$ is the expected per-capita contribution of individuals of class j to the abundance of class i at the next time-step.

dim

The number of classes of the model.

I.1.b.  Irreducibility and related quantities

irreducible

Whether the projection matrix A is irreducible, meaning for any two classes i and j, there exists a directed path from i to j in the life-cycle graph. Equivalently, a non-negative matrix A is irreducible if and only if for all i, j there exists k such that the (i, j)-th entry of ${\mathbf{A}}^k$ is positive.

See is_irreducible() from the module mathtools for a function that can be used on any NumPy array and for details about the implementation.

irreducible_components

The decomposition of the projection matrix A into irreducible components, as a couple (submatrices, classes) where submatrix[k] is the submatrix of A that corresponds to the k-th irreducible component and classes[k][i] is the class of A that corresponds to the i-th class of that component. Example:

>>> m = MPM(A = "1 0 0 0; 2 0 3 0; 0 4 0 0; 0 0 5 0")
>>> m.A
array([[1, 0, 0, 0],
       [2, 0, 3, 0],
       [0, 4, 0, 0],
       [0, 0, 5, 0]])
>>> m.irreducible_components
([array([[1]]), array([[0, 3], [4, 0]])], [[0], [1, 2]])

Note that classes are indexed from 0 to n - 1, where n = dim(A), and that classes that do not belong to any directed cycle (as is the case above with class 3) are not part of any irreducible component.

See strongly_connected_components() from the module mathtools for a function that can be used on any NumPy array and for details about the implementation. See also normal_form() for an alternative representation of the decomposition that retains the transitions between the classes.

normal_form

A normal form of the projection matrix A is obtained by reordering its classes so as to get a projection matrix of the form

$$\begin{array}{r}\left(\begin{array}{cccc}{\mathbf{A}}_{11}& \mathbf{0}& \cdots & \mathbf{0}\\ {\mathbf{A}}_{21}& {\mathbf{A}}_{22}& \ddots & ⋮\\ ⋮& ⋮& \ddots & \mathbf{0}\\ {\mathbf{A}}_{2m}& {\mathbf{A}}_{2m}& \cdots & {\mathbf{A}}_{mm}\end{array}\right),\end{array}$$

where each of the diagonal blocks ${\mathbf{A}}_{11},\dots ,{\mathbf{A}}_{mm}$ is either an irreducible matrix or a null matrix – see e.g. Chapter XIII, §4 of [Gant59].

The new matrix is returned with a NumPy array permutation representing the reordering of the classes: permutation[i] is the index of the row/column of A that corresponds to the i-th row/column in the normal form. Indexing starts at 0. Example:

>>> m = MPM(A = "1 0 2 0; 3 0 0 4; 5 0 0 0; 0 6 7 0")
>>> m.A
array([[1, 0, 2, 0],
       [3, 0, 0, 4],
       [5, 0, 0, 0],
       [0, 6, 7, 0]])
>>> m.normal_form
(array([[1, 2, 0, 0],
        [5, 0, 0, 0],
        [3, 0, 0, 4],
        [0, 7, 6, 0]]), array([0, 2, 1, 3]))

Note that the normal form is not unique, as we can reorder the classes inside each irreducible component to get a different one (if some of the blocks under the diagonal are null, we may also be able to permutate the irreducible blocks). This implementation tries to preserve the original order of the classes. Also note that several authors (e.g, [Varg00] and [HoJo13]) define the normal form to be block upper-triangular, rather than block lower-triangular. However, this is at odds with the usage in matrix population models, where projection matrices are typically lower-triangular (e.g, Leslie matrices).

See irreducible_components for a function that extracts the submatrices corresponding to the irreducible components.

quasi_irreducible

Whether the projection matrix A is quasi-irreducible, i.e. whether it has only one irreducible component whose dominant eigenvalue is equal to λ.

Unlike irreducibility, quasi-irreducibility is not a standard notion. However, it is useful to treat models with post-reproductive classes or metapopulation models with uni-directional migration: indeed, the vast majority of matrix population are quasi-irreducible, and many of the results for irreducible matrices are easily adapted to quasi-irreducible ones. See quasi_primitive for more information.

I.1.c.  Aperiodicity and related quantities

periods

The vector of periods of the projection matrix A. The period of the class i is the greatest common divisor of the lengths of the cycles goings through the corresponding vertex in the life-cycle graph (if there is no such cycle, then it is None). See also index_of_imprimitivity for the greatest common divisor of the periods the classes.

See periods() from the mathtools module for a function that can be used on any NumPy array and for details about the implementation.

index_of_imprimitivity

The index of imprimitivity of the projection matrix A, also known as its period.

When A is irreducible, the coefficient of imprimitivity is both the greatest common divisor of the lengths of the cycles in the life-cycle graph and the number of eigenvalues whose modulus is equal to the asymptotic growth rate λ. It corresponds to the period of the asymptotic oscillations of the population structure.

When A is reducible, the coefficient of imprimitivity is defined as the least common multiple of the periods of its dominant components (see the documentation of the function index_of_imprimitivity() from the module mathtools for details). It corresponds to the longest possible period of the asymptotic oscillations: if we start from any¹ population vector with individuals in every classes, there will be asymptotic oscillations whose period is the coefficient of imprimitivity; but some initial population vectors can produce asymptotic oscillations with a period that is a divisor of the index of imprimitivity.

When the index of imprimitivity is equal to 1, the population structure can never exhibit asymptotic oscillations and A is said to be aperiodic.

aperiodic

Whether the projection matrix A is aperiodic, i.e. whether the index of imprimitivity is equal to 1. Aperiodicity is a necessary condition for the convergence of the population structure: otherwise, it can exhibit asymptotic oscillations (see index_of_imprimitivity for details).

I.1.d.  Primitivity and quasi-primitivity

primitive

Whether the projection matrix A is primitive, meaning that there exists an integer k such that ${\mathbf{A}}^k$ is positive. Equivalently, a non-negative matrix is primitive if and only if it is irreducible and aperiodic.

Primitivity is a central notion of the Perron-Frobenius theory. It ensures that, as $t\to \mathrm{\infty }$,

$$\frac{{\mathbf{A}}^t}{{\lambda }^t}\to \mathbf{w}\mathbf{v},$$

with λ the dominant eigenvalue of A, v and w the corresponding left and right eigenvectors (scaled so that vw = 1), and wv the matrix $(w_i{v}_j)$. As a result, for any¹ initial population vector $\mathbf{n}(0)$, $\mathbf{n}(t)/{\lambda }^t\to c\mathbf{w}$, where $c=\mathbf{vn}(0)$. Note that these conclusions also hold if A is quasi_primitive.

See is_primitive() from the module mathtools for a function that can be used on any NumPy array and for details about the implementation.

quasi_primitive

Whether the projection matrix A is quasi-primitive, i.e. quasi-irreducible and aperiodic. Equivalently, a non-negative matrix is quasi-primitive if and only if it has a single eigenvalue of maximum modulus and this eigenvalue has algebraic multiplicity 1.

Quasi-primitivity is not a standard notion, but it is a useful relaxation of primitivity when working with matrix population models. Indeed, like primitivity it ensures that ${\mathbf{A}}^t/{\lambda }^t\to \mathbf{w}\mathbf{v}$ as $t\to \mathrm{\infty }$. The only difference with the primitive case is that the dominant eigenvectors v and w are allowed to have entries that are equal to zero (corresponding to “sink” classes for v and to “source” ones for w).

Further generalizations are possible – see e.g. [Roth81] – but they are more technical and typically not useful for matrix population models, where the structure of matrices is subject to biological constraints. See Section 9.3 of [Hogb06] for an overview of Perron-Frobenius theory for reducible matrices.

I.2.  Eigen-elements and related quantities

I.2.a.  Asymptotic dynamics

lmbd

The asymptotic growth rate λ. For large times t, the population size will grow (or decrease, if $\lambda <1$) like ${\lambda }^t$ times a constant that depends on the initial composition of the population.

Mathematically, λ is the spectral radius of the projection matrix A, i.e. the maximum of the modulus of its eigenvalues. Note that the spectral radius of a non-negative matrix is always an eigenvalue of that matrix (see e.g. Section 8.3 of [Meye00]). However, there can be several eigenvalues of maximum modulus; in fact, λ is the only eigenvalue of maximum modulus if and only if A is aperiodic – see Fact 3.(a) in Section 9.3 of [Hogb06]. When λ is the only eigenvalue of maximum modulus, we refer to it as the dominant eigenvalue.

w

The stable class distribution vector w. If the projection matrix A is primitive, then for any¹ initial population structure $\mathbf{n}(0)$ the population vector n will satisfy

$$\frac{\mathbf{n}(t)}{\Vert \mathbf{n}(t){\Vert }_1}\to \mathbf{w}\text{as}t\to \mathrm{\infty },$$

where $\Vert \mathbf{n}(t){\Vert }_1=\sum _i|n_i(t)|$ denotes the total population size. This convergence will also hold if A is quasi-primitive and we start with at least one individual in the dominant component (or one of its sources); see the documentation of quasi_primitive for details.

Mathematically, w is the right eigenvector associated to the asymptotic growth rate λ. If A is not quasi-primitive, then w is not well-defined (the asymptotic population structure may either depend on the initial one, due to reducibility; or not converge, due to periodicity). In that case, a vector of nan will be returned.

v

The reproductive values vector v, scaled so that $\mathbf{vw}=\sum _i{v}_i{w}_i=1$. The term “reproductive value” is justified by the fact that, for quasi-primitive projection matrices, as $t\to \mathrm{\infty }$ the population size is asymptotically equivalent to $c{\lambda }^t$, where

$$c=\mathbf{v}\mathbf{n}(0)=\sum _i{v}_i{n}_i(0).$$

Thus, $v_i$ is the relative contribution of an individual of class i to the composition of the population in the distant future.

Mathematically, v is the left eigenvector associated to the asymptotic growth rate λ. If A is not quasi-primitive, then v is not well-defined. In that case, a vector of nan will be returned.

Note that even though NumPy does not make a distinction between row vectors and column vectors (both correspond to 1-D arrays), in our mathematical expressions we consider v to be a row vector, This explains why we write vw for the scalar $\sum _i{v}_i{w}_i$ and wv for the matrix $(v_i{w}_j)$.

I.2.b.  Second order dynamics

damping_ratio

The damping ratio of the projection matrix A, defined as the growth rate λ divided by the maximum of the modulus of the subdominant eigenvalues, that is,

$$\frac{\lambda }{\tau }\text{where}\tau =max\{|\mu |:\mu \text{ is an eigenvalue with }|\mu |<\lambda \},$$

where $max\varnothing =0$ – i.e. when there is no eigenvalue $\mu$ such that $0<|\mu |<\lambda$, the damping ratio is infinite and inf is returned.

The damping ratio corresponds to the rate of convergence to the steady regime. Indeed, for primitive projection matrices, for any $r>\tau$,

$${\mathbf{A}}^t={\lambda }^t\mathbf{wv}+O(r^t)$$

see e.g. Theorem 1.2 in [Sene06] for a slightly more precise statement.

Note that the name “damping ratio” is now well-established in the literature on matrix population models (see e.g. [Casw00]), but that it is slightly unfortunate since $\lambda /\tau$ does not correspond to the notion of damping ratio in classical physics, which also depends on the period of the second-order oscillations (see second_order_period). The damping ratio from matrix population models is actually a geometric rate of convergence; in all likeliness the confusion comes from the interpretation of the word ratio (which, in physics, refers to the ratio of the actual damping to the critical one).

second_order_period

The period of the second-order oscillations,

$$\mathcal{P}=\left|\frac{2\pi }{\arg {\lambda }_2}\right|,$$

where ${\lambda }_2$ and its complex conjugate ${\lambda }_2^{\ast }$ are the only eigenvalues with largest modulus strictly less than λ (note that these two eigenvalues give the same $\mathcal{P}$ so it does not matter which one we choose in the definition). When ${\lambda }_2={\lambda }_2^{\ast }\ge 0$, then there are no oscillations associated to the eigenvalue of second-largest modulus and nan is returned. Note that there could still be oscillations associated to eigenvalues of smaller modulus.

In the unlikely case where there are more than two eigenvalues of second-largest modulus, then the corresponding oscillations are periodic if and only if the corresponding periods ${\mathcal{P}}_1,\dots ,{\mathcal{P}}_q$ are commensurable, i.e. if there exists positive integers $k_1,\dots ,k_q$ such that $k_1{\mathcal{P}}_1=\cdots =k_q{\mathcal{P}}_q$. The period is then the smallest possible such common value of $k_i{\mathcal{P}}_i$. Of course, it does not make sense to ask if real numbers are commensurable when they are represented by floating point numbers. Thus, in that case we test if there are small integers such that the relation above holds up to numerical precision; otherwise a CannotCompute error is raised.

The values that we compute may differ from those calculated by other software. Whenever we identified discrepancies, they were due to the fact that other software chose their arguments in $[0,2\pi ]$ or failed to take into account the possibility of more than two eigenvalues of second largest modulus, which are errors of implementation.

I.2.c.  Other eigen-elements

eigenvalues

The eigenvalues of the projection matrix A, in decreasing order of their modulus and repeated according to their algebraic multiplicities.

See right_eigenvectors and left_eigenvectors for the corresponding eigenvectors, and eigen_elements() from the module mathtools for details about the implementation.

right_eigenvectors

Right eigenvectors of the projection matrix A, that is, each ${\mathbf{w}}^{(k)}$ = right_eigenvectors[k] is a solution of

$$\mathbf{A}{\mathbf{w}}^{(k)}={\lambda }^{(k)}{\mathbf{w}}^{(k)},$$

where ${\lambda }^{(k)}$ = eigenvalues[k]. Note that there are always some small errors (typically of the order of the machine epsilon, although there are no strict guarantees on that) in the obtention of those solutions:

>>> m = MPM(A = "0.1 2; 0.4 0.8")
>>> lmbd0, w0 = m.eigenvalues[0], m.right_eigenvectors[0]
>>> m.A @ w0
array([0.85212626, 0.55834237])
>>> lmbd0 * w0
array([0.85212626, 0.55834237])
>>> m.A @ w0 - lmbd0 * w0
array([-1.11022302e-16,  0.00000000e+00])

Each vector is scaled so that

$$\Vert {\mathbf{w}}^{(k)}{\Vert }_1:=\sum _i|w_i^{(k)}|=1,$$

and is chosen to be non-negative when possible.

Importantly, note that unless A is diagonalizable, some of its eigenvectors will be linearly dependent (that includes the possibility of repeated or null eigenvectors). Thus, although right_eigenvectors will always return exactly dim eigenvectors (because the eigenvalues are repeated according to their algebraic multiplicities in eigenvalues), this should not be interpreted as meaning that A “has dim eigenvectors” (an abusive formulation that implicitly refers to linearly independent eigenvectors). If needed, functions to compute the geometric multiplicities and to test linear independence are provided in the module mathtools.

left_eigenvectors

Left eigenvectors of the projection matrix A, that is, each ${\mathbf{v}}^{(k)}$ = left_eigenvectors[k] is a solution of

$${\mathbf{v}}^{(k)}\mathbf{A}={\lambda }^{(k)}{\mathbf{v}}^{(k)},$$

where ${\lambda }^{(k)}$ = eigenvalues[k].

Note that some authors define left eigenvectors as column vectors, and will sometimes introduce them as “the right eigenvectors of the conjugate transpose matrix ${\mathbf{A}}^{\ast }$”. In the context of matrices with real entries, all definitions are equivalent up to transposition, since ${\mathbf{A}}^{\ast }={\mathbf{A}}^{\mathrm{\top }}$ and $({\mathbf{A}}^{\mathrm{\top }}\mathbf{x}{)}^{\mathrm{\top }}=\mathbf{x}\mathbf{A}$.

Although in our mathematical expressions we consider left eigenvectors as row vectors, NumPy does not make a distinction between row vectors and column vectors and all vectors in this library are represented using NumPy 1-D arrays.

The left eigenvectors returned here are all scaled so that

$$\Vert {\mathbf{v}}^{(k)}{\Vert }_1:=\sum _i|v_i^{(k)}|=1,$$

however note that a different scaling is used for the vector of reproductive values v.

Finally, the comments on the accuracy of the determination of the eigenvectors made in the documentation of right_eigenvectors also apply here.

I.2.d.  Derivatives (sensitivities and elasticities)

sensitivities

The matrix of sensitivities of the growth rate λ to the entries of the projection matrix A, that is,

$$s_{\lambda }(a_{ij})=\frac{\partial \lambda }{\partial a_{ij}}.$$

If A is primitive, then

$$s_{\lambda }(a_{ij})=v_i{w}_j,$$

where the vector of reproductive values v and the stable class distribution w are such that vw = 1. In fact, this formula is valid whenever A is quasi-irreducible; however for this v and w have to be defined as the only non-negative eigenvectors associated to λ, which cannot necessarily be interpreted as the reproductive values and stable distribution.

If A is not quasi-irreducible, then λ is not a differentiable function of its entries. In that case, a warning is issued and a matrix of nan is returned.

The standard reference for the formula $s_{\lambda }(a_{ij})=v_i{w}_j$ is [Casw78], which was the first time it appeared in the biological literature. However, this paper does not contain a mathematical proof. The first proof seems to have been published in the earlier paper [Vahr76]. As of June 2021, this paper was cited only 6 times.

See Chapter 3 of [KiNe19] for a more general yet very accessible treatment.

elasticities

The matrix of elasticities of the growth rate λ to the entries of the projection matrix A, that is,

$$e_{\lambda }(a_{ij})=\frac{a_{ij}}{\lambda }\frac{\partial \lambda }{\partial a_{ij}}=\frac{\partial \log \lambda }{\partial \log a_{ij}}.$$

Contrary to the sensitivities, which measure the absolute change in λ in response to an additive perturbation of the entries of A, the elasticities measure the relative change in λ in response to a multiplicative perturbation of the entries of A.

The elasticities also have a simple interpretation as the asymptotic frequencies of the transitions of the life cycle that we encounter as we follow the lineage of a focal individual back in time; see [BALM17] for more on this.

If A is primitive, then

$$e_{\lambda }(a_{ij})=v_i{a}_{ij}{w}_j/\lambda ,$$

where the vector of reproductive values v and the stable class distribution w are such that vw = 1. See the documentation of sensitivities for details about the validity of this formula and for references.

I.3.  The genealogical Markov chain

I.3.a  Transition matrix and stationary distribution

P

The transition of matrix of the genealogical Markov chain, $\mathbf{P}=(p_{ij})$ where

$$p_{ij}=\frac{a_{ij}{w}_j}{\lambda w_i},$$

with λ the asymptotic growth rate and w the stable distribution.

The genealogical Markov chain describes the sequence of classes that we encounter as we go “up” the genealogy of the population by following the lineage of a focal individual backwards in time: $p_{ij}$ is the probability that an individual in class i at time t was in class j at time t-1 (or had its mother in class j, if the individual was not born yet at time t-1).

The genealogical Markov chain was first introduced in [Deme74] to define population entropy, but has recently found its use in a variety of problems. See [BALM17] for a complete introduction.

See also Ps and Pf for the survival and fertility components of P, and pi for its stationary distribution.

pi

The vector of class reproductive values ${\pi }_i=v_i{w}_i$. This is also the stationary distribution of the Markov chains corresponding to P, see [BALM17].

I.3.b  Entropy rate and other statistics

entropy_rate

The entropy rate of the genealogical Markov chain P, that is

$$H=-\sum _{i,j}{\pi }_i{p}_{ij}{\log }_2{p}_{ij}.$$

This measures the information density of the trajectories of the genealogical chain, and can therefore be seen as a measure of the complexity of the life cycle (although, when comparing life cycles with different number of classes, remember that H takes values in $[0,{\log }_{2}n]$, where n is the number of classes). Mathematically, if we let $X_1,X_2,\dots$ be a trajectory of the genealogical chain and $H(X_1,\dots ,X_t)$ denote the joint entropy of $X_1,\dots ,X_t$, then

$$H=\underset{t\to \mathrm{\infty }}{lim}\frac{1}{t}H(X_1,\dots ,X_t)$$

The first use of the entropy rate of P in matrix population models is due to Demetrius in [Deme74], where it is refered to as the population entropy. The same terminology is used in [Tulj82], but Demetrius later changed it²; nowadays, $H$ is often referred to as the evolutionary entropy, the term population entropy being reserved to the quantity $S=H{T}_a$, where $T_a$ is the generation time.

hitting_times

(coming in version 0.2.0)

kemeny_constant

(coming in version 0.2.0)

II.  Descriptors of the survival and fertility matrices

II.1.  The A = S + F decomposition

II.1.a  Matrices

split

Whether the A = S + F decomposition of the projection matrix into survival probabilities vs fertilities is available. Note that many descriptors of matrix population models require this decomposition in order to be computed.

S

The survival matrix $\mathbf{S}=(s_{ij})$, where $s_{ij}$ is the probability that an individual of class j survives and goes to class i at the next time-step. See also F and split.

F

The fertility matrix $\mathbf{F}=(f_{ij})$, where $f_{ij}$ is the expected number of offspring of class i of an individual in class j. See also S and split.

mixed_transitions

Whether the model has “mixed transitions”, i.e. transitions in the life-cycle that can correspond to both survival and reproduction. In other words, whether there are overlapping entries between the matrices S and F.

fundamental_matrix

The fundamental matrix associated to S, that is,

$$\mathbf{N}=(\mathbf{I}-\mathbf{S}{)}^{-1}.$$

The (i, j)-th entry of N is the expected number of times that an individual from class j will visit the class i during its remaining lifetime.

See e.g. Chapter III of [KeSn76] for more on the mathematical properties of N.

Ps

Survival component of the genealogical matrix P: $p_s(i,j)$ is the probability that an individual in class i at time t was alive in class j at time t-1.

Pf

Fertility component of the genealogical matrix P: $p_f(i,j)$ is the probability that an individual in class i has just been born to a mother from class j.

fundamental_matrix_Ps

The fundamental matrix associated to ${\mathbf{P}}_s$, that is

$$(\mathbf{I}-{\mathbf{P}}_s{)}^{-1}.$$

The (i, j)-th entry of this matrix is the expected number of times that an individual from class i has been in class j in the past.

II.1.b  Leslie and Usher models

leslie

Whether the model is a Leslie model, i.e. an age-classified model with a projection matrix is of the form

$$\begin{array}{r}\mathbf{A}=\left(\begin{array}{ccccc}{f}_1& f_2& \cdots & f_{n-1}& f_n\\ s_1& 0& \cdots & 0& 0\\ 0& s_2& \ddots & & ⋮\\ ⋮& & \ddots & 0& 0\\ 0& 0& \cdots & s_{n-1}& s_n\end{array}\right)\end{array}$$

where $s_i$ is the probability of surviving from class i to class i + 1 and $f_i$ is the fertility of class i. The last class regroups all individuals older than a certain age.

Note that here the classes have to be labeled by increasing age, otherwise the model will not be considered a Leslie model, even if the life-cycle graph is isomorphic to a Leslie model. You can therefore rely on the fact that any model such that leslie == True is in the exact form above. To test whether a model is a “relabeled” Leslie model, see relabeled_leslie.

Also note that for a model to be considered a Leslie model, the survival and fertility matrices have to have been provided when creating the model. Otherwise there simply is no way to know whether this is really a Leslie model and this could lead to incorrect results.

>>> m0 = MPM(A = "0.5 1; 0.3 0.8")
>>> m0.leslie
Traceback (most recent call last):
  ...
matpopmod.utils.NotAvailable: A = S + F decomposition required.
>>> m1 = MPM(A = "0.5 1; 0.3 0.8", F = "0.5 1; 0 0")
>>> m1.leslie
True
>>> m1.R0
2.0
>>> m2 = MPM(A = "0.5 1; 0.3 0.8", F = "0 1; 0 0")
>>> m2.leslie
False
>>> m2.R0 # Different from the corresponding Leslie model m1. 
3.0000000000000004

usher

Whether the model is an Usher model (also sometimes called “Lefkovitch model” or “standard size-based model”), i.e. has a projection matrix of the form

$$\begin{array}{r}\mathbf{A}=\left(\begin{array}{ccccc}{f}_1+r_1& f_2& \cdots & f_{n-1}& f_n\\ g_1& r_2& \cdots & 0& 0\\ 0& g_2& \ddots & & ⋮\\ ⋮& & \ddots & r_{n-1}& 0\\ 0& 0& \cdots & g_{n-1}& r_n\end{array}\right)\end{array}$$

where $f_i$ is the fertility of class i, $g_i$ the probability of “growing” from class i to class i + 1, and $r_i$ the probability of surviving and remaining in class i. Of course the term “growing” is an image, since we have no way to know whether the classes are actually based on size.

Note that, mathematically, Usher models are a generalization of Leslie models. Thus, every model such that usher == True is also such that leslie == True.

See also relabeled_usher for whether the model can be obtained by reordering the classes of an Usher model.

relabeled_leslie

(coming in version 0.2.0)

relabeled_usher

(coming in version 0.2.0)

II.1.c  Newborn and reproductive classes

newborn_classes

The vector whose i-th entry is equal to:

• $0$ if no individual in class i is a newborn, i.e. if $f_{ij}=0$ for all j.

• $1$ if every individual in class i is a newborn, i.e. if $s_{ij}=0$ for all j (and there exists j such that $f_{ij}>0$).

• nan otherwise. In that case, class i may or may not be considered newborn class, since some but not all of its individuals are newborns.

If you want to use a different convention, you can use NumPy’s nan_to_num() function:

>>> m = MPM(S = "0.1 0; 0.4 0.8", F = "0 2; 0 0")
>>> m.newborn_classes
array([nan,  0.])
>>> numpy.nan_to_num(m.newborn_classes, nan = 0)
array([0., 0.])
>>> numpy.nan_to_num(m.newborn_classes, nan = 1)
array([1., 0.])

(this requires NumPy ≥1.17)

See also proportion_newborns for the proportion of newborns in each stage, and nu and class_of_birth for related quantities.

unique_newborn_class

Whether the model has a single, well-defined newborn class, i.e. whether the two following conditions are met: (1) all individuals are born in the same class and (2) all individuals in that class are newborns.

reproductive_classes

The vector whose i-th entry is equal to $1$ if i is a reproductive class and $0$ otherwise.

postreproductive_classes

The vector whose i-th entry is equal to $1$ if i is a post-reproductive class and $0$ otherwise.

Note that if i is post-reproductive then its reproductive value $v_i$ is equal to zero, but that the converse is not true: in a meta-population model, the classes of a sink patch can have $v_i=0$ without being post-reproductive.

proba_maturity

(coming in version 0.2.0)

II.1.d  Steady state distributions

proportion_newborns

The vector whose i-th entry is the proportion of the individuals of class i that have just been born, assuming the population is at its stable class distribution. That is,

$$\sum _j\frac{f_{ij}{w}_j}{\lambda w_i}.$$

If the proportion of newborns in the class is not defined (e.g, because it depends on the composition of the population and there is no stable class distribution), the corresponding entry is set to nan.

nu

The vector $\mathit{\nu }=({\nu }_i)$ giving the fraction of offspring produced in a given year that are born in class i, assuming the population is at its stable class distribution. That is,

$$\mathit{\nu }=\frac{\mathbf{Fw}}{\Vert \mathbf{Fw}{\Vert }_1},$$

where $\Vert \cdot {\Vert }_1$ denotes the $L_1$-norm (sum of the absolute values of the components). If there is no stable class distribution but all individuals are born in class k, then the vector $\mathit{\nu }=({\nu }_i)$, where ${\nu }_i=1$ if and only if i = k, is returned.

class_of_birth

The matrix $\mathbf{B}=(b_{ij})$ whose (i, j)-th entry gives the probability that an individual currently in class j was born in class i, that is,

$$\begin{array}{rl}{b}_{ij}& =\mathbb{P}(\text{born in }i\mid \text{in }j\text{ at }t=0)\\ & =\sum _{t\ge 0}\mathbb{P}(\text{born at }-t|\text{in }i\text{ at }-t)\mathbb{P}(\text{in }i\text{ at }-t|\text{in }j\text{ at }t=0)\\ & =\sum _{t\ge 0}\sum _k{p}_f(i,k)p_s^{(t)}(j,i),\end{array}$$

with ${\mathbf{P}}_f=(p_f(i,j))$ and ${\mathbf{P}}_s^t=(p_s^{(t)}(i,j))$ the genealogical matrices. In matrix form,

$$\mathbf{B}=(\mathbf{I}-{\mathbf{P}}_s^{\mathrm{\top }}{)}^{-1}\otimes {\mathbf{P}}_f\mathbf{E}$$

where E is a square matrix of ones and $\otimes$ denotes the entry-wise product.

class_of_death

(coming in version 0.2.0)

II.1.d  Fertility and survival excesses

fertility_excess

The fertility excess (non-standard terminology), defined as the constant c > 0 such that $\rho (\mathbf{S}+c^{-1}\mathbf{F})=1$, where $\rho$ denotes the spectral radius (i.e. the asymptotic growth rate λ, in the case of a primitive matrix). In other words, the fertility excess is the factor by which all fertilities must be rescaled to get a stationary projection matrix.

survival_excess

The survival excess (non-standard terminology), defined as the constant c > 0 such that $\rho (c^{-1}\mathbf{S}+\mathbf{F})=1$, where $\rho$ denotes the spectral radius (i.e. the asymptotic growth rate λ, in the case of a primitive matrix). In other words, the survival excess, if it exists, is the factor by which all fertilities must be rescaled to get a stationary projection matrix.

When λ > 1, if the spectral radius of F is greater than 1, then c is not defined and NaN will be returned. When λ < 1, the matrix $c^{-1}\mathbf{S}$ can have columns that sum to more than 1. This means that the original model is such that it is impossible to get a sationary projection matrix by increasing survival by a common factor, even though c is mathematically well-defined.

theta

(coming in version 0.2.0)

II.2.  Biological descriptors

II.2.a  Reproductive output

G

The next generation matrix

$$\mathbf{G}=\mathbf{F}(\mathbf{I}-\mathbf{S}{)}^{-1},$$

whose (i, j)-th entry is the expected number of offspring of class i that an individual of class j will produce in its remaining lifetime. See R0 for the dominant eigenvalue of this matrix, a commonly used measure of the net reproductive rate.

w_G

Cohort stable distribution (dominant right-eigenvector of G, normalized so that its entries sum to 1).

v_G

Cohort reproductive values (dominant left-eigenvector of G, normalized so that $v_c\cdot w_c=1$).

total_reproductive_output

The matrix $\mathbf{R}=(r_{ij})$ such that $r_{ij}$ is the expected number of offspring of class i that an individual currently in class j is expected to produce over the course of its life. This includes offspring produced before the present:

$$\begin{array}{rl}\mathbf{R}& =\underset{\text{before}}{\underbrace{\sum _{t\ge 1}\mathbf{F}({\mathbf{P}}_s^{\mathrm{\top }}{)}^t}}+\underset{\text{after}}{\underbrace{\sum _{t\ge 0}\mathbf{F}{\mathbf{S}}^t}}\\ & =\mathbf{F}(\mathbf{I}-{\mathbf{P}}_s^{\mathrm{\top }}{)}^{-1}{\mathbf{P}}_s^{\mathrm{\top }}+\mathbf{G}.\end{array}$$

where ${\mathbf{P}}_s$ is the survival component of the genealogical matrix P.

proba_repro

The vector whose i-th entry is the probability that an individual currently in class i will reproduce during its remaining lifetime. Use nu @ proba_repro to get the proportion of newborns that reproduce before dying in the stable population.

See life_expectancy_repro for details on how this is calculated.

II.2.b  Net reproductive rate

R0

The measure of the net reproductive rate (also called net reproduction number) introduced by Cushing and Zhou in [CuZh94], namely the dominant eigenvalue of the next generation matrix $\mathbf{G}=\mathbf{F}(\mathbf{I}-\mathbf{S}{)}^{-1}$.

Although this is a relevant measure of the net reproductive rate, one should note that if there are several newborn classes then this does not necessarily correspond to the expected number of offspring that newborn will produce over the course of its life, in the usual sense. This is because there are two, non-equivalent ways to choose a newborn:

• According to the dominant eigenvector of G, as done here.

• Uniformly among the individuals born in a given year, as done in cohort_R0.

These two measures can differ significantly:

>>> mpm.examples.thalia_democratica.R0
4.001730883164397
>>> mpm.examples.thalia_democratica.cohort_R0
5.222416673531082

cohort_R0

The net reproductive rate, as the average number of offspring that a newborn is expected to produce over its lifetime in the stable population. This is given by

$$R_0^{\ast }=\Vert \mathbf{G}\mathit{\nu }{\Vert }_1$$

with G the next-generation matrix and $\mathit{\nu }=({\nu }_i)$, where ${\nu }_i$ is the fraction of newborns that are born in class i when the population is at its stable distribution (see nu).

II.2.c  Generation time

T_a

The generation time, as the mean age of mothers in the stable population (that is, when there is a birth in the stable population we record the age of the mother).

As suggested by [CoEl92], all births are not counted equally: instead, they are weighted by the reproductive value of the newborns. This yields a mathematically more tractable measure. In particular, [BiLe15] showed that $T_a$ can be interpreted as the expected time between two reproductive events along the ancestral lineage of an individual, and is given by

$$T_a=\frac{\lambda }{\mathbf{vFw}}.$$

See [Elln18] for a detailed discussion.

T_R0

The $R_0$ generation time, i.e. the time it takes for the population to grow by a factor of its net reproductive rate $R_0$:

$$T_{R_0}=\frac{\log R_0}{\log \lambda }.$$

When $\lambda \to 1$, $T_{R_0}\to T_a$, the average age of mothers in the stable population; see [Elln18]. Thus we can extend $T_{R_0}$ by continuity at $\lambda =1$.

For coherence with the literature and other libraries, the dominant eigenvalue of the next generation matrix G is used to compute the net reproductive rate $R_0$. However, the interpretation of this quantity is not as straightforward as is usually assumed, see R0 for details.

T_G

A variant of the classic “cohort generation time” mu1. Here births are weighted by the cohort reproductive values v_G instead of equally, as suggested by [StTC14]. This yields the simple formula

$$T_G={\mathbf{v}}_G(\mathbf{I}-\mathbf{S}{)}^{-1}{\mathbf{w}}_G,$$

as shown by [Elln18]. However the interpretation of this measure is subject to the same difficulties as that of mu1, and this should not be used as a measure of the typical age at which an individual is expected to reproduce.

mu1

The mean age of the mothers when considering all offspring produced by a cohort. This is sometimes referred to as the cohort generation time and is given by the formula

$${\mu }_1=\frac{\mathbf{e}\mathbf{F}(\mathbf{I}-\mathbf{S}{)}^{-2}\mathbf{F}\mathbf{w}}{\mathbf{e}\mathbf{F}(\mathbf{I}-\mathbf{S}{)}^{-1}\mathbf{F}\mathbf{w}},$$

where e is a vector of ones. See [CoEl92] and [Elln18].

Despite the name “mean age of mothers in a cohort”, this is not a good measure of the typical age at which an individual is expected to reproduce. In particular, ${\mu }_1$ can be greater than the life expectancy conditional on reproduction:

>>> astrocaryum = mpm.examples.astrocaryum_mexicanum
>>> astrocaryum.mu1
275.1594139265308
>>> astrocaryum.life_expectancy_repro
232.18934046629715

See [Bien19] for a detailed explanation of the problems associated with ${\mu }_1$, and mean_age_repro() for an alternative measure of the typical age at reproduction.

mean_age_repro(n=1000, reproduction='poisson', ini=None, target_err=0.01, return_err=False)

Note: unlike most descriptors, this function will return a different value every time it is called. This is because the numerical estimation is based on individual-based simulations, not an analytical formula. The last computed value is stored in self._LAST_COMPUTED_MEAN_AGE_REPRO.

Be aware that the computing time can be very long, especially when few newborns get to reproduce.

Numerical estimate of the mean age at which a typical mother produces offspring, using individual-based simulations. Formally, this is the expected age at offspring production, conditional on reproduction; see [Bien19].

n	The number of replicates.
reproduction	The distribution of the number of offspring. Use "bernoulli" for monotocous species and "poisson" (the default) otherwise.
ini	The birth class of the focal individual, as a vector of weights. The default is to use nu, the stable distribution of newborns.
target_err	if None, this parameter will be ignored and n replicates will be used. Otherwise, the number of replicates will start at n and double that number until err < target_err * mean (see below).
return_err	Whether to return the number of replicates used and $\mathsf{err}={\mathrm{t}}_{n-1}\hat{\sigma }/\sqrt{n}$, where ${\mathrm{t}}_{n-1}$ is the 0.975 quantile of the Student distribution with n-1 degrees of freedom. Thus, err is the half-width of the 95% confidence interval for the mean of normally distributed variables. If return_err is set, then a tuple (mean, reps, err) is returned; otherwise only mean is returned.

II.2.d  Life expectancy

Important remark: an unavoidable problem when using matrix population models to compute age-specific quantities is that the exactitude of the formulas depend on the specific assumptions made when building the model.

For instance, in pre-breeding models, “newborns” are $1-\epsilon$ years old, with $\epsilon$ small, whereas in post-breeding models they are $\epsilon$ years old. Similarly, if an individual is alive at time t and not alive at time t + 1, then we know that it died between times t and t + 1; but we have no way of knowing when.

As a result, when interpreting quantities such as the life expectancy or the mean age in a class, one should remember that “age” is counted in number of observations of the individuals, rather than absolute time. One may also want to add/subtract a correction term, e.g, based on knowledge on when deaths are most likely to occur in the projection interval.

survivorship(t)

The survivorship function $\ell$, which gives the probability that a newborn reaches age t in the stable population:

$$\ell (t)={\mathbf{e}\mathbf{S}}^t\mathit{\nu },$$

where S is the survival matrix; e is a vector of ones; and $\mathit{\nu }$ is the distribution of newborns in the stable population (see nu).

class_survivorship(t)

The class-structured survivorship, which gives the probability that an individual survives at least t time steps, as a function of the current class of that individual. That is,

$${\ell }_i(t)=\mathbb{P}(T_i\ge t),$$

where $T_i$ is the remaining lifespan of an individual from class i. For instance, ${\ell }_i(1)$ is the probability that the individual survives to the next time step. We have:

$$\ell (t)={\mathbf{e}\mathbf{S}}^t,$$

where S is the survival matrix and e is a vector of ones.

life_expectancy

The expected number of projection intervals that an individual born in the stable population lives. The convention used here is to return 1 for annuals. For instance,

>>> m = MPM(S = "0", F = "1.2")
>>> m.life_expectancy
1.0

Note that whether this corresponds to the expected age at death depends on the type of census (pre-breeding, post-breeding, birth-flow…), as well as on when deaths occur in the projection interval. It is also possible to get different results for models that describe the same population, if the models are build differently:

>>> eg = matpopmod.examples
>>> eg.passerine_prebreeding.life_expectancy
1.7
>>> eg.passerine_postbreeding.life_expectancy
1.3399999999999999

This is normal: individuals who die shortly after their birth are not taken into account in the pre-breeding model, which biases the life expectancy upwards.

With our convention, the life expectancy is given by

$$\sum _{t\ge 0}\ell (t)=\mathbf{e}(\mathbf{I}-\mathbf{S}{)}^{-1}\mathit{\nu },$$

where $\ell$ is the survivorship, e is a vector of ones and $\mathit{\nu }$ is the distribution of newborns in the stable population (see nu).

remaining_life_expectancy

The remaining life expectancy of an individual, as a function of its current class. See life_expectancy for a discussion of what this quantity represents. It is given by

$$\sum _{t\ge 0}{\ell }_i(t)=\mathbf{e}(\mathbf{I}-\mathbf{S}{)}^{-1},$$

where ${\ell }_i$ is the survivorship of class i and e is a vector of ones.

conditional_life_expectancy

(coming in version 0.2.0)

life_expectancy_repro

The life expectancy conditional on reproducing, assuming that the fertilities are Poisson-distributed. See life_expectancy for a discussion of what we mean by “life expectancy”.

To compute this quantity, we build a Markov chain tracking the movement of individuals in the life-cycle, with two absorbing states: corresponding to dead individuals that reproduced and dead individuals that did not reproduce. This requires knowing the probability of reproducing in each class of the model, which explains that we need to assume that the fertilities are Poisson-distributed.

Once we have this extended Markov chain, we condition it on absorption in the “died after reproducing” state (see e.g. Section 11.1.2 of [KeCa05], or Chapter III of [KeSn76] for a more detailed treatment) and compute the expected time to absorption when starting from the distribution of newborns in the stable population (see nu).

remaining_life_expectancy_repro

The remaining life expectancy of an individual, as a function of its current class and conditional on this individual reproducing in the future. See life_expectancy_repro for details on how this is calculated.

variance_lifespan

(coming in version 0.2.0)

variance_remaining_lifespan

(coming in version 0.2.0)

II.2.e  Age structure

mean_age_class

The vector whose i-th component is the mean age in class i, when the population is at its stable structure. It is given by

$$\mathbf{y}=((\mathbf{I}-{\lambda }^{-1}\mathbf{S}{)}^{-2}\mathit{\nu })\oslash ((\mathbf{I}-{\lambda }^{-1}\mathbf{S}{)}^{-1}\mathit{\nu }),$$

where $\mathit{\nu }$ is the distribution of newborns in the stable population and $\oslash$ denotes the entry-wise division – see e.g. Equation (23) in [CoEl92].

For post-breeding models, one may want to substract 1 from this quantity; here individuals in newborn classes are 1 year old.

mean_age_population

The mean age in the stable population, given by $\mathbf{y}\cdot \mathbf{w}$, where y is the vector whose i-th entry is the mean age in class i (see mean_age_class).

For post-breeding models, one may want to substract 1 from this quantity.

variance_age_class

(coming in version 0.2.0)

variance_age_population

(coming in version 0.2.0)

age_specific_fertility(t)

The matrix $\mathbf{\Phi }(t)$, whose (i, j)-th entry gives the expected number of class-i offspring produced at age t by an individual in class j at age 0. It is given by

$$\mathbf{\Phi }\mathbf{(}\mathbf{t}\mathbf{)}=\mathbf{F}\tilde{\mathbf{S}}(t),$$

where the (i, j)-th entry ${\tilde{s}}_{ij}(t)$ of $\tilde{\mathbf{S}}(t)=({\tilde{s}}_{ij}(t))$ is the probability that an individual in class j at time 0 is in class i at time t, conditional on having survived to that time:

$${\tilde{s}}_{ij}(t)=\frac{{\mathbf{S}}^t(i,j)}{\sum _k{\mathbf{S}}^t(k,j)}.$$

See e.g. Section 11.3.2 of [KeCa05].

mean_age_maturity

(coming in version 0.2.0)

mean_age_first_repro

(coming in version 0.2.0)

II.2.f  Entropies

See entropy_rate above for the entropy rate H of the genealogical Markov chain.

population_entropy

The population entropy, defined by Demetrius as²

$$S=H\cdot T_a$$

where H is the entropy_rate of the model and $T_a$ is the generation time, as given by the mean age of mothers in the stable population.

For Leslie models, S has several interpretations; it corresponds to both:

• The Shannon entropy of the difference of age between mothers and daughters in the stable population.

• The Shannon entropy of a random trajectory in the genealogy of the population, which we refer to as the genealogical entropy.

Moreover, for those models S has the following simple expression:

$$S=-\sum _i{\varphi }_i{\lambda }^{-i}{\log }_2{\varphi }_i{\lambda }^{-i},$$

where ${\varphi }_i$ is the age-specific net fertility and $\lambda$ is the asymptotic growth rate – see [Deme74], but note that in this article Demetrius refers to H, not S, as the entropy of the population.

For general projection matrices, there is no clear interpretation of S as an entropy in the usual sense.

birth_entropy

(coming in version 0.2.0)

genealogical_entropy

(coming in version 0.2.0)

lifetable_entropy

The life-table entropy, also known as Keyfitz’s entropy to distinguish it from the population entropy introduced by Demetrius. Introduced by Keyfitz in [Keyf77], is defined as

$$-\frac{1}{L}\sum _{t\ge 0}\ell (t){\log }_2\ell (t),$$

where $\ell$ is the survivorship, that is, the probability of reaching age t and $L=\sum _{t\ge 0}\ell (t)$ is the life expectancy.

Note that despite its accepted name, the life-table entropy is not an entropy in the usual sense. This makes its interpretation challenging – all the more so given that in the literature it is often conflated with Demetrius’ population entropy (see e.g. Section 4.3.1 of [KeCa05]).

III.  Functions of the initial population

III.1  Distance to the stable population

keyfitz_delta(n0)

Keyfitz’s $\mathrm{\Delta }$, which corresponds to the total variation distance between the stable class distribution w and the initial class distribution $\mathbf{x}=\mathbf{n}(0)/\Vert \mathbf{n}(0){\Vert }_1$. It is given by

$$\mathrm{\Delta }=\frac{1}{2}\Vert \mathbf{x}-\mathbf{w}{\Vert }_1=\frac{1}{2}\sum _i|x_i-w_i|.$$

The factor 1/2 ensures that $\mathrm{\Delta }\in [0,1]$.

cohen_D1(n0)

Cohen’s $D_1$ index of cumulative distance, introduced in [Cohe79]:

$$D_1=\Vert \sum _{t\ge 0}(\mathbf{n}(t)/{\lambda }^t-\mathbf{vn}(0)\mathbf{w}){\Vert }_1.$$

This corresponds to accumulating the difference between $\mathbf{n}(t)/{\lambda }^t$ and its limit along the trajectory of the population, then taking sum of the absolute values of the entries of the resulting vector.

What makes $D_1$ noteworth is that when the projection matrix A is (quasi-)primitive we have

$$D_1=((\mathbf{I}+\mathbf{wv}-{\lambda }^{-1}\mathbf{A})-\mathbf{wv})\mathbf{n}(0),$$

which makes it straightforward to compute $D_1$ in practice. However, one of the major drawbacks of $D_1$ is that when $n_i(t)/{\lambda }^t$ fluctuates around $\mathbf{vn}(0)w_i$, the terms of the time series tend to cancel out. For instance, consider the following situation:

>>> dp = mpm.examples.dipsacus_sylvestris
>>> dp.w
array([0.63771991, 0.26395418, 0.01215485,
       0.06931315, 0.01224112, 0.00461679])
>>> n0 = [0.2635, 0.6420, 0.0157, 0.0483, 0.0013, 0.0292]
>>> dp.cohen_D1(n0)
0.037129817800633416

Here, $D_1$ is small even though the initial population structure is very different from the stable one, and the fluctuations of $\mathbf{n}(t)/{\lambda }^t$ around $\mathbf{vn}(0)\mathbf{w}$ are initially much larger than $D_1$:

>>> traj = dp.trajectory(n0, t_max=5)
>>> numpy.sum(traj.rescaled_Y - (dp.v @ n0) * dp.w, axis=1)
array([-1.43193488,  2.16076214, -0.48105167,
       -0.9132861 ,  1.02463032, -0.31244266])

This makes the relevance of $D_1$ questionable.

population_momentum(n0, auto=True)

Population momentum is the name given by Keyfitz in [Keyf71] to the following phenomenon: when a population undergoes a demographic transition, it can keep growing even after its asymptotic growth rate has become equal to 1. The reason for this is that the population is not at its stable structure, and that the increase / decrease that we see corresponds to the transient regime.

The standard way to quantify this is to use the following quantity, often referred to as the momentum of the population:

$$M=\underset{t\to +\mathrm{\infty }}{lim}\frac{\Vert \mathbf{n}(t){\Vert }_1}{\Vert \mathbf{n}(0){\Vert }_1},$$

where t = 0 is the time at which the asymptotic growth rate of the population becomes equal to 1. For (quasi-)primitive, stationary projection matrices, this limit always exists and is given by

$$M=\frac{\mathbf{vn}(0)}{\Vert \mathbf{n}(0){\Vert }_1},$$

where v is the vector of reproductive values of the projection matrix, rescaled so that vw = 1.

Computing the population momentum requires knowing two things: (1) the population vector at the time t = 0 of the transitions and (2) the projection matrix after the transition; specifying the projection matrix before the transition is not sufficient. Therefore, this function makes the following assumptions:

• If λ = 1, we assume that the projection matrix of the model describes the dynamics of the population after the demographic transition, and we return $\mathbf{vn}(0)/\Vert \mathbf{n}(0){\Vert }_1$.

• If λ ≠ 1, we do not know what projection matrix should be used after the transition, so NotAvailable is raised. However, if auto is True, then instead of raising an error we assume that the projection matrix after the transition is

  $${\mathbf{A}}^{\prime }=\mathbf{S}+c^{-1}\mathbf{F}$$

  where c is the fertility_excess, i.e. the constant such that $\lambda (\mathbf{S}+c^{-1}\mathbf{F})=1$. This corresponds to Keyfitz’s original idea of assuming that all fertilities get reduced by the same factor during the demographic transition.

III.2  Trajectories

trajectory(n0, t_max, step=1)

The trajectory of the population vector associated to the initial condition n0, up to time t = t_max. This is the solution of

$$\mathbf{n}(t+1)=\mathbf{A}\mathbf{n}(t),t=0,\dots ,t_{max}.$$

The population vector will only be computed for multiples of the optional argument step. For instance, use step = 3 to get the population vector at times $t=0,3,6,\dots$

Returns a Trajectory object. Trajectories have a set of methods that facilitate their study and can be plotted using the function trajectory() from the module plot.

stochastic_trajectories(n0, t_max, reps, reproduction='poisson', Z=None)

Returns a list of reps independent stochastic trajectories corresponding to the matrix population model. Each trajectory starts at time t = 0 from the integer-valued population vector n0, and ends at time t = t_max.

In ecological terms, those stochastic trajectories account for demographic stochasticity – that is, the randomness stemming from the variation in the realizations of the demographic rates of the individuals (as opposed to randomness in the environment affecting the demographic rates population-wide).

Formally, the trajectories are realizations of a multitype Galton–Watson process whose types correspond to the classes of the model. At each time-step, each class-j individual x:

• produces a random number $F_{ij}(x)$ of class-i individuals, where $\mathbb{E}(F_{ij}(x))=f_{ij}$, the (i, j)-th entry of the fertility matrix F;

• survives and goes to class i with probability $s_{ij}$, the (i, j)-th entry of the survival matrix S.

What happens to individuals is independent of what happens to other individuals or of anything that happened in the past. This crucial assumption is the cornerstone of the link between multitype Galton–Watson processes and matrix population models. It ensures that $\mathbf{n}(t)$, the expected value of the (integer-valued) vector tracking the number of individuals in each class of the stochastic process, satisfies the linear equation

$$\mathbf{n}(t+1)=\mathbf{A}\mathbf{n}(t)$$

with A = S + F is the projection matrix of the matrix population model.

Each Galton–Watson process corresponds to a single matrix population model, but there are an infinity of Galton–Watson processes for a given matrix population model. To fully specify a Galton–Watson process, one needs to provide the joint law of the family $(S_{ij},F_{ij}{)}_i$, for every j. Two standard sets of assumptions can be selected via the argument reproduction:

• "poisson" — for all j, $(F_{ij}{)}_i$ are independent Poisson variables with respective parameters $f_{ij}$, and are independent of the family $(S_{ij}{)}_i$. If the survival and fertility matrices are not available, instead of the model described above we simply take the number of class-i individuals that a class-j individual leaves at the next time-step to follow a $\mathrm{Poisson}(a_{ij})$ distribution (instead of $\mathrm{Bernoulli}(s_{ij})+\mathrm{Poisson}(f_{ij})$ in the more detailed model).

  This is the assumption that we recommend using in the absence of additional information.

• "bernoulli" — for all j, $(F_{ij}{)}_i$ are incompatible Bernoulli variables that are independent of $(S_{ij}{)}_i$. In other words, each class-j individual produces an offspring with probability $\sum _k{f}_{kj}\le 1$, and that offspring is produced in class i with probability $f_{ij}/\sum _k{f}_{kj}$. If the survival and fertility matrices are not available or if $\sum _k{f}_{kj}>1$, an error will be raised.

  This assumption is what we recommend using for monotocous species.

Alternatively, one can use a custom Galton–Watson process by specifying the argument Z. In that case, Z should be a random Python function such that the dynamics of the population vector $N(t)$ is defined by:

$$N(t+1)=Z(\mathbf{F},\mathbf{S},N(t))$$

if S and F are known, and $N(t+1)=Z(\mathbf{A},\mathrm{None},N(t))$ otherwise.

Because this function focuses on efficiency, it is not possible to specify $(S_{ij},F_{ij}{)}_i$ only; the user has to use them to work out how to sample Z. See the module ibm to generate stochastic trajectories directly from $(S_{ij},F_{ij}{)}_i$.

Footnotes

1. When we say that the population satisfies $\mathbf{n}(t)\sim \mathbf{vn}(0){\lambda }^t\mathbf{w}$ for “any” initial population vector, in fact this is true if and only if $\mathbf{vn}(0)\ne 0$. For instance, if $\mathbf{n}(0)$ were colinear with another right-eigenvector of the projection matrix – say ${\mathbf{w}}^{\prime }$ with associated eigenvalue ${\lambda }^{\prime }$ – then we would have $\mathbf{n}(t)\sim \mathbf{n}(0)({\lambda }^{\prime }{)}^t$.

   However, the set of vectors x such that vx = 0 is a vector space of dimension n - 1 in the n-dimensional space of population vectors. Thus, when choosing an initial population vector according to any “reasonable” probability distribution, the probability of getting a vector such that $\mathbf{vn}(0)=0$ is equal to 0, just as there is a probability zero of falling exactly on a given line when throwing a point in the plane. As a result, the case $\mathbf{vn}(0)=0$ can be considered a degenerate case with no biological relevance.

2. The terminology regarding the quantities H and S has not been very consistent; in fact Demetrius himself has used both population entropy and evolutionary entropy to refer to either of the two. The convention we use here corresponds to a distinction explicitly made by Demetrius in a correspondence with FB, where he referred to H as the evolutionary entropy and to S as the population entropy. However, we favor the term entropy rate over evolutionary entropy because it is more explicit.

   While the definition and interpretation of S are very clear and natural in the context of Leslie models, this is not the case for general matrix population models. Based on the same correspondence with Demetrius, we adopt the point of view that S should be defined so as to ensure that $S=H/T_a$. This, however, means that we have to renounce to the probabilistic interpretations that S has in the context of Leslie models.