Ihsen Hedhli

Gabriele Moser

Josiane Zerubia

Sebastiano B Serpico

NEW CASCADE MODEL FOR HIERARCHICAL JOINT CLASSIFICATION OF MULTITEMPORAL, MULTIRESOLUTION AND MULTISENSOR REMOTE SENSING DATA

Keywords: Image time series, Multitemporal classification, hierarchical multiresolution Markov random fields

. An especially novel element of the proposed approach is the use of multiple quad-trees in cascade, each associated with each new available image at different dates, with the aim to characterize the temporal correlations associated with distinct images in the input time series. Experimental results are shown with multitemporal and multiresolution Pléiades data.

INTRODUCTION

The capabilities to monitor the Earth surface play primary roles from multiple social, economic, and human viewpoints. Nowadays, a wide variety of remote sensing images are available, that convey a huge potential for such applications, as they allow a spatially distributed and temporally repetitive view of the monitored area at the desired spatial scales. In this framework, accurate and time-efficient classification methods are especially important tools to support rapid and reliable assessment of the ground changes. Given the huge amount and variety of data available currently from lastgeneration very-high resolution (VHR) satellite missions, We would like to thank the French Space Agency (CNES) for providing the data and for partial financial support.

(such as Pléiades, COSMO-SkyMed, or WorldView-2), the main difficulty is to develop a classifier that can take benefit of multiband, multiresolution, multidate, and possibly multisensor input imagery. The proposed method addresses the problem of multitemporal image classification and allows input data collected at multiple resolutions to be effectively exploited for classification purposes without resampling. The approach consists of a supervised Bayesian classifier that combines a joint class-conditional statistical model for pixelwise information and a hierarchical Markov random field (MRF) for spatio-temporal and multiresolution contextual information [2] [3]. This paper is organized as follows: In Section 2, we focus on the new hierarchical method for classification in multidate imagery by detailing the mathematical modeling used in this classifier. First results of the use of this new hierarchical model in time series classification are presented in Section 3. Finally, we conclude and highlight some future works in Section 4.

THE PROPOSED MULTIDATE HIERARCHICAL MODEL

Given an input time series of remote sensing images acquired at multiple spatial resolutions, a multiscale and multitemporal model is proposed to fuse the related spatial, temporal, and multiresolution information. Two Bayesian approaches can generally be adopted for this purpose. The 'cascade' classifier (e.g., [START_REF] Swain | Bayesian classification in a time-varying environment[END_REF]) removes the coupling between the spectral and temporal dimensions and classifies each image in the input series on the basis of itself and of the previous images. The 'mutual' approach classifies each image on the basis of the previous and the subsequent images in the series [5][6]. The two approaches are basically complementary in terms of applicability: online processing can be feasible within the cascade domain, whereas the mutual domain can be generally more appropriate for batch processing. The cascade approach is adopted in the method proposed in this paper. This model is a hierarchical spatiotemporal and multiresolution MRF integrated in a quad-tree structure.

Notations and definitions

Here, we define the notations that are used throughout this paper. defines a site (pixel) and belongs to a finite set . describes the neighborhood of s and is the set of neighborhoods. The couple defines a graph . We shall focus on a specific graph, namely a tree. The set of sites is then hierarchically partitioned, i.e.:

