Galton-Watson process and bayesian inference: A turnkey method for the viability study of small populations

TL;DR

This paper introduces a Bayesian inference method combined with the Galton-Watson process to accurately predict extinction times and assess viability in small populations, enhancing conservation strategies.

Contribution

It presents a novel, non-computational Bayesian approach integrating the Galton-Watson model for small population viability analysis, surpassing classical models in flexibility and information gain.

Findings

Effective prediction of extinction times in small populations

Compatibility with classical matrix population models

Application to real-life conservation scenarios

Abstract

1 Sharp prediction of extinction times is needed in biodiversity monitoring and conservation management. 2 The Galton-Watson process is a classical stochastic model for describing population dynamics. Its evolution is like the matrix population model where offspring numbers are random. Extinction probability, extinction time, abundance are well known and given by explicit formulas. In contrast with the deterministic model, it can be applied to small populations. 3 Parameters of this model can be estimated through the Bayesian inference framework. This enables to consider non-arbitrary scenarios. 4 We show how coupling Bayesian inference with the Galton-Watson model provides several features: i) a flexible modelling approach with easily understandable parameters ii) compatibility with the classical matrix population model (Leslie type model) iii) A non-computational approach which then leads to more information with less computing iv) a non-arbitrary choice for scenarios, parameters... It can be seen to go one step further than the classical matrix population model for the viability problem. 5 To illustrate these features, we provide analysis details for two examples whose one of which is a real life example.

Figures

Tables

Table 1. Table 1. Synthetic data of Section  5 , N ​ ( t ) 𝑁 𝑡 N(t) represents the number of individuals at time t 𝑡 t and 𝐧 = ( n 1 , 1 ​ ( k , t ) ) 0 ≤ k ≤ 4 , 0 ≤ t ≤ 4 𝐧 subscript subscript 𝑛 1 1 𝑘 𝑡 formulae-sequence 0 𝑘 4 0 𝑡 4 \mathbf{n}=(n_{1,1}(k,t))_{0\leq k\leq 4,0\leq t\leq 4} the learning data.

Times $t$	0	1	2	3	4	5
$N$	100	75	59	43	33	22
$n_{1,1}(0,t)$	47	30	29	20	18
$n_{1,1}(1,t)$	39	37	23	17	11
$n_{1,1}(2,t)$	8	4	2	3	2
$n_{1,1}(3,t)$	4	2	4	2	1
$n_{1,1}(4,t)$	2	2	1	1	1

Table 2. Table 2. Posterior for all parameters p i , j ​ ( k ) subscript 𝑝 𝑖 𝑗 𝑘 p_{i,j}(k) and the vector valued p 4 , 0 subscript 𝑝 4 0 p_{4,0} for the bears population example. We use the data [ 10 ] and the non-informative prior.

Parameter	Posterior	Mean	IC 90%
$p_{0,1}\left(1\right)$	Beta (18,4)	0.8181	[0.67,0.93]
$p_{1,2}\left(1\right)$	Beta(17,1)	0.9444	[0.84,1]
$p_{2,3}\left(1\right)$	Beta(13,3)	0.8125	[0.64,0.94]
$p_{3,4}\left(1\right)$	Beta(13,1)	0.9286	[0.79,1]
$p_{4,4}\left(1\right)$	Beta(73,4)	0.9480	[0.90,0.98]
$p_{4,0}$	Dir(71,9,7,1)
$p_{4,0}\left(0\right)$	Beta(71,17)	0.8068	[0.73,0.87]
$p_{4,0}\left(1\right)$	Beta(9,79)	0.1023	[0.05,0.16]
$p_{4,0}\left(2\right)$	Beta(7,81)	0.0795	[0.04,0.13]
$p_{4,0}\left(3\right)$	Beta(1,87)	0.0114	[0,0.03]

Equations

$$\mathbb{E}[N(t)|\mathbf{n}]=\mathbb{E}[M^t\cdot N(0)|\mathbf{n}].$$

$$\mathbb{P}(\lambda >1|\mathbf{n}).$$

$$U(t):t\mapsto \mathbb{E}[{\lambda }^t|\mathbf{n}]\times \frac{{\sum }_{j=1}^K{N}_j(0)v_j}{{\min }_{1\le j\le K}{v}_j}.$$

$$\mathbb{P}(\text{Ext})=\mathbb{E}\left[s_1^{N_1(0)}\times \cdots \times s_K^{N_K(0)}|\mathbf{n}\right]$$

$$N(t+1)=A\cdot N(t),$$

$$N_j(t+1)=\sum _{i=1}^K{A}_{i,j}{N}_i(t).$$

$$N_1(t+1)=A_{1,1}{N}_1(t)\Rightarrow N(t)=A_{1,1}^t{N}_1(0),$$

$$\lim_{t\to\infty}N_{1}(t)=\left\{\begin{array}[]{ll}0&\text{if }A_{1,1}<1\\ N_{1}(0)&\mbox{if }A_{1,1}=1\\ +\infty&\mbox{if }A_{1,1}>1.\end{array}\right.$$

$$T_{\gamma }=\frac{\log (\gamma )-\log (N_1(0))}{\log (A_{1,1})},$$

$$N_j(t+1)=\sum _{i=1}^K\sum _{l=1}^{N_i(t)}{\xi }_{i,j,l,t}.$$

$$A_{i,j}(t)=\frac{1}{N_i(t)}\sum _{l=1}^{N_i(t)}{\xi }_{i,j,l,t},$$

$$N(t+1)=A(t)\cdot N(t).$$

$$p_{i,j}\sim \mathbf{Dir}\left({\alpha }_{i,j}(0),...,{\alpha }_{i,j}({\kappa }_{i,j})\right),$$

$$\mathbb{P}$$

$$=\prod _{l,t\in \mathbb{N}}\prod _{i,j=1}^K\sum _{k\in A_{i,j,l,t}}{p}_{i,j}(k),$$

$$M_{i,j}=\sum _{k=0}^{{\kappa }_{i,j}}k{p}_{i,j}(k).$$

$$uM=\lambda M,Mv=\lambda v$$

$$\sum _{i=1}^K{u}_i{v}_i=1,\sum _{i=1}^K{u}_i=1$$

$$\mathbb{P}(T_{\text{Ext}}>t|N_1(0),\dots ,N_K(0),(p_{i,j)})$$

$$\displaystyle=\$$

$$\displaystyle\leq\$$

$$\mathbb{P}$$

$$=s_1^{N_1(0)}\times \cdots \times s_K^{N_K(0)},$$

$${\phi }_i(s_1,...,s_K)=\sum _{i_1,...,i_k}{s}_1^{i_1}{p}_{i,1}^{i_1}\dots s_K^{i_K}{p}_{i,K}^{i_K}.$$

$$\mathbb{P}$$

$$=s_i.$$

$$\mathcal{L}(p_{i,j}|\mathbf{n})$$

$$=$$

$$\dots ,{\alpha }_{i,j}({\kappa }_{i,j})+\sum _{t=0}^T{n}_{i,j}({\kappa }_{i,j},t)).$$

