domain: sde stat

text:
Edith Gabriel
email: edith.gabriel@univ-avignon.fr

Florent Bonneu

Pascal Monestiez

Joël Chadoeuf

Adapted kriging to predict the intensity of partially observed point process data

Keywords:

We consider a stationary and isotropic spatial point process whose realization is observed within a large window. We assume it to be driven by a stationary random field U. In order to predict the local intensity of the point process, λ(x|U), we propose to define the first-and second-order characteristics of a random field, defined as the regularized counting process, from the ones of the point process and to interpolate the local intensity by using a kriging adapted to the regularized process.

Introduction

In many projects the study window is too large to extensively map the intensity of the point process of interest since observation methods may be available at a much smaller scale only. That is for instance the case when studying the spatial repartition of a bird species at a national scale, while the observations are made in windows of few hectares. The intensity must then be estimated from data issued out of samples spread in the study window, and hence, from a partial realization of the point process in this window.

In the following, we consider a stationary and isotropic point process, Φ, which we assume to be driven by a stationary random field, U. We define the local intensity of Φ by its intensity conditional to the random field U. We denote it λ(x|U). A simple example of such a process Φ is the Thomas process which is a Poisson cluster process where the cluster centers (parents) are assumed to be Poisson and the offsprings are normally distributed around the parent point. This process is stationary and the local intensity corresponds to the intensity of the inhomogeneous Poisson process of offsprings, i.e. the conditional intensity given the parent process. We will refer to the estimation of the local intensity when we want to the know it at point locations lying in the observation window of the point process, and to its prediction when point locations are outside the observation window.

Usually, when estimating a non-constant intensity, we observe the full point pattern within a window and we want to know its local changes over a given mesh. This issue has been addressed in several ways: kernel smoothing, see [START_REF] Silverman | Density Estimation for Statistics and Data Analysis[END_REF] and [START_REF] Guan | On consistent nonparametric intensity estimation for inhomogeneous spatial point processes[END_REF] in presence of covariates, and van Lieshout (2012) for a general class of weight function estimators that encompasses both kernel and tessellation based estimators; or parametric methods; see for instance [START_REF] Illian | Statistical Analysis and Modelling of Spatial Point Patterns[END_REF] for a review. A recurrent and remaining question in these approaches is which bandwidth/mesh should we use? This has been addressed by using cross-validation [START_REF] Härdle | Smoothing Techniques, with Implementation in S[END_REF] or double kernel [START_REF] Devroye | The double kernel method in density estimation[END_REF].

In contrast to the previous methods which look at the intensity changes inside the observation window, our main interest lies in predicting the intensity outside the observation window, all the more when it is not connected as it frequently happens when sampling in plant ecology. To predict the intensity we could use [START_REF] Tscheschel | Statistical reconstruction of random point patterns[END_REF]'s reconstruction method based on the first-and second-order characteristics of the point process. Once the empirical point pattern predicted within a given window, one can get the intensity by kernel smoothing. As it is a simulation-based method, it requires long computation times, especially when the prediction window is large and/or the point process is complex. As alternative method, few authors model the point pattern by a point process with the intensity driven by a stationary random field. In [START_REF] Diggle | Model-Based Geostatistics[END_REF] and [START_REF] Diggle | Spatial and spatio-temporal log-gaussian Cox processes: Extending the geostatistical paradigm[END_REF], the approach is heavily based on a complete modeling and considers a log-Gaussian model. The parameter estimation, the intensity estimation and its prediction outside the observation window are obtained using a Bayesian framework. The method developed in [START_REF] Monestiez | Geostatistical modelling of spatial distribution of balaenoptera physalus in the northwestern mediterranean sea from sparse count data and heterogeneous observation efforts[END_REF] and [START_REF] Bellier | Reducing the uncertainty of wildlife population abundance: model-based versus design-based estimates[END_REF] is close to classical geostatistics. Basically, it consists of counting the number of points within some grid cells, computing the related empirical variogram and theoretically relating it to the one obtained from the random field driving the intensity. Then, the variogram is fitted and kriging is used to predict the intensity. Its advantage is that the estimation is only based on its first-and secondorder moments so that the model does not need to be fully specified. While this approach requires less hypotheses, the model remains constrained within the class of Cox processes. Moreover, the mesh size is arbitrary defined. van Lieshout and [START_REF] Van Lieshout | Extrapolating and interpolating spatial patterns[END_REF] developed, for a wider class of parametric models, a Bayesian approach for extrapolating and interpolating clustered point patterns.

Here, we propose to interpolate the local intensity by an adapted kriging, where the kriging weights depend on the local structure of the point process. Hence, our method uses all the data to locally predict at a given point, which it is not the case of most of kernel methods. It also uses the information at a fine scale of the point process, which it is not the case in geostatistical approaches. Furthermore, it does not require a specific model but only (an estimation of) the first-and second-order characteristics of the point process.

In Section 2 we define the regularized process as a random field of point counts on grid cells and we link up the mean and variogram of this random field to the intensity and pair correlation function of the point process. The kriging weights, the related interpolator and its properties are presented in Section 3 as well as the optimal mesh of the interpolation grid. In Section 4 we use our kriging interpolator to estimate and predict the intensity of Montagu's Harriers' nest locations in a region of France. In Section 5, we discuss the influence of the mesh and the rate and shape of unobserved areas on the statistical properties of our kriging interpolator from numerical results.

Linking up characteristics of two theories

About geostatistics

For any real valued random field Z (x), x ∈ R 2 , the first-order characteristic is the mean value function: E [Z (x)] = m(x) and the second-order characteristics are classically described in

) - 3
geostatistics [START_REF] Matheron | Traité de Géostatistique Appliquée: Mémoires du Bureau de Recherches Géologiques et Minières[END_REF][START_REF] Gabriel | Spatial Statistics[END_REF] by the (semi)-variogram, i.e. the mean squared difference at lag

r: γ (r) = 1 2 E  (Z(x) -Z (x + r)) 2 
. For a stationary and isotropic random field, we have

E [Z (x)] = m, γ (r) = σ 2 -Cov(Z(x), Z (x + h)), (1)
where σ 2 is the field variance and Cov(Z(x), Z (x + r)) is the auto-covariance of the random field.

We can interpolate the value Z (x o ) at the unsampled location x o by using the best linear unbiased predictor, so-called kriging interpolator:  Z (x o ) = µ T z, where z = {Z(x i )} i=1,...,n is the observation vector of the random field and µ is the n-vector of weights. In the case of ordinary kriging [START_REF] Cressie | Statistics for Spatial Data[END_REF], which will be of interest here since the mean value of the random field will be unknown, we have

µ = C -1 C o + 1 -1 T C -1 C o 1 T C -1 1 C -1 1, (2)
where

C =  Cov  Z (x i ), Z (x j )  i,j=1,...,n is the covariance matrix between the observations, C o =  Cov  Z (x i ), Z (x o )  i=1,.
..,n is the covariance vector between the observations and Z (x o ) and 1 is the n-vector of 1 (see e.g. [START_REF] Cressie | Statistics for Spatial Data[END_REF][START_REF] Wackernagel | Multivariate Geostatistics: An Introduction with Applications[END_REF].

About point processes

Let Φ be a stationary and isotropic point process defined in R 2 and B a Borel set centered at 0. Following the notations in [START_REF] Chiu | Stochastic Geometry and Its Applications[END_REF], a realization of Φ within a window S obs will be denoted by Φ S obs and the random counting measure for a Borel set B by Φ(B).

The first-and second-order characteristics of Φ are described through its intensity λ and the Ripley's K -function or the pair correlation function g:

λ = E [Φ(S obs )] ν(S obs ) , (3)
K * (r) = 1 λ E [Φ(b(0, r)) -1|0 ∈ Φ] ,
(4)

g(r) = 1 2π r
∂K * (r) ∂r ,

(5)

where ν(S obs ) is the area of S obs and b(0, r) is the disc centered at 0, with radius r. The intensity λ is thus the expected number of points per unit area, λK * (r) is the mean number of points in a circle of radius r centered at a typical point of the point process, whereas g(r) measures how K * changes with r. See for instance [START_REF] Chiu | Stochastic Geometry and Its Applications[END_REF] for a review about the theory of point processes.

Lemma 2.1. Let Φ be a point process with intensity λ and B, D two Borel sets. Then,

1. If ν(B), ν(D) → 0, then P [{Φ(B) = 1} ∩ {Φ(D) = 1}] = λν(B ∩ D) + λ 2  B×D g(x -y) dx dy + o (ν(B ∪ D)), 2. E  Φ 2 (B)  = λν(B) + λ 2  B×B g(x -y) dx dy, 3. Var(Φ(B)) = λν(B) + λ 2  B×B (g(x -y) -1) dx dy, 4. If B ∩ D = ∅, then Cov (Φ(B), Φ(D)) = λ 2  B×D (g(x -y) -1) dx dy.
The proof of 1. is given in Appendix A. Decompositions 2.-4. can be found in [START_REF] Chiu | Stochastic Geometry and Its Applications[END_REF].

Linking up

In our context, data are defined as informative point locations (the realization of the point process Φ) while the geostatistical calculations (kriging) need to be carried out over the values of a random field Z observed at several sampling locations, grid cell centers for example. Thus, we must regularize our process over a compact set [START_REF] Zhang | Scale in Spatial Information and Analysis[END_REF]. This consists in defining Z (x) by the count of the point process over the grid cell B centered at x ∈ R 2 i.e. Z (x) = Φ(B ⊕ x). Such a random field is of interest in our case as we want to estimate a non-constant intensity, classically defined by

lim ν(B)→0 E  Φ(B⊕x) ν(B)  .
From the first-and second-order moments defined in the previous sections, we can link up the characteristics of the point process Φ to the ones of the random field of point counts Z . Because of the stationary assumption it can also be related to the auto-covariance function (Eq. ( 1)), thus in the following we shall consider the latter.

Proposition 2.2. For the count random field defined by Φ(B ⊕ x), where B is a given Borel set, we have:

1. m = λν(B), 2. For B and D two regularization blocks, B = D + r, B D = B \ D, D B = D \ B, 2γ (r) = λ (ν(B D ) + ν(D B )) + λ 2  B D ×B D g(x -y) dx dy +  D B ×D B g(x -y) dx dy -2  B D ×D B g(x -y) dx dy  . 3. If B = D + r, then for ν(B) = ν(D) → 0 Cov (Φ(B), Φ(D)) ≈ λν(B)  I {B=D} + λν(B)  g(r) -1   . ( 6
)
The proof of Proposition 2.2 is straightforward from Lemma 2.1 and from the approximation

P [{Φ(B) = 1} ∩ {Φ(D) = 1}] ≈ λ 2 ν(B)ν(D)g(r) (see Appendix B).

Adapted kriging for point processes

In what follows, we consider that the stationary and isotropic point process, Φ, is driven by a stationary random field, U. We want to interpolate the local intensity of Φ, λ(x|U), x ∈ S ⊂ R 2 , i.e. its intensity conditional to the random field U, from its realization within an observation window S obs . Hence, we use the relation between point processes and geostatistics (Section 2.3) and approximate the point process by the counting process within a grid of elementary cell B.

For sake of clarity, in the following we denote by S the region of interest so that S unobs define the complementary of S obs within S. We consider a regular grid superimposed on S with a square-mesh.

We denote by B an elementary square centered at 0, B i = x i ⊕ B the elementary square centered at x i such that B i ∩ B j = ∅, and n (resp. n obs ) the number of grid cell centers lying in S (resp. S obs ).

Defining the interpolator

According to the classical geostatistical method defined in Section 2, the kriging interpolator of the local intensity at x o ∈ S, λ(x o |U), should be written as

µ T  λ(x 1 |U), . . . , λ(x n obs |U)  ,
for some well-chosen kriging weights µ where x i , i = 1, . . . , n obs correspond to data sample locations, i.e. here to the cell centers of S obs . Note that in our case we cannot observe the local intensity at x i , thus we can estimate it by

Φ(B i ) ν(B)
. Furthermore because of the cell-point relation, we cannot have an exact interpolation of the local intensity.

Proposition 3.1. Given the elementary square B, the interpolator at x o defined by

 λ(x o |U) =  x i ∈S obs µ i Φ(B i ) ν(B) , ( 7
)

) - 5 where µ = (µ 1 , . . . , µ n obs ) = C -1 C o + 1-1 T C -1 C o 1 T C -1 1 C -1 1, is the best linear unbiased predictor (BLUP) of Φ(B o ) ν(B)
and the asymptotically BLUP of λ(x o |U).