where R corresponds to the coarsest resolution (the root of the tree), corresponds to the finest resolution (the leaves of the tree) and indicates the subset of nodes associated with the level of the tree ( . In the considered structure a parent-child relationship can be defined: for each site s in { } a unique parent and several children are well defined. Sites at the root level do not own any parent and sites at the finest level do not own any children. We also define as the set including and its descendants. For the specific case in which owns four children, the tree structure is named quad-tree. The aim of classification is to estimate a set of hidden labels (e.g., land cover class labels) given a set of observations (e.g., satellite data) attached to the sites (pixels). X and Y are considered as random processes, the restriction of (resp. Y) to the level n is { } (resp. { }). Some extra hypotheses are needed to ensure that is Markovian with respect to scale, i.e.

( | ) | (1)
Where and are scales and is a realization of at scale .

Multitemporal MPM model

An especially novel element of the proposed approach is the use of multiple quad-trees in cascade, each associated with the set of images available at each observation date y in the considered time series. Specifically, for each date, the input images are inserted in the quad-tree structure on the basis of their resolutions, while missing levels of the tree are filled in with wavelet transforms of the images embedded in finerresolution levels [START_REF] Voisin | Supervised classification of multisensor and multiresolution remote sensing images with a hierarchical copulabased approach[END_REF]. This approach is aimed at both exploiting multiscale information, which is known to play a crucial role in very high resolution image analysis, and to support input images acquired at different resolutions. The choice of a quad-tree structure relies on the good properties of this configuration such as the causality in scale, under Markovian assumption, which allows the use of a non-iterative algorithm and to integrate an exact estimator of the maximizer of posterior marginal (MPM) [START_REF] Marroquin | Probabilistic solution of ill-posed problems in computational vision[END_REF] [START_REF] Laferte | Discrete Markov modeling and inference on the quad-tree[END_REF]. The aim is to maximize recursively the posterior marginal at each site :

̂ | (2)
i.e., the probability distribution of the label of each site, conditioned to the whole set of available observations in the quad-tree. Specifically, a novel formulation of MPM is proposed for the resulting multitemporal quad-tree. The posterior marginal ( | ) of the label of each spatiotemporal site is expressed as a function not only of the posterior marginal ( | ) of the parent node in the corresponding quad-tree but also of the posterior marginal ( | ), of the parent node in the quad-tree associated with the previous date, with the aim to characterize the temporal correlations associated, at different scales, with distinct images in the input time series:

The steps defined above involve the pixel-wise classconditional PDFs of the image data at each node of each quad-tree, given the corresponding (satellite or wavelet) features and the transition probabilities between consecutive scales and consecutive dates.

Likelihood term

Given a training set for each input date, for each class, scale, acquisition time, and (satellite or wavelet) feature, we model the corresponding class-conditional marginal PDF using a finite Gaussian mixture [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. The use of finite mixtures instead of single PDFs offers the possibility to consider heterogeneous PDFs, usually reflecting the contributions of the different materials present in each class. Such class heterogeneity is relevant since we deal with VHR images. The parameters of the mixture model are estimated through the stochastic expectation maximization (SEM) algorithm [START_REF] Celeux | Stochastic versions of the EM algorithm: an experimental study in the mixture case[END_REF], which is an iterative stochastic parameter estimation algorithm developed for problems characterized by data incompleteness and approaching, under suitable assumptions, maximum likelihood estimates. Given the resulting marginal conditional PDF estimates, joint conditional PDF estimates are obtained through copula functions. According to Sklar's theorem [START_REF] Nelsen | An introduction to copulas[END_REF], an arbitrary joint PDF can be expressed in terms of the corresponding marginal PDFs and of a copula function. In the proposed method, we wish to determine this specific joint PDF for each scale level of each quad-tree in the cascade model given the aforementioned marginal distributions. Specifically, the approach proposed in [START_REF] Krylov | Supervised high resolution dual polarization SAR image classification by finite mixtures and copulas[END_REF] is used for classconditional PDF estimation, that is based on a dictionary of parametric copula models. The advantage of using copulas, over the choice of a specific parametric model (e.g., multivariate Gaussian or elliptically contoured distributions), is that they enable modeling the dependence structure of any type of features. The use of copulas in the joint PDF modeling becomes especially relevant when the dependence between the source features is strong [START_REF] Voisin | Supervised classification of multisensor and multiresolution remote sensing images with a hierarchical copulabased approach[END_REF] [START_REF] Krylov | Supervised high resolution dual polarization SAR image classification by finite mixtures and copulas[END_REF].

Transition probabilities

The transition probabilities between consecutive scales and consecutive dates determine the properties of the hierarchical MRF because they formalize the causality of the statistical interactions involved. As such, they need to be carefully defined.

Transition between scales

The fundamental hypothesis in the described hierarchical MRF is to consider the Markovianity on scale (see Eq. ( 1)). To define the transition probability ( | ) we use the formulation introduced by Bouman and Shapiro [START_REF] Bouman | A multiscale random field model for Bayesian image segmentation[END_REF], for all sites and all scales

( | ) { (7)
With the parameter ⁄ where M is the number of considered classes.

Transition between consecutive observation times

In the proposed method, there is two kinds of probabilities that involve time: The temporal transition in the same scale ( | ), estimated using a specific formulation of the expectation-maximization (EM) algorithm [START_REF] Bruzzone | A Neural-Statistical Approach to Multitemporal and Multisource Remote-Sensing Image Classification[END_REF], and the child-parent transition probability ( | ). We use here a new formulation of the transition probability model proposed by Bouman and Shapiro, which favors identity between children and parents (in current and previous dates), all other transitions being unlikely:

( | ) { (8)
With the parameters ⁄ and ⁄ ⁄ .

EXPERIMENTAL RESULTS

Preliminary experiments have been performed with a time series of panchromatic and multispectral Pléiades images acquired over Port-au-Prince (Haiti). For the sake of brevity, only results associated with a pair of images acquired on 2011 and 2013 are shown in this paper (see Figure 2). The finest resolution of the multiresolution pyramid (level 0) is set equal to the finest resolution of the input images (0.5-m panchromatic). Co-registered 2-m multispectral images are integrated in level 2 of the pyramid. Level 1 is filled in through Haar wavelet decomposition of the panchromatic image [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. Preliminary experiments suggested the Haar transform to be especially effective in the application within the proposed method. Five land cover classes have been considered: urban (red), water (blue), vegetation (green), soil (yellow) and containers (mauve). A preliminary visual analysis of the resulting classification maps has suggested that the proposed hierarchical method leads to accurate results, especially as compared to separate hierarchical classification at each individual date (e.g., [START_REF] Swain | Bayesian classification in a time-varying environment[END_REF]). These results suggest the effectiveness of the proposed multitemporal hierarchical model in fusing the temporal, spatial, and multiresolution information associated with the input data.

In particular, the main source of misclassification in the single-date results is the confusion between the "urban" and "soil" classes (e.g., see highlighted regions in Figure 2(c), (d), and (e)); this misclassification is reduced in the multitemporal classification obtained by the proposed method thanks to the modeling of the temporal relationships among the input multiresolution data.

CONCLUSION

The proposed method allows performing joint classification on multidate, multiband, and multiresolution imagery. It combines a joint statistical modeling of the considered input images and a hierarchical MRF model for spatio-temporal and multiresolution contextual information, leading to a statistical supervised classification approach. Moreover, one main advantage of the proposed classifier is that it can be extended to the use not only of optical data, but also of synthetic aperture radar (e.g. COSMO-SkyMed) or multisensor data. The extension to the multisensor case will be a major direction of further research.

Figure 2 :Figure 1 :

 21 Figure 2: (a) Panchromatic image of Port au Prince (Pléiades) in 2011 © CNES, (b) Panchromatic image of Port au Prince (Pléiades) in 2013 © CNES, (c) Classification map for 2011, obtained using the multispectral and panchromatic images acquired in 2011, (d) Classification map for 2013, obtained through the method in [7] using only the multispectral and panchromatic images acquired in 2013, (e) Classification map for 2013, obtained through the proposed cascade method using the multispectral and panchromatic images acquired in 2011 and 2013.

Where bold type highlights the role of the recursively computed posterior marginals. These equations involve two conditional independence assumptions: (i) the label depends only on the data of the site and their descendants, so ( | ) ( | ); and (ii) the label at parent is independent of the label at parent at the previous date, when conditioned to the data , so | | . The multidate MPM algorithm runs in two passes on a quad tree, referred to as "bottom-up" and "top-down" passes, similar to the classical MPM in a single date (see Figure 1). The aim is to calculate recursively the posterior marginal ( | ) as well as the probabilities ( | ) are made available. This is achieved by a preliminary bottom-up pass based on:

The first factor ( | ) corresponds to the childparent transition probability; ( ) is the prior probability computed using a top down pass via:

) is the temporal transition in the same scale. An upward recursion allows computing the posterior marginal ( | ) using the formulation introduced by Laferte [START_REF] Melgni | A Markov Random Field Approach to Spatio-Temporal Contextual Image Classification[END_REF] for the single date case:

( | )


JOS Patrick Cattiaux ♠

É R Le Ón ♣

C L Émentine Prieur

Patrick Cattiaux
email: cattiaux@math.univ-toulouse.fr

José Le Ón

Clémentine Prieur
email: clementine.prieur@imag.fr

RECURSIVE ESTIMATION FOR STOCHASTIC DAMPING HAMILTONIAN SYSTEMS

Keywords: MSC 2010 : 62G07, 60H10, 60F05 hypoelliptic diffusion, nonparametric estimation, recursive estimation, partial observations

In this paper, we complete our previous works ( Cattiaux, Leon and Prieur (2014-a), Cattiaux, Leon and Prieur (2014-b), Cattiaux, Leon and Prieur (2014-c)) on the (non-parametric) estimation of the characteristics (invariant density, drift term, variance term) of some ergodic hamiltonian systems, under partial observations. More precisely, we introduce recursive estimators using the full strength of the ergodic behavior. We compare the theoretical results obtained with these estimators to the results we proved for the ones we have introduced previously.

INTRODUCTION.

Let Z t := (X t , Y t ) ∈ R 2d , t ≥ 0 be governed by the following Ito stochastic differential equation :

dX t = Y t dt dY t = σ dW t -(c(X t , Y t )Y t + ∇V (X t ))dt.
(1.1)

Each component Y i (1 ≤ i ≤ d)
is the velocity of a particle i with position X i . Function c is called the damping force and V the potential, σ is some symmetric positive definite constant matrix and W is a standard Brownian motion. We shall assume that c and V are regular enough for the existence and uniqueness of a non explosive solution of (1.1). We shall also assume that the process is ergodic with a unique invariant probability measure µ, and that the convergence in the ergodic theorem is quick enough. Some sufficient conditions will be discussed below. These models are important due to their physical relevance. They have a long history. We refer to [START_REF] Wu | Large and moderate deviations and exponential convergence for stochastic damping Hamiltonian systems[END_REF] and Cattiaux et al. (2014-a) for a detailed bibliography. We have chosen the terminology "damping Hamiltonian systems" in reference to Wu. Such systems are also called "kinetic diffusions" by several authors.

Once the probabilistic picture is well understood, it is particularly relevant to build statistical tools for these models. In a series of previous papers, we have built several estimators for the characteristics of such processes: non parametric estimators for the invariant density in Cattiaux, Leon and Prieur (2014-a), for the drift term in Cattiaux, Leon and Prieur (2014-b), and various estimators for the noise or the diffusion term in Cattiaux, Leon and Prieur (2014-c). The most important fact is that all these estimators are built under partial observations, i.e. we assume that the position only can be observed, the velocity is then estimated. This is in accordance with the concrete situations where these processes can be used as a pertinent model. The second common fact in these papers is that we assume high frequency observations, i.e. the process X t is observed at times kh n (k ∈ R + ) for some parameter h n going to 0 as n goes to infinity (but of course nh n → +∞).

In the present paper, we propose recursive estimators. Recursive estimation is of course particularly well adapted to the ergodic situation, where one can observe the process on a long time interval. To be more precise, denote by p s the invariant density (see the next section for its existence) and by b the drift term b(x, y) = -c(x, y)y -∇V (x) .

We introduce the two natural estimations of Y and d Y ,

∆ 1 (X t , h) = 1 h (X t+h -X t ) and h ∆ 2 (X t , h) = 1 h (X t+h + X t-h -2X t ) . (1.2)
We define the following three estimators pr,n (x, y)

:= 1 n n i=1 1 b d 1,i b d 2,i K x -X i b 1,i , y -∆ 1 (X i , h i ) b 2,i , (1.3)
for the invariant density, pr,n (x, y) br,n (x, y)

:= 1 n n i=1 1 b d 1,i b d 2,i K x -X i b 1,i , y -∆ 1 (X i , h i ) b 2,i ∆ 2 (X i+2h i , h i ) , (1.4)
for the drift term, and finally σr,n := 3 2n

n i=1 1 h i ∆ 2 (X i+h i , h i ) >< ∆ 2 (X i+h i , h i ) , (1.5)
for the (constant) diffusion matrix, where A = x >< x denotes the symmetric matrix with entries A j,k = x j x k , x j denoting the jth coordinate of x.

The recursion formulae are thus pr,n (x, y) = n -1 n pr,(n-1) (x, y)

+ 1 n b d 1,n b d 2,n K x -X n b 1,n , y -∆ 1 (X n , h n ) b 2,n ,
pr,n (x, y) br,n (x, y) = n -1 n pr,(n-1) (x, y) br,(n-1) (x, y)

+ + 1 n b d 1,n b d 2,n K x -X n b 1,n , y -∆ 1 (X n , h n ) b 2,n ∆ 2 (X n+2hn , h n ) , σr,n = n -1 n σr,(n-1) + 3 2nh n ∆ 2 (X n+hn , h n ) >< ∆ 2 (X n+hn , h n ) .
We intend to give sufficient conditions on the discretization steps h n , bandwidth b 1,n and b 2,n and kernels K, for these estimators to be consistent and provide some confidence intervals. In the whole paper K will be a smooth kernel with compact support, satisfying K(x, y)dxdy = 1 .

It is worth noticing that the observations will be only locally high frequency, i.e. what is needed is the joint observation of (X n , X n+hn , X n+2hn , X n+3hn ) and we will choose h n going to 0 and assume that 3h n < 1 for all n. In the case of complete observation, we only need to observe the position process at integer times.

Of course, all these estimators are the "recursive" form of the estimators we have introduced in our previous works.

Recursive estimation for the density of stationary sequences has a long history too, starting presumably with [START_REF] Wegman | Remarks on some recursive estimators of a probability density[END_REF]. One can mention among others a series of papers by E. [START_REF] Masry | Recursive probability density estimation for weakly dependent stationary processes[END_REF], [START_REF] Masry | Almost sure convergence of recursive density estimators for stationary mixing processes[END_REF], [START_REF] Masry | Nonparametric estimation of conditional probability density and expectation of stationary processes; strong consistency and rates[END_REF], [START_REF] Masry | Strong consistency and rates for recursive probability density estimators of stationary processes[END_REF]) and more recently [START_REF] Liang | Asymptotic normality of recursive density estimates under some dependence assumptions[END_REF] and [START_REF] Amiri | Sur une famille paramétrique d'estimateurs séquentiels de la densité pour un processus fortement mélangeant[END_REF]. We refer to [START_REF] Amiri | Estimateurs fonctionnels récursifs et leurs applications à la prévision[END_REF] for a more complete bibliography. Though these papers do not study explicitly diffusion processes, the results therein can be applied to the case of complete observation of some stationary diffusion processes. However we decided to give here complete new proofs for two major reasons: first the proofs of the Central Limit Theorem given in these papers, based on an intricate block decomposition can and will be simplified, second we are dealing with the partial observation case. In addition, some key assumption in all these papers is not clearly satisfied for our model, though it is presumably true. A precise description and comparison will be given later (see Section 5.2 for further discussion).

The plan of the paper is the following: after recalling some properties of the model, we study recursive estimators for the invariant density and the drift under complete observation. Then, as in Cattiaux et al. (2014-a) and Cattiaux et al. (2014-b), we show that looking at estimators with partial observation, introduce a small perturbation of the fully observed case, at least for well chosen h n . Finally we study the diffusion coefficient in the spirit of the final section of Cattiaux et al. (2014-c). Large parts of the proofs being similar to what we have done in our previous works, we only indicate the modifications or simplifications in the present setting.

The main results of the paper can be summarized as follows: choosing h n = n -β for some β > 1/2 as close as 1/2 we want, one can find kernels K and bandwidths b n,j such that the recursive estimators simultaneously converge to the invariant density, drift term and diffusion term. In addition we are able to find confidence intervals for all these quantities, of respective size n -θp for θ p < 1/2, n -θ b for θ b < 1/4 and n -1/2 .

The model and its properties.

We are obliged to recall some facts on the model. A more detailed discussion is contained in Cattiaux et al. (2014-a).

We shall first give some results about non explosion and long time behaviour. In a sense, coercivity can be seen in this context as some exponential decay to equilibrium.

Let | • | denote the euclidean norm in R d .
Let us first introduce some sets of assumptions: Hypothesis H 1 :

(i) the potential V is lower bounded, smooth over R d , V and its first three partial derivatives have polynomial growth at infinity and there exists v > 0 such that

+∞ ≥ lim inf |x|→+∞ x.∇V (x) |x| ≥ v > 0 ,
the latter being often called "drift condition", (ii) the damping coefficient c(x, y) is smooth and bounded, its first two derivatives are bounded by some polynomial, and there exist c, L > 0 so that c s (x, y) ≥ cId > 0, ∀(|x| > L, y ∈ R d ), where c s (x, y) is the symmetrization of the matrix c(x, y), given by 1 2 (c ij (x, y) + c ji (x, y)) 1≤i,j≤d , (iii) σ is uniformly elliptic, i.e. here σ -1 exists. These conditions ensure that there is no explosion, and that the process is positive recurrent with a unique invariant probability measure µ. In some of our previous papers we assumed that σ = Id, but of course multiplying Y by σ -1 we are immediately reduced (for the probabilistic properties) to the case σ = Id. Furthermore µ admits some exponential moment, hence polynomial moments of any order. We will denote by P t f (z) = E z (f (Z t )) which is well defined for all bounded function f , P t extends as a contraction semi-group on L p (µ) for all 1 ≤ p ≤ +∞.

Another key feature is that the process is actually α-mixing, i.e. Proposition 2.1. There exist some constants C > 0 and ρ < 1 such that:

∀ g, f ∈ L ∞ (µ) , ∀ t ≥ 0, |Cov µ (f (Z t ), g(Z 0 ))| ≤ C ρ t/2 g -gdµ ∞ f -f dµ ∞ .
(2.2)

i.e., in the stationary regime, (Z t , t ≥ 0) is α-mixing with exponential rate.

Actually, using the semi-group property and the Riesz-Thorin interpolation theorem, one can extend the previous mixing condition in the following way: if g, f ∈ L q (µ) for some q ≥ 2, then for t > 0,

|Cov µ (f (Z t ), g(Z 0 ))| ≤ C ρ t(q-2)/2q g -gdµ q f -f dµ q .
(2.3)

As explained in section 2.2 of Cattiaux et al. (2014-a), the infinitesimal generator L is hypo-elliptic, which implies that µ(dz) = p s (z) dz for some smooth function p s . At the same time, for all t > 0 the law of the process has a density p(t, z, .) with respect to Lebesgue measure, which is smooth, and satisfies for any t 0 > 0 and for all compact subsets A and A ′ , sup

z∈A , z ′ ∈A ′ sup t≥t 0 |p(t, z, z ′ )| = Q(A, A ′ ) < +∞ .
(2.4)

One can relax the C ∞ assumption on the coefficients into a C k assumption, for a large enough k, but this is irrelevant. Furthermore it can be shown that p s is everywhere positive.

One can relax some assumptions and still have the same conclusions: Notice that one can relax the repelling strength of the potential, and obtain, no more exponential but sub-exponential or polynomial decay (see the discussion in [START_REF] Bakry | Rate of convergence for ergodic continuous Markov processes : Lyapunov versus Poincaré[END_REF]).

Hypothesis H 2 : (a)
From now on in the whole paper we will assume that Hypothesis H 1 (or H 2 ) is fulfilled. In all the proofs of the paper C denotes some constant which may vary from line to line. We also use the expression in the stationary regime when the process is stationary, i.e. when Z 0 is distributed according to µ.

3. Estimation of the invariant density and of the drift term in the stationary regime: complete observation case.

To begin with, we look at the recursive estimators under complete observation. This is not our main goal, but as in Cattiaux et al. (2014-a), it is a first step for our study.

3.1. Invariant density. Introduce pr,n (x, y) := 1 n n i=1 1 b d 1,i b d 2,i K x -X i b 1,i , y -Y i b 2,i . (3.1)
We have the following first result Lemma 3.2. Assume that b 1,n and b 2,n are non-increasing sequences of real numbers (say between 0 and 1), going to 0 as n → +∞. Then for all (x, y), as n → +∞,

E µ (p r,n (x, y)) → p s (x, y) . If in addition nb d 1,n b d n,2 → +∞ as n → +∞, then Var µ (p r,n (x, y)) ≤ C n b d 1,n b d n,2
, hence pr,n (x, y) converges to p s (x, y) in L 2 (µ) norm, as n → +∞.

Proof. A simple change of variables shows that

E µ 1 b d 1,i b d 2,i K x -X i b 1,i , y -Y i b 2,i = K(u, v) p s (x -ub 1,i , y -vb 2,i ) du dv
so that using the smoothness of p s and the fact that K has compact support, it follows that

E µ 1 b d 1,i b d 2,i K x -X i b 1,i , y -Y i b 2,i → p s (x, y)
as i goes to infinity. The first result then follows from the Cesaro rule.

For the variance we have the decomposition Var

µ (p r,n (x, y)) = V 1 (n) + 2V 2 (n) with V 1 (n) = 1 n 2 n i=1 1 b 2d 1,i b 2d 2,i Var µ K x -X 0 b 1,i , y -Y 0 b 2,i (3.3) = 1 n 2 n i=1 1 b d 1,i b d 2,i K 2 (u, v) p s (x -ub 1,i , y -vb 2,i ) du dv - 1 n 2 n i=1 1 b 2d 1,i b 2d 2,i b d 1,i b d 2,i K(u, v)p s (x -ub 1,i , y -vb 2,i )dudv 2 ≤ C n b d 1,n b d 2,n where C = p s L ∞ (suppK+1) K 2 (u, v)dudv, and
V 2 (n) = 1 n 2 n i=1 j>i 1 b d 1,i b d 2,i 1 b d 1,j b d 2,j (3.4) Cov µ K x -X 0 b 1,i , y -Y 0 b 2,i , K x -X j-i b 1,j , y -Y j-i b 2,j .
In order to bound the previous quantity, we have two options.

On one hand we may write

E µ K x -X 0 b 1,i , y -Y 0 b 2,i K x -X j-i b 1,j , y -Y j-i b 2,j = b d 1,i b d 2,i b d 1,j b d 2,j K(u, v) K(u ′ , v ′ ) p s (x -ub 1,i , y -vb 2,i ) p(j -i, (x -ub 1,i , y -vb 2,i ), (x -u ′ b 1,j , y -v ′ b 2,j )) dudvdu ′ dv ′ so that, using (2.4), 1 b d 1,i b d 2,i 1 b d 1,j b d 2,j E µ K x -X 0 b 1,i , y -Y 0 b 2,i K x -X j-i b 1,j , y -Y j-i b 2,j ≤ C .
At the same time,

E µ K x -X 0 b 1,i , y -Y 0 b 2,i = b d 1,i b d 2,i K(u, v) p s (x -ub 1,i , y -vb 2,i ) dudv ≤ C b d 1,i b d 2,i so that 1 b d 1,i b d 2,i 1 b d 1,j b d 2,j Cov µ K x -X 0 b 1,i , y -Y 0 b 2,i , K x -X j-i b 1,j , y -Y j-i b 2,j ≤ C .
On the other hand, we may use (2.3) in order to get for any 2 ≤ q ≤ +∞,

Cov µ K x -X 0 b 1,i , y -Y 0 b 2,i , K x -X j-i b 1,j , y -Y j-i b 2,j ≤ C b d/q 1,i b d/q 2,i b d/q 1,j b d/q 2,j ρ (j-i)(q-2)/2q .
Using roughly the monotonicity of the sequences b .,i , we thus get

V 2 (n) ≤ C n 2 n i=1 j>i min   1, ρ (j-i)(q-2)/2q b 2d(1-1 q ) 1,n b 2d(1-1 q ) 2,n   ≤ C n 2 n k=1 (n -k) min   1, ρ k(q-2)/2q b 2d(1-1 q ) 1,n b 2d(1-1 q ) 2,n   .
The sum splits into two terms, depending on which term reaches the minimal value. But roughly, for all 1 ≤ k n ≤ n, we have, with q = +∞ for simplicity, a bound of the form

V 2 (n) ≤ C k n n + ρ kn/2 n b 2d 1,n b 2d 2,n = o 1 n b d 1,n b d 2,n provided ln(1/b 1,n ) + ln(1/b 2,n ) ≪ k n ≪ 1 b d 1,n b d 2,n
, hence the result.

The previous proof gives us the feeling of the right normalization for a Central Limit Theorem. Let us state the result Theorem 3.5. Assume that H 1 or H 2 is fulfilled. Assume in addition that the following assumption

H K is satisfied (H K ) there exists m ∈ N * such that for all multi-index (m 1 , m 2 ) ∈ N d × N d such that m ≥ |m 1 | + |m 2 | ≥ 1 , u m 1 v m 2 K(u, v)dudv = 0 . Let b 1,n and b 2,n be non increasing sequences of real numbers satisfying (1a) b 1,n and b 2,n go to 0 as n → +∞, (2a) nb d 1,n b d 2,n → +∞ as n → +∞, (3a) there exists Σ > 0 such that lim n→+∞ 1 n n i=1 b d 1,n b d 2,n b d 1,i b d 2,i = Σ. (4a) nb d 1,n b d 2,n 1 n n i=1 max(b 1,i , b 2,i ) m+1 → 0 as n → +∞.
Then, in the stationary regime

nb d 1,n b d 2,n (p r,n (x, y) -p s (x, y)) D ----→ n→+∞ N 0, Σ p s (x, y) K 2 (s, t)dsdt .
Remark 3.6. Typical allowed bandwidths are b j,n = n -α j for j = 1, 2 with α j > 0, d(α 1 + α 2 ) < 1 and then choosing m such that 1 < d(α 1 + α 2 ) + 2(m + 1) min(α 1 , α 2 ). In this case Σ = 1/(1 + d(α 1 + α 2 )). These conditions are satisfied in particular once

1 2(d + m + 1) < min(α 1 , α 2 ) ≤ max(α 1 , α 2 ) < 1 2d ,
in which case the rate we obtain is at most n τ with τ < m+1 2(d+m+1) . ♦ Remark 3.7. It is always possible to build some kernel K satisfying H K for a given m. First we may take K(z) = Π 2d j=1 K 1 (z j ) and choose K 1 satisfying H K in dimension 1. To build K 1 choose for instance some even function ψ which is C ∞ and compactly supported in [-1, +1], positive on its support. Then consider the space of polynomials R(u), equipped with the L 2 (ψ(u)du) norm, and its completion H(ψ). It is not difficult to see that, for all k ∈ N * , the family 1, u, u 2 , ..., u k is an independent family of vectors for the euclidean norm induced by the L 2 (ψ(u)du) norm. We may thus find some polynomial M with degree m which is orthogonal to all the u j for j = 1, ..., m and satisfies < M,

1 >= 1 with < •, • > the usual scalar product in L 2 (ψ(u)du). It remains to choose K 1 = M ψ.
The previous shows that we may choose m as large as we want.

n b d 1,n b d 2,n V 1 (n) . But n b d 1,n b d 2,n V 1 (n) = V 11 (n) + V 12 (n) with V 11 (n) = 1 n n i=1 b d 1,n b d 2,n b d 1,i b d 2,i p s (x, y) K 2 (u, v)dudv , while V 12 (n) = 1 n n i=1 b d 1,n b d 2,n b d 1,i b d 2,i K 2 (u, v)(p s (x -ub 1,i , y -vb 2,i ) -p s (x, y)) dudv .
Using the regularity of p s (here the local Lipschitz property) we see that

V 12 (n) ≤ C b 1,n b 2,n n n i=1 b d-1 1,n b d-1 2,n b d-1 1,i b d-1 2,i ≤ C b 1,n b 2,n
thanks to the monotonicity of the sequences, and the latter goes to 0. Finally

lim n→+∞ n b d 1,n b d 2,n Var µ (p r,n (x, y)) = Σ p s (x, y) K 2 (u, v)dudv . (3.8)
To prove the Central Limit Theorem it remains to use the Lindeberg method. As shown in [START_REF] Bardet | Dependent Lindeberg Central Limit Theorem and some applications[END_REF] Theorem 1, the proof amounts to check two conditions. Define

U i,n = b d 1,n b d 2,n √ n b d 1,i b d 2,i K x -X i b 1,i , y -Y i b 2,i
then what is required is (1) there exists some δ > 0 such that

lim n→+∞ A n := lim n→+∞ n i=1 E µ (|U i,n | 2+2δ ) = 0 ,
(2) for all a ∈ R,

lim n→+∞ T n := lim n→+∞ n i=1
Cov µ e ia(U 1,n +...+U i-1,n ) , e iaU i,n = 0 .

We have

A n = 1 n 1+δ n i=1 b d(1+δ) 1,n b d(1+δ) 2,n b 2d(1+δ) 1,i b 2d(1+δ) 2,i E µ K 2+2δ x -X 0 b 1,i , y -Y 0 b 2,i ≤ C n 1+δ n i=1 b d(1+δ) 1,n b d(1+δ) 2,n b d(1+2δ) 1,i b d(1+2δ) 2,i ≤ C n δ b dδ 1,n b dδ 2,n 1 n n i=1 b d 1,n b d 2,n b d 1,i b d 2,i
so that, under our assumptions, A n → 0 as n → +∞.

As in [START_REF] Bardet | Dependent Lindeberg Central Limit Theorem and some applications[END_REF], we write

Cov µ e ia(U 1,n +...+U i-1,n ) , e iaU i,n = i-1 j=1
Cov µ e ia(U 1,n +...+U j,n )e ia(U 1,n +...+U j-1,n ) , e iaU i,n

(3.9) (as usual empty sums are equal to 0). Consider a random variable U * i,n independent from (U 1,n , . . . , U j-1,n ), with same distribution as U i,n . Then, following [START_REF] Bardet | Dependent Lindeberg Central Limit Theorem and some applications[END_REF], we note that each element in the sum of the right hand term of Equality (3.9) is equal to Cov µ e ia(U 1,n +...+U j,n )e ia(U 1,n +...+U j-1,n ) , e iaU i,ne iaU * i,n which is bounded, arguing as in the proof of Lemma 3.2, by

C(a) min E µ (|U j,n ||U i,n |), ρ i-j/2 b d 1,n b d 2,n n b d 1,j b d 2,j b d 1,i b d 2,i ≤ C(a) b d 1,n b d 2,n n min 1, ρ i-j/2 b d 1,i b d 2,i b d 1,j b d 2,j . T n ≤ C(a) n i=1 i-1 j=1 b d 1,n b d 2,n n min 1, ρ i-j/2 b d 1,i b d 2,i b d 1,j b d 2,j ≤ C(a) n k=1 (n -k) min b d 1,n b d 2,n n , ρ k/2 nb d 1,n b d 2,n ≤ C(a)(k(n) b d 1,n b d 2,n + ρ k(n) b d 1,n b d 2,n ) so that choosing again 1 b d 1,n b d 2,n ≫ k(n) ≫ ln(1/b 1,n ) + ln(1/b 2,n )
we obtain that T n → 0.

Remark 3.10. As we have seen H K is only required to control the bias term. Lemma 3.2 is essentially contained in Amiri (2010) chapter 2, where a still more complicated estimator is studied. For the Central Limit Theorem, usual proofs in the recursive framework are using intricate block decompositions. Here we follow the beautiful and now classical Lindeberg method (using Rio's decomposition). The calculations are done in [START_REF] Bardet | Dependent Lindeberg Central Limit Theorem and some applications[END_REF]

) := 1 n n i=1 1 b d 1,i b d 2,i K x -X i b 1,i , y -Y i b 2,i ∆ 2 (X i+2h i , h i ) . (3.11)
As for the density, using Cesaro's rule, it is easy to see that Hr,n is an asymptotically unbiased estimator of b(x, y)p s (x, y), provided the bandwidths b .,n and h n go to 0 at infinity. Next it is enough to closely follow the proof of Proposition 4.2 in Cattiaux et al. (2014-b) keeping in mind that if the times ih n are replaced by i, this does not introduce any modification since what is important is the size of the windows defining ∆ j (X, h i ).

Since some things will differ (in particular the normalization of the kernels), we will give a proof (indicating the differences with

) n h n b d 1,n b d 2,n → +∞ , (3b) there exists θ > 0 such that 1 n n i=1 hn b d 1,n b d 2,n h i b d 1,i b d 2,i → θ, (4b) nb d 1,n b d 2,n h n 1 n n i=1 max(b 1,i , b 2,i ) m+1 → 0, (5b) there exists some r > 1 such that nb d 1,n b d 2,n h n 1 n n i=1 h i b d 1,i b d 2,i 1 r -1 → 0 , 2b
) there exists some p > 1 such that

hnb d 1,n b d 2,n n n i=1 h 2 i (b d 1,i b d 2,i ) 1 p -2 → 0, Then, in the stationary regime nb d 1,n b d 2,n h n Hr,n (x, y) -b(x, y)p s (x, y) D ----→ n→+∞ σN 0, θ 3 p s (x, y)( K 2 (s, t)dsdt) Id .
Remark 3.13. Assume that h n is such that n i=1 h 2 i < +∞. Then for (6) to be satisfied, it suffices that

hn(b d 1,n b d 2,n ) 1/p n b d 1,n b d 2,n
→ 0, which is automatic, thanks to (2). This situation is not the most interesting. Now assume that b j,n = n -α j and h n = n -β for some positive α j and β. It is not difficult to see that, choosing p as close to 1 as necessary, ( 6) is again automatically satisfied, whatever the positive value of β. For (1) to ( 5) to be satisfied, it is then enough to assume

β + d(α 1 + α 2 ) < 1 < β + d(α 1 + α 2 ) + 2 min(β, (m + 1) min(α 1 , α 2 )) ,
and in this case we have θ = 1/(1 + β + d(α 1 + α 2 )). As we have seen in Remark 3.7, we may choose m as large as we want, so that the previous conditions amounts to

β + d(α 1 + α 2 ) < 1 < 3β + d(α 1 + α 2 ) .
(3.14)

Hence we may choose the α j 's as small as possible and β larger than but close to 1/3 to get the quasi optimal rate n 1/3 . ♦ Proof. We give the main modifications w.r. Define

I i := i+3h i i+2h i t i+2h i b(X s , Y s )dsdt + i+2h i i+h i i+2h i t b(X s , Y s )dsdt , and
W i := i+3h i i+2h i (W s -W i+2h i )ds + i+2h i i+h i (W i+2h i -W s )ds .
The vector nb d 1,n b d 2,n h n Hr,n (x, y) can be decomposed in two terms: the one driving the bias in the central limit theorem

S n,1 (x, y) := nb d 1,n b d 2,n h n 1 n n i=1 1 b d 1,i b d 2,i K x -X i b 1,i , y -Y i b 2,i 1 h 2 i I i ,
and the one driving the variance

S n,2 (x, y) := nb d 1,n b d 2,n h n 1 n n i=1 1 b d 1,i b d 2,i K x -X i b 1,i , y -Y i b 2,i 1 h 2 i σ W i .
Notice that ES n,2 (x, y) = 0. We thus have

I 2n = E µ S n,1 (x, y) -n b d 1,n b d 2,n h n b(x, y)p s (x, y) , (3.15) while I 1n = (S n,1 -E µ S n,1 (x, y)) + S n,2 (x, y) .
(3.16)

First step: Showing that I 2n → 0.

We define

P i := i+3h i i+2h i t i+2h i (P s b(X i , Y i ) -b(X i , Y i ))dsdt + i+2h i i+h i i+2h i t (P s b(X i , Y i ) -b(X i , Y i ))dsdt .
Thanks to stationarity, it holds

I 2n = nb d 1,n b d 2,n h n E µ 1 n n i=1 1 b d 1,i b d 2,i K x -X 0 b 1,i , y -Y 0 b 2,i 1 h 2 i P 0 + + nb d 1,n b d 2,n h n E µ 1 n n i=1 1 b d 1,i b d 2,i K x -X 0 b 1,i , y -Y 0 b 2,i b(X 0 , Y 0 ) -b(x, y)p s (x, y) .
The second summand in the above expression can be treated as for the density and goes to 0 as soon as

nb d 1,n b d 2,n h n 1 n n i=1 max(b 1,i , b 2,i ) m+1 → 0 .
The first summand can be decomposed into a sum

nb d 1,n b d 2,n h n n i=1 (A 1i + A 2i ) with A 1i = 1 h 2 i 3h i 2h i t 2h i (P s b(u, v) -b(u, v)) 1 b d 1,i b d 2,i K x -u b 1,i , y -v b 2,i ds dt µ(du, dv) ,
A 2i being similar just changing

3h i 2h i t 2h i into 2h i h i 2h i t
. We thus only study A 1i for which we have (|.| denoting the norm in R d ) as in Cattiaux et al. (2014-b)

|A 1i | ≤ CM p h i 1 b d 1,i b d 2,i K r x -u b 1,i , y -v b 2,i µ(du, dv) 1/r ≤ CM p h i b d 1,i b d 2,i 1 r -1 , with p > 1 , r < +∞ ∈ N * such that 1 p + 1 r = 1.
Accordingly the first summand goes to zero as n tends to infinity, provided

nb d 1,n b d 2,n h n 1 n n i=1 h i b d 1,i b d 2,i 1 r -1 → 0 .
Second step: Study of S n,2

We now consider the term driving the variance. As in the proof of Proposition 4.2 in Cattiaux et al. (2014-b) one can study the behavior of the characteristic function Φ n (t) of S n,2 , recall that

S n,2 (x, y) := nb d 1,n b d 2,n h n 1 n n i=1 1 b d 1,i b d 2,i K x -X i b 1,i , y -Y i b 2,i σ W i h 2 i . and that ES n,2 (x, y) = 0. If < •, • > denotes the usual scalar product in R d , introduce the random function f n (t) (t ∈ R d ) f n (t) = e -1 6n |σt| 2 b d 1,n b d 2,n hn n k=1 1 b 2d 1,k b 2d 2,k h k K 2 x-X k b 1,k , y-Y k b 2,k
. One can follow the method used in Cattiaux et al. (2014-b) to show that the convergence of Φ n amounts to the convergence in Probability of f n . Now define

Z n = b d 1,n b d 2,n h n 1 n n k=1 1 b 2d 1,k b 2d 2,k h k K 2 x -X k b 1,k , y -Y k b 2,k . Let A = θ p s (x, y) K 2 (u, v)dudv ,
where

θ = lim n→+∞ 1 n n k=1 b d 1,n b d 2,n h n b d 1,k b d 2,k h k . Then E µ (|Z n -A|) ≤ 1 n n k=1 b d 1,n b d 2,n h n b d 1,k b d 2,k h k K 2 (u, v) |p s (x -ub 1,k , y -vb 2,k ) -p s (x, y)|dudv + 1 nθ n k=1 b d 1,n b d 2,n h n b d 1,k b d 2,k h k -1 A ≤ 1 n n k=1 K 2 (u, v) |p s (x -ub 1,k , y -vb 2,k ) -p s (x, y)|dudv + 1 nθ n k=1 b d 1,n b d 2,n h n b d 1,k b d 2,k h k -1 A ,
using the monotonicity properties of all sequences. This proves, using Cesaro's rule again that E µ (|Z n -A|) → 0. Thus Z n → θ p s (x, y) K 2 (u, v)dudv, in L 1 (µ). Using the bounded convergence theorem, we deduce that f n (t)

P → e -θ|σt| 2 6 ps(x,y) K 2 (u,v)dudv , so that we get that

S n,2 (x, y) D ----→ n→+∞ σ N 0, θ 3 p s (x, y)( K 2 (s, t)dsdt) Id .
Third step: showing that S n,1 (x, y)) -E µ (S n,1 (x, y)) → 0 Let us denote by P k i the k th coordinate of the vector P i, . Defining

Γ k (i, x, y, X, Y ) = 1 b d 1,i b d 2,i K( x -X b 1,i , y -Y b 2,i ) 1 h 2 i P k i ,
we have

S k n,1 = nh n b d 1,n b d 2,n 1 n n i=1 Γ k (i, x, y, X i , Y i ), so that n b d 1n b d 2n hn Var µ (S k n,1 ) = n i=1 Var µ (Γ k (i, x, y, X i , Y i )) + i =l Cov µ (Γ k (i, x, y, X i , Y i ), Γ k (l, x, y, X l , Y l )) .
To bound the above expression we first bound as we did for the first step

E µ ((Γ k (i, x, y, X i , Y i )) 2 ) ≤ ≤ C h 4 i b 2d 1,i b 2d 2,i E µ K 2 x -X 0 b 1,i , y -Y 0 b 2,i h i 0 (h i -s) (P s b(X 0 , Y 0 ) -b(X 0 , Y 0 )) ds 2 ≤ C h i b 2d 1,i b 2d 2,i h i 0 K 2 x -u b 1,i , y -v b 2,i (P s b(u, v) -b(u, v)) 2 ds dµ ≤ C h i b 2d 1,i b 2d 2,i h i 0 K 2 x -u b 1,i , y -v b 2,i P s b(u, v) -b(u, v) s 2 dµ ds
We may argue as in the first step, this time using Hölder inequality for some conjugate pair (p, q) to conclude that

E µ ((Γ k (i, x, y, X i , Y i )) 2 ) ≤ C h 2 i (b 1,i b 2,i ) (d/p)-2d . (3.17) It follows h n b d 1,n b d 2,n n n i=1 Var µ (Γ k (i, x, y, X i , Y i )) ≤ C h n b d 1,n b d 2,n n n i=1 h 2 i (b d 1,i b d 2,i ) 1 p -2 ,
the latter going to 0 thanks to (6).

Let us now compute the covariances. One has thanks to stationarity and mixing inequality (2.2)

(remark that b d 1,i b d 2,i Γ k (i, .) is bounded) b d 1,n b d 2,n h n n i =l Cov µ Γ k (i, x, y, X i , Y i ), Γ k (l, x, y, X l , Y l ) ≤ ≤ C b d 1,n b d 2,n h n n n i=1 i-1 j=1 min 1, ρ i-j/2 b d 1,i b d 2,i b d 1,j b d 2,j ≤ Ch n n-1 j=1 min b d 1,n b d 2,n , ρ j/2 b d 1,n b d 2,n ≤ C h n k n b d 1,n b d 2,n , ρ kn/2 b d 1,n b d 2,n
which goes to 0 for an appropriate choice of k n .

Putting together the previous two theorems, we may deduce exactly as in Cattiaux et al. (2014-b) proof of theorem 4.1 the following Corollary 3.18. Assume that all the assumptions of Theorem 3.5 and Theorem 3.12 are satisfied. Then, defining br,n = Hr,n /p r,n we have in the stationary regime

nb d 1,n b d 2,n h n br,n (x, y) -b(x, y) D ----→ n→+∞ σN 0, θ 3p s (x, y) ( K 2 (s, t)dsdt) Id .
Remark 3.19. In the setting of Remark 3.13 the quasi optimal rate is still n 1/3 . ♦ 4. Estimation of the invariant density and of the drift term in the stationary regime: partial observation case.

We turn to the study of the estimators pr,n and br,n defined in (1.3) and (1.4).

Estimation of the invariant density.

If we can only observe the position process, we have to consider instead of Y i its natural approximation

X i+h i -X i h i . Recall our estimator pr,n (x, y) := 1 n n i=1 1 b d 1,i b d 2,i K x -X i b 1,i , y - X i+h i -X i h i b 2,i .
Theorem 4.1. Assume that H 1 (or H 2 ) and H K are fulfilled. Assume in addition that conditions (1a) up to (6a) are satisfied, where (1a) up to (4a) are given in Theorem 3 and

(5a) lim n→+∞ nb d 1,n b d 2,n 1 n n i=1 √ h i b d+1 2,i = 0 , (6a) for some p > 1, lim n→+∞ nb d 1,n b d 2,n 1 n n i=1 h i b d+1 2,i b d-(d/p) 1,i = 0 .
Then in the stationary regime

nb d 1,n b d 2,n (p r,n (x, y) -p s (x, y)) D ----→ n→+∞ N 0, Σ p s (x, y) K 2 (s, t)dsdt .
Proof. We want to evaluate the expectation of the difference :

nb d 1,n b d 2,n (p s (x, y) -ps (x, y)) (4.2) = b d 1,n b d 2,n √ n n i=1 1 b d 1,i b d 2,i K x -X i b 1,i , y - X i+h i -X i h i b 2,i -K x -X i b 1,i , y -Y i b 2,i , = b d 1,n b d 2,n √ n n i=1 1 b d 1,i b d 2,i A i .
We will closely follow Step 2 in the proof of Theorem 3.3 of Cattiaux et al. (2014-a).

Introduce M u = 1 h u i (Y s -Y i )ds, defined for i ≤ u ≤ i + h i .

Then we may write

A i = -1 h i b 2,i i+h i i ∇ y K x -X i b 1,i , y -Y i -M u b 2,i .(Y u -Y i )du .
Arguing as in the referred proof we can show that for all conjugate pair (p, q) with 1 < p < +∞,

E µ (|A i |) ≤ C h i b 2,i h 3/2 i b d 1,i + b d/p 1,i h 2 i . Hence E µ nb d 1,n b d 2,n |p r,n (x, y) -pr,n (x, y)| ≤ b d 1,n b d 2,n √ n n i=1 1 b d 1,i b d 2,i E µ (|A i |) ≤ C b d 1,n b d 2,n √ n n i=1 h 3/2 i b d 1,i + b d/p 1,i h 2 i b d 1,i b d+1 2,i h i .
It follows that nb d 1,n b d 2,n (p r,n (x, y)pr,n (x, y)) goes to 0 in L 1 as soon as (5a) and (6a) are satisfied. The proof follows from Slutsky's lemma and Theorem 3.

Remark 4.3. Assume that b j,n = n -α j and h n = n -β . For (1a) up to (4a) to be satisfied, we have seen in Remark 3.6 that the condition

d(α 1 + α 2 ) < 1 < d(α 1 + α 2 ) + 2(m + 1) min(α 1 , α 2 ) (4.4)
is sufficient and in Remark 3.7 that we can choose m as large as we want, so that the previous amounts to d(α 1 + α 2 ) < 1.

Since we can choose p as close to 1 as we want, we easily see that (6a) is automatically satisfied as soon as (5a) is satisfied. The latter implies

1 + α 2 (d + 2) < β + dα 1 . (4.5)
It follows that, if we are only interested in estimating the invariant density, we may take β large enough for the latter to be satisfied. But of course we are interested in the simultaneous estimation of both the invariant density and the drift term, so that (4.5) will become some restrictive condition, we shall study in more details at the end of the next subsection. ♦ 4.2. Estimation of the drift.

Introduce

Ĥr,n (x, y)

:= 1 n n i=1 1 b d 1,i b d 2,i K x -X i b 1,i , y - X i+h i -X i h i b 2,i ∆ 2 (X i+2h i , h i ) . (4.6)
We shall first prove the analogue of Theorem 3.12.

Theorem 4.7. Assume that H 1 (or H 2 ) and H K are fulfilled. Assume in addition that conditions (1b) up to (8b) are satisfied, where (1b) up to (6b) are given in Theorem 3.12 and

(7b) for some ε ≥ 1, lim n→+∞ nh n b d 1,n b d 2,n 1 n n i=1 √ h i b d+1 2,i b dε/(1+ε) 1,i = 0 , ( 8b
) for some ε > 0 and some r > 1 + ε,

lim n→+∞ nh n b d 1,n b d 2,n 1 n n i=1 h i b d+1 2,i b d(r-1)/r 1,i = 0 .
Then in the stationary regime

nh n b d 1,n b d 2,n Ĥr,n (x, y) -p s (x, y)b(x, y) D ----→ n→+∞ σN 0, θ 3 p s (x, y)( K 2 (s, t)dsdt) Id .
Proof. Starting from Theorem 3.12, as in the previous subsection, it remains to consider D n defined by

D n := nh n b d 1,n b d 2,n 1 n n i=1 1 b d 1,i b d 2,i A i ∆ 2 (X i+2h i , h i ) ,
where

A i = K x -X i b 1,i , y -Y i b 2,i -K x -X i b 1,i , y - X i+h i -X i h i b 2,i ,
as in the previous subsection.

As before we decompose ∆ 2 (X i+2h i , h i ) = 1 h 2 i (I i + M i ), I i being the bounded variation part and M i the martingale part. More precisely we define

M i = σ i+3h i i+2h i (W s -W i+2h i )ds + i+2h i i+h i (W i+2h i -W s )ds , I i = i+3h i i+2h i t i+2h i g(X s , Y s )dsdt + i+2h i i+h i i+2h i t g(X s , Y s )dsdt .
We want now to control

M n = E µ 1 n n i=1 1 h 2 i b d 1,i b d 2,i |A i ||I i | and N n = E µ 1 n n i=1 1 h 2 i b d 1,i b d 2,i A i M i .
For the first term, we use Hölder inequality for some ε > 0, (1+ε) .

E µ (|A i I i |) ≤ E µ |A i | 1+ε 1/(1+ε) E µ |I i | (1+ε)/ε ε/
Similarly to what we have done in the previous proof we may show on one hand that

E µ (|A i | 1+ε ) ≤ C(ε) h i b 1+ε 2,i (h (3+ε)/2 i b d 1,i + h 2+ε i b d/p 1,i ) . (4.8)
On the other hand, according to the proof in Cattiaux et al. (2014-b) (see formula (7.7) therein), one has

E µ |I i | (1+ε)/ε ≤ C(ε, u, v) h 2(1+ 1 ε ) i . Finally, E µ (|A i I i |) ≤ C(ε, u, v) h 1+2ε 1+ε i b 2,i (h (3+ε)/2 i b d 1,i + h 2+ε i b d/p 1,i ) 1 1+ε ≤ C(ε, u, v) h 5/2 i b d/(1+ε) 1,i b 2,i + h 3 i b d/p(1+ε) 1,i b 2,i .
Hence nh n b d 1,n b d 2,n M n → 0 according to (7b) and (8b).

It remains now to bound N n .

To this end we compute the expectation of the square. But the terms A i M i are centered and uncorrelated. It follows (similarly to Cattiaux et al. (2014-b)) that

N 2 n ≤ C n 2 n i=1 1 h 4 i b 2d 1,i b 2d 2,i h 3 i E µ (|A i | 2 ) ≤ C n 2 n i=1 1 h 4 i b 2d 1,i b 2d 2,i h 3 i 1 h i b 2 2,i (h 2 i b d 1,i + h 3 i b d/p 1,i )
where we have used (4.8) with ε = 1, and finally

N n ≤ C n n i=1 1 b d/2 1,i b d+1 2,i + n i=1 √ h i b d+1 2,i b d(1-(1/2p)) 1,i . Again nh n b d 1,n b d 2,n
N n → 0 according to (7b) and (8b).

Remark 4.9.

If lim n→+∞ nh n b d 1,n b d 2,n 1 n n i=1 h κ 0 i b κ 2 2,i b κ 1 1,i
= 0 for some positive κ j 's, then the same holds for all η j 's with η 0 ≥ κ 0 , η 1 ≤ κ 1 and η 2 ≤ κ 2 . Hence it is easy to see that (7b) and (8b) amounts to (7-8-b') there exists r > 1 2 such that lim

n→+∞ nh n b d 1,n b d 2,n 1 n n i=1 √ h i b d+1 2,i b rd 1,i = 0 .
Now assume as usual that b j,n = n -α j and h n = n -β for some positive α j and β.

(7-8-b') is satisfied as soon as 1 + (d + 2)α 2 < 2β . In particular, for all the conditions in Theorem 4.7 and Theorem 4.1 to be satisfied it is sufficient that both conditions (3.14) and 1 + (d + 2)α 2 < β + β ∧ dα 1 are fulfilled.

Unfortunately, as for our estimators in Cattiaux et al. (2014-b), both conditions are incompatible. If we do not want to multiply the observations, we have to choose the same discretization step h n in both Ĥr,n and pr,n . But as we did in Cattiaux et al. (2014-b) Theorem 4.1 and proposition 4.4, we may choose different bandwidths for both estimators, denoted by b p and b H for respectively the density estimator and the drift estimator. The required conditions are then, on one hand

d(α p 1 + α p 2 ) < 1 and 1 + α p 2 (d + 2) < β + β ∧ dα p 1 , (4.10)
and, on the other hand

β + d(α H 1 + α H 2 ) < 1 < 3β + d(α H 1 + α H 2 ) and 1 + α H 2 (d + 2) < 2β + dα H 1 . (4.11) We also want β + d(α H 1 + α H 2
) to be as small as possible. But it is easily seen that this quantity has to be larger than 1/2. We will describe a typical situation in the next Corollary.

♦

According to the previous discussion we have Corollary 4.12. Assume that H 1 (or H 2 ) and H K are fulfilled. Assume in addition that b H j,n = n -α H j (resp. b p j,n = n -α p j ) and h n = n -β for some positive α j 's and β. Denote by Ĥr,n (resp. pr,n ) the estimators built with the bandwidths b H j,n (resp. b p j,n ) and the discretization step h n . Finally assume that all the parameters satisfy (4.10) and (4.11) and that m min(α H 1 , α H 2 , α p 1 , α p 2 ) > 1 . Then, defining br,n = Ĥr,n /p r,n we have in the stationary regime

nb d 1,n b d 2,n h n br,n (x, y) -b(x, y) D ----→ n→+∞ σN 0, θ 3p s (x, y) ( K 2 (s, t)dsdt) Id , where θ = 1/(1 + β + d(α H 1 + α H 2 )).
Remark 4.13. In the previous Corollary we may choose any β > 1/2, then dα p 1 > 1/2, then d max(α H 1 , α H 2 , α p 2 ) sufficiently small and finally m large enough for all the assumptions of the previous Corollary to be satisfied. It yields a quasi optimal rate of order n 1/4 .

Notice that if we want to estimate simultaneously p s and b, we may use two estimators of p s , say pp r,n and pH r,n , both with h n = n -β for some β > 1/2, the first one with bandwidths α p j , the second one with α H j and look at (p H r,n , ( Ĥr,n /p p r,n )) simultaneously. ♦

5.

Additional results and comments on the previous sections.

5.1. Non stationary case. This subsection will be very short and can be reduced to the following sentence : "All the results of the two previous sections are still true if the initial distribution is either a Dirac mass or is absolutely continuous with respect to µ." Indeed the situation here is very simple: if we replace n i=1 by n i=ln with (l n /n) → 0, we obtain exactly the same results as the first l n terms go to 0 starting from any initial distribution. Next we have to compare the law of Z ln to µ. But thanks to geometric ergodicity (or exponential mixing see Proposition 2.1), it is known that the total variation distance between both is less than Cρ ln provided the initial distribution satisfies the conditions we have mentioned (see e.g. [START_REF] Bakry | Rate of convergence for ergodic continuous Markov processes : Lyapunov versus Poincaré[END_REF]). The conclusion is then straightforward using Markov property and l n = c log(n).

5.2.

Comparison with existing results. The rate of convergence for the recursive estimators are much better than the ones we obtain in Cattiaux et al. (2014-a) and Cattiaux et al. (2014-b) (see e.g Proposition 4.4 in Cattiaux et al. (2014-b)). For the density the main reason is that we do not need the small time bounds for the transition kernel as in Cattiaux et al. (2014-a) (see [START_REF] Konakov | Explicit parametrix and local limit theorems for some degenerate diffusion processes[END_REF], [START_REF] Cattiaux | Calcul stochastique et opérateurs dégénérés du second ordre. I Résolvantes, théorème de Hörmander et applications[END_REF], [START_REF] Cattiaux | Hypoellipticité et hypoellipticité partielle pour les diffusions avec une condition frontière[END_REF]), for the unnormalized drift estimator Ĥr,n the gain is not as important. In comparison with the existing literature on recursive estimation for the density (see e.g. Amiri (2010) for a review), the results in Section 3 are similar. Of course in many quoted references the authors also studied almost sure results in the spirit of the log-log law. Again see Amiri ( 2010) for a review.

The assumption made in all these references is that the transition kernel p s,t (., .) satisfies sup |t-s|≥1 p s,t (., .)p s ⊗ p s ∞ < +∞ .

In our situation, we do not know about such uniform bound, though it is presumably true.

Estimation of the variance.

We turn to the estimation of σ. To this end we merely follow Section 4 in Cattiaux et al. (2014-c), though here we only consider the σ constant case, which introduces various simplifications. Recall that σr,n := 3 2n

n i=1 1 h i ∆ 2 (X i+h i , h i ) >< ∆ 2 (X i+h i , h i ) .
Theorem 6.1. Assume that H 1 (or H 2 ) and H K are fulfilled. Also assume that σ is symmetric (i.e. we take the symmetric square root of σ * σ). Assume in addition that, as

n → +∞, 1 √ n n i=1 h i → 0 . (6.2)
Then, in the stationary regime or in the situation of subsection 5.1

√ n (σ r,n -σ 2 ) D ----→ n→+∞ σN d,d σ ,
where N d,d is a symmetric gaussian random matrix whose entries N k,l d,d are N (0, 1 + δ k,l ), all entries for k ≥ l being independent. Remark 6.3. Condition (6.2) is satisfied when h n = n -β with β > 1/2. ♦ Proof. In the sequel the indices (k, l) denote respectively the index of the row and of the column of a matrix.

(SS i + SV i + V V i )
where

SS k,l i = 1 h 3 i i+2h i i (h i -|u -i -h i |) u i (h i -|s -i -h i |)((σdW s ) k (σdW u ) l + (σdW s ) l (σdW u ) k ) SV k,l i = 1 h 3 i i+2h i i (h i -|u -i -h i |) u i (h i -|s -i -h i |)((σdW s ) k b l (Z u )du + (σdW s ) l b k (Z u )du) + 1 h 3 i i+2h i i (h i -|u -i -h i |) u i (h i -|s -i -h i |)((σdW u ) k b l (Z s )ds + (σdW u ) l b k (Z s )ds) V V k,l i = 1 h 3 i i+2h i i (h i -|u -i -h i |) u i (h i -|s -i -h i |)(b k (Z s )ds b l (Z u )du + b l (Z s )ds b k (Z u )du) .
It is immediate that the random matrices (M i = 3 2 σ -1 SS i σ -1 ) i are independent, with centered entries such that Var µ (M k,l i ) = 1 + δ k,l . Hence according to the usual Central Limit Theorem, It remains to prove that the remaining terms go to 0 in L 1 . First, it is not hard to see that

E µ (|V V k,l i |) ≤ C h i so that 1 √ n n i=1 E µ (|V V k,l i |
) → 0 according to (6.2). Next, to control SV i we are led to control two types of terms: either martingale terms

N k,l i = 1 h 3 i i+2h i i (h i -|u -i -h i |) β l i (u) (σdW u ) k
or bounded variation terms

B k,l i = 1 h 3 i i+2h i i (h i -|u -i -h i |) γ l i (u) b k (Z u )du .
The martingale terms are centered and to prove the convergence of their normalized sum to 0 in L 1 (µ), it is enough to prove the convergence to 0 of the variance. But for i = j the covariance Cov µ (N k,l i , N k,l j ) = 0 so that the variance

Var µ 1 √ n n i=1 N k,l i ≤ C n n i=1 1 h 6 i i+2h i i (h i -|u -i -h i |) 2 E µ ((β l i (u)) 2 ) du ≤ C n n i=1 h i
goes to 0 provided h n goes to 0, hence we get convergence to 0 in L 1 . Here we have used

E µ i+2h i i θ(u)du 2 ≤ C h i E µ i+2h i i θ 2 (u)du .
B k,l i,1 = 1 h 3 i i+2h i i (h i -|u -i -h i |) γ l i (u) b k (Z i )du ,
are thus centered and satisfy Cov µ (B k,l i,1 , B k,l j,1 ) = 0 for i = j. As before it follows

Var µ 1 √ n n i=1 B k,l i,1 ≤ C n n i=1 1 h 5 i i+2h i i (h i -|u -i -h i |) 2 E µ ((γ l i (u)b k (Z i )) 2 ) du . But E µ ((γ l i (u)b k (Z i )) 2 ) = E µ (b k (Z i )) 2 u i (h i -|s -i -h i |) 2 c(σ)
B k,l i,2 = 1 h 3 i i+2h i i (h i -|u -i -h i |) γ l i (u) (b k (Z u ) -b k (Z i ))du .
But, according to what we have done before

E µ (γ l i (u) (b k (Z u ) -b k (Z i ))) ≤ (E µ ((γ l i (u)) 2 )) 1/2 (E µ ((b k (Z u ) -b k (Z i )) 2 )) 1/2 ≤ Ch 3/2 i h 1/2 i . It follows that 1 √ n E µ n i=1 |B k,l i,2 | ≤ 1 √ n n i=1 h i ,
which goes to 0 thanks to (6.2).

♦

  Proof. First, the bias term B n = n b d 1,n b d 2,n (E µ (p r,n (x, y))p s (x, y)) can be treated exactly as in Cattiaux et al. (2014-a) section 5 step 3. It is for this term that the conditions on K and (4a) are required. Next we show that n b d 1,n b d 2,n Var µ (p r,n (x, y)) converges as n → +∞. According to the proof of the previous Lemma, we see that this convergence amounts to the convergence of

  t. the proof of Proposition 4.2 in Cattiaux et al. (2014-b). We decompose S n := n b d 1,n b d 2,n h n Hr,n (x, y)b(x, y)p s (x, y) = n b d 1,n b d 2,n h n Hn (x, y) -E Hr,n (x, y) + E Hr,n (x, y)b(x, y)p s (x, y) =: I 1n + I 2n .

For

  the bounded variation term we use the same trick as inCattiaux et al. (2014-c) writing b k (Z u ) = (b k (Z u )b k (Z i )) + b k (Z i ). The terms

  Cattiaux et al. (2014-b)) of the following statement Theorem 3.12. Assume that H 1 or H 2 , and H K are fulfilled. Let h n , b 1,n and b 2,n be non increasing sequences of real numbers satisfying, as n → +∞, (1b) h n < 1/3 for all n, h n , b 1,n and b 2,n go to 0, (

  Using Ito's formula we haveh 2 (∆ 2 (X t , h) >< ∆ 2 (X t , h)) k,l = (X i+h i , h i ) >< ∆ 2 (X i+h i , h i ) -

						t+h
				=	t-h	(h -|u -t|) (H k u (σdW u ) l + H l u (σdW u ) k ) +	(6.4)
						t+h
					+	t-h	(h -|u -t|) (H k u b l (Z u ) + H l u b k (Z u ))du
						t+h
					+	t-h	(h -|u -t|) 2 σ 2 k,l du ,
	with						u
					H u =
							t-h
							(6.5)
	Now we can write
	σr,n -σ 2 = =	3 2n	n i=1	1 h i	i h 2 ∆ 2 1	i+2h i
	=	3 2n	n i=1

(h -|s -t|) (σdW s + b(Z s )ds) .

As explained in

Cattiaux et al. (2014-c)

formula (4.5) applied in our situation, it holds for all j ∈ N and s between th and t + h,

E µ (|H s | 2j ) ≤ C(j) (s -(th)) j h 2j . i (h i -|s -(i + h i )|) 2 σ 2 ds


Ahmad Mheich

Mahmoud Hassan

Olivier Dufor

Mohamad Khalil

Claude Berrou

Fabrice Wendling

Spatiotemporal Analysis of Brain Functional Connectivity

Keywords: High resolution EEG, PLV, K-means, functional connectivity dynamics

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

A key challenge in neuroscience is to identify distinctive networks underlying specific brain functions, from data provided by neuroimaging techniques, which can be either structural (diffusion imaging, DTI) or functional (electroencephalography, -EEG magnetoencephalography -MEG-or even functional magnetic resonance imaging fMRI-). These neuroimaging techniques can be used to identify brain networks involved in normal brain functions(behavioral response to stimulus, learning , memory) as well as in neurological disorders like epilepsy, Alzheimer, Schizophrenia [START_REF] Stam | Graph theoretical analysis of magnetoencephalographic functional connectivity in Alzheimer's disease[END_REF].

In this context, fMRI has significantly evolved over the past two decades and is now commonly used to characterize brain connectivity. However, the short duration of cognitive processes (~500ms for the picture naming, for example) requires the use of techniques that have a very high time resolution (on the order of ms), which is not the case of fMRI (~1s), although it has the advantage of having excellent spatial resolution.

Several studies have indicated that the use of electroencephalography (EEG) -which provides s c a l p s i g n a l s a t excellent temporal resolution (sub-millisecond) -combined with a p p r o p r i a t e s i g n a l p r o c e s s i n g t e c h n i q u e s c a n b r i n g relevant information about normal networks activated during cognitive activity [START_REF] Rodriguez | Perception's shadow: long-distance synchronization of human brain activity[END_REF] or about disturbed networks associated for instance with tumors [START_REF] Bartolomei | How do brain tumors alter functional connectivity? A magnetoencephalography study[END_REF].

This excellent temporal resolution allows us to analyze the dynamic properties of the cognitive process, a challenging issue so far addressed in a few studies only. In [START_REF] Brunet | Spatiotemporal analysis of multichannel EEG: CARTOOL[END_REF][START_REF] Murray | Topographic ERP analyses: a step-by-step tutorial review[END_REF], proposed algorithms were based on the Event Related Potentials (ERPs) amplitude. However, these algorithms do not take into account the connectivity between signals (electrode space) or brain regions (source space). Other methods rely on the decomposition of EEG measurements into topographic and time-frequency elements. These methods based on frequency analysis, show that a large set of data can be represented by a small number of topographic distributions [START_REF] Koenig | Brain connectivity at different time-scales measured with EEG[END_REF].Again, these methods ignore the connectivity aspects.

Regarding the results based on the connectivity analysis, most of reported methods make use of a constant time window to track the EEG dynamic connectivity. This window is chosen empirically or based on a priori information about the analyzed task [START_REF] Rodriguez | Perception's shadow: long-distance synchronization of human brain activity[END_REF]. A few attempts have been recently reported to avoid this constraint [START_REF] Dimitriadis | Tracking brain dynamics via time-dependent network analysis[END_REF][START_REF] Allen | Tracking whole-brain connectivity dynamics in the resting state[END_REF].However, most of the proposed algorithms are not adequate to track changes over very short durations (as in the case of responses evoked by visual stimuli).

In this paper, we propose a new algorithm to track the brain functional connectivity dynamics during picture naming task. The proposed algorithm is based on the K-means clustering of the connectivity graphs obtained by the trial by trial Phase Locking Values (PLV) method.

II. MATERIALS AND METTHODS

A. DATA

Six subjects were shown pictures on a screen using E-Prime 2.0 software (Psychology Software Tools, Pittsburgh, PA). They were asked to name the displayed objects. The 148 images were selected from a database of 400 pictures standardized for French [START_REF] Alario | A set of 400 pictures standardized for French: norms for name agreement, image agreement, familiarity, visual complexity, image variability, and age of acquisition[END_REF] and were used during two sessions (about eight minutes each) of 74 stimuli. The brain activity was recorded using an hr-EEG system (256 electrodes, EGI, Electrical Geodesic Inc.). EEG signals were collected at a 1 kHz sampling frequency and were bandpass filtered between 3 and 45 Hz. Each trial was visually inspected, and epochs contaminated by eye blinking, movements or any other noise source were rejected and excluded from the analysis performed using the EEGLAB open source toolbox [START_REF] Delorme | EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis[END_REF]. This study was approved by the National Ethics Committee for the Protection of Persons (CPP), (conneXion study, agreement number 2012-A01227-36, promoter: Rennes University Hospital). We excluded the electrodes located on the face as well as the few electrodes showing too high impedance. Overall, in our hr-EEG montage, 180 (over 256) electrodes were retained as providing excellent quality signals over all subjects.

B. Functional connectivity measure

Functional connectivity is defined as the temporal correlation (wide sense) among the activity of different neuronal assemblies. Several methods have been proposed to characterize the brain functional connectivity. In this paper we used the phase locking value (PLV) method proposed by Lachaux et al. [START_REF] Lachaux | Measuring phase synchrony in brain signals[END_REF]to measure synchronization on a trial by trial basis. For each channel pair, x and y, at time t, and for all the trials (n=1…N), PLV is defined as: [START_REF] Stam | Graph theoretical analysis of magnetoencephalographic functional connectivity in Alzheimer's disease[END_REF] Where and are the unwrapped phases of signals x and y respectively. We applied a normalization procedure so that the values were compared with the 200ms baseline preceding the presentation of the image. Let and be the mean and standard deviation computed from a 200ms pre-stimulus baseline. The normalized PLV is then computed as . The connectivity matrices are computed in the low gamma frequency band (30 Hz -45 Hz). A thresholding procedure is then applied on the functional connectivity values in order to retain a small fraction (10%) of the strongest functional connections, as discussed in [START_REF] Sporns | Structure and function of complex brain networks[END_REF].

C. Segmentation algorithm

The proposed algorithm can be summarized as follows:

1-Compute the adjacency matrices for each subject using PLV to obtain T connectivity graphs. The instantaneous averaged connectivity PLV (t) is computed over all the subjects. 2-Select randomly K graphs (K varies from 3 to-10) from the T graphs. To avoid the choice of very close graphs , an additional constraint has been added which consists in rejecting the K graphs if the time interval between two graphs in any pairs is less than 30 ms. 3-Calculate the spatial correlation (sC) between each of the K and T graphs. sC is computed based on PLV(t). For two pairwise PLV P and Q, at instant t, sC(t) is defined as:

(

) 2
Where M is the number of combinations (180x(180-1)/2). 4-From these spatial correlation values, the Global Explained Variance (GEV) is calculated as defined in [START_REF] Murray | Topographic ERP analyses: a step-by-step tutorial review[END_REF]. While GEV is unstable (i.e. doesn't reach its highest value over the iterations), the template maps are redefined averaging all the graphs yielding to the same cluster. When reaching the highest GEV, K+1 template maps are then selected and all the above procedure is repeated. To choose the optimal number of template graphs, we use a method based on the Cross Validation (CV) criterion as introduced in [START_REF] Murray | Topographic ERP analyses: a step-by-step tutorial review[END_REF] which is the ratio between GEV and the degrees of freedom for a given set of graphs. Its minimum gives the optimal number of segments. Finally, an additional constraint is used to prevent the small clusters, a given connectivity graph must "survive" for at least 30ms. 5-The result obtained in the averaged data is then compared with the moment-by-moment connectivity graphs of individual subjects PLVs. Each time point is labelled according to the graph with which it best correlates, yielding a measure of graph presence. This procedure referred to as 'fitting' allows for establishing how well a cluster map explains individual activity and its duration.

III. RESULTS

To start with, after the random selection of the template graphs, we compute the spatial correlation between these template maps and all the graphs. An example of the sC curves (at each instant) is shown in Fig. 1 for 3 random different instants. The figure shows how sC is able to detect the spatial correlation between graphs. We can notice also that around the highest correlation peak, the correlation values are higher than temporary distant graphs. This step is followed by the computation of GEV for different template graphs. In our study the number of clusters varies from 3 to 10 clusters. The lowest CV value is obtained for a number of clusters equal to 6.

We then choose this number of clusters to identify the spatiotemporal behavior of the connectivity graphs obtained with the algorithm. Results are exemplified in Fig. 2. The figure shows the segmentation of ERP signals (180 channels) into 6 clusters. The first cluster corresponds to the period from 0(onset) to 116 ms. In this cluster, we can observe a network located in the occipital lobe. A shorter graph is then observed between117 and 153ms with also strong connections at the occipital and temporal lobe. The 3 rd cluster is located at 154 -190 ms where strong connections appear in the temporal and parietal lobes. A network appears in the cluster 191-316ms with the presence of connectivity in frontal regions, followed by a cluster (317-480 ms) with dense connectivity in frontal and occipital regions. The networks (481-620 ms) then become denser in temporal left and frontal regions.

Table 1 shows the results obtained after the 'fitting' process which consists in testing the inter-subject variability. The map presence values are presented in column, the templates (T3, T6) have high values of maps presence between subjects (73%, 78%) and low values in the template (T1, T4). Preliminary results were presented to demonstrate that the proposed algorithm can be used to track the dynamics of the brain functional connectivity. Very interestingly, results are qualitatively consistent with the state of the art of the analyzed task regarding the activation starting in the occipital lobes during picture recognition followed by activation in the temporo-frontal r e g i o n s [START_REF] Salmelin | Dynamics of brain activation during picture naming[END_REF].One o f t h e c hallenge faced by our study is to specify the maximal number of clusters. The increase of the maximal number of clusters increases the computation cost. In our study we used 10 clusters. A classical and still unsolved difficult question relates to the setting of threshold values applied on the connectivity matrices. In our study, the chosen threshold value was 10%, [START_REF] Sporns | Structure and function of complex brain networks[END_REF]. Other approaches can be explored like those based on surrogate data [START_REF] Rodriguez | Perception's shadow: long-distance synchronization of human brain activity[END_REF], although requiring an even higher computation time. Finally, this algorithm is a good start to then compare the dynamics of the brain connectivity for two different cognitive tasks such as the difference between spelling [START_REF] Hassan | Graph-based analysis of brain connectivity during spelling task[END_REF] and naming task which is our ongoing work. Also the algorithm will be applied on EEG source connectivity graphs to track brain functional connectivity dynamics at cortical source level.

V. CONCLUSION

In this paper we proposed a new algorithm to track the functional brain connectivity dynamics and segment hr-EEG signals recorded during picture recognition and naming task using K-means c l u s t e r i n g o f t h e c o n n e c t i v i ty graphs. This algorithm shows good performance to demonstrate the stability of the functional brain connectivity overs some short periods of time and to segment the cognitive process into functional connectivity microstates.

Fig1.

  Fig1.Spatial correlation (sC) of 3 templates (t=120, t=300, t=450) with all the others graphs (T).

Fig 2 .

 2 Fig 2.Event Related Potentials (180 electrodes) for the picture naming task and temporal distribution of the functional connectivity graphs revealed by spatio-temporal segmentation algorithm detailed in the paper. Connectivity graphs and their specific time windows are showed by each rectangle.

ACKNOWLEDGEMENTS

This work has received support from European Research Council (ERC) and from French government granted to the CominLabs excellence laboratory and managed by the National Research Agency in the "Investing for the Future" program under reference Nb. ANR-10-LABX-07-01.

Author:

Ahmad Mheich Institute:

LTSI, University of Rennes1 Street:

Compus de Beaulieu City:

Rennes Country:

France Mail:

ahmad.mheich@univ-rennes1.fr