$$\mathbb{P}(Ext|(n)),$$

$$\mathbb{E}[M_{i,j}|\mathbf{n}]=\frac{{\sum }_{k=1}^{{\kappa }_{i,j}}k\left({\alpha }_{i,j}(0)+{\sum }_{t=0}^T{n}_{i,j}(k,t)\right)}{{\sum }_{k=0}^{{\kappa }_{i,j}}\left({\alpha }_{i,j}(0)+{\sum }_{t=0}^T{n}_{i,j}(k,t)\right)}.$$

$$p_R=(0.8,0.1,0.05,0.05).$$

$$p_{1,1}^{\star }(0)=(1-p_S)p_R(0)=0.48$$

$$p_{1,1}^{\star }(1)=p_S{p}_R(0)+(1-p_S)p_R(1)=0.38$$

Full text

2 The Galton-Watson process is a classical stochastic model for describing population dynamics. Its evolution is like the matrix population model where offspring numbers are random. Extinction probability, extinction time, abundance are well known and given by explicit formulas. In contrast with the deterministic model, it can be applied to small populations.

3 Parameters of this model can be estimated through the Bayesian inference framework. This enables to consider non-arbitrary scenarios.

4 We show how coupling Bayesian inference with the Galton-Watson model provides several features: i) a flexible modelling approach with easily understandable parameters ii) compatibility with the classical matrix population model (Leslie type model) iii) A non-computational approach which then leads to more information with less computing iv) a non-arbitrary choice for scenarios, parameters… It can be seen to go one step further than the classical matrix population model for the viability problem.

5 To illustrate these features, we provide analysis details for two examples whose one of which is a real life example.

1. Introduction

Population viability analysis (PVA) [7, 45, 8, 40, 38] aims to manage or predict population development from some ecological, genetic or demographic assumptions. It generally relies on modelling purposes. The simulations are frequently used to work on these questions even if some works highlight the use of Mathematical results for these problems as [13, 25, 30, 39, 19] among many others.

Habitat degradation is pointed out as the major cause of extinction for many threatened species around the world. Yet, for species that are exploited for food or other purposes (e.g., exploited fish populations, bushmeat etc…), or purposely destroyed (e. g., large predators), extinction often occurs in intact habitats. For these cases, the assumption that population dynamics takes place in a fixed habitat seems reasonable. Then, assuming a given environment and a short period of time without brutal variation of the environment, are we able to propose a model to support decision making for management?

In this paper, we will illustrate how to answer some questions emanating from conservation biology such as:

• (I)

What are the odds on saving a population? How probable is extinction? For how long can the population survive?

• (II)

How sharply are numbers expected to fall in the short term?

• (III)

How can we proceed to avoid extinction [14, 15], to maintain a population size [43] or wipe out some population [6, 33, 23]. For instance, if we are able to add some individuals in a population, how many re-introductions are needed to reach a threshold of variability for the population with high probability?

These simple questions, when the environment is assumed unchanged/stable are not currently being resolved successfully and we propose a reformulation of classical models that leads to understand the population viability with or without possible actions.

To illustrate this, Cairns, Ross and Taimre in [9] underlined that the usual class of diffusion models are simple but lead to inaccurate predictions of critical values such as the expected time to extinction. They advocated that a more appropriate model may be a discrete-state Markov process describing the actual number of individuals in the population like the birth–death processes. These models are not used in practice because they are difficult to work with, from both analytical and computational points of view.

Starting from the famous Matrix population [4, 32, 31], we will introduce the simplest birth-death process of Galton Watson which leads asymptotically to exponential growth and does not take into account some density-dependence dynamics. This assumption is relevant if we address population dynamics in populations at low density, whose dynamics is generally not subject to density dependence factors. On the other hand, these populations may exhibit small numbers, which may preclude the use of deterministic models. We will show step by step how to use it for extinction prediction.

Although extensive literature exists on the frequentist estimation for the Galton-Watson processes (see for instance [1, 12, 18, 34] and references therein), we choose a Bayesian reasoning [41]. In contrast with a frequentist approach, it has the advantage that a prior with expert knowledge or previous analyses can be introduced in the models. This may be especially useful for small data sets. Indeed, data sets for rare and elusive species, which are often the target of PVA, are generally scarce and incomplete. Moreover, a posterior allows a precise quantification of the uncertainty due to the parameter estimations. This gives posterior probability estimates or credible intervals for the various parameters allowing to build natural scenarios without arbitrary choice. This is of particular interest for practitioners.

Besides, our approach is based on a remarkable relationship between the Galton-Watson processes and Bayesian inference. As we will see, the Galton-Watson processes have multinomial transitions and verify the Markov property, the law of the parameters conditioned according to the historical paths (namely the posterior) can have an explicit form. Indeed, choosing a classical conjugate prior of the multinomial distribution enables recovery of this distribution as a posterior (with different parameters) after one generation and then the Markov property enables to generalize this property for all generations.

This allows to use classical formulas on Galton-Watson processes and integrate it under the posterior to predict quantitatively the behaviour of the population; that is the growth, extinction, relative abundance…

This simple and powerful expression of the posterior is one of the main strengths of our approach. Indeed, this remarkable property seems to be observed only in [37] in the mono-type case (i.e. $K=1$), and, to the best of our knowledge, has never been applied in practice with a mono-type nor with many types. We will show how this leads to a simple and powerful way to analyse the viability of the population.

Outline: In the two next Sections, Section 2 and Section 3, we present the Galton-Watson processes and our stastical model. Section 3 is more technical and contains several formulas and properties of such processes. Section 4 states the posterior on which we will apply the properties developed in Section 3. Section 5 is devoted to two pedagogical examples. the first is based on synthetic data and is concerned with extinction. The second concerns the French bear population and survival, which apart from illustrating our method, is of great interest in is own right. Finally Section 6 gives some useful simple extensions.

Box 0. Notation

•

$K$ is the number of types.

•

$N(t)=(N_{1}(t),\dots,N_{K}(t))^{\intercal}$ designes the state vectors (sample size of each type) at time $t$.

•

$p_{i,j}(k)$ is the probability that an individual of type $i$ gives $k$ individuals of type $j$ after one unit of time. It is supposed to be a random number.

•

$\kappa_{i,j}$ is the maximal number of offspring of type $j$ from an individual of type $i$.

•

$(\alpha_{i,j}(k))$ design the prior knowledge in terms of sample size of individuals in type $i$ which gives $k$ individuals of type $j$.

•

$\xi_{i,j,l,t}$ is the (random) number of offspring of type $j$ of the $l$-th individual of type $j$.

•

$M=(M_{i,j})=\left(\sum_{k=0}^{\kappa_{i,j}}kp_{i,j}(k)\right)$ is the random mean matrix.

•

$\lambda$ is the principal eigenvalue of $M$.