The weights depend on

• the covariance matrix C = λν(B)I + λ 2 ν 2 (B)(G -1), where G =  g ij  i,j=1,...,n obs , with g ij = 1 ν 2 (B)
 B×B g(x i -x j +u-v) du dv, and I is the n obs ×n obs -identity matrix,

• the covariance vector

C o = λν(B)I x o + λ 2 ν 2 (B)(G o -1),
where G o = {g io } i=1,...,n obs , and I x o is the n obs -vector with zero values and one term equals to one where x o = x i (which only happens in estimation).

When the distance between x o and x i , for all i, is larger than the range of interaction, the predicted value tends to λ if x o ̸ ∈ S obs and the estimated value tends to 0 if x o ∈ S obs .

Proof. At the scale of B, the kriging weights such that  λ(

x o |U) is a BLUP of Φ(B o )
ν(B) are given by the ordinary kriging equations [START_REF] Cressie | Statistics for Spatial Data[END_REF][START_REF] Chilès | Geostatistics: Modeling Spatial Uncertainty[END_REF].

At a finer scale we have that

E  Φ(B⊕x) ν(B)
 tends to λ(x|U) when ν(B) tends to 0 (as λ(x|U) is assumed to be continuous). Thus we propose to interpolate λ(x o |U) by using  λ(

x o |U) =  x i ∈S obs µ i Φ(B i ) ν(B) , with the constraint  n obs i=1 µ i = 1. Minimizing the error variance Var   λ(x o |U) -λ(x o |U) 
under this constraint and using Eq. ( 6) lead to the following kriging weights:

µ = ν(B)C -1  C o + 1 -ν(B)1 T C -1  C o 1 T C -1 1 C -1 1, where  C o =  Cov  Φ(B i ), λ  x o |U   i=1,...,n obs . To get Cov  Φ(B i ), λ  x o |U  
, note that, for λ Φ B (x|U) denoting the local intensity of the point process Φ given its realization in B:

E  Φ(B o )|Φ B i  =  B o λ Φ B i (x|U) dx = λ Φ B i (x o |U)ν(B) + o(ν(B)) and λ Φ B i (x o |U) = E  λ Φ S obs (x o |U)|Φ B i 
. Thus, we have

• for x i ̸ = x o , E [Φ(B o )Φ(B i )] = E  Φ(B i )E  Φ(B o )|Φ B i  = E  Φ(B i )ν(B)λ Φ B i (x o |U)  = E  Φ(B i )ν(B)E  λ Φ S obs (x o |U)  |Φ B i  = ν(B)E  E  Φ(B i )λ Φ S obs (x o |U)  |Φ B i  = ν(B)E  Φ(B i )λ Φ S obs (x o |U)  which leads to E [Φ(B i )λ(x o |U)] = 1 ν(B) E [Φ(B o )Φ(B i )]
as in our case the local intensity is conditional to the realization of the process Φ in S obs .

• for

x i = x o , E  Φ 2 (B o )  = E  E  Φ 2 (B o )|Φ S obs  = E  Φ(B o )E  Φ(B o )|Φ S obs  = E [Φ(B o )ν(B)λ(x o |U)] = ν(B)E [Φ(B o )λ(x o |U)] which leads to E [Φ(B o )λ(x o |U)] = 1 ν(B) E  Φ 2 (B o )  . E. Gabriel et al. / Spatial Statistics ( ) - Consequently,  C o = 1 ν(B) C o and we get µ = C -1 C o + 1 -1 T C -1 C o 1 T C -1 1 C -1 1.
Interpolating Φ(B o )/ν(B) or λ(x o |U) leads to the same kriging weights.

Finally,

E   λ(x o |U)  = E   x i ∈S obs µ i Φ(B i ) ν(B)  = E  Φ(B o ) ν(B)  ----→ ν(B)→0 λ(x o |U)
shows that  λ(x o |U) is an asymptotically unbiased predictor of λ(x o |U).

Properties of the interpolator

In order to develop the variance of the kriging interpolator, we use the following Neumann series (see e.g. [START_REF] Petersen | The Matrix Cookbook[END_REF] to invert the covariance matrix C , which holds when λν(B)

tends to 0:

C -1 = 1 λν(B) [I + λν(B)J λ ] , (8)
where a generic element of the matrix J λ is given by

J λ [i, j] = ∞  k=1 (-1) k λ k-1  g(x i , x l 1 ) -1   g(x l k-1 , x j ) -1  ×  S k-1 obs k-2  m=1 (g(x l m , x l m+1 ) -1) dx l 1 . . . dx l k-1 . Proposition 3.2. When λν(B) tends to 0, the variance of  λ(x o |U) is Var   λ(x o |U)  = λ ν(B) + 2λ 2 I T x o J λ I x o + 2λ 3 ν(B)I T x o J λ (G o -1) + λ 3 ν 2 (B)(G o -1) T (G o -1) + λ 4 ν 3 (B)(G o -1) T J λ (G o -1) + 1 -  1 + λν(B)1 T J λ I x o + λν(B)1 T (G o -1) + λ 2 ν 2 (B)1 T J λ (G o -1)  2 ν(S obs ) λ + ν 2 (B)1 T J λ 1 . ( 9
)
In estimation, i.e. when x o ∈ S obs , we get the following approximation,

Var   λ(x o |U)  ≈ λ ν(B) . ( 10
)
In prediction, i.e. when x o ̸ ∈ S obs , the variance reduces to

Var   λ(x o |U)  = λ 3 ν 2 (B)(G o -1) T (G o -1) + λ 4 ν 3 (B)(G o -1) T J λ (G o -1) + 1 -  λν(B)1 T (G o -1) + λ 2 ν 2 (B)1 T J λ (G o -1)  2 ν(S obs ) λ + ν 2 (B)1 T J λ 1 . ( 11
)
Proof. From e.g. [START_REF] Cressie | Statistics for Spatial Data[END_REF], the variance of the predictor is given by

Var   λ(x o |U)  = Var   x i ∈S obs µ i Φ(B i ) ν(B)  = 1 ν 2 (B) µ T C -1 µ = 1 ν 2 (B)  C T o C -1 C o + 1 -(1 T C -1 C o ) 2 1 T C -1 1  . ) - 7
• When estimating the local intensity, i.e. for x o lying in the observation window, we have

C o = λν(B)I x o + λ 2 ν 2 (B)(G o -1).
Thus, from Eq. ( 8):

C T o C -1 C o = λν(B)  1 + λν(B)  J λ (x o , x o ) + 2I T x o (G o -1)  λ 2 ν 2 (B)  2I T x o J λ (G o -1) + (G o -1) T (G o -1)  λ 3 ν 3 (B)(G o -1) T J λ (G o -1)  , where J λ (y, z) =  ∞ k=1 (-1) k λ k-1  S k-1 obs (g(y, x l 1 ) -1)  k-2 m=1 (g(x l m , x l m+1 ) -1)(g(x l k-1 , z) - 1) dx l 1 . . . dx l k-1 , 1 T C -1 C o = 1 + λν(B)  1 T J λ I x o + 1 T (G o -1) + λν(B)1 T J λ (G o -1)  ,
and

1 T C -1 1 = 1 λν(B)  n obs + λν(B)1 T J λ 1  = ν(S obs ) λν 2 (B) + 1 T J λ 1.
Then, if ν(B) is very small,

C T o C -1 C o ν 2 (B) varies in λ ν(B) and 1-(1 T C -1 C o ) 2 ν 2 (B)1 T C -1 1 in λ ν(S obs )
. Thus, we get Eq. ( 10).

• When predicting the local intensity, i.e. for x o outside the observation window, we have

C o = λ 2 ν 2 (B)(G o -1). Thus, from C T o C -1 C o = λ 3 ν 3 (B)(G o -1) T (G o -1) + λ 4 ν 4 (B)(G o -1) T J λ (G o -1)
and

1 T C -1 C o = λν(B)1 T (G o -1) + λ 2 ν 2 (B)1 T J λ (G o -1)
we get Eq. (11).

Defining an optimal mesh size

The Integrated Mean Squared Error of  λ(x|U) is defined as

IMSE   λ(x|U)  =  S   λ(x|U) -E[  λ(x|U)]  2 + Var   λ(x|U)   dx.
When estimating the local intensity, this leads to the following approximation:

IMSE   λ(x|U)  ≈ √ ν(B) 12  S obs ∥∇λ(x|U)∥ 2 dx + λν(S obs ) ν(B) . ( 12
)
We propose to find the optimal mesh of the estimation grid by minimizing IMSE   λ(x|U) 

(see Appendix C), and we get:

ν opt (B) =  12λν(S obs )  S obs ∥∇λ(x|U)∥ 2 dx . ( 13
)
Note that because the optimal mesh depends on the inverse of squared L 2 -norm of the gradient of the local intensity, it decreases for clustered point patterns. Conversely, it increases for regular point patterns.

In practice the optimal mesh can be approximated by estimating the gradient of the intensity over a fine grid (see Appendix D). When predicting the local intensity, the smaller the mesh, the better. Computation time is the only limit.

Real case study

In this section we estimate and predict the intensity of Montagu's Harriers' nest locations in the Zone Atelier ''Plaine & Val de Sèvre'' 1 (Fig. 1), a NATURA2000 site in France of 450 km 2 , designated for its remarkable diversity of bird species. Dots in Fig. 1 represent the exhaustive collection of Montagu's Harriers' nest locations. The area in the center of the Zone Atelier delineates the administrative boundaries of the commune Saint-Martin-de-Bernegoue, which will be used for prediction.

Estimation of the pair correlation function

