Minimal quasi-stationary distribution approximation for a birth and death process

In a first part, we prove a Lyapunov-type criterion for the $\xi\_1$-positive recurrence of absorbed birth and death processes and provide new results on the domain of attraction of the minimal quasi-stationary distribution. In a second part, we study the ergodicity and the convergence of a Fleming-Viot type particle system whose particles evolve independently as a birth and death process and jump on each others when they hit $0$. Our main result is that the sequence of empirical stationary distributions of the particle system converges to the minimal quasi-stationary distribution of the birth and death process.


Introduction
Let X be a stable birth and death process on N = {0, 1, 2, . ..} absorbed when it hits 0. The minimal quasi-stationary distribution (or Yaglom limit) of X, when it exists, is the unique probability measure ρ on N * = {1, 2, . ..} such that ρ(•) = lim t→∞ P x (X t ∈ • | t < T 0 ) , for all x ∈ N * , where T 0 = inf{t ≥ 0, X t = 0} is the absorption time of X.The probability measure ρ is called a quasi-stationary distribution because it is stationary for the conditioned process, in the sense that ρ = P ρ (X t ∈ • | t < T 0 ), for all t ≥ 0.
These notions and important references on the subject are recalled with more details in Section 2, with important definitions and well known results on quasi-stationary distributions.We also provide a new Lyapunov-type criterion ensuring that a probability measure µ belongs to the domain of attraction of the minimal quasi-stationary distribution, which means that (1.1) These results are illustrated with several examples.
We use these new results in Section 3 to extend existing studies on the long time and high number of particles limit of a Fleming-Viot type particle system.The particles of this system evolve as independent copies of the birth and death process X, but they undergo rebirths when they hit 0 instead of being trapped at the origin.In particular, the number of particles that are in N * remains constant as time goes on.Our main result is a sufficient criterion ensuring that the empirical stationary distribution of the particle system exists and converges to the minimal quasi-stationary distribution of the underlying birth and death process.
We conclude the paper in Section 4, providing a numerical study of the speed of convergence of the Fleming-Viot empirical stationary distribution expectation to the minimal quasi-stationary distribution for a linear birth and death process and a logistic birth and death process.This numerical results suggest that the bias of the approximation is surprisingly small for linear birth and death processes and even smaller for logistic birth an death processes.
2 Quasi-stationary distributions for birth and death processes Let (X t ) t≥0 be a birth and death process on N = {0, 1, 2, . ..} with birth rates (b i ) i≥0 and death rates (d i ) i≥0 .We assume that b i > 0 and d i > 0 for any i ≥ 1 and b 0 = d 0 = 0.The stochastic process X is a N-valued pure jump process whose only absorption point is 0 and whose transition rates from any point i ≥ 1 are given by i → i + 1 with rate b i , i → i − 1 with rate d i , i → j with rate 0, if j / ∈ {i − 1, i + 1}.
Such processes are extensively studied because of their conceptual simplicity and pertinence as demographic models.It is well known (see for instance [20,Theorem 10 and Proposition 12]) that X is stable, conservative and hits 0 in finite time almost surely (for any initial distribution) if and only if The divergence of this series will be assumed along the whole paper.In particular, for any probability measure µ on N, the law of the process with initial distribution µ is well defined.We denote it by P µ (or by P x if µ = δ x with x ∈ N) and the associated expectation by E µ (or by , we thus have where, for any subset F ⊂ N, M 1 (F ) denotes the set of probability measures on F .
A quasi-stationary distribution for X is a probability measure ρ on N * = {1, 2, . ..} such that The probability measure ρ is thus stationary for the conditioned process (and, as a matter of fact, was called a stationary distribution in the seminal work [7]).The property "ρ is a quasi-stationary distribution for X" is directly related to the long time behaviour of X conditioned to not being absorbed.Indeed (see for instance [22] or [20]), a probability measure ρ is a quasi-stationary distribution if and only if there exists µ ∈ M 1 (N * ) such that We refer the reader to [25,20,10] and references therein for an account on classical results concerning quasi-stationary distributions for different models.
For a given quasi-stationary distribution ρ, the set of probability measures µ such that (2.2) holds is called the domain of attraction of ρ.It is nonempty since it contains at least ρ and may contains an infinite number of elements.In particular, when the limit in (2.2) exists for any µ = δ x , x ∈ N * , and doesn't depend on the initial position x, then ρ is called the Yaglom limit or the minimal quasi-stationary distribution.Thus the minimal quasi-stationary distribution, when it exists, is the unique quasi-stationary distribution whose domain of attraction contains {δ x , x ∈ N * }.From a demographical point of view, the study of the minimal quasi-stationary distribution of a birth and death process aims at answering the following question: knowing that a population isn't extinct after a long time t, what is the probability that its size is equal to n at time t?
One of the oldest and most understood question for quasi-stationary distributions of birth and death processes concerns their existence and uniqueness.Indeed, van Doorn [22] gave the following picture of the situation: a birth and death process can have no quasi-stationary distribution, one unique quasi-stationary distribution or an infinity (in fact a continuum) of quasistationary distributions.In order to determine whether a birth and death process has 0, one or an infinity of quasi-stationary distributions, one define inductively the sequence of polynomials (Q n (x)) n≥0 for all x ∈ R by (2.3) As recalled in [22, eq. (2.13)], one can uniquely define the non-negative number ξ 1 satisfying Also, the useful quantity can be easily computed (see [1,Section 8.1]), since, for any z ≥ 1, The following theorem answers the question of existence and uniqueness of a QSD for birth and death processes.
3. If S = +∞ and ξ 1 > 0, then there is a continuum of QSDs, given by the one parameter family (ρ x ) 0<x≤ξ 1 : and the minimal quasi-stationary distribution is given by ρ ξ 1 .
Remark 1. Theorem 2.1 gives a complete description of the set of quasistationary distributions for a birth and death process but is not well suited for the numerical computation of the Yaglom limit of a given birth and death process.Indeed, the polynomials Q n have in most cases quickly growing coefficients, so that the value of ξ 1 cannot be easily obtained by numerical computation.
Theorem 2.1 is quite remarkable since it describes completely the possible outcomes of the existence and uniqueness problem for quasi-stationary distributions.However, it only partially answers the crucial problem of finding the domain of attraction of the existing quasi-stationary distributions and in particular of the minimal quasi-stationary distribution.The following theorem answers the problem when there exists a unique quasi-stationary distribution.
Theorem 2.2 (Martínez, San Martín, Villemonais 2013 [19]).Let X be a birth and death process such that Then there exists γ ∈ [0, 1) such that, for any probability measure µ on N * , where • T V denotes the total variation norm and ρ is the unique quasistationary distribution of the process.In particular, the domain of attraction of the unique quasi-stationary distribution is the whole set M 1 (N * ) of probability measures on N * .
A weaker form of Theorem 2.2 has also been proved in [28] but the strong form (with uniform convergence in total variation norm) is necessary to derive the results of the next section.A generalized version of Theorem 2.2 has been rencently derived in [8], with complementary results on the so-called Q-process (the process conditioned to never being absorbed).
The case where there exists an infinity of quasi-stationary distributions is trickier and can be partially solved, as we will show, when the birth and death process is ξ 1 -positive recurrent.
Definition The birth and death process X is said to be ξ 1 -positive recurrent if ξ 1 > 0 in Theorem 2.1 and if, for some i ∈ {1, 2, . ..} and hence for all i ∈ {1, 2, . ..}, we have In the following theorem, we provide a new Lyapounov-type criterion ensuring the ξ 1 -positive recurrence of a birth and death process.As will be shown in the examples below, this criterion can be checked on a wide variety of examples and has its own interest in the domain of ξ 1 -classification for birth and death processes (see [17] and [24] for an account on this area).
Theorem 2.3.Let X be a birth and death process with infinitesimal generator L. We assume that there exists C > 0, λ 1 > d 1 and φ : N → R + such that φ(i) goes to infinity when i → ∞ and Then X admits a quasi-stationary distribution and the birth and death process X is ξ 1 -positive recurrent.
In the next theorem, we assume that the process is ξ 1 -positive recurrent and we exhibit a subset of the domain of attraction for the minimal quasistationary distribution.
Theorem 2.4.Let X be a ξ 1 -positive recurrent birth and death process with infinitesimal generator L. Then the domain of attraction of the minimal quasi-stationary distribution of X contains the set D defined by Assume moreover that there exists C > 0, λ 1 > ξ 1 and φ : N → R + such that φ(i) goes to infinity when i → ∞ and Then the domain of attraction of the minimal quasi-stationary distribution of X contains the set D φ defined by As it will be shown in the proof, we have D φ ⊂ D for all function φ satisfying the assumptions of Theorem 2.4.However, Q • (ξ 1 ) cannot be computed explicitly but in few situations.As a consequence, we won't be able to use the first criterion to determine whether a probability distribution µ belongs or not to the domain of attraction of the minimal quasi-stationary distribution.
On the contrary, we will be able to give explicit functions φ satisfying the Lyapunov criterion of our theorem for a wide range of situations.
Note that, since d 1 ≥ ξ 1 , our results immediately imply that, if the process X fulfils the assumptions of Theorem 2.3 with a Lyapunov function φ, then the process is ξ 1 -positive recurrent and the domain of attraction of its minimal quasi-stationary distribution contains D φ .This consequence is used in the following examples.
Example 1.We consider the case where b i = b i a and d i = d i a for all i ≥ 1, where b < d are two positive constants and a > 0 is fixed.Now, defining φ(0) = 0 and Since i a → ∞ when i → ∞, we immediately deduce that there exists C > 0 and λ 1 > d 1 such that φ satisfies Lφ ≤ −λ 1 φ + C. Now Theorem 2.3 implies that the process is ξ 1 -positive recurrent and Theorem 2.4 implies that the domain of attraction of the minimal quasi-stationary distribution contains Example 2. We consider now the case where the birth and death rates are constant for all i ≥ 2, that is b i = b > 0 and where b < d are positive constants.We assume that and the value of b 1 > 0 can be chosen arbitrarily.Using the same function as in the previous example, that is In particular, there exists C > 0 and λ 1 > d 1 such that φ satisfies Lφ ≤ −λ 1 φ + C. Once again, we deduce from Theorem 2.3 that the process is ξ 1positive recurrent, which was already known in this case (see [23, eq.(6.6)]).We also deduce the following new result from Theorem 2.4: the domain of attraction of the minimal quasi-stationary distribution contains the set Example 3. In the two previous examples, the birth and death rates are nondecreasing and proportional to each other.This is coincidental and is only useful to get straightforward calculations.The aim of the present example is to illustrate this on a particular case without monotony nor proportionality between the birth and death rates: we choose b i = | sin(iπ/2)|i + 1 and we get, for all i ≥ 2, As above, we deduce that the process is ξ 1 -positive recurrent and that the domain of attraction of the minimal quasi-stationary distribution contains the set of probability measures defined by The end of this section is dedicated to the proof of Theorems 2.3 and 2.4.
Lemma 2.5.We assume that there exists λ 1 > ξ 1 , C > 0 and φ : N → R + such that φ(i) goes to infinity when i → ∞ and Then there exists a constant i 0 ≥ 1 and α > 0 such that φ(i) For all i ≥ i 0 , we have φ (i) ≥ 1 and Hence, since Q i 0 (ξ 1 ) > 0 (see [22, eq.(2.13)]) and replacing φ by Q i 0 (ξ 1 )φ , we can assume without loss of generality that Our aim is now to prove that, for all i ≥ i 0 , φ(i) ≥ Q i (ξ 1 ).Because of the changes made with the function φ, this implies the inequality claimed in the lemma with the constant α = 1/Q i 0 (ξ 1 ).Assume the contrary, which means that there exists i This is feasible, since x → Q i (x) is a polynomial function of x and is then continuous in x for any fixed i. Then the function ϕ x : N → R + defined by Let us now prove that this inequality extends to any j ≥ i 1 + 1.Indeed, using the equality Lϕ x = −xϕ x and the inequality Lφ ≤ −λ 1 φ, we have, for all j > i 1 + 1, and We deduce that and thus φ(j) < ϕ x (j) for all j ≥ i 1 + 1.
Lemma 2.6.Let φ : N → R + be a function such that φ(0) = 0, φ(1) = 1 and Lφ ≤ 0. Then Proof.Under our assumption, (φ(X t )) t≥0 is a super martingale.As a consequence, for all i ≥ 1, But T 1 ≤ T 0 < ∞ almost surely and X T 1 = 1 almost surely, so that Proof of Theorem 2.3.The main difficulty is to prove that the minimal quasistationary distribution ρ ξ 1 for X exists and that Once this is proved, Lemma 2.5 implies that which is a sufficient condition for X to be ξ 1 -positive recurrent (see [23,Theorem 5.2]).Let us prove that (2.5) holds.For any M ≥ 1, let us denote by (P M t ) t≥0 the semi-group of the process X M evolving in {0, 1, . . ., M } and defined as where T M = inf{t ≥ 0, X t = M }.We also define φ M by φ M (M ) = 0 and φ M (i) = φ(i), i ∈ {0, 1, . . ., M − 1}.Now, denoting by L M the generator of the stopped process X M and setting ϕ(i) = 1 i≥1 , we thus have Hence, using Kolmogorov equations for the finite state space continuous time Markov chain X M , we deduce that .
This implies that, for any t ≥ 0, we have But X M 0 = 1 under P 1 , so that Now, by dominated convergence, we have By monotone convergence, we also deduce that We finally deduce that, for all t ≥ 0, The first consequence of this inequality is that the process X t conditioned on the event X t = 0 does not diverge to infinity.As a consequence, ξ 1 > 0 and there exists a minimal quasi-stationary distribution ρ ξ 1 for X (see [22,Theorem 4.1]).In particular, X t conditioned on X t = 0 converges in law to ρ ξ 1 .Hence, we deduce from the above inequality that, for any K ≥ 0, By monotone convergence, we obtain by letting K tend to ∞ that Proof of Theorem 2.4.Let X be a ξ 1 -positive recurrent birth and death process with minimal quasi-stationary distribution ρ ξ 1 .We prove that the domain of attraction of ρ ξ 1 contains the set of probability measures Once this is proved, the second assertion of Theorem 2.4 follows immediately from Lemma 2.5.

Approximation of the minimal quasi-stationary distribution
This section is devoted to the study of the ergodicity and the convergence of a Fleming-Viot type particle system.
Fix N ≥ 2 and let us describe precisely the dynamics of this system with N particles, which we denote by (X 1 , X 2 , . . ., X N ).The process starts at a position (X 1 0 , X 2 0 , . . ., X N 0 ) ∈ (N * ) N and evolves as follows: -the particles X i , i = 1, . . ., N , evolve as independent copies of the birth and death process X until one of them hits 0; this hitting time is denoted by τ 1 ; -then the (unique) particle hitting 0 at time τ 1 jumps instantaneously on the position of a particle chosen uniformly among the N − 1 remaining ones; this operation is called a rebirth; -because of this rebirth, the N particles lie in N * at time τ 1 ; then the N particles evolve as independent copies of X and so on.
We denote by As a consequence, the particle system (X 1 t , X 2 t , . . ., X N t ) t≥0 is well defined for any time t ≥ 0 in an incremental way, rebirth after rebirth (see Figure 1 for an illustration of this construction with N = 2 particles).This Fleming-Viot type system has been introduced by Burdzy, Holyst, Ingermann and March in [5] and studied in [6], [12], [27], [13] for multidimensional diffusion processes.The study of this system when the underlying Markov process X is a continuous time Markov chain in a countable state space has been initiated in [11] and followed by [4], [2], [14], [3] and [9].We also refer the reader to [15], where general considerations on the link between the study of such systems and front propagation problems are considered.
We emphasize that, because of the rebirth mechanism, the particle system (X 1 , X 2 , . . ., X N ) evolves in (N * ) N .For any t ≥ 0, we denote by µ N t the empirical distribution of (X 1 , X 2 , . . ., X N ) at time t, defined by where M 1 (N * ) is the set of probability measures on N * .A general convergence result obtained in [26] ensures that, if µ N 0 → µ 0 , then The generality of this result does not extend to the long time behaviour of the particle system, which is the subject of the present study.We provide a sufficient criterion ensuring that the process (µ N t ) t≥0 is ergodic.Denoting by X N its empirical stationary distribution (a random measure whose law is the stationary distribution of µ N ), our criterion also implies that where ρ is the minimal quasi-stationary distribution of the birth and death process X.Our result applies (1) to birth and death processes with a unique quasi-stationary distribution (such as logistic birth and death processes) and (2) to birth and death processes with a minimal quasi-stationary distribution satisfying an explicit Lyapunov condition (fulfilled for instance by linear birth and death processes).These two different conditions are summarized in Assumptions H1 and H2 below.
Assumption H1.There exist a function φ : N → R + and two constants Assumption H2.The birth and death process X admits a unique quasistationary distribution (S < +∞).
Theorem 3.1.Assume that Assumption H1 or Assumption H2 is satisfied.Then, for any N > λ 1 λ 1 −d 1 under H1 and any N ≥ 2 under H2, the measure process (µ N t ) t≥0 is ergodic, which means that there exists a random measure If H1 holds, then Moreover, if Assumption H1 or H2 is satisfied, then where ρ is the minimal quasi-stationary distribution of X.
1. Assumption H1 is the Lyapunov criterion which is used in Theorem 2.3 to ensure ξ 1 -positivity (and hence the existence of a minimal quasistationary distribution).This assumption also implies that the conditions of Theorem 2.4, where we determine a subset of the domain of attraction of the quasi-stationary distribution, are also satisfied.For instance, the birth and death processes of Examples 1, 2 and 3 in the previous section satisfy Assumption H1.
2. Assumption H2 is satisfied for processes that come fast from infinity to compact sets, as the logistic birth and death process (where b i = b i and d i = d i + c i(i − 1) for all i ≥ 1 with b, c, d > 0).Note that, in this particular example, an easy calculation shows that Assumption H1 is also satisfied with φ(i) = 2 i .However, this assumption is useful for any situation where it is easy to check that S < ∞, but difficult to find an explicit Lyapunov function satisfying Assumption H1.In particular, we cannot apply Theorem 2.3 on ξ 1 -positivity and, in fact, it is known that the pure drift birth and death process is not ξ 1 -positive recurrent (see [23]).As a consequence, the additional difficulty is not a technical one and the following proof cannot work in the pure drift situation.We emphasize that Theorem 3.1 for pure drift birth and death processes remains an open problem.See for instance [3] and the numerical investigation in [18] for more details.
Since the proof of Theorem 3.1 differs whether one assumes H1 or H2, it is split in two different subsections : in Subsection 3.1, we prove the theorem under Assumption H1 and, in Subsection 3.2, we prove the result under assumption H2.

Proof under Assumption H1: exponential ergodicity via a Foster-Lyapunov criterion
Step 1. Proof of the exponential ergodicity by a Forster-Lyapunov criterion We define the function where φ is the Lyapunov function of Assumption H1.Fix N ≥ 2 and let us express the infinitesimal generator L N of the empirical process (µ N t ) t≥0 applied to f at a point µ ∈ M 1 (N * ) given by where (x 1 , . . ., x N ) ∈ (N * ) N .In order to shorten the notations, we introduce, for any y ∈ N * , the probability measure For a fixed N > λ 1 λ 1 −d 1 and any constant k > 0, the set of probability measures µ is irreducible (this is an easy consequence of the irreducibility of the birth and death process X).Thus, using the Foster Lyapunov criterion of [21, Theorem 6.1, p.536] (see also [16,Proposition 1.4] for a simplified account on the subject), we deduce that the process µ N is exponentially ergodic and, denoting by X N a random measure distributed following its stationary distribution, we also have This concludes the proof of the first part of Theorem 3.1.
Step 2. Convergence to the minimal QSD Since φ(i) goes to infinity when i → ∞, we deduce from (3.2) that the family of random measures (X N ) N is tight.In particular, the family admits at least one limiting random probability measure X , which means that X N converges in law to X , up to a subsequence.
Let µ N t be the random position at time t of the particle system with initial (random) distribution X N .On the one hand, the stationarity of X N implies that µ N t ∼ X N for all t ≥ 0, and thus On the other hand, the general convergence result of [26] implies that As an immediate consequence Using Theorem 2.4, we deduce that X belongs to the domaine of attraction of the minimal QSD ρ almost surely, that is Thus the random measure X converges in law to the deterministic measure ρ, which implies that X = ρ almost surely.
In particular, ρ is the unique limiting probability measure of the family (X N ) N , which ends the proof of Theorem 3.1 under Assumption H1.

Proof under Assumption H2: exponential ergodicity by a Dobrushin coefficient argument
Fix N ≥ 2 and let us prove that the process is exponentially ergodic.Under assumption (H2), it is well known (see for instance [19]) that the process X comes back in finite time from infinity to 1, which means that Since the particles of a Fleming-Viot type system are independent up to the first rebirth time, we deduce that This implies that the FV process is exponentially ergodic.
Let us now denotes by X N the empirical stationary distribution of the system (X 1 , . . ., X N ), for each N ≥ 2. Theorem 2.2 implies that there exists γ > 0 such that, for any t ≥ t , any initial distribution µ 0 and any function f : But, for any t ≥ 0, [26] implies that As a consequence, In particular, for any > 0, there exists t and N such that But µ N t converges in law to X N , so that This inequality being true for any > 0, this concludes the proof of Theorem 3.1 under Assumption (H2).
4 Numerical simulation of the Fleming-Viot type particle system In this section, we present numerical simulations of the Fleming-Viot particle system studied in Section 3. Namely, we focus on the distance in total variation norm between the expectation of the empirical stationary distribution (i.e.E(X N )) and the minimal quasi-stationary distribution of the underlying Markov process X, when N goes to infinity.This means that we aim at studying the bias of the approximation method.
We start with the linear birth and death process case in Subsection 4.1.This is one of the rare situation where explicit computation of the minimal quasi-stationary distribution can be performed (see for instance [20]).In Subsection 4.2, we provide the results of numerical simulations in the logistic birth and death case.

The linear birth and death case
We assume in this section that b i = i and d i = 2i for all i ≥ 0. This is a sub-case of Example 1 and thus one can apply Theorem 3.1: the empirical stationary distribution of the process X N exists and converges in law, when the number N of particles goes to infinity, to the minimal quasi-stationary distribution ρ of the process, which is known to be given by (see [20]) The results of the numerical estimations of E(X N ) − ρ T V for different values of N (from 2 to 10 4 ) are reproduced on Table 1.One interesting point is the confirmation that E(X N ) is a biased estimator of ρ.A second interesting point is that the bias decreases quickly when N increases.Up to our knowledge, there exists today no theoretical justification of this fact, despite its practical implications.Indeed, one drawback of the speed of the the sequence of rebirths times.Since the rate at which rebirths occur is uniformly bounded above by N d 1 , lim n→∞ τ n = +∞ almost surely.

Figure 1 :
Figure 1: One path of a Fleming-Viot system with two particles.

Remark 3 .
The pure drift birth and death process (b i = b and d i = d for all i ≥ 1, where b < d are two positive constants) does not satisfy Assumption H1 nor Assumption H2.Note that this process is the same as in Example 2 but does not satisfy ( √ d − √ b) 2 > d 1 .

Figure 2 :
Figure 2: Estimated value of the minimal quasi-stationary distribution ρ(n) for a logistic birth and death process.