•

$u,v$ are respectively the left and right eigenvector of $M$.

•

$T_{\text{Ext}}$ is the (random) extinction time. Namely the duration after the present time $T$ to have extinction.

•

$s_{i}$ is the probability of extinction starting from one individual with type $i$.

•

$\mathbf{n}=(n_{i,j}(k,t))_{1\leq i,j\leq K,k\geq 0,T\geq t\geq 0}$ are the data. they correspond to the number of individuals of type $i$ with $k$ offspring of type $j$ at time $t$ for all types and times over $\{0,...,T\}$.

•

$\alpha$ is a statistical threshold. Classicaly $\alpha=0.05$ or $0.1$.

Box 1. Pathways to PVA through Bayesian-Galton-Watson approach

Here, we summarize some steps for carrying out a viability study using data. Several details are given in Appendix 8.2.

(1)

Previous knowledge: Initialize your hyper parameters $(\alpha_{i,j}(k))$.

(2)

Compute the following quantity of interest:

•

Short time evolution: number of individuals:

$\mathbb{E}[N(t)\ |\ \mathbf{n}]=\mathbb{E}[M^{t}\cdot N(0)\ |\ \mathbf{n}].$

•

Probability of viability:

$\mathbb{P}(\lambda>1\ |\ \mathbf{n}).$

If $\lambda<1$ then reintroduction is useless, one has to act on the environment.

•

Time to extinction: if $\lambda<1$, extinction is certain and one can draw

$U(t):t\mapsto\mathbb{E}[\lambda^{t}\ |\ \mathbf{n}]\times\frac{\sum_{j=1}^{K}N_{j}(0)v_{j}}{\min_{1\leq j\leq K}v_{j}}.$

Indeed when $U(t)$ is below a threshold $\alpha$ (as $1$ or $5\%$) we know that after this time there is more than $1-\alpha$ that all individuals are extinct; that is $T_{\text{Ext}}\leq t$. An upper bound for $T_{\text{Ext}}$ is given in Section 8.2.

•

Probability of extinction: if $\lambda>1$ then the population has the capacity to survive but extinction can occur if the population is too small:

$\mathbb{P}(\text{Ext})=\mathbb{E}\left[s_{1}^{N_{1}(0)}\times\cdots\times s_{K}^{N_{K}(0)}\ |\ \mathbf{n}\right]$

•

How to plan reintroduction: draw the laws of $s_{1}$,$\dots$, $s_{K}$ to see the survival potential of each type and therefore choose to reintroduce the more efficient one.

2. Galton Watson process: a Stochastic Matrix Population Model

One of the most famous and powerful mathematical settings for PVA is the Matrix population model. They were introduced by [4, 32, 31] and they have become a classical tool to predict the evolving size of a population as illustrated in the books [13, 46]. In this model, a population is classified into discrete types and its evolution is in discrete time. The state of the population, at a certain time $t\in\mathbb{N}=\{0,1,...\}$, is given by a (column) vector, $N(t)=(N_{1}(t),\dots,N_{K}(t))^{\intercal}$, where $K$ is the number of types (and the exponent $\intercal$ designs the usual transposed vector) . This vector evolves from time $t$ to $t+1$ with the following matrix product:

[TABLE]

that is, for every type $j$,

[TABLE]

Matrix $A$ is called the population projection matrix and the entry $A_{i,j}$ tells how many individuals of type $i$ appear per individual of type $j$. Classically the types correspond to the age (or the location), and $A$ then codes the development (or the migration), the birth and the death.

This classical model makes it possible to understand population variation and resulting extinction or viability. Indeed, let us consider for instance the one-type case $K=1$. We then have

[TABLE]

for all $t\geq 0$, which entails

[TABLE]

In particular when $A_{1,1}>1$ the population is viable and when $A_{1,1}<1$ the population goes into extinction. This result can easily be generalized with $K$ types using the principal eigenvalue of $A$. Namely if $\lambda$ is the largest eigenvalue of $A$ (in modulus) and $N$ the associated eigenvector then $N(t)\underset{t\to\infty}{\sim}\lambda^{t}N$. The vector $N$ then relates the asymptotic relative abundance although $\lambda$ represents the growth rate. See for instance [13] for details.

In case of extinction, a central question in PVA is to have an estimate of the extinction time. Nevertheless, any prediction of the extinction time is impossible since population size is a deterministic real number. Namely, we have $N_{1}(t)>0$ for all $t\geq 0$ as soon as $N_{1}(0)>0$. When $K=1$ and $A_{1,1}<1$, or more generally $\lambda<1$, extinction means $\lim_{t\to\infty}N(t)=0$ but at each time there are some individuals even if this number can be very small. Hence, mathematically the extinction time cannot be defined as the time $T_{\text{Ext}}$ that $N_{1}(T_{\text{Ext}})=0$ for the matrix population model. The usual setting to quantify time to extinction in this setting is therefore the so-called quasi-extinction time [24, 28, 35, 17]. That is, we fix some arbitrary threshold $\gamma$ under which the number of individuals will be insufficient to ensure persistence of the population, and the associated quasi-extinction time $T_{\gamma}$ is defined as the moment the population size reaches the level $\gamma$. This threshold has to be chosen in order to take into account demographic stochasticity, Allee effect, etc. When $K=1$, Equation (2) gives

[TABLE]

and we see that the problem is that this time crucially depends on the arbitrary choice of $\gamma$.

Another drawback to considering real numbers to model the population size is the unit choice. Indeed does $N_{1}(t)=1$ mean that there is $1,10,10^{2},10^{3},...$ individuals? Then how can we use this type of model for small population sizes as in [10, 21] or [19, Section 5] for instance? Note that, as deterministic models are proved to be realistic only in large populations, this unit choice has to be large. But how many individuals is enough to be large depends on the model and its parameters. Even if deterministic models can give suitable results for only a hundred individuals, a unit choice of $10^{3}$ (or even $10^{6}$) may not be enough to avoid atto-fox type problem [11] (that is having $10^{-18}$ individuals able to avoid extinction).

As pointed out in [25, 30], to overcome these problems, it is more natural to consider the population size as a discrete number (namely an integer) as this is the case in Galton-Watson processes.

Multi-type Galton-Watson processes, were introduced in [5, 47], and can be seen as a matrix population model with offspring randomness (but with non-random environment). In this approach, at each time $t$, each individual $l$ of type $i$ gives a discrete random number $\xi_{i,j,l,t}\in\{0,1,2,\dots\}$ of child type $j$; namely formula (1) became

[TABLE]

Of course setting

[TABLE]

we recover the form of equation (1) that is

[TABLE]

With this formalism, $A(\cdot)$ is a random matrix changing at each time. However it is not the same as in [46, Chapter 3] because this sequence of matrices is not a sequence of independent and identically distributed matrices. To see how demography or environment play a role in this matrix notation and to see the importance of considering demographic stochasticity (even in larger populations), we can see [19] and references therein.