The pair correlation is estimated as defined in [START_REF] Stoyan | Fractals, Random Shapes and Point Fields: Methods of Geometrical Statistics[END_REF]:

 g(r) = 1 2π r  ξ ∈Φ S obs ̸ =  ζ ∈Φ S obs k h (r -∥ξ -ζ ∥) prop  S obs ∩ S obs,ξ -ζ 
where k h is the Epanechnikov kernel with bandwidth h, the optimal Stoyan's bandwidth equals to 0.15/ √ Φ(S obs )/ν(S obs ) and prop  S obs ∩ S obs,ξ -ζ  is the proportion of translations of (ξ , ζ ) which have both ξ and ζ inside S obs . Fig. 2(a) shows the pair correlation function estimated from either all data point locations (solid line) or only the ones outside the boundaries of Saint-Martin-de-Bernegoue (dashed line). These estimates are characteristic of a Thomas cluster process with an infinite range of correlation, see [START_REF] Illian | Statistical Analysis and Modelling of Spatial Point Patterns[END_REF].

Intensity estimation

For our kriging estimator, the optimal mesh is obtained by minimizing the IMSE. Usual nonparametric estimation methods also require to preliminary set the smoothing parameter and this parameter is chosen as an optimal value minimizing a specific criterion (typically mean square error, integrated bias, asymptotic mean square error). In our case, we have an explicit formula of the optimal mesh (Eq. ( 13)), which depends on the unknown terms λ and λ(x|U). If λ = Φ(S obs )/ν(S obs ) appears to be a natural candidate to estimate λ, the challenging goal is to estimate  S obs ∥∇λ(x|U)∥ 2 dx.

Based on simulation experiments (Appendix D), we consider a Gaussian kernel [START_REF] Silverman | Density Estimation for Statistics and Data Analysis[END_REF], with a bandwidth minimizing the mean-square error criterion defined by [START_REF] Diggle | A kernel method for smoothing point process data[END_REF], to get a good 1 http://www.za.plainevalsevre.cnrs.fr/. the difference between our estimation and the one obtained by Gaussian kernel smoothing at the same grid resolution. Our kriging interpolator gives higher values of the local intensity (in blue in Fig. 2(c)) close to aggregated observation points than the kernel estimator, while the maximum value may be higher for the later. This illustrates that our method may be particularly relevant for point patterns strongly aggregated at a small scale.

Intensity prediction

In order to apply our kriging predictor to the real dataset, we consider an unobserved window S unobs defined by the administrative boundaries of the commune Saint-Martin-de-Bernegoue in the center of the 'Zone Atelier' (Fig. 1). Thus, we remove the points in this area (red dots in Fig. 3) and use the remaining nest locations (blue dots in Fig. 3) to predict the local intensity within S unobs . We consider a grid of size 100 × 100 over S to make the prediction. The estimated pair correlation function is plotted in Fig. 2(a) (dashed line). The result, zoomed in Fig. 3, shows that the kriging predictor is able to reproduce the second-order structure of the point process. In particular, it reproduces clusters as soon as there are points close enough to the boundary of the unobserved area. This will be further illustrated and discussed in the next section. Note that at distances greater than the range of the pair correlation function, the method can only provide a constant intensity estimate.

Illustrative simulation experiments

Objectives

Now, we focus on the kriging predictor and explore its accuracy through simulation experiments, varying rate and shape of the observation window S obs . To measure the quality of prediction, we compute the mean bias (MB) and the mean square error of prediction (MSEP):

MB = 1 n obs  x∈Φ S obs 1 n sim n sim  k=1   λ k (x|U) -λ k (x|U)  , MSEP = 1 n obs  x∈U 1 n sim n sim  k=1   λ k (x|U) -λ k (x|U)  2 ,
where λ k and  λ k correspond respectively to the intensity and its predictor on the kth simulation and n sim is the number of simulations. We also compute the coefficient of determination R 2 of the regression between the predicted values of the local intensity and the theoretical ones.

Experimental design

Throughout our experimental study, in order to simplify the analysis of the two parameters of interest (rate and shape of S obs ), we decide to simulate all point patterns from a single spatial point process model. We consider in the sequel a Thomas process, for which we have explicit formulas of the intensity, pair correlation function and others characteristics (see Illian et al., 2008, p. 377):

λ(x|U) =  ξ ∈Φ S obs µ 2π σ 2 exp  - ∥x -ξ ∥ 2 2σ 2



, for all x ∈ S, ( 14) Several windows of interest S i , with i = 1, . . . , 24, are considered (Fig. 4), corresponding to different observation rates (83%, 66%, 50%, 33% and 17%). The unobserved windows S unobs are defined by the union of bands, with varying width (in gray in Fig. 4).

g(r) = 1 + 1 4π κσ 2 exp  - r 2 4σ 2  , for r ≥ 0.
Because the weights in our kriging interpolator depend on the pair correlation function, in our experiment we compare results arising from the theoretical pair correlation function, and from its estimate defined in Section 4.1. In order to estimate the pair correlation function from similar numbers of points in each window S i , we first simulated point patterns within a larger window (the initial one extended on the right side), so that the area of observation zones equals to one. The pair correlation function is then estimated from this first pattern and the prediction is made on its restriction to the initial unit square.

Results

The mean bias and the mean square error of prediction are presented in Fig. 5, with theoretical values of the pair correlation function (solid lines) or an estimate (dotted line). It shows that the mean bias has no effect on the MSEP. With theoretical values of the pair correlation function, the mean bias is close to zero whatever the width of the bands defining S unobs , which numerically reveals the unbiasedness statistical property of our predictor. When the pair correlation function is estimated, λ k is under-estimated and the discrepancy is higher when the observation rate decreases than when the width of the unobserved bands increases.

At a given observation rate, the MSEP increases when the width of the unobserved bands increases. Indeed, the geometry of our windows of interest implies that wider the unobserved bands, less numerous they are. Consequently, for some simulated patterns, cluster points can completely fall within an unobserved band, what damages the quality of prediction. At a given value of the unobserved band, we obviously see a slight increase of the MSEP when the observation rate decreases.

E. Gabriel et al. / Spatial Statistics (

) -Fig. 5. Mean Bias (left) and Mean Square Error of Prediction (right) of our kriging predictor related to the width of the unobserved bands, and to the observation rates: 83% in cyan, 66% in blue, 50% in green, 33% in red and 17% in black. The lines correspond to a linear approximation of the MSEP values when g is known (solid lines) or estimated (dotted lines). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

We first illustrate the influence of the estimation of the pair correlation function onto the accuracy of prediction on a single simulation. The simulated pattern, the associated theoretical intensity and the observation window, with a rate of 50% of observed areas, are represented in Fig. 6(a)-(c) respectively. The theoretical (dotted line) and estimated (solid line) pair correlation function are given in Fig. 6(d).

Fig. 6(e) and (f) illustrate the theoretical local intensity in S obs and the prediction in S unobs on a 96 × 96 grid, using the true (e) and the estimated (f) pair correlation function. In the first case, the prediction is relatively smooth and gives accurate results. In the second case, the prediction is more noisy, but recover the same blocks with high intensity values. In both cases, the method correctly predict the clusters when there are observations close to the unobserved bands. That is the case for all clusters located at the right hand side of the vertical line x = 0.25. When the full cluster falls in the unobserved band, as the ones located at the left hand side of the vertical line x = 0.25 the method fails in predicting the cluster. We plotted (Fig. 7) the boxplot of the coefficient of determination resulting from 100 simple linear regressions between the predicted values of the local intensity and the theoretical ones, for different grid size (24 × 24, 48 × 48 and 96 × 96), when the pair correlation function is estimated (Fig. 7(b))

or not (Fig. 7(a)). These results are related to an observation rate of 50%, according to the window configuration highlighted Fig. 4. We obviously see that the goodness of prediction increases when the grid resolution increases and when the pair correlation function is known. We considered a 96 × 96 grid as it is a trade off between computation times and a small mesh, allowing a good description of the intensity variations due to clusters. We obtained, in the worst case where the pair correlation function is estimated, that the coefficient of determination R 2 is around 0.8 (median).

Discussion

Our kriging method introduced to estimate and/or predict the local intensity of a stationary and isotropic point process has a large number of advantages, particularly in prediction. Taking into account the spatial structure of the point pattern allows to perform the intensity estimation for point processes highly aggregated at fine scale. In the prediction framework, our kriging method is innovative for the interpolation of the local intensity and presents good statistical properties (unbiasedness, low variance. . . ) when the pair correlation function is known. When it is estimated, the quality of our interpolator is slightly reduced but our results can be improved by better taking into account the double estimation of the pair correlation function and the local intensity on the same pattern.

This prediction method is less time consuming than the reconstruction methods and appears a promising way in prediction of intensity of a spatial point pattern. Note that existing prediction methods are constrained within a class of point processes (van Lieshout and Baddeley, 2001, in particular Cox processes [START_REF] Monestiez | Geostatistical modelling of spatial distribution of balaenoptera physalus in the northwestern mediterranean sea from sparse count data and heterogeneous observation efforts[END_REF][START_REF] Diggle | Model-Based Geostatistics[END_REF][START_REF] Bellier | Reducing the uncertainty of wildlife population abundance: model-based versus design-based estimates[END_REF][START_REF] Diggle | Spatial and spatio-temporal log-gaussian Cox processes: Extending the geostatistical paradigm[END_REF], making any comparison with our method very restrictive relatively to its broad scope of applications. That is for instance the case of any point process obtained by a weak dependent process (e.g. Thomas, Markov) with a parameter driven by a stationary random field at a larger scale (e.g. Cox), but not only.

Relaxing the stationary assumption implies to make further assumptions. The formalism should be quite similar to the one of this paper, but with some confounding effects as the ones observed when using the same point pattern to estimate both a spatially varying intensity and second-order characteristics (Diggle et al., 2007;[START_REF] Gabriel | Estimating second-order characteristics of inhomogeneous spatio-temporal point processes: influence of edge correction methods and intensity estimates[END_REF]. For instance, if we consider a non stationary Cox process, we cannot disentangle the first-order non-stationarity to the second-order non-stationarity. One could thus allocate the effects at different scales.

E. Gabriel et al. / Spatial Statistics (

) -

In our simulations and application, we used the R function solve, based on the LU factorization, to compute the inverse of the covariance matrix C . This matrix is of dimension the square of the number of cells of the grid superimposed on the observation window. Thus, it can quickly become heavy to inverse. In such cases, estimating the matrix C -1 using Eq. ( 8) would be somewhat cumbersome. Thus, we propose instead to inverse the covariance matrix numerically. Several approximations could be used, depending mainly on the width of B with respect to λ and the curvature of the pair correlation function:

1. if the diameter of B is large, the covariance between two tiles at distance r is equal to λν(B)I r=o + λ 2  B×D (g(r + xy) -1) dx dy. It can be approximated numerically by computing for example the integral on a fine grid. Then the finer the grid, the smaller the difference between the exact values and the approximations, but the computing time cost can become prohibitive.

2. when the diameter of B becomes small, the integral can be approximated by ν(B) 2 g(r) so that the covariance is approximated by λν(B)I r=o + ν(B) 2 λ 2 g(r), 3. when the diameter becomes very small, C may be approximated by λν(B)I, a situation seldom met in practice, since it needs a tile B small enough to neglect point dependence.

Approximation (2) will thus be the most reasonable one, needing only a B small enough to consider that g(r + x) is almost constant for x ∈ B, but avoiding too small B leading to large matrix inversion time.

Our estimation is roughly pixellated compared to kernel methods, but it does not oversmooth the intensity of highly aggregated point processes. We could take the benefit of the two approaches to get smoother estimations. Our on-going work consists in regularizing the counting process by a kernel and in defining a kriging estimator for the related random field. Our optimal grid could then be used to define an optimal bandwidth, thus eliminating the Poisson aspect of classical kernels.

Our method provides good predictions in areas at small distances of data locations. From the definition of the kriging predictor, at distances larger than the range of interaction, it only provide a constant mean value. To improve it and make it more relevant in practice, we could consider further information provided by covariates. From our application point of view, wheat field mapping could be of interest as Montagu's Harriers nest in there. From a methodological point of view, including covariates would imply that we should either consider external drift kriging (or any other universal kriging) rather than ordinary kriging; or spatial regression.

Finally, our kriging predictor depends on the count data in the grid cells, B i , and not on exact data locations in B i . Thus we can further consider count datasets, as it is often the case in biodiversity measures, e.g. plant species abundance. The exact position of each plant is rarely given, but we know its abundance per small unit areas. So, once the pair correlation function is estimated from the point data subset, one can apply our method to interpolate the intensity. much better values of the optimal grid size. While λ(x|U) may be roughly estimated, the integral of its gradient is sufficiently well approximated to obtain good results. sizes (in number of pixels) obtained from different size of the estimation grid: 100 × 100 (solid line), 200 × 200 (dashed line) and 500 × 500 (dotted line). The theoretical grid size is in black and the ones derived from the kernel smoothing and the k-nearest neighbors method are in red and green respectively. This figure is related to the Thomas process with parameters κ = 10, µ = 50, and we get similar results from the other set of parameters. It appears that the k-nearest neighbors based method is very sensitive to the size of the estimation grid and tends to over-estimate the optimal grid size. The kernel based method under-estimates the optimal grid size when the estimation grid is not fine enough and when the scale of clustering in very small.

From this simulation study, we recommend the kernel based method on a 200 × 200 grid to first estimate ∇λ(x|U) and then compute the optimal mesh or equivalently the optimal grid size.

Fig. 1 .

 1 Fig. 1. Montagu's Harriers nest locations in the Zone Atelier ''Plaine & Val de Sèvre''.

Fig. 2 .

 2 Fig. 2. (a) Estimation of the pair correlation function from all points (solid line) and without those lying in Saint-Martin-de Bernegoue (dashed line). Estimation of Montagu's Harriers nest locations using our kriging interpolator on the optimal grid (b) and using a Gaussian kernel on a 128 × 128 grid (d). (c) Difference between our estimator and the Gaussian kernel smoothing at the same resolution: dark gray (resp. blue) indicates higher values of the Gaussian kernel (resp. our kriging) estimates. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) approximation of the gradient of λ(x|U) on a 200 × 200 grid. This methodology applied to the real dataset leads to a value of ν opt (B) equals to 23.19 hectares, which corresponds to a grid of 64 × 53 cells. Fig. 2(b) shows the kriging estimate on the optimal grid. Fig. 2(d) represents an estimate obtained by a Gaussian kernel, with a bandwidth selected as previously mentioned, on a 128 × 128 grid (default of the spatstat function 'density.ppp', Baddeley and Turner, 2005). Fig. 2(c) illustrates

Fig. 3 .

 3 Fig. 3. Prediction within the commune of Saint-Martin-de-Bernegoue (data location is in blue; the red ones are removed for the prediction). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .

 4 Fig. 4. Windows of interest with S obs (resp. S unobs ) the union of white (resp. gray) bands.

Fig. 6 .Fig. 7 .

 67 Fig. 6. (a) Simulated pattern from a Thomas process (parent points in red and children points in black). (b) Theoretical intensity from the simulated pattern. (c) Observed window (light gray) and unobserved window (red) with an observation rate of 50%. (d) Theoretical (dotted line) and estimated (solid line) pair correlation function g. (e) and (f) Theoretical intensity in the observed window and predicted local intensity in the unobserved window obtained with the true pair correlation function (e) and an estimate (f). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

  =y∈Φ S I B (x)I D (y)

Fig. D. 8 .

 8 Fig. D.8.Mean of optimal grid sizes computed from 100 Thomas processes with parameters κ = 10, µ = 50 and from the k-nearest neighbor based method (green), the kernel based method (red), both evaluated on a 100 × 100 grid (solid line), a 200 × 200 grid (dashed line) and a 500 × 500 grid (dotted line), compared to the theoretical value (black). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

  Fig. D.8 shows the optimal grid

(1)

(2) follows from Lemma 2.1:



(3) follows from the approximation

Appendix C. Proof of Eqs. (12) and (13)

Let B a square centered at 0 of area ν(B) = b 2 . We denote by ∇λ(x) the gradient vector

By the following Taylor expansion around the origin

we obtain that:

Appendix D. Optimal mesh in practice

The optimal mesh of the estimation grid, ν opt (B), depends on the unknown terms λ and ∇λ(x|U). In this section, we compare different methods which could be used to compute ν opt (B) in practice.

We simulated point patterns in the unit square from Thomas point processes with different set of parameters: (κ, µ) ∈ {(10, 50); (22,23); (50, 10)} and σ ∈ {0.001, 0.0025, 0.005, 0.01, 0.025, 0.05}. Then, ν opt (B) is computed as follows. First we estimate/compute the intensity on a N × N grid. Second, we deduce its gradient from the rate of change between the estimated intensity and its one-cell translated value. Third, we compute ν opt (B) and its related grid size (in number of pixels).

We used different methods to estimate the intensity. The counting method consists in estimating the local intensity in each pixel

The global kernel smoothing method is based on a Gaussian kernel estimator with global bandwidth and without border correction,

where w is the density function of a standard normal distribution. The k-nearest neighbors method is an adaptive nonparametric estimation so that  λ(x|U) = 1/(π d k (x) 2 ), with d k (x) the distance of x to its k-nearest neighbor. Note that we also compute the theoretical value of the intensity from Eq. ( 14).

We considered three N × N-grids, with N ∈ {100, 200, 500}. We compared the distributions of the optimal grid sizes obtained from 100 realizations of each process and from the different methods, to the theoretical ones. We select the method which provides the more accurate results and the less sensitivity to the different scenario. If the counting method works well when the pattern is strongly aggregated at very small scale (σ ≤ 0.005), it requires a fine grid and is inaccurate for other scales of clustering. Thus in the following it will be no longer considered. The other methods provide globally


domain: math info

text:
Bertrand Cloez
email: bertrand.cloez@supagro-inra.fr

Marie-Noémie Thai

FLEMING-VIOT PROCESSES : TWO EXPLICIT EXAMPLES

Keywords: AMS 2000 Mathematical Subject Classification: 60K35, 60B10, 37A25 Fleming-Viot process -quasi-stationary distributions -coupling -Wasserstein distance -chaos propagation -commutation relation

The purpose of this paper is to extend the investigation of the Fleming-Viot process in discrete space started in a previous work to two specific examples. The first one corresponds to a random walk on the complete graph. Due to its geometry, we establish several explicit and optimal formulas for the Fleming-Viot process (invariant distribution, correlations, spectral gap).

The second example corresponds to a Markov chain in a two state space. In this case, the study of the Fleming-Viot particle system is reduced to the study of birth and death process with quadratic rates.

INTRODUCTION

In discrete space, the Fleming-Viot particle system has been studied by many authors [START_REF] Asselah | Quasistationary distributions and Fleming-Viot processes in finite spaces[END_REF][START_REF] Asselah | Fleming-Viot selects the minimal quasi-stationary distribution: The Galton-Watson case[END_REF][START_REF] Asselah | A note on the rightmost particle in a Fleming-Viot process[END_REF][START_REF] Ferrari | Quasi stationary distributions and Fleming-Viot processes in countable spaces[END_REF][START_REF] Marić | Fleming-Viot particle system driven by a random walk on $\mathbb{N}[END_REF]. Such a system is a mean field particle system described in the following way: we consider N copies of an absorbed Markov chain and, instead of being absorbed, one chain jumps randomly on the state of another one. It is well known that, when the number of copies tends to infinity, the empirical measure converges to the law of the initial chain conditioned not to be absorbed, see for instance [START_REF] Cloez | Quantitative results for the Fleming-Viot particle system and quasi-stationary distributions in discrete space[END_REF][START_REF] Del Moral | A Moran particle system approximation of Feynman-Kac formulae[END_REF][START_REF] Groisman | Simulation of quasi-stationary distributions on countable spaces[END_REF][START_REF] Villemonais | General approximation method for the distribution of Markov processes conditioned not to be killed[END_REF]. Convergence to equilibrium as time goes to infinity is less known. In [START_REF] Rousset | On the control of an interacting particle estimation of schrödinger ground states[END_REF][START_REF] Villemonais | Interacting particle systems and Yaglom limit approximation of diffusions with unbounded drift[END_REF], this question is addressed for some models. Nevertheless, to our knowledge, there are few results on the expression of the invariant distribution or on the explicit rates of convergence. This paper is concerned with studying two specific models for which the invariant distribution is explicit and extend the investigation started by Cloez and Thai [START_REF] Cloez | Quantitative results for the Fleming-Viot particle system and quasi-stationary distributions in discrete space[END_REF]. Let Q = (Q i,j ; i, j ∈ F * ∪ {0}) be the transition rates matrix of an irreducible and positive recurrent continuous time Markov process (X t ) t≥0 on a countable state space F = F * ∪ {0}. We think of 0 as an absorbing state. Let µ be the initial law of (X t ) t≥0 and let µT t be its law at time t conditioned on non absorption up to time t. That is defined, for all non-negative function f on F * , by µT t f = µP t f µP t 1 {0} c = y∈F * P t f (y)µ(y)

y∈F * P t 1 {0} c (y)µ(y)
, where (P t ) t≥0 is the semigroup associated with the transition matrix Q and we use the convention f (0) = 0. For every x ∈ F * , k ∈ F * and non-negative function f on F * , we also set T t f (x) = δ x T t f and µT t (k) = µT t 1 {k} , ∀t ≥ 0.

A quasi-stationary distribution (QSD) for Q is a probability measure ν qs on F * satisfying, for every t ≥ 0, ν qs T t = ν qs .

The particle system we are focusing on was initially introduced in [START_REF] Del Moral | On the stability of measure valued processes with applications to filtering[END_REF][START_REF] Del Moral | A Moran particle system approximation of Feynman-Kac formulae[END_REF] for approximating the conditioned semigroup (T t ) t≥0 and the QSD ν qs . It is convenient to think of particles as being indistinguishable, and to consider the occupation number η with, for k ∈ F * , η(k) = η (N ) (k) representing the number of particles at site k. The configuration (η t ) t≥0 is a Markov process with state space E = E (N ) defined by

E = η : F * → N | i∈F * η(i) = N .
Applying its generator to a bounded function f gives

Lf (η) = L (N ) f (η) = i∈F * η(i)   j∈F * (f (T i→j η) -f (η)) Q i,j + Q i,0 η(j) N -1   , (1)
for every η ∈ E, where, if η(i) = 0, the configuration T i→j η is defined by T i→j η(i) = η(i) -1, T i→j η(j) = η(j) + 1, and

T i→j η(k) = η(k) k / ∈ {i, j}.
The present paper is a continuation of [START_REF] Cloez | Quantitative results for the Fleming-Viot particle system and quasi-stationary distributions in discrete space[END_REF] in which the following limits are studied and quantified:

µ N t t→∞ N →∞ μT t t→∞ ! ! µ N ∞ N →∞ } } ν qs
where µ N is the associated empirical distribution of the particle system defined, for η ∈ E, by