Consistently, as in matrix populations, one can recover the asymptotic type distributions and the increasing rate of the population through decomposition of a matrix, but one can go further and find sharp results on extinction. The credit that can be given to these predictions can also be estimated. That is, we can give a precise confidence interval for the extinction risk, the extinction time, $etc.$ This again is one of the main differences with matrix population models. The proof of these results can be found in [2, 29, 26, 27, 36]. Some of them have already been applied for real-life examples; see for instance [14] or [25, page 208]. However, besides using the mathematical properties of Galton-Watson processes, we also propose in this work an efficient and consistent statistical approach to using it. Indeed, as noticed in [39] "Conducting a population viability analysis involves the steps of choosing an appropriate model, fitting the model to data, and using the fitted model to predict the extinction risk ", the approach that we describe aims to cover all these PVA steps from some life tables or more generally from demographic data.

To conclude, one should keep in mind, that as with matrix population models, the Galton-Watson process crucially depends on the two following strong assumptions:

1. (1)

interaction between individuals has no effect on the population dynamics; 2. (2)

parameters of one individual depend on this individual only through its type. In particular, all individuals of the same type have the same distribution of offspring and this does not vary in time.

3. Our model: a Galton-Watson process embedded in a Bayesian framework

To define a multi-type Galton-Watson process, we need to introduce some probabilities $p_{i,j}$ on $\{0,1,...\}$. They will represent the progeny of type $j$ from an individual of type $i$. More precisely, an individual oftype $i$ will give, after one unit of time, $k$ individuals of type $j$, with probability $p_{i,j}(k)$. We assume the existence of a maximal number of offspring $\kappa_{i,j}$; that is $p_{i,j}(k)=0$, for all $k>\kappa_{i,j}$. We have then $K^{2}$ probabilities $p_{1,1},p_{1,2},...,p_{K,K}$ which are represented in a vector form $p_{i,j}=(p_{i,j}(0),...,p_{i,j}(\kappa_{i,j}))$.

For real life problems, these (fundamental) parameters $p_{i,j}$ are unknown. In a frequentist framework, these parameters of interest would have been assumed to be unknown, but fixed. Namely, it would be assumed that in the population, for all $i,j$ there is only one true probability distribution $p_{i,j}$. In the Bayesian view of subjective probability, all of these unknowns are treated as uncertain and therefore should be described by a probability distribution over probability distribution.

We will then suppose that $(p_{i,j})_{1\leq i,j\leq K}$ is distributed according to some probability distribution, usually called the prior. All through this paper, this prior will be the convolution of Dirichlet distribution [41, Section A.8 p. 521]; namely $p_{i,j}$ are independent from each other and probability-valued random variables whose law is given for every $i,j$ by

[TABLE]

for some sequence $(\alpha_{i,j}(k))_{k=0,...,\kappa_{i,j}}$. This is a classical law on discrete probability whose expression is given in Appendix 8.1. Parameters $(\alpha_{i,j})$ are called hyper parameters and correspond to the information that you want to incorporate in your estimation. They are fixed by the users and can take into account previous studies or expert knowledge.

Among all the choices, two main strategies can be chosen.

The first one corresponds to the non-informative choice. If we have no information on the ecological system, we wish to decrease the effects of prior outcomes. We then take $p_{i,j}$ as drawn uniformly at random over all probability vectors. This corresponds to the choice $\alpha_{i,j}(0)=\dots=\alpha_{i,j}(\kappa_{i,j})=1$. Namely, we consider as prior, the uniform law on the simplex of probability measure on a discrete space. This choice is motivated by the minimization of the entropy of these distributions (See [41, Section 3.2.3] and [42, Section 6]). Roughly, it is the more random choice or the less informative one.

Another choice can consist of putting information in the hyperparameters $(\alpha_{i,j})$. Indeed, as Dirichlet distribution is uni-modal (at least for large parameters) then one can easily choose $(\alpha_{i,j})$ so that $\mathbb{E}[p_{i,j}(k)]$ is equal to a presupposed valued $m_{k}$ plus or minus a fixed error estimate $\sigma$ (uniform in $k$). Values of $m_{k}$ and $\sigma$ can for instance be taken from a previous study where data are not available but just the estimator with a confidence interval as in [15, 14]. Then the expert choice is driven by external data and the weight of the expertize is automatic. Finally $(\alpha_{i,j})$ can incorporate belief which does not come from any data. These parameters also represent the balance between data and expertize.

Details for initializing $(\alpha_{i,j})$ in practice are given in Appendix 8.2.

Forthwith, a Galton-Watson process $N=(N(t))_{t=1,2,...}$ is a random vector-valued sequence satisfying equation (3). In this equation, conditionally on $(p_{i,j})$, random variables $\xi_{i,j,l,t}$ are supposed to be independent and $\xi_{i,j,l,t}$ are distributed according to $p_{i,j}$ on $\{0,1,...,\kappa_{i,j}\}$. Namely,

[TABLE]

for all sets $A_{i,j,l,t}\subset\mathbb{N}$. This is a discrete-time model. Although there is a continuous-time version of the Galton-Watson processes and our result should hold for continuous times, we restrict ourself to a discrete-time setting because of the general nature of the data. Indeed data is often collected at punctual moments. Number $N(t)$ then represents the number of individuals at $t$ years, $t$ days or more relevant choices of times.

An important quantity associated with this model the random mean matrix $M$ whose entries are

[TABLE]

It is a random matrix since $p_{i,j}$ are supposed random (Bayesian framework). This matrix is the counterpart to the population projection matrix in matrix population model [13, 46]. When $p_{i,j}(0)=1$ for $j\notin\{0,i,i+1\}$ that is one can just get old or reproduce then $M_{i,j}=0$ for $j\notin\{0,i,i+1\}$ and we recover the classical Leslie matrix. Under general assumptions [26, Section 2.3.1], the evolution of $N(t)$ when $t$ becomes large only depends on the largest eigenvalue $\lambda$ of this matrix and the two associated eigenvectors $u,v$ (with positive coordinates). That are those verifying

[TABLE]

and

[TABLE]

On the event $\lambda<1$ then the process goes to extinction; that is there exists a finite random time $T_{\text{Ext}}$ such that $N(t)=0$ for all $t\geq T_{\text{Ext}}$ (and $N(T_{\text{Ext}}-1)>0$). Even if one can estimate $T_{\text{Ext}}$ by simulation, sharp analytic bounds have been proved. For instance,

[TABLE]

See for instance [26, Equation (5.54)]. This bound allows to find an upper bound of the extinction time with high probability (see Box 1 and Section 8.2). To our knowledge, there was no equivalent lower bound in great generality. However, [26, Box 5.2 p. 119] gives such a bound when $K=1$. This allows to have a lower bound on the extinction time. In Appendix 8.2, we generalize this result for $K>1$ and show how use it and (6) to bound the extinction time.