µ N t = 1 N k∈F * η(k)δ {k} .
For countable space F , the ergodicity of the Fleming-Viot process is not guaranteed. In [START_REF] Cloez | Quantitative results for the Fleming-Viot particle system and quasi-stationary distributions in discrete space[END_REF],

Cloez and Thai show that under some conditions, the particle system converges exponentially fast to equilibrium for a suitable Wasserstein coupling distance. Let us recall the different distances given by the autors. For η, η ∈ E, let d be the distance defined by

d(η, η ) = 1 2 j∈F |η(j) -η (j)|, (2)
and for any two probability measures µ and µ on E, let W d (µ, µ ) be the Wasserstein coupling distance between these two laws defined by

W d (µ, µ ) = inf X∼µ X∼µ E d(X, X) , (3)
where the infimum runs over all the couples of random variables with marginal laws µ and µ .

Theorem 1.1 (Theorem 1.1 of [8]). Let λ = inf i,i ∈F *   Q i,i + Q i ,i + j =i,i Q i,j ∧ Q i ,j   and for i ∈ F * , p 0 (i) = Q i,0 . If ρ = λ -(sup(p 0 ) -inf(p 0 )
) then for any processes (η t ) t>0 and (η t ) t>0 generated by [START_REF] Ané | Sur les inégalités de Sobolev logarithmiques, volume 10 of Panoramas et Synthèses [Panoramas and Syntheses[END_REF], and for any t ≥ 0, we have

W d (Law(η t ), Law(η t )) ≤ e -ρt W d (Law(η 0 ), Law(η 0 )).
In particular, if ρ > 0 then there exists a unique invariant distribution ν N satisfying for every t ≥ 0,

W d (Law(η t ), ν N ) ≤ e -ρt W d (Law(η 0 ), ν N ).
This theorem gives the existence and uniqueness of the invariant distribution, but this one is not explicit. Actually, there are few models for which an explicit formula of the invariant distribution is given. This observation is one of the motivations of the present paper. So, it is interesting to consider and study a model satisfying the previous point : the example of random walk on a complete graph. An interesting point of the complete graph approach is that it permits to reduce the difficulties of the Fleming-Viot to the interaction. Due to its simple geometry, several explicit formulas are obtained such as the invariant distribution, the correlations and the spectral gap. It seems to be new in the context of Fleming-Viot particle systems. A second model for which the invariant distribution is explicit is the two point case, the study of the particle system is then reduced to a birth and death process with quadratic rates. The bound obtained in Theorem 1.1 is not optimal. Nevertheless, the coupling introduced in [START_REF] Cloez | Quantitative results for the Fleming-Viot particle system and quasi-stationary distributions in discrete space[END_REF] in order to prove Theorem 1.1, permits us to obtain the spectral gap as rate of convergence. Moreover, we show that the spectral gap of the Fleming-Viot process is always bounded from below by a positive constant not depending on the number of particles.

The remainder of the paper is as follows. Section 2 is dedicated to the study of random walk on the complete graph and Section 3 to that of the two point case.

COMPLETE GRAPH DYNAMICS

In all this section, we study the example of a random walk on the complete graph. Let us fix K ∈ N * , p > 0 and N ∈ N * , the dynamics of this example is as follows: we consider a model with N particles and K +1 vertices 0, 1, . . . , K. The N particles move on the K vertices 1, . . . , K uniformly at random and jump to 0 with rate p. When a particle reaches the node 0, it jumps instantaneously over another particle chosen uniformly at random. This particle system corresponds to the model previously cited with parameters

Q i,j = 1 K , ∀i, j ∈ F * = {1, . . . , K}, i = j and Q i,0 = p, ∀i ∈ F * .
The generator of the associated Fleming-Viot process is then given by

Lf (η) = K i=1 η(i)   K j=1 (f (T i→j η) -f (η)) 1 K + p η(j) N -1   , (4)
for every function f and η ∈ E. A process generated by ( 4) is an instance of inclusion processes studied in [START_REF] Giardinà | Correlation inequalities for interacting particle systems with duality[END_REF][START_REF] Grosskinsky | Condensation in the inclusion process and related models[END_REF][START_REF] Grosskinsky | Dynamics of condensation in the symmetric inclusion process[END_REF]. It is then related to models of heat conduction. One main point of [START_REF] Giardinà | Correlation inequalities for interacting particle systems with duality[END_REF][START_REF] Grosskinsky | Condensation in the inclusion process and related models[END_REF] is a criterion ensuring the existence and reversibility of an invariant distribution for the inclusion processes. In particular, they give an explicit formula of the invariant distribution of a process generated by (4) and we give this expression in Subsection 2.3. They also study different scaling limits which seem to be irrelevant for our problems. Another application of this example comes from population genetics. Indeed, this model can also be referred as neutral evolution, see for instance [START_REF] Etheridge | Some mathematical models from population genetics[END_REF][START_REF] Watterson | Reversibility and the age of an allele. I. Moran's infinitely many neutral alleles model[END_REF]. More precisely, consider N individuals possessing one type in F * = {1, . . . , K} at time t. Each pair of individuals interacts at rate p. Upon an interacting event, one individual dies and the other one reproduces. In addition, every individual changes its type (mutates) at rate 1 and chooses uniformly at random a new type in F * . The measure µ N t gives the proportions of types. The kind of mutation we consider here is often referred as parent-independent or the house-of-cards model. In all this section, for any probability measure µ on E, we set in a classical manner

E µ [•] = F * E x [•]µ(dx) and P µ = E µ [1 • ];
similarly Cov µ and Var µ are defined with respect to E µ .

2.1. The associated killed process. We define the process (X t ) t≥0 by setting

X t = Z t if t < τ 0 if t ≥ τ,
where τ is an exponential variable with mean 1/p and (Z t ) t≥0 is the classical complete graph random walk (i.e. without extinction) on {1, . . . , K}. We have, for any bounded function f ,

T t f (x) = E [f (X t ) | X 0 = x, X t = 0] , t ≥ 0, x ∈ F * .
The conditional distribution of X t is simply given by the distribution of Z t :

P(X t = i | X t = 0) = P(Z t = i).
The study of (Z t ) t≥0 is trivial. Indeed, it converges exponentially fast to the uniform distribution π K on {1, . . . , K}. We deduce that for all t ≥ 0 and all initial distribution µ,

d TV (µT t , π K ) = K i=1 |P µ (X t = i | τ > t) -π K (i)| ≤ e -t .
Thus in this case, the conditional distribution of X converges exponentially fast to the Yaglom limit π K .

2.2.

Correlations at fixed time. The special form of L, defined at (4), makes the calculation of the two-particle correlations at fixed time easy.

Theorem 2.1 (Two-particle correlations). For all k, l ∈ {1, . . . , K}, k = l and any probability measure µ on E, we have for all t ≥ 0

Cov µ (η t (k), η t (l)) = E µ [η 0 (k)η 0 (l)] e - 2K(N -1+p) K(N -1) t + -N + 1 + 2pN K(N -1 + 2p) (E µ [η 0 (k)] + E µ [η 0 (l)])e -t -E µ [η 0 (k)] E µ [η 0 (l)] e -2t + -N 2 (p + 1) + N K 2 (N -1 + p) .
Remark 2.2 (Limit t → +∞). By the previous theorem, we find for any probability measure µ

lim t→+∞ Cov µ (η t (k), η t (l)) = -N 2 (p + 1) + N K 2 (N -1 + p) = Cov(η(k), η(l)),
where η is distributed according to the invariant distribution; it exists since the state space is finite, see the next section. Remark 2.3 (Limit N → +∞). If Cov µ (η 0 (k), η 0 (l)) = 0 then for all k, l ∈ {1, . . . , K}, k = l and any probability measure µ, we have

Cov µ η t (k) N , η t (l) N ∼ N e -2t Cov µ η 0 (k) N , η 0 (l) N ,
where

u N ∼ N v N iff lim N →+∞ u N v N = 1.
Proof of Theorem 2.1. For k, l ∈ {1, .., K}, let ψ k,l be the function η → η(k)η(l). Applying the generator (4) to ψ k,l we obtain

Lψ k,l (η) = - 2K(N -1 + p) K(N -1) η(k)η(l) + N -1 K (η(k) + η(l)).
So, for all t ≥ 0,

Lψ k,l (η t ) = - 2K(N -1 + p) K(N -1) η t (k)η t (l) + N -1 K (η t (k) + η t (l)).
Using Kolmogorov's equation, we have

∂ t E µ (η t (k)η t (l)) = - 2K(N -1 + p) K(N -1) E µ (η t (k)η t (l)) + N -1 K (E µ (η t (k)) + E µ (η t (l))). (5) Now if ϕ k (η) = η(k) then Lϕ k (η) = N K -η(k).
We deduce that, for every t ≥ 0,

∂ t E µ (η t (k)) = N K -E µ (η t (k)) and E µ (η t (k)) = E µ (η 0 (k))e -t + N K .
Solving equation ( 5) ends the proof.

2.3.

Properties of the invariant measure. As (η t ) t≥0 is an irreducible Markov chain on a finite state space, it is straightforward that it admits a unique invariant measure. In fact, this invariant distribution is reversible and we know its expression.

Theorem 2.4 (Invariant distribution). The process (η t ) t≥0 admits a unique invariant and reversible measure ν N , which is defined, for every η ∈ E, by

ν N ({η}) = Z -1 K i=1 η(i)-1 j=0 N -1 + Kpj j + 1 ,
where Z is a normalizing constant.

This result is a slight generalisation of [START_REF] Giardinà | Correlation inequalities for interacting particle systems with duality[END_REF]Section 4] and [15, Theorem 2.1].

Proof. A measure ν is reversible if and only if it satisfies the following balance equation

ν({η})C(η, ξ) = ν({ξ})C(ξ, η) (6)
where

ξ = T i→j η and C(η, ξ) = L1 ξ (η) = η(i)(K -1 + pη(j)(N -1) -1 ).
Due to the geometry of the complete graph, it is natural to consider that ν has the following form

ν({η}) = 1 Z K i=1 l(η(i)),
where l : {0, . . . , N } → [0, 1] is a function and Z is a normalizing constant. From (6), we have

l(η(i))l(η(j))η(i)(N -1 + Kpη(j)) = l(η(i) -1)l(η(j) + 1)(η(j) + 1)(N -1 + Kp(η(i) -1)),
for all η ∈ E and i, j ∈ {1, . . . K}. Hence,

l(n) l(n -1) n N -1 + Kp(n -1) = l(m) l(m -1) m N -1 + Kp(m -1) = u,
for every m, n ∈ {1, . . . , N } and some u ∈ R.

({η}) = K i=1   u η(i) η(i)-1 j=0 N -1 + Kpi i + 1 l(0)   = l(0) K u N K i=1 η(i)-1 j=0 N -1 + Kpj j + 1 ,
and Z = 1/(l(0) K u N ).

In particular, we have directly Corollary 2.5 (Invariant distribution when p = 1/K). If p = 1/K then the process (η t ) t≥0 admits a unique invariant and reversible measure ν N , which is defined, for every η ∈ E, by

ν N ({η}) = Z -1 K i=1 N -2 + η(i) N -2 ,
where Z is a normalizing constant given by

Z = (K + 1)N -K -1 KN -K -1 .
Corollary 2.6 (Marginal laws when p = 1/K). If p = 1/K then for all i ∈ {1, . . . , K} we have

P ν N (η(i) = x) = 1 Z N -2 + x N -2 KN -K -x (K -1)N -K ,
Proof. Firstly let us recall the Vandermonde binomial convolution type formula: let n, n 1 , . . . , n p be some non-negative integers satisfying

p i=1 n i = n, we have r -1 n -1 = r 1 +•••+rp=r p j=1 r j -1 n j -1 .
The proof is based on the power series decomposition of z

→ (z/(1 -z)) n = p i=1 (z/(1 -z)) n i .
Using this formula, we find

P ν N (η(i) = x) = x∈E 1 P ν N (η = (x 1 , . . . , x i-1 , x, x i+1 . . . , x K )) = 1 Z N -2 + x N -2 x∈E 1 i-1 l=1 K l=i+1 N -2 + x l N -2 = 1 Z N -2 + x N -2 (K -1)(N -1) + N -x -1 (K -1)(N -1) -1 ,
where

E 1 = {x = (x 1 , . . . , x i-1 , x i+1 . . . , x K )|x 1 + • • • + x i-1 + x i+1 • • • + x K = N -x}.
We are now able to express the particle correlations under this invariant measure.

Theorem 2.7 (Correlation estimates). For all i = j ∈ {1, . . . , K}, we have Proof. Let η be a random variable with law ν N . As η(1), . . . , η(K) are identically distributed and

|Cov ν N (η(i)/N, η(j)/N )| ∼ N p + 1 K 2 N ,
K i=1 η(i) = N we have Cov ν N (η(i)/N, η(j)/N ) = - Var ν N (η(i)/N ) K -1 .
Using the results of Section 2.4, we have

L(η(i) 2 ) = η(i) 2 -2 - 2p N -1 + η(i) 2N K + 2pN N -1 + K -2 K + N K .
Using the fact that L(η(i) 2 )dν N = 0 and η(i)dν N = N K , we deduce that

η(i) 2 dν N = N [(2N + K -2)(N -1) + 2KN p + K(N -1)] 2K 2 (N -1 + p) .
Finally,

Var ν N (η(i)) = η(i) 2 dν N - η(i)dν N 2 = N (K -1)(N p + N -1) K 2 (N -1 + p) ,
and thus, for i = j,

|Cov ν N (η(i)/N, η(j)/N )| ∼ N p + 1 K 2 N .
Remark 2.8 (Proof through coalescence methods). Maybe we can use properties of Kingman's coalescent type process (which is a dual process) to recover some of our results (as for instance the previous correlation estimates). Indeed, after an interacting event, all individuals evolve independtly and it is enough to look when the first mutation happens (backwards in time) on one of the genealogical tree branches. Nevertheless, we prefer to use another approach based on Markovian techniques.

Remark 2.9 (Number of sites). Theorem 2.7 gives the rate of the decay of correlations with respect to the number of particles, but we also have a rate with respect to the number of sites K.

For instance when p = 1/K and if η is distributed under the invariant measure, then

|Cov ν N (η(i)/N, η(j)/N )| ∼ K 1 K(K -1)N .
The previous theorem shows that the occupation numbers of two distinct sites become noncorrelated when the number of particles increases. In fact, Theorem 2.7 leads to a propagation of chaos:

Corollary 2.10 (Convergence to the QSD). We have

E ν N d TV (µ N , π K ) ≤ K(p + 1) N ,
where π K is the uniform measure on {1, . . . , K}.

Proof. By the Cauchy-Schwarz inequality, we have

E ν N η(k) N - 1 K ≤ E ν N η(k) N - 1 K 2 1 2 = Var ν N η(k) N 1/2 ≤ (K -1)(p + 1) K 2 N .
Summing over {1, . . . , K} ends the proof.

E ν N d TV (µ N , π K ) ≤ C N θ .
All constants are explicit and we have θ < 1/2. The last corollary then gives a better bound. To our knowledge, it is the first time that this rate of convergence is obtained for the Fleming-Viot process in discrete space. With spectral arguments, this type of result was obtained for diffusion processes in [START_REF] Rousset | On the control of an interacting particle estimation of schrödinger ground states[END_REF]. This bound is achieved because of the absence of bias term. Indeed,

∀k ∈ F * , E ν N [µ N (k)] = 1 K = π K (k).
The bad term in [8, Theorem 1.2] comes from, with the notations of its proof, the estimation of |u k (t) -v k (t)| and Gronwall Lemma.

Remark 2.11 (Parameters depending on N ). A nice application of explicit rates of convergence is to consider parameters depending on N . For instance, we can now consider that p = p N depends on N , this does not change neither the conditioned semi-goup nor the QSD but this changes the dynamics of our interacting-particle system. The last corollary gives that if lim

N →∞
p N /N = 0 then the empirical measure converges to the uniform measure.

2.4.

Long time behavior and spectral analysis of the generator. In this subsection, we point out the optimality of Theorem 1.1 in this special case. Conditions in Theorem 1.1, which seems to be a bit strong, are tight in the complete graph dynamics. In that case, λ = ρ = 1 and the bound obtained is optimal in terms of contraction. Moreover, the obtained rate is exactly the spectral gap.

Corollary 2.12 (Wasserstein contraction). For any processes (η t ) t>0 and (η t ) t>0 generated by (4), and for any t ≥ 0, we have

W d (Law(η t ), Law(η t )) ≤ e -t W d (Law(η 0 ), Law(η 0 )).
In particular, when (η 0 ) follows the invariant distribution ν N associated to (4), we get for every

t ≥ 0 W d (Law(η t ), ν N ) ≤ e -t W d (Law(η 0 ), ν N ).
In particular, if λ 1 is the smallest positive eigenvalue of -L, defined at (4), then we have

1 = ρ ≤ λ 1 .
Indeed, on the one hand, let us recall that, as the invariant measure is reversible, λ 1 is the largest constant such that lim

t→+∞ e 2λt R t f -ν N (f ) 2 L 2 (ν N ) = 0, (7)
for every λ < λ 1 and f ∈ L 2 (ν N ), where (R t ) t≥0 is the semi-group generated by L. See for instance [START_REF] Bakry | L'hypercontractivité et son utilisation en théorie des semigroupes[END_REF][START_REF] Saloff-Coste | Lectures on finite Markov chains[END_REF]. On the other hand, if λ < 1 then, by Theorem 1.1, we have

e 2λt R t f -ν N (f ) 2 L 2 (ν N ) = e 2λt E ((δ η R t )f -(ν N R t )f ) 2 ν N (dη) ≤ 2e 2λt f 2 ∞ E W d (δ η R t , ν N R t ) 2 ν N (dη) ≤ 2e 2(λ-1)t f 2 ∞ E W d (δ η , ν N ) 2 ν N (dη),
and then (7) holds. Now, the constant functions are trivially eigenvectors of L associated with the eigenvalue 0, and if, for k ∈ {1, . . . , K}, l ≥ 1 we set ϕ

(l) k : η → η(k) l then the function ϕ (1) k satisfies Lϕ (1) k = N/K -ϕ (1)
k .

In particular ϕ

k -N/K is an eigenvector and 1 is an eigenvalue of -L. This gives λ 1 ≤ 1 and finally λ 1 = 1 is the smallest eigenvalue of -L. By the reversibility, we have a Poincaré (or spectral gap) inequality

∀t ≥ 0, R t f -ν N (f ) 2 L 2 (ν N ) ≤ e -2t f -ν N (f ) 2 L 2 (ν N ) .
Remark 2.13 (Complete graph random walk). If (a i ) 1≤i≤K is a sequence such that

K i=1 a i = 0 then the function K i=1 ϕ (1)
i is an eigenvector of L. However, if L is the generator of the classical complete graph random walk, La = -a and then a is also an eigenvector of L with the same eigenvalue.

Let us finally give the following result on the spectrum of L:

Lemma 2.14 (Spectrum of -L). The spectrum of -L is included in

K i=1 λ l i | l 1 , . . . , l K ∈ {0, . . . , N } , where ∀l ∈ {0, . . . , N }, λ l = l + l(l -1)p N -1 .
Proof. For every k ∈ {1, . . . , K} and l ∈ {0, . . . , N }, we have

Lϕ (l) k (η) = -λ l ϕ (l) k (η) + Q l-1 (η)
, where Q l-1 is a polynomial whose degree is less than l -1. A straightforward recurrence shows that whether there exists or not a polynomial function ψ (l) k , whose degree is l, satisfying Lψ

(l) k = -λ l ψ (l) k (namely ψ (l)
k is an eigenvector of L). Indeed, it is possible to have ψ (l) k = 0 since the polynomial functions are not linearly independent (F is finite). More generally, for all l 1 , . . . , l K ∈ {1, . . . , N }, there exists a polynomial Q with K variables, whose degree with respect to the i th variable is strictly less than l i , such that the function φ : η

→ K i=1 η(k i ) l i + Q(η) satisfies Lφ = -λφ where λ = K i=1 λ l i .
Again, provided that φ = 0, φ is an eigenvector and λ an eigenvalue of -L. Finally, as the state space is finite, using multivariate Lagrange polynomial, we can prove that every function is polynomial and thus we capture all the eigenvalues.

Remark 2.15 (Cardinal of E). As card(F

* ) = K, we have card(E) = N + K -1 K -1 = (N + K -1)! N !(K -1)! .
In particular, the number of eigenvalues is finite and less than card(E).

Remark 2.16 (Marginals). For each k, the random process (η t (k)) t≥0 , which is a marginal of a process generated by (4), is a Markov process on N N = {0, . . . , N } generated by for every function f on N N and x ∈ N N . We can express the spectrum of this generator. Indeed, let ϕ l : x → x l , for every l ≥ 0. The family (ϕ l ) 0≤l≤N is linearly independent as can be checked with a Vandermonde determinant. This family generates the L 2 -space associated to the invariant measure since this space has a dimension equal to N + 1. Now, similarly to the proof of the previous lemma, we can prove the existence of N + 1 polynomials, which are eigenvectors and linearly independent, whose eigenvalues are λ 0 , λ 1 , . . . , λ N .