On the event $\lambda>1$, extinction occurs with positive probability (but not equal to $1$). In fact, if Ext denotes the extinction event then we have

[TABLE]

where $s=(s_{1},...,s_{K})$ is the unique solution of $s=\varphi(s)$ in the simplex $\triangle=\{s\in[0,1]^{K}\ |\ \sum_{i=1}^{K}s_{i}=1\}$. Function $\varphi=(\varphi_{i})_{1\leq i\leq K}$ is the generating function associated with $p$; for $i\in\{1,...,K\}$, it is defined by

[TABLE]

In particular, when we start with only one individual with type $i$ (that is $N_{i}(0)=1$ and $N_{j}(0)=0$ for $j\neq i$) then

[TABLE]

Additionally, Process $(N(t))$ is a Markov process. That means that for a current time $T$, the future evolution of the population only depends on the distribution of the population at time $T$ and not on periods of time in the past (conditionally on $(p_{i,j})$ of course). As a consequence, Formula (6) Formula (7) hold when replacing $N(0)$ by $N(T)$ and conditioning on $N(T)>0$. Know that some individuals remain at time $T$, we have to take into account this information.

Finally, in addition to the spectral analysis, all results of a matrix population model (sensitivity analysis,…) can be applied using the matrix $M$ for the study of the Galton-Watson process. Also, for short time, the Galton-Watson process evolves in the same way as this associated deterministic model up to error around $\sqrt{N(t)}$.

All of these mathematical results can be found in [2, 29, 26, 27, 36].

4. Main result: explicit expression of the posterior

We assume that our demographic data has the form of life tables $(n_{i,j}(k,t))_{1\leq i,j\leq K,k,t\geq 0}$; where $n_{i,j}(k,t)$ corresponds to the number of individuals at time $t$ of type $i$ which gives $k$ individuals of type $j$ at time $t+1$. As for example in Table 1. We assume we know the population over the time interval $\{0,\dots,T\}$; namely our data is $\mathbf{n}=(n_{i,j}(k,t))_{1\leq i,j\leq K,k\geq 0,T\geq t\geq 0}$.

Due to the well known conjugation property between multinomial and Dirichlet distribution (see for instance [41, 42]), and the form of the Galton-Watson transitions, the form of the posterior is simple. More precisely, we can see this for every $i,j$, the posterior of the parameter $p_{i,j}$ conditioned on the data $\mathbf{n}$ are

[TABLE]

The proof of this formula is given in Appendix 8.1.

This posterior turns into an extensive summary of the information available on the parameter, integrating simultaneously modelling choices and available observations. To answer our main questions, we must summarize this object into the principal objects of interest (risk of extinction, extinction time, abundance, …). Box 1 and Appendix 8.2 gives a typical "to-do" list in order to answer some main questions in PVA. We directly use the posterior to calculate precisely the probability of some events. For instance, one can calculate

[TABLE]

directly integrating Equation (7) under the Dirichlet distribution given in Equation (9). This has the advantage of not being conditioned on one $(p_{i,j})$ because we take a mean value of all $(p_{i,j})$ combinations, using our posterior. This perfectly describes the prediction randomness coming from the demographic stochasticity and the statistical estimation.

One possible drawback is that all quantities of interest (abundance…) are random objects and may be difficult to understand or describe. Using the full posterior can then be difficult for practitioners and it can be helpful to sum up the information. We describe here two possible ways of aggregating the posterior according to the objectives of the outcomes: an estimator approach and a scenario approach.

Firstly, as explained in detail in [41, Section 4], the posterior mean, the posterior mode, the posterior median, $etc.$ can be used to build an estimator of the input of interest. If a choice must be made among the types of estimators above, there is no way of selecting a best estimator without using an error-type criterion; see [41, Section 4] for details. Note that some of these types of estimator can be easily computed; for instance

[TABLE]

The main drawback of this approach is that it loses the error estimation which is naturally embedded in the posterior. However, as soon as these values are fixed one can directly calculate $\lambda$ (estimation of viability), probability of extinction $s_{1},s_{2}$ …

Finally, an intermediate approach between the full posterior approach and the punctual estimation consists of building scenarios from the posterior. In PVA, in order to make some decisions, it is (more or less) usual to build some scenarios for the future; see for instance [14, 15]. Nevertheless, in these papers, the variation of the main estimation to build the difference between the scenarios is not so clear. A natural choice for making scenarios that take into account the data and the error (due to the variance) that we make in the estimation is to consider different quantiles of the posterior. For instance, from Table 2, we can choose for $p_{0,1}(1)$, the values $0.67$, $0.82$, $0.93$ for scenarios that underestimate, estimate or overestimate the mortality rate (note that there is no symmetry).

5. Examples

In this section, we give two examples. The first one based on simulated data aims to show in a pedagogical way our method and compare it with previous study when population goes into extinction. The second example is a real life example. It is concerned with the French Pyrenean brown bear whose viability has a major interest from an ecological and political point of view.

Synthetic data: comparison of methods

Let us consider a simple population model going to extinction. We assume there is only one type $K=1$ and the maximum number of offspring is $\kappa_{1,1}=3$. We assume that survival and reproduction events are independent. All individuals survive with probability $p_{S}=0.4$ and have $k$ offspring with probability $p_{R}(k)$ with $p_{R}$ given by the vector

[TABLE]

Namely each individual gives birth to $1$ child with probability $0.1$, $2$ with probability $0.05$ and $3$ with probability $0.05$. This is a Galton-Watson process with non-random parameters. With our previous notation, we have

[TABLE]

The expected mean matrix $M^{\star}=\sum_{k=0}^{4}kp^{\star}_{1,1}(k)=0.75$ is just a number and $\lambda^{\star}=M^{\star}=0.75<1$ and then extinction is certain. We note $M^{\star},\lambda^{\star}$, $p_{1,1}^{\star}$, $etc.$ for the true parameters instead of $M,\lambda,p_{1,1},$ $etc.$ because we keep this notation for the random variables distributed according to our prior which will serve to estimate $M^{\star},\lambda^{\star}$, $p_{1,1}^{\star}$, $etc.$ Let us begin by simulating a population for $5$ years to create a learning data set, that is our knowledge on the population. This is resumed in Table 1. Focus on this table, we recover that

[TABLE]

For instance, at time [math] there are $100$ individuals of whom $47$ die without offspring, $39$ survive without offspring or die with one offspring $etc.$ There are then $75$ individuals at time $1$ and so on.

Continuing the simulation we see that the population evolves as follows

[TABLE]

[TABLE]

[TABLE]

Therefore $T_{\text{Ext}}^{\star}=7$ in this sample. The aim of what follows is to give bounds on this quantity using only Table 1.

We choose a non-informative prior $\alpha_{1,1}=1$. We then find that, conditionally on data of Table 1, $p_{1,1}$ is distributed according to a Dirichlet law with parameter $(145,128,20,14,8)$. Indeed it is the sum of each line of Table 1 plus $\alpha_{1,1}=1$. We then find

[TABLE]

That is we find that the posterior mean estimator is relatively close to the true value $M^{\star}=0.75$ and that the population is not viable with a probability very close to $1$.

Now using the method developed in Appendix 8.2 with thresholds $\alpha=0.05$ we find that the extinction time should be between $3$ and $31$; which is what happens because $T_{\text{Ext}}=7$.

We compare our results with two other methods.

The first one is based on a diffusion model [22]. Using the method and notation of Foley [22], we find that the mean and variance of the growth rate are respectively $r_{d}=-0.3028$ and $v_{r}=0.0041$ in log scale. We have $e^{r_{d}}=0.7387$ which is relatively close to $M^{\star}=0.75$. We then obtain a good parameter estimation. However to estimate the time to extinction, we have to choose a carrying capacity charge. At best, we obtain $36548$ which is far from the real extinction time. This could be expected because we are not in the setting of [22] which is concerned with a more stable population (not decreasing one) whose extinction arises with rare events (designed as environmental stochasticity by Foley).

We also use a naive approach based on linear regression. Starting from Equation (2), we can fit the trend line:

[TABLE]

and then look at the time spent given by the upper and lower confidence lines (at $90\%$ threshold) to reach [math]. The lower bound is $7$ and the upper bound is $9$. These bounds are correct because they include the true value $T_{\text{Ext}}=7$.

Taking only one simulation has only a qualitative or pedagogical interest. Let us repeat the previous steps $1000$ times to see how many times the extinction is well predicted. With our method the extinction time is in our interval in around $93\%$ of simulations instead of $49\%$ from the naive regression. To end this example, we plot the density of the extinction through these $1000$ simulations in Figure 1 and Figure 2. In Figure 1, we add the normal law with the same first two moments. We see that the distribution of the extinction time is not Gaussian at all: it is not symmetric and has a heavier tail. It is therefore not suitable to search for some values and build a confidence interval through Gaussian quantiles. In Figure 2, we can illustrate the law of our lower and upper bound.

Real data: French Pyrenean brown bears

The Pyrenean brown bear (Ursus arctos) population is considered as one of the most seriously threatened with extinction in Western Europe. In the 90’s and early 2000s, the reinforcement of the population by the introduction of a few individuals from Slovenia allowed to avoid population extinction. Yet, the population remains fragile, and the introduction of new individuals gives rise to tumultuous debates in France. Let us give here a short study on the powerful properties of our approach (and the difference with the preceding example).

The data we used comes from [10]. It corresponds to the exhaustive supervision of all the population between $2005$ and $2016$ of the French Pyrenean brown bears. This population is split into two different isolated subpopulations. In the western part of the Pyrenees, before 2018, there were only two males (of which only one is indigenous), therefore we will focus on the subpopulation living in the central part of the Pyrenees.

We consider the same structure model as [15]. Namely, we only consider the evolution of the number of females and we consider a population structured by age, with $K=5$ classes. For classes $i\in\{1,\dots,4\},$ this corresponds to bears whose age is $i-1$ years. For $i=5$ this corresponds to the bears whose age is greater than or equal to $4$ years. The life of a bear is modelled as follows: during its first 4 age stages, it can either die or pass to the next age. When it is in the last fully developed stage, it can die, survive and reproduce. See [15] for biological motivations.

We assume that surviving and reproducing are two independent events even if generally it is assumed that reproduction is possible only when the female parent survive. Both choices are modelization preferences that marginally add a bias to the prediction. It can have an importance mostly for the interpretation of the estimated parameters but we are not doing this here. Indeed, as we estimate our paramaters with this assumption, it decreases the values of reproducing parameters (because we will not see offspring when parents die) but it will increase the offspring predictions (because it is possible to have offspring when parents die). These two facts are balanced because they are assumed for the estimation and for the prediction. In the same way, we assume that females can reproduce each year although it is an exceptional event. This changes nothing because it will naturally divide by $2$ the reproduction rate.

Let us now follow the steps of Box 1.

We start by choosing the information that we want to take into account. Begin by considering the non-informative setting, that is, for all $i\in\{1,...,4\}$ and $k\in\{0,1\}$

[TABLE]

and for $k\in\{0,...,3\}$

[TABLE]

Using Formula (9), the posterior law is easily calculable from the data of [10] and is then given in Table 2. The Beta distribution can be seen as a particular case of Dirichlet distribution.

However, we only have 11 years of data and one would think that this is not sufficient to estimate all the parameters we need. Let us show how we can integrate the statistical estimation of [15]. The authors used data coming from different articles with not necessarily the full description of the estimation procedure. If we look at the infantile mortality rate $p_{1,2}(0)$ then it was based on [48]. Their estimation is $0.4$ for a population of $150$ bear cubs. The resulting estimated variance is then $1.96\sqrt{0.4\times 0.6/150}=0.0784$. To take into account this study, we can then choose

[TABLE]

Using these hyper parameters and the data, we find that the posterior law is a Beta$(18.2417,3.7855)$. This is very close to the non-informative estimation Beta$(18,4)$ given in Table 2. Hence, even if the sample size data seems limited,there is enough information to limit the sensitivity of the prior. Note that we can nevertheless force $(\alpha_{i,j})$ to be large enough to have an impact.

We then restrict ourselves to the non-informative prior because it gives similar results and several incomes are unknown as for instance for choosing $(\alpha_{1,5})$.

We can now see the evolution of the population over short time spans. To illustrate this, let us represent for each year $t$, the prediction we give from the previous years $2005,\dots,t-1$. We represent these successive predictions in Figure 3. Each boxplot represents the distribution of the prediction. The mean estimator given by Formula (10) was added (diamond symbol). To compare with the genuine evolution of the population, we add the total number of female bears at each year (red circles). In this figure, we can for instance see the learning step of our algorithm.

Let us now focus on $\lambda$; that is the viability problem. For this model the principal eigenvalue $\lambda$ is then the solution of the equation

[TABLE]

which has nevertheless no explicit solution (this equation comes from the fact that $\lambda$ is just the root of the characteristic polynomial). However, numerically, using the algorithm given in Appendix 8.2 one can find

[TABLE]

Hence, with high probability the process is super-critical and the reproduction capacity seems to guarantee viability. However if the population is too small, the survival is not guaranteed. Again, surviving probabilities $s$ in Equation (7) have no explicit solutions but are functions of $p_{i,j}$ and one can use the algorithm of Appendix 8.2 to find that, with the current population size, the probability of extinction can be rounded to [math]. It seems that the size of the population and the reproduction capacity is large enough to guaranty survival without reintroduction in the central part of the Pyrenees. Let us emphasize again that this result is only valid under our main assumptions. That includes no change in the environmental conditions including hunting or other human interference. Note that we also assumed no change in the basic demographic parameters that could be due to inbreeding depression. In fact, the current population is highly inbreed [10].