Gf (x) = (N -x) 1 K + px N -1 (f (x + 1) -f (x)) + x K -1 K + p(N -x) N -1 (f (x -1) -f (x)),

THE TWO POINT SPACE

In all this section we denote by p 0 the function i ∈ F * → Q i,0 . We consider a Markov chain defined on the states {0, 1, 2} where 0 is the absorbing state. Its infinitesimal generator G is defined by

G =   0 0 0 p 0 (1) -a -p 0 (1) a p 0 (2) b -b -p 0 (b),  
where a, b > 0, p 0 (1), p 0 (2) ≥ 0 and p 0 (1) + p 0 (2) > 0. The generator of the Fleming-Viot process with N particles applied to bounded functions f :

E → R reads Lf (η) = η(1) a + p 0 (1) η(2) N -1 (f (T 1→2 η) -f (η)) + η(2) b + p 0 (2) η(1) N -1 (f (T 2→1 η) -f (η)). (8)
3.1. The associated killed process. The long time behavior of the conditionned process is related to the eigenvalues and eigenvectors of the matrix:

M = -a -p 0 (1) a b -b -p 0 (2) .
Indeed see [18, section 3.1]. Its eigenvalues are given by

λ + = -(a + b + p 0 (1) + p 0 (2)) + (a -b + p 0 (1) -p 0 (2)) 2 + 4ab 2 , λ -= -(a + b + p 0 (1) + p 0 (2)) -(a -b + p 0 (1) -p 0 (2)) 2 + 4ab 2 ,
and the corresponding eigenvectors are respectively given by

v + = a -A + √ A 2 + 4ab and v -= a -A - √ A 2 + 4ab , where A = a -b + p 0 (1) -p 0 (2). Also set ν = v + /(v + (1) + v + (2)).
From these properties, we deduce that Lemma 3.1 (Convergence to the QSD). There exists a constant C > 0 such that for every initial distribution µ, we have ∀t ≥ 0, d TV (µT t , ν) ≤ Ce -(λ + -λ -)t .

Proof. See [START_REF] Méléard | Quasi-stationary distributions and population processes[END_REF]Theorem 7] and [START_REF] Méléard | Quasi-stationary distributions and population processes[END_REF]Remark 3].

Note that

λ + -λ -= (a + b) 2 + 2(a -b)(p 0 (1) -p 0 ( 2

3.2.

Explicit formula of the invariant distribution. Firstly note that, as

∀η ∈ E, η(1) + η(2) = N,
each marginal of (η t ) t≥0 is a Markov process: Lemma 3.2 (Markovian marginals). The random process (η t (1)) t≥0 , which is a marginal of a process generated by [START_REF] Cloez | Quantitative results for the Fleming-Viot particle system and quasi-stationary distributions in discrete space[END_REF], is a Markov process generated by G defined by

Gf (n) = b n (f (n + 1) -f (n)) + d n (f (n -1) -f (n)), (9)
for any function f and n ∈ N N = {0, . . . , N }, where

b n = (N -n) b + p 0 (2) n N -1 and d n = n a + p 0 (1) N -n N -1 .
Proof. For every η ∈ E, we have η = (η(1), N -η( 1)) thus the Markov property and the generator are easily deducible from the properties of (η t ) t≥0 .

From this result and the already known results on birth and death processes [START_REF] Chafaï | Intertwining and commutation relations for birth-death processes[END_REF][START_REF] Chen | From Markov chains to non-equilibrium particle systems[END_REF], we deduce that (η t (1)) t≥0 admits an invariant and reversible distribution π given by

π(n) = u 0 n k=1 b k-1 d k and u -1 0 = 1 + N k=1 b 0 • • • b k-1 d 1 • • • d k ,
for every n ∈ N N . This gives

π(n) = u 0 N n n k=1 b(N -1) + (k -1)p 0 (2) a(N -1) + (N -k)p 0 (1) , and u
-1 0 = 1 + N k=1 b(N -1) + kp 0 (2) a(N -1) + kp 0 (1) .
Similarly, as η t (2) = N -η t (1), the process (η t (2)) t≥0 is a Markov process whose invariant distribution is also easily calculable. The invariant law of (η t ) t≥0 , is then given by

ν N ((r 1 , r 2 )) = π ({r 1 }) , ∀(r 1 , r 2 ) ∈ E.
Note that if p 0 is not constant then we can not find a basis of orthogonal polynomials in the L 2 space associated to ν N . It is then very difficult to express the spectral gap or the decay rate of the correlations.

3.3. Rate of convergence. Applying Theorem 1.1, in this special case, we find:

Corollary 3.3 (Wasserstein contraction).
For any processes (η t ) t>0 and (η t ) t>0 generated by [START_REF] Cloez | Quantitative results for the Fleming-Viot particle system and quasi-stationary distributions in discrete space[END_REF], and for any t ≥ 0, we have

W d (Law(η t ), Law(η t )) ≤ e -ρt W d (Law(η 0 ), Law(η 0 )), where ρ = a + b -(sup(p 0 ) -inf(p 0 )).
In particular, when (η 0 ) follows the invariant distribution ν N of (8), we get for every t > 0

W d (Law(η t ), ν N ) ≤ e -ρt W d (Law(η 0 ), ν N ).
This result is not optimal. Nevertheless, the error does not come from the coupling choice of [START_REF] Cloez | Quantitative results for the Fleming-Viot particle system and quasi-stationary distributions in discrete space[END_REF] but it comes from how the distance is estimated. Indeed, this coupling induces a coupling between two processes generated by G defined by [START_REF] Del Moral | On the stability of measure valued processes with applications to filtering[END_REF]. More precisely, let L = L Q + L p be the generator of the coupling introduced in the proof of [START_REF] Cloez | Quantitative results for the Fleming-Viot particle system and quasi-stationary distributions in discrete space[END_REF]Theorem 1.1] in this special case. We set G = G Q +G p , where for any n, n ∈ N N and f on E × E, and

L Q f ((n, N -n), (n , N -n )) = G Q ϕ f (n, n ), L p f ((n, N -n), (n , N -n )) = G p ϕ f (n, n ),
ϕ f (n, n ) = f ((n, N -n), (n , N -n )). It satisfies, for any function f on N N and n > n two elements of N N , G Q f (n, n ) = na f (n -1, n -1) -f (n, n ) + (N -n )b f (n + 1, n + 1) -f (n, n ) + (n -n)b f (n + 1, n ) -f (n, n ) + (n -n)a f (n, n -1) -f (n, n ) ,
and

G p f (n, n ) = p 0 (1) n(N -n ) N -1 f (n -1, n -1) -f (n, n ) + p 0 (2) n(N -n ) N -1 f (n + 1, n + 1) -f (n, n ) + p 0 (1) n(n -n) N -1 f (n -1, n ) -f (n, n ) + p 0 (2) (N -n )(n -n) N -1 f (n, n + 1) -f (n, n ) + p 0 (2) n(n -n) N -1 f (n + 1, n ) -f (n, n ) + p 0 (1) (N -n )(n -n) N -1 f (n, n -1) -f (n, n ) .
Now, for any sequence of positive numbers (u k ) k∈{0,...,N -1} , we introduce the distance δ u defined by

δ u (n, n ) = n -1 k=n u k ,
for every n, n ∈ N N such that n > n. For all n ∈ N N \{N }, we have Gδ u (n, n + 1) ≤ -λ u δ u (n, n + 1) where

λ u = min k∈{0,...,N -1} d k+1 -d k u k-1 u k + b k -b k+1 u k+1 u k ,
and thus, by linearity, Gδ u (n, n ) ≤ -λ u δ u (n, n ), for every n, n ∈ N N . This implies that for any processes (X t ) t≥0 and (X t ) t≥0 generated by G , and for any t ≥ 0, W δu (Law(X t ), Law(X t )) ≤ e -λut W δu (Law(X 0 ), Law(X 0 )).

Note that, for every n, n ∈ N N , we have

min(u)d((n, N -n), (n , N -n )) ≤ δ u (n, n ) ≤ max(u)d((n, N -n), (n , N -n )),
and then for any processes (η t ) t≥0 and (η t ) t≥0 generated by [START_REF] Cloez | Quantitative results for the Fleming-Viot particle system and quasi-stationary distributions in discrete space[END_REF], and for any t ≥ 0, we have

W d (Law(η t ), Law(η t )) ≤ max(u) min(u) e -λut W d (Law(η 0 ), Law(η 0 )).
Finally, using [START_REF] Chen | From Markov chains to non-equilibrium particle systems[END_REF]Theorem 9.25], there exists a positive sequence v such that λ v = max u λ u > 0 is the spectral gap of the birth and death process (η t (1)) t≥0 . These parameters depend on N and so we should write the previous inequality as

W d (Law(η t ), Law(η t )) ≤ C(N )e -λ N t W d (Law(η 0 ), Law(η 0 )), (10)
where C(N ) and λ N are two constants depending on N . In conclusion, the coupling introduced in Theorem 1.1 gives the optimal rate of convergence but we are not able to express a precise expression of λ N and C(N ). Nevertheless, in the section that follows, we will prove that, whatever the value of the parameters, the spectral gap is always bounded from below by a positive constant not depending on N .

3.4.

A lower bound for the spectral gap. In this subsection, we study the evolution of (λ N ) N ≥0 . Calculating λ N for small value of N (it is the eigenvalue of a small matrix) and some different parameters show that, in general, this sequence is not monotone and seems to converge to λ + -λ -.

We are not able to prove this, but as it is trivial that for all N ≥ 0, λ N > 0, we can hope that it is bounded from below. The aim of this section is to prove this fact. Firstly, using similar arguments of subsection 2.4, we have λ N ≥ ρ, for every N ≥ 0. This result does not give us information in the case ρ ≤ 0. However, we can use Hardy's inequalities [START_REF] Ané | Sur les inégalités de Sobolev logarithmiques, volume 10 of Panoramas et Synthèses [Panoramas and Syntheses[END_REF]Chapter 6] and mimic some arguments of [START_REF] Miclo | An example of application of discrete hardy's inequalities. Markov Process[END_REF] to obtain:

Theorem 3.4 (A lower bound for the spectral gap). If ρ ≤ 0 then there exists c > 0 such that

∀N ≥ 0, λ N > c.
The rest of this subsection aims to prove this result. Hardy's inequalities are mainly based on the estimation of the quantities B N,+ and B N,-defined for every i ∈ N by

B N,+ (i) = max x>i   x y=i+1 1 π(y)d y   π([x, N ]), (11)
and

B N,-(i) = max x<i i-1 y=x 1 π(y)b y π([1, x]).
We recall that π = π N is the invariant distribution defined in Subsection 3.2 and jumps rates b and d also depend on N . More precisely, [START_REF] Miclo | An example of application of discrete hardy's inequalities. Markov Process[END_REF]Proposition 3] shows that if one wants to get a "good" lower bound of the spectral gap, one only needs to guess an "adequate choice" of i and to apply the estimate

λ N ≥ 1 4 max{B N,+ (i), B N,-(i)}
.

So, we have to find an upper bound for these two quantities. Before to give it, let us prove that the invariant distribution π is unimodal. Indeed, it will help us to choose an appropriate i.

Lemma 3.5 (Unimodality of π). The sequence (π(i + 1)/π(i)) i≥0 is decreasing.

Proof of Lemma 3.5. For all i ∈ {1, . . . , N }, we set

g(i) = π(i + 1) π(i) = (N -i)(b(N -1) + ip 0 (2)) (i + 1)((a + p 0 (1))(N -1) -ip 0 (1))
.

It follows that

g(i + 1) -g(i) = Λ N (i) (i + 1)((a + p 0 (1))(N -1) -ip 0 (1))(i + 2)((a + p 0 (1))(N -1) -(i + 1)p 0 (1))
where

Λ N (i) = (N -i -1)(b(N -1) + (i + 1)p 0 (2))(i + 1)((a + p 0 (1))(N -1) -ip 0 (1)) -(N -i)(b(N -1) + ip 0 (2))(i + 2)((a + p 0 (1))(N -1) -(i + 1)p 0 (1)) = -[b(N -1) -p 0 (2)] [(N + 1) (a(N -1) -p 0 (1)) + p 0 (1)(N -i)(N -i -1)] -p 0 (2) i 2 + 3i + 2 (a(N -1) -p 0 (1)) ≤ 0.
We deduce the result. Proof of Theorem 3.4. Without less of generality, we assume that p 0 (1) ≥ p 0 (2) and we recall that ρ ≤ 0. We would like to know where π reaches its maximum i * since it will be a good candidate to estimate B N,+ (i * ) and B N,-(i * ). From the previous lemma, to find it, we look when π(i + 1)/π(i) is close to one. We have, for all i ∈ {1, . . . , N },

π(i + 1) π(i) = b i d i+1 = 1 + (p 0 (1) -p 0 (2))(i -i 1 )(i -i 2 ) (i + 1) ((a + p 0 (1))(N -1) -ip 0 (1)) , (12)
where i 1 and i 2 are the two real numbers given by

i 1 = N (a + b + p 0 (1) -p 0 (2)) -(a + b + 2p 0 (1)) - √ ∆ 2(p 0 (1) -p 0 (2)) and i 2 = N (a + b + p 0 (1) -p 0 (2)) -(a + b + 2p 0 (1)) + √ ∆ 2(p 0 (1) -p 0 (2))
,

where

∆ = [N (a + b + p 0 (1) -p 0 (2)) -(a + b + 2p 0 (1))] 2 -4(N -1)(bN -a -p 0 (1))(p 0 (1) -p 0 (2)).
In particular, 1 ≤ i 1 ≤ N ≤ i 2 . Furthermore, if . denotes the integer part then

π( i 1 + 2) π( i 1 + 1) ≤ 1 ≤ π( i 1 + 1) π( i 1 ) .
Let us define m N = i 1 + 1 and l N = 2( √ N + 1). Using a telescopic product, we have

π(m N + l N ) π(m N ) = π(m N + l N - √ N -1) π(m N ) √ N +1 j=1 π(m N + l N -j + 1) π(m N + l N -j) ,
Using Lemma 3.5 and the previous calculus, we have that the sequences (π(i)) i≥m N and (π(i + 1)/π(i)) i≥0 are decreasing and then

π(m N + l N ) π(m N ) ≤ π(m N + l N - √ N ) π(m N + l N - √ N -1) √ N +1
. Now using [START_REF] Ferrari | Quasi stationary distributions and Fleming-Viot processes in countable spaces[END_REF] and some equivalents, there exists a constant δ 1 > 0 (not depending on N ) such that

π(m N + l N - √ N ) π(m N + l N - √ N -1) ≤ 1 - δ 1 √ N .
Using the fact that 1 -x ≤ e -x for all x ≥ 0, we finally obtain π(m N + l N )/π(m N ) ≤ e -δ 1 . Similar arguments entail the existence of δ 2 > 0 (also not depending on N ) such that π(m Nl N )/π(m N ) ≤ e -δ 2 . In conclusion, using Lemma 3.5, we have shown that for all i ≥ m N and j ≤ m N , the following inequalities holds:

π(i + l N ) ≤ e -δ 1 π(i) and π(j -l N ) ≤ e -δ 2 π(j).
We are now armed to evaluate B N,+ (m N ) defined in [START_REF] Etheridge | Some mathematical models from population genetics[END_REF]. Firstly, using the expressions of the death rate d and m N , there exist γ > 0 (not depending on N ) and N 0 ≥ 0 such that for all N ≥ N 0 and all i ≥ m N + 1, d i ≥ γN . Let us fix x ≥ m N + 1, using that (π(i)) i≥m N is decreasing, we

x y=m N +1 1 π(y) = {i,k|m N +1≤k-il N ≤x} 1 π(k -il N ) ≤ {i,k|m N +1≤k-il N ≤x} e -δ 1 i π(k) ≤ 1 1 -e -δ 1 x k=x-l N +1 1 π(k) ≤ l N π(x) 1 1 -e -δ 1 .
Similarly, we have

π([x, N ]) = {k,i|x≤k+il N ≤N } 1 {x+il N ≤N } Π N (k + il N ) ≤ l N π(x) 1 -e -δ 1 .
Using these three estimates, we deduce that, for every N ≥ N 0 , 3.5. Simulation and evolution of the spectral gap (of the Fleming-Viot process). As stated in Lemma 3.1, the rate of convergence λ of the conditioned semi-group to the QSD is explicit. For a fixed N ≥ 1, the Fleming-Viot process is a reversible Markov chain (for this example) on a finite space, it then converges to equilibrium at rate λ N , where λ N is the eigenvalue of its generator closer to 0. Theorem 3.4 shows that the rate of convergence to equilibrium of the Fleming-Viot process is uniformly bounded. Nevertheless, there is some natural questions such as the convergence of λ N to λ or the monotonicity of this sequence. Also, do we have λ N ≤ λ or λ N ≥ λ? When there are few particles, the generator of the Fleming-Viot process is a square matrix of size 2N and one can explicitly calculate its spectrum (with the help of a computer). In figure 1, each graphic represents, with different parameters, the evolution of the spectral gap of the conditioned process λ = λ + -λ -(detailed in Lemma 3.1), of the Fleming Viot particle system λ N and of the upper-bound ρ in Theorem 1.1, with respect to the number of particles N . Graphics (i) and (ii) illustrate that when osc(p 0 ) = p 0 (2) -p 0 (1) remains constant for two parameter choices then the Fleming-Viot dynamics are different although the conditioned semigroups are the same. More p 0 (1) is large more jumps there are; these graphics seem to show that the rate of convergence is dragged down by the interactions. Graphic (iii) shows that λ N is neither increasing nor decreasing and neither upper nor lower than λ. Graphic (iv) shows that when osc(p 0 ) is small then all curves are close.

B N,+ (m N ) ≤ 1 γN l N 1 -e -δ 1
In any case, it seems that λ N converges to λ but this point remains an open question. where C(N ) is defined in [START_REF] Del Moral | A Moran particle system approximation of Feynman-Kac formulae[END_REF]. Even if Theorem 3.4 gives an estimate of λ N , C(N ) is not (completely) explicit and we do not know if the right-hand side of the previous expression tends to 0 as N tends to infinity. This example shows the difficulty of finding explicit and optimal rates of the convergence towards equilibrium and the decay of correlations.

  Comment citer ce document : Cloez, B., Thai, M.-N. (2016). Fleming-Viot processes: two explicit examples. Alea-Latin American Journal of Probability and Mathematical Statistics, 13 (1), 337-356.

  Finally, Comment citer ce document : Cloez, B., Thai, M.-N. (2016). Fleming-Viot processes: two explicit examples. Alea-Latin American Journal of Probability and Mathematical Statistics, 13 (1), 337-356.

  ν

  Comment citer ce document : Cloez, B., Thai, M.-N. (2016). Fleming-Viot processes: two explicit examples. Alea-Latin American Journal of Probability and Mathematical Statistics, 13 (1), 337-356.

  Comment citer ce document : Cloez, B., Thai, M.-N. (2016). Fleming-Viot processes: two explicit examples. Alea-Latin American Journal of Probability and Mathematical Statistics, 13 (1), 337-356.

  )) + (p 0 (1) -p 0 (2))2 and ρ = a + b -(sup(p 0 ) -inf(p 0 )) defined in Theorem 1.1. We have then λ + -λ -> ρ when sup(p 0 ) > inf(p 0 ). In particular Theorem 1.1 seems not optimal. Comment citer ce document : Cloez, B., Thai, M.-N. (2016). Fleming-Viot processes: two explicit examples. Alea-Latin American Journal of Probability and Mathematical Statistics, 13 (1), 337-356.

  Comment citer ce document : Cloez, B., Thai, M.-N. (2016). Fleming-Viot processes: two explicit examples. Alea-Latin American Journal of Probability and Mathematical Statistics, 13 (1), 337-356.

  Comment citer ce document : Cloez, B., Thai, M.-N. (2016). Fleming-Viot processes: two explicit examples. Alea-Latin American Journal of Probability and Mathematical Statistics, 13 (1), 337-356.

  Comment citer ce document : Cloez, B., Thai, M.-N. (2016). Fleming-Viot processes: two explicit examples. Alea-Latin American Journal of Probability and Mathematical Statistics, 13 (1), 337-356.

  have

FIGURE 1 .

 1 FIGURE 1. Evolution of the spectral gap with respect to the number of particles. Details are described in Subsection 3.5

3. 6 . 2 1 2 N - 1 .

 6221 Correlations. Using [8, Theorem 2.6], we have Corollary 3.6 (Correlations). If (η t ) t≥0 is a process generated by (8) then we have for all t ≥ 0,Cov(η t (k)/N, η t (l)/N ) ≤ 2 N -e -2ρt ρ N (a ∨ b) + sup(p 0 ) N If ρ ≤ 0,the right-hand side of the previous inequality explodes as t tends to infinity whereas these correlations are bounded by 1. Nevertheless, using Theorem [8, Theorem 2.6], Remark [8, Remark 2.7] and Inequality (10), we can prove that there exists two constants C (N ), depending on N , and K, which does not depend on N , such that sup t≥0 Cov(η t (k)/N, η t (l)/N ) ≤ C (N ) = KC(N ) N λ N ,

  BERTRAND CLOEZ AND MARIE-NOÉMIE THAIResult [8, Theorem 1.2] and its corollaries states that there exist C, θ > 0 such that

Comment citer ce document : Cloez, B., Thai, M.-N. (2016). Fleming-Viot processes: two explicit examples. Alea-Latin American Journal of Probability and Mathematical Statistics, 13 (1), 337-356.

  Comment citer ce document : Cloez, B., Thai, M.-N. (2016). Fleming-Viot processes: two explicit examples. Alea-Latin American Journal of Probability and Mathematical Statistics, 13 (1), 337-356.