However one can still see in Figure 4 that it is more beneficial to reintroduce older bears than younger ones. In particular, in 2018, two females were reintroduced in the western part of the Pyrenees. Assuming similar evolution, there is around $15\%$ risk of extinction for this sub-population (of this two males and 2 reintroduced females). Additionally, if these females are pregnant and then if they survive and give birth each to one female cub then this risk will be divided by $2$ (around $8\%$). Closely related to theses questions, using Formula (7), we can calculate the effective population size; namely the number of reproducing individuals we need to have an extinction risk lower than a defined threshold. Here, to have an extinction probability lower than $0.05$ then we need at least $5$ reproducing females.

6. Extension to different types of dataset

To apply the present method, we need to know exactly the sum over the time of the number of individuals passing from type $i$ to $j$ at each time. This type of information is almost as difficult to collect as the complete genealogical tree of each individual. Even if several such datasets exist and to the best of our knowledge, few or no previous methods exhaustively use complete information. It is not always possible to have access to this type of full dataset [10, 14, 21]. We explain some extensions to generalize our approach in this section. This was not done simply to clarify.

In the synthetic example of Section 5, we have to know about the whole life of each individual (survival and reproduction) but in general we only know the number of survivors and the number of offspring, then we can write

[TABLE]

as seen written in this example. However, suppose that $p_{S}$ and $p_{R}$ are both Dirichlet distributed and independent, this is not the same modelling framework and with this modelling choice, we can study survival and reproduction separately.

Another generalisation comes from the second example of the French Pyrenean brown bears. In general we do not know the sex of the bear cub before until a few years after their birth. Consequently, we do not know $n_{4,1}(\cdot,t)$ for the last year, $t$. However, one can slightly change the Bayesian model to take into account this missing information. Indeed, one can add some structure and assume that

[TABLE]

where $q$ is the law of the number of spring for a female bear (male and female bear cubs) and $p$ is the probability that a cub is a female. Using this (natural) representation then we keep the conjugation property, and all the posterior remains explicit. For the dataset [10], this leads, for instance, to $0.583$ as posterior mean estimator for the sex ratio. This is slightly more favourable than the estimation of $0.5$ used in [15]. Using the number of males in the population and supposing that the survival rates are similar can reduce the variance of our estimation; this is a compromise between the data and expertise.

The assumption that we have made, which can be restrictive, is that surviving and reproducing are independent, but for prediction in contrast with parameter interpretation, this is not a real problem. Nevertheless, one can extend our setting to model the necessity of survival to make reproduction possible. The assumption that all $p_{i,j}$ are independent from each other can be removed. A careful reading of the proof of Formula (9) shows that we can extend it for the following setting: for every $k_{1},...,k_{K}$, we replace all $p_{i,j}(k)$ by $p_{i}(k_{1},...,k_{K})$ which is the probability that an individual with type $i$ has $k_{1}$ offspring of type $1$ and $k_{2}$ offspring of type $2$ $etc.$ It remains to assume that all $p_{i}$ are independent and keep a Dirichlet law (on $\prod_{j=1}^{K}\{1,..,\kappa_{i,j}\}$ instead of $\{1,...,\kappa_{i,j}\}$) as prior to have the same type of approach.

Finally, In many examples coming from ecology, the number of offspring may be very large and it is not possible to fix any value $\kappa_{i,j}$ for the maximum numbers of offspring. It is therefore not possible to directly use our approach unless we take a large $\kappa_{i,j}$ (this requires too much information). In order to keep the same idea for our approach, we can suppose that offspring are given by a Poisson law $\mathcal{P}(m)$ (which arises naturally as $\kappa_{i,j}\to\infty$ as a limit of our model under general assumptions [3]). Also the naturally associated prior will be the Gamma law. It also verifies properties (conjugation, non-informative, expertise, $etc$.) developed in Section 3 and we also have an explicit posterior distribution.

7. Discussion

We have shown how the Galton-Watson process with Bayesian inference is a simple but powerful setting for PVA analysis. On top of giving the same information as the matrix population model, it easily gives precise results for probability of extinction, time to extinction, $etc$. It is perfectly adapted for small populations and then perfectly completes previous studies, see [9, 19, 24, 28, 39]. Deterministic models including differential equation and matrix products are well suited to large population (but sometimes more than $10^{6}$ individuals [11]). Stochastic differential equations or hybrid models (including piecewise deterministic models) are more appropriate at a mesoscopic scale; at least for a thousand individuals. Few models are concerned with smaller population size although these populations are more concerned by extinction.

For conservation purposes, there is more and more population monitoring and there is more and more exhaustive information to understand population evolution. Nevertheless, classical statistical methods do not use this exhaustive information to calibrate demographic model. The Bayesian setting we present offers the advantage of taking into account all the information of these types of dataset and can be adapted for missing data or less informative data.

The combination between Galton-Watson processes and Bayesian, and known mathematical results [2, 26, 27, 36] implies simple formulas and there is no need for large simulations to answer simple questions. Results are found without using long computational algorithms with unknown errors. Consequently, it is easily scalable to big data: in contrast with ABC methods [20] or other MCMC methods [44], there is no need for a large number of simulations; this permits us to use this method for large datasets. One can easily use this model as part of a larger model to motivate control of the population or others strategies. Furthermore, there is no doubt about the convergence of this type of algorithm and the associated errors due to such computations.

For small data sets, our approach enables to add expert knowledge and gives precise confidence intervals or error estimation for the prediction. This takes into account statistical errors and the demographic stochasticity. However, the environmental stochasticity is not directly included in this estimation even if the Bayesian learning picks up some part of it with the data.

This model is founded on several assumptions and in general, some of them can be questionable such as the absence of carrying capacity, the independence of individuals and the static environment but they are very natural considering a small population during short time period. However, genetic aspects that can be of importance for small populations is not taken into account. Considering the example of the bears, we do not speak about the male reintroduction which can have a real interest to avoid inbreeding depression. We focus on demographic aspects and numbers. Including the inbreeding issue in our setting may be feasible and could provide interesting insights [16, 27].

8. Appendix

8.1. Proof of the expression of the Dirichlet distribution

First we recall that

[TABLE]

means that for all $a_{1},...a_{\kappa_{i,j}},b_{1},...,b_{\kappa_{i,j}}$, setting $A=[a_{1},b_{1}]\times\cdots\times[a_{\kappa_{i,j}},b_{\kappa_{i,j}}]$, we have

[TABLE]

where $\Gamma$ is the classical Gamma function. Then, by definition of the model we have the following prior:

[TABLE]

Let us now detail the likelihood of the model using multinomial laws. Let us denote by $N_{i,j}(t,k)$ the number of individuals with type $i$ having $j$ descendants of type $i$ at time $t$. Namely

[TABLE]

We have for all $t\geq 0$

[TABLE]

where $\mathbf{M}(n,q_{0},...,q_{\kappa})$ denotes the usual multinomial law with $n$ trials and $q_{0},...,q_{m}$ events probabilities. Namely $M\sim\mathbf{M}(n,q_{0},...,q_{\kappa})$ means that for all $m_{0},\dots,m_{\kappa}$ we have

[TABLE]

when $\sum_{i=0}^{\kappa}m_{i}=n$ either it is equal to [math]. Using these two explicit formulas (multinomial and Dirichlet densities) and the Bayes Theorem, it is easy to check the classical conjugation result:

[TABLE]

Further details can be found in [41, 42]. It remains to prove that it holds true for any times through a chain rule type argument. The Markov property (derived by the independence assumption) gives that

[TABLE]

Using this formula, the expressions of multinomial and Dirichlet densities and Bayes Theorem again finalise the proof of Equation (9). Note that $\mathbf{n}$ that represents the data is a realisation of the random variables $\left(N_{i,j}(k,t)\right)_{0\leq k\leq\kappa_{i,j},0\leq t\leq T}$.

8.2. Pathways to PVA through the Bayesian-Galton-Watson approach: the expanded version

Here we give details for each step of Box 1.

Previous knowledge

Initializes $(\alpha_{i,j}(k))$ as follows:

1. i)

If one cannot pass from state $i$ to state $j$, set

[TABLE]

or equivalently, in a certain sense,

[TABLE] 2. ii)

If you have no information, set

[TABLE] 3. iii)

If we use previous estimates and we expect that $p_{i,j}=m_{i,j}(k)\pm\epsilon_{i,j}(k)$ and $\epsilon_{i,j}(k)$ is centred and has variance of magnitude

[TABLE]

then set

[TABLE]

if $\sigma^{2}_{i,j}<1$ else $\alpha_{i,j}(k)=1$. 4. iv)

If we want to impose an expert opinion. The experts think that

[TABLE]

and that their opinion is as important as $M$ data then we can choose

[TABLE]

We can see in case iii), that this choice enables to match the expectation to the variance. More precisely, to match the expectation, we need

[TABLE]

with $S_{i,j}=\sum_{k=0}^{\kappa_{i,j}}\alpha_{i,j}(k)$. We have then $\alpha_{i,j}(k)=S_{i,j}m_{i,j}(k)$. To match the variance, we have

[TABLE]

hence

[TABLE]

Of course we need $\sigma_{i,j}^{2}<1$. However, if this is not the case, the noise is too large and the information on $p_{i,j}$ is therefore less informative than no information.

Calculus of quantities of interest

For $K$ ($K=1$ or $2$), one can directly give an expression for almost all quantities through some integral calculus. However, in general, it is difficult to express them. The better way to calculate this quantity is perhaps the following simple Monte Carlo algorithms.

Due to the central limit theorem the error of this simple algorithm can be estimated through an estimator of the variance but in general one cannot fix a number of simulations before doing the calculation. However for the probability of viability and the probability of extinction, the variance can be bounded by $1/4$ and then the error of this simple algorithm is lower than $(4\sqrt{n_{\text{prec}}})^{-1}$. For an error lower than $0.5\%$ we can take $n_{\text{prec}}=2500$.

Below, simulate $p_{i,j}$ means drawing a random variable with law

[TABLE]

Short time evolution:

Do $n_{\text{prec}}$ times the steps

• •

Simulate the $K^{2}$ probability valued random variables $p_{i,j}$.

• •

Calculate $M^{n}\cdot X_{0}$ with $M$ defined by (4).

Then do a mean value of these proportions.

Probability of viability

Do $n_{\text{prec}}$ times the steps

• •

Simulate the $K^{2}$ probability valued random variables $p_{i,j}$.

• •

Take for $\lambda$ the largest eigenvalue of the $M$ defined by (4).

Then count the proportion of $\lambda$ being larger than $1$.

Time to extinction: To bound the extinction time $T_{\text{Ext}}$ between two values

[TABLE]

with probability $1-\alpha$, for some fixed threshold $\alpha$, do $n_{\text{prec}}$ times the steps

• •

Simulate the $K^{2}$ probability valued random variables $p_{i,j}$.

• •

Calculate the function

[TABLE]

and

[TABLE]

where

[TABLE]

where $v$ is given in Equation (5).

Then do a mean value of functions $U$ and $L$ and choose $(T_{-},T_{+})$ such that

[TABLE]

These two functions are based on Equation (6) and Lemma 8.1. Note that these bounds on extinction time only hold when $\lambda<1$. To avoid some computational problems due to the fact that $\mathbb{P}(\lambda<1)>0$ (even if it is very small), it can be useful, when calculating the extinction time conditioned on the event $\{\lambda<1\}$. In the previous algorithm, this conditioning translates into doing a mean value only for $\lambda$ satisfying $\lambda<1$.

Probability of extinction:

Do $n_{\text{prec}}$ times the steps

• •

Simulate the $K^{2}$ probability valued random variables $p_{i,j}$.

• •

Find the (vector) solution $s$ of equation $\varphi(s)=s$, with $\varphi$ is the generating function defined in (8). To do this step one can use any classical optimization algorithm.

• •

Calculate $s_{1}^{N_{1}(T)}\times\cdots\times s_{K}^{N_{K}(T)}$.

Then do a mean of these quantities.

How to plan reintroduction: Do $n_{\text{prec}}$ times the steps

• •

Simulate the $K^{2}$ probability valued random variables $p_{i,j}$.

• •

Find the (vector) solution $s$ of equation $\varphi(s)=s$, with $\varphi$ is the generating function defined in (8). To do this step one can use any classical optimization algorithm.

Then do a histogram or smoothed density through classical kernel density methods for instance.

8.3. Lower bound for the extinction time

Lemma 8.1**.**

We have

[TABLE]

where

[TABLE]

Moreover, the right-hand side of equation (11) is equivalent to

[TABLE]

Proof.

The proof is inspired from [26, Box 5.2 p.119] We have

[TABLE]

where

[TABLE]

Yet, we have $\mathbb{E}[X(t)]=\lambda^{t}\sum_{j=1}^{K}v_{j}N_{j}(0)$, and it remains to study the second moment. Recall

[TABLE]

Setting $\mathbf{N}(t)=\sum_{j}N_{j}(t)$ for the population size and reindexing the population, we find

[TABLE]

where $i(l)$ is the type of the individual $l$ and $\zeta$ its offspring. We have

[TABLE]

In the third line we conditioned on $\mathbf{N}(t)$ and use the independance property. In the fourth line, we use the notation

[TABLE]

and then under our assumptions

[TABLE]

if $j\neq j^{\prime}$ and either

[TABLE]

Then

[TABLE]

and also setting $\overline{\sigma_{j}^{2}}=\sup_{i}\sigma^{2}_{i,j}$

[TABLE]

Finally setting $\Xi=\sum_{j=1}^{K}v_{j}^{2}\overline{\sigma_{j}^{2}}/\min(v)$, we obtain

[TABLE]

where we used an iteration argument in the third line. To conclude

[TABLE]

∎