Christophe Denis
email: christophe.denis@parisdescartes.fr

Classification in postural style based on stochastic process modeling

Keywords: Change point estimation, Multiclass classification, Cross-validation, Postural maintenance, Stochastic process modeling

We address the statistical challenge of classifying subjects as hemiplegic, vestibular or normal based on complex trajectories obtained through two experimental protocols which were designed to evaluate potential deficits in postural control. The classification procedure involves a dimension reduction step where the complex trajectories are summarized by finite-dimensional summary measures based on a stochastic process model for a real-valued trajectory. This allows us to retrieve from the trajectories information relative to their temporal dynamic. A leave-one-out evaluation yields a 79% performance of correct classification for a total of n = 70 subjects, with 22 hemiplegic (31%), 16 vestibular (23%) and 32 normal (46%) subjects.

Introduction

This article contributes to the study of postural maintenance. Posture is fundamental for physical activity. A deficit in postural maintenance often results in falling, which is particularly hazardous in elderly people. The main objective of the research in postural maintenance is to adapt protocols for functional rehabilitation for people who display deficits in maintaining posture. We focus here on the issue of classifying subjects in terms of postural maintenance. A cohort of 70 subjects has been followed at the center for the study of sensorimotor functioning (CESEM) of the University Paris Descartes. The subjects who did not exhibit deficit in postural maintenance were labelled as normal. The others were hemiplegic and vestibular subjects and labelled accordingly. A hemiplegic subject suffers from a disorder of his proprioceptive system, which pertains to the sense of position, location, orientation of the body and its parts. A vestibular subject suffers from a deficit of his vestibular system, which is the system composed by the inner ear and the vestibular nerve and contributes to the sense of balance.

Each subject completed two experimental protocols designed to evaluate his/her potential deficits in maintaining posture. Each protocol is divided into three phases: a first phase of 15 seconds with no postural perturbation, a second phase of 35 seconds with postural perturbation followed by a last phase of 20 seconds without postural perturbation. During each protocol, measurements of the center-of-pressure of each foot are performed at discrete times, which results in temporal trajectories. Our objective is to classify subjects as hemiplegic, vestibular, and normal based on these trajectories. This is a significantly more difficult extension of the problem considered by [START_REF] Chambaz | Classification in postural style[END_REF], where only hemiplegic and normal subjects were to be classified (particularly because classifying into three classes is more difficult than classifying into two classes), and thus another step in the direction of clustering subjects in terms of their postural style. We refer to the bibliography in [START_REF] Chambaz | Classification in postural style[END_REF] for a review on the analysis of postural control.

Our present study is related to the topic of functional data classification. There is a sizeable literature dedicated to this topic. A variety of methods have been proposed, relying for instance on linear discriminant analysis [START_REF] James | Functionnal linear discriminant analysis for irragularly sampled curves[END_REF], principal component analysis [START_REF] Hall | A functionnal data-analytic approach to signal discrimination[END_REF], or a functional data version of the nearest-neighbors classification rule [START_REF] Biau | Functional classification in hilbert spaces[END_REF]. We refer to [START_REF] Ramsay | Functional Data Analysis[END_REF] for a general introduction to functional data analysis. Our approach here involves a dimension reduction step of our complex trajectories based on change points estimation and inference on a stochastic process model. This allows us to better use the data at hand than in [START_REF] Chambaz | Classification in postural style[END_REF] in the sense that we manage to retrieve from the trajectories information relative to their temporal dynamic. In contrast, [START_REF] Chambaz | Classification in postural style[END_REF] retrieve static information from the trajectories in the sense that the dimension reduction step relies on comparisons of basic statistics (such as the mean value of a segment of the trajectory, see Section 5.2) computed on arbitrarily chosen time intervals starting or ending where the perturbation phase starts or ends.

The data set and its modeling are introduced in Section 2. We present our inference procedure in Section 3. We carry it out on real and simulated data, and summarize the results in the latter section. The classification procedure is presented in Section 4. The results of its application to the real dataset are exposed and commented on in Section 5, which concludes on some perspectives for future research.

Data and modeling

The dataset at hand is described in Section 2.1. We introduce and motivate its modeling in Section 2.2.

Original dataset

The dataset was collected at the center for the study of sensorimotor functioning (CESEM) of the University Paris Descartes. It is composed of a cohort of n = 70 subjects. Among the 70 subjects, 22 are hemiplegic (due to cerebrovascular accidents), 16 are vestibular, while the 32 remaining subjects are identified as normal based on an initial medical evaluation.

Each subject completes two protocols which are designed to evaluate potential deficits in maintaining posture. These protocols have been identified as the most informative among four similar protocols for classifying hemiplegic versus normal subjects in the earlier study [START_REF] Chambaz | Classification in postural style[END_REF]. Both protocols are divided into three phases: a first phase of 15 seconds with no postural perturbation, a second phase of 35 seconds with postural perturbations (either some muscular perturbations or a combination of muscular and visual perturbations), and a third phase of 20 seconds with no postural perturbation. We expect that the subject's postural sway change around the beginning and the end of the perturbation phase (around 15 and 50 seconds).

For each protocol, the center-of-maximal-pressure exerted by each foot on a forceplatform is recorded at equispaced discrete times. Thus each protocol results in a trajectory (L i , R i ) i≤N , where (L i , R i ) is the observation at time t i = iδ for i = 1, . . . , N = 2800 and δ = 0.025 seconds. For each i = 1, . . . , N , 1: Specifics of the two protocols considered in this study. A protocol is divided into three phases: a first phase with no postural perturbation is followed by a second phase with perturbations, which is itself followed by a third phase without perturbation. Different kind of perturbations are considered. In protocol 1, one perturbs the processing by the brain of proprioceptive information (through muscular stimulation). In protocol 2, one perturbs the processing by the brain of both visual information (subjects must close their eyes) and proprioceptive information (through muscular stimulation).

L i = (L 1 i , L 2 i ) ∈ R 2 and R i = (R 1 i , R 2 i ) ∈
Following [START_REF] Chambaz | Classification in postural style[END_REF], we derive from (L i , R i ) i≤N the one-dimensional trajectory C 1:N (sometimes we will denote a n-tuple (x 1 , . . . , x n ) by x 1:n ) characterized by

C i = L i + R i 2 -γ 2 ,
where γ is defined as the componentwise median value of ( 1 2 (L i + R i )) i≤400 over the ten first seconds of the protocol. The process C 1:N provides a relevant description of the sway of the body during the course of the protocol. In Figure 1 we display the trajectories C 1:N corresponding to protocol 1 (left plot) and protocol 2 (right plot) as completed by a single hemiplegic subject. Figure 1 confirms the intuition that the subject's postural sway is not necessarily instantaneously affected by the start and end of perturbations.

Data modeling

We model the trajectory C 1:N as the observation at discrete times of a stochastic process (C(t)) t∈ [δ,N δ] characterized by a stochastic differential equation.

We model the effects of the perturbations by possible changes in the volatility and drift functions at two unknown change points. In order to account for the fact that subjects react differently, we also assume that the change points differ among subjects. We denote by (T 1 , T 2 ) the change points of a trajectory, with

T 0 = τ 0 δ = δ < T 1 = τ 1 δ < T 2 = τ 2 δ < T 3 = 70 (τ 0 = 1,
τ 1 and τ 2 are integers). We assume that there exist two functions a, b, and two parameters (φ 1 , φ 2 , φ 3 ) and (σ 1 , σ 2 , σ 3 ) such that, for all t

∈ [T k-1 , T k [ and k = 1, 2, 3, dC(t) = a(C(t), φ k )dt + b(C(t), σ k )dW (t) C(T k-1 ) = C τ k-1 ,
where (W (t)) t∈[δ,N δ] is a standard Wiener process.

We now specify the parametric forms of the volatility and drift functions a and b. Because Figure 1 indicates that the variance of C 1:N is not constant over time and because, in addition, the process C 1:N takes positive values, we decide to rely on the classical Cox-Ingersoll-Ross (CIR) process [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equation[END_REF], by setting a(

x, φ = (λ, µ)) = λ(µ -x) and b(x, σ) = σ √
x (for all x ≥ 0). In summary, we model the trajectory C 1:N as the observation at discrete times of the stochastic process (

C(t)) t∈[δ,N δ] characterized, for t ∈ [T k-1 , T k [ (k = 1, 2, 3), by    dC(t) = λ k (µ k -C(t))dt + σ k C(t)dW (t) C(T k-1 ) = C τ k-1 , (1)
where the initial condition

C(T k-1 ) is the observation C τ k-1 at time T k-1 , λ k , µ k and σ k are positive constants. It is known that if 2µ k λ k /σ 2
k ≥ 1 then the process remains positive and admits a stationary distribution. This distribution is a Gamma distribution with shape parameter 2µ k λ k /σ 2 k and scale parameter σ 2 k /2λ k [START_REF] Kloeden | Numerical Solution of Stochastic Differential Equation[END_REF]. Therefore, the complete parameter writes as (τ

1 , τ 2 , θ 1 , θ 2 , θ 3 ) with θ k = (λ k , µ k , σ k ) for k = 1, 2, 3.
We assume that the changes necessarily affect the drift parameter through a change in µ, or put in other words that µ 1 = µ 2 and µ 2 = µ 3 . Changes can also affect the other parameters.

Inference

In this section we address the inference of the parameter (τ 1 , τ 2 , θ 1 , θ 2 , θ 3 ). We first describe the estimating procedure in Section 3.1, then we present its results on the real dataset in Section 3.2, and finally we summarize a simulation study performed to evaluate its performances in Section 3.3.

A two-step estimating procedure

We define a two-step estimating procedure: we first estimate the sequence of change points, then we estimate (θ 1 , θ 2 , θ 3 ) on each interval characterized by the previously estimated change points. A great variety of methods have been proposed for the estimation of change points (Basseville and [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF][START_REF] Bai | Least squares estimation of a shift in linear processes[END_REF][START_REF] Cuenod | Parameter estimation and change-point detection from Dynamic Contrast Enhanced MRI data using stochastic differential equations[END_REF][START_REF] Lavielle | Using penalized contrasts for the change-point problem[END_REF]Bardet et al., 2012, among many others). We choose to adopt a popular approach originally proposed by [START_REF] Bai | Least squares estimation of a shift in linear processes[END_REF]. The approach relies on a least-squares criterion and aims at detecting change points which affect the mean of a linear process. The estimator (τ 1 , τ2 ) of the sequence of change points (τ 1 , τ 2 ) is defined as:

(τ 1 , τ2 ) = arg min (τ 1 ,τ 2 ) 1 N 3 k=1 τ k -1 i=τ k-1 C i -Ck 2 ,
where Ck is the arithmetic mean of C τ k-1 :τ k -1 . This makes sense because if (C(t)) t∈[T k-1 ,T k [ reaches a stationary regime then the random variables C τ k , . . . , C τ k -1 are identically distributed from a stationary distribution whose mean parameter is µ k .

Set k ∈ {1, 2, 3}. On each interval [ Tk-1 , Tk [= [τ 1 δ, τ2 δ[, we estimate θ k by minimizing a contrast function which is based on the log-likelihood of the approximated discrete-times Euler-Marumaya's scheme [START_REF] Kessler | Estimation of an ergodic diffusion from discrete observations[END_REF], for instance). The latter scheme with step size δ guarantees that, for i = τk-1 , . . . , τk -1,

C i+1 ≈ (1 -δλ k )C i + δλ k µ k + σ k √ δ C i η i+1 , (2)
where (η i ) τk-1 <i≤τ k is a sequence of independent random variables with standard Gaussian distribution. By (2), it is convenient to consider the following equivalent parametrization:

θ k = (θ 1,k , θ 2,k , θ 3,k ), with θ 1,k = (1 -δλ k ), θ 2,k = δλ k µ k and θ 3,k = σ k √ δ.
The estimator θk of parameter θ k is defined as: θk = arg min

θ k ∈R 3 + τk -1 i=τ k-1 (C i+1 -θ 1,k C i -θ 2,k ) 2 C i θ 2 3,k + (τ k -τk-1 ) log(θ 2 3,k ),
which actually yields closed-form expressions:

θ1,k = (τ k -τk-1 ) C i+1 C i -C i+1 1 C i (τ k -τk-1 ) 2 -C i 1 C i , θ2,k = C i+1 -θ1,k C i (τ k -τk-1 ) , θ3,k = (C i+1 -θ1,k C i -θ2,k ) 2 ) C i τk -τk-1
(the sums in the above expressions range over [τ k-1 , τk -1]).

Under mild conditions, and if the process reaches the stationary regime, then (τ 1 , τ2 ) consistently estimates (τ 1 , τ 2 ) (see [START_REF] Lavielle | Detection of multiple changes in a sequence of dependent variables[END_REF]Lavielle and Ludeña, 2000, for instance). Furthermore, if the true change points (τ 1 , τ 2 ) are known then, under another set of mild conditions, the estimators ( θ1 , θ2 , θ3 ) consistently estimate (θ 1 , θ 2 , θ 3 ) (see [START_REF] Kessler | Estimation of an ergodic diffusion from discrete observations[END_REF], Theorem 1). To the best of our knowledge, there is no satisfactory result in the literature regarding the joint estimation of the change points (τ 1 , τ 2 ) and the parameter (θ 1 , θ 2 , θ 3 ).

Application to the real dataset

We undertake a simulation study of the properties of the two-step estimating procedure presented in the previous section, and summarize its results in the next section. Because we characterize our simulation scheme based on the results obtained when applying the latter procedure to the real dataset, we first present them here.

For each subject and protocol, we estimate (τ 1 , τ 2 , θ 1 , θ 2 , θ 3 ) from the corresponding (real) observed trajectory. The results pertaining to the estimation of (τ 1 , τ 2 ) are summarized in Table 2 (the mean and standard deviation of the estimates (τ 1 δ, τ2 δ) computed over each group of subjects are provided) and illustrated in Figure 2. Results pertaining to the estimation of (θ 1 , θ 2 , θ 3 ) are summarized in Table 5 (the mean and standard deviation of the estimates ( θ1 , θ2 , θ3 ) computed over each group of subjects are provided).

protocol 1 protocol 2 hemiplegic subjects 18.4 (7.4) 43.3 (12.2) 16.8 (1.0) 44.5 (12.4) vestibular subjects 20.9 (5.4) 49.9 (6.5) 19.5 (6.9) 50.3 (5.3) normal subjects 21.4 (8.2) 47.5 (10.7)

22.4 (8.1) 51.4 (2.7)

Table 2: For each subject and protocol, we estimate the change points (τ 1 , τ 2 ). For each group of subjects, we compute over the group the mean and standard deviation (given between parentheses) of the estimates (τ 1 δ, τ2 δ).

Three features of Table 2 are worth commenting on. First, one notes that there is no significant difference across groups of subjects (however, the means of τ1 δ and τ2 δ are slightly shifted to the left in the group of hemiplegic subjects relative to the two other groups). Second, the means of τ1 δ are close to the time of start of perturbations for all groups and protocols. As for the means of τ2 δ, they are close to the time of end for vestibular and normal subjects and both protocols. In hemiplegic subjects, one notes that the standard deviations are quite large and that the mean of τ2 δ is shifted to the left relative to the time of end of perturbations. This is due to the fact that for each protocol, 30% of hemiplegic subjects feature an estimator τ2 δ close to 30 seconds. Third, judging by the standard deviations, all hemiplegic subjects tend to react similarly by adjusting quickly to the perturbations undergone in protocol 2. Likewise, all normal subjects tend to react similarly by adjusting quickly to the end of perturbations undergone in protocol 2. In protocol 1, the large standard deviation associated to the mean value of τ2 δ computed over the group of normal subjects reflects the fact that 20% of these subjects feature an estimator τ2 δ close to 30 seconds.

Regarding Table 5, we consider in turn the parameters θ 1,k (k = 1, 2, 3), θ 2,k (k = 1, 2, 3) and θ 3,k (k = 1, 2, 3). First, the estimates θ1,k (k = 1, 2, 3) behave similarly across groups of subjects, protocols and intervals [ Tk-1 , Tk [ (k = 1, 2, 3). On the contrary, the estimates of θ2,k (k = 1, 2, 3) behave quite differently across groups of subjects, protocols and intervals Tk [ (k = 1,2,3). This is a promising feature for the sake of classifying subjects by group, which is our main problem at stake. As for the estimates θ3,k (k = 1, 2, 3), for any given group of subjects and protocol, they behave quite similarly. However, for protocol 2, it seems that the estimates θ3,k (k = 1, 2, 3) in normal subjects behave differently from their counterpart in hemiplegic or vestibular subjects. This is another promising feature.

In conclusion, and for the sake of illustrating our modeling and two-step estimating procedure, (i) we arbitrarily choose a hemiplegic subject and a normal subject;

(ii) for each of them, we retrieve the estimator (τ 1 , τ2 , θ1 , θ2 , θ3 ) of (τ 1 , τ 2 , θ 1 , θ 2 , θ 3 ) based on the trajectory obtained under protocol 2;

(iii) for each of them, we simulate a trajectory under (τ 1 , τ2 , θ1 , θ2 , θ3 ) (we refer to Section 3.3 for the specifics of the simulation);

(iv) we represent on the same plots the real and simulated trajectories, see Figure 3.

It appears that the real and simulated trajectories are relatively resembling.

Simulation study

We carry out a simulation study to evaluate the performances of the two-step estimating procedure presented in Section 3.1. We directly simulate the trajectory C 1:N . Specifically, (i) we set (T 1 , T 2 ) = (τ 1 δ, τ 2 δ) = (15, 50);

(ii) we rely on the Euler scheme (2) with step size δ/10 to approximate the sampling of (C(t)) t∈ [δ,N δ] from the distribution characterized by (1) where, for each k ∈ {1, 2, 3}, θ k equals the mean of its 22 estimates based on the 22 real trajectories associated to the 22 hemiplegic subjects and protocol 1 (see Section 3.2 and Table 5);

(iii) we conclude by sub-sampling (C(t)) t∈ [δ,N δ] to derive C 1:N .

We sample B = 100 independent copies of C 1:N . For each copy we estimate the parameters (τ 1 , τ 2 , θ 1 , θ 2 , θ 3 ). The means and standard deviations of the estimated change points and parameters computed over the 100 independent replications are reported in Table 3.

Three comments on Table 3 are in order. First, regarding the estimation of (τ 1 , τ 2 ), we note that the means of the estimated change points are very close to their respective true values. Moreover, the standard deviations are small. Second, regarding the estimation of (θ 1,k , θ 3,k ) (k = 1, 2, 3), we note that the means of the estimates of θ 1,k and θ 3,k are very close to their respective true values and that the standard deviations are small. Third, regarding the estimation of θ 2,k (k = 1, 2, 3), we emphasize that the means are quite apart Table 3: For each of the B = 100 independently simulated datasets, we derive the estimates (τ 1 δ, τ2 δ, θ1 , θ2 , θ3 ). We report here the true values and means and standard deviations (between parentheses) computed over the B = 100 replications.

k = 1 k = 2 k = 3 θ1,k 0.
from their respective true values. Moreover, the standard deviations are not small. Overall this indicates a poorer estimation of θ 2,k (k = 1, 2, 3) than of (θ 1,k , θ 3,k ) (k = 1, 2, 3). This is probably due to the fact that the time intervals [T k-1 , T k [ (k = 1, 2, 3) are relatively narrow given the value of δ.

Classification

In Section 4.2 we describe our classification procedure of subjects as hemiplegic, vestibular or normal based on their trajectories obtained under the two protocols. It is built upon the previous section. Indeed, it does not rely on the trajectories themselves but rather on their finite-dimensional summary measures whose definition, given in Section 4.1, depends on the results of our two-step estimation procedure.

Summary measures

Most classification procedures based on trajectories involve a preliminary step of dimension reduction where the high-dimensional trajectories are summarized into a low-dimensional summary measure [START_REF] Biau | Functional classification in hilbert spaces[END_REF][START_REF] Ramsay | Functional Data Analysis[END_REF]. Here we build a tailored finite-dimensional summary measure of every trajectory which relies on the estimates (τ 1 , τ2 , θ1 , θ2 , θ3 ) derived from it by applying our two-step estimating procedure. We argue in Section 3.2 that, overall, τ1 , τ2 , θ2,k and θ3,k (k = 1, 2, 3) may be relevant for the sake of classifying subjects as hemiplegic, vestibular or normal. Because the ratio θ2,k /(1 -θ1,k ) = μk is easier to interpret than θ2,k (and since θ1,k varies very little across subjects and protocols) and because θ3,k = σk √ δ, we choose to define our finite-dimensional summary measure as (τ 1 , τ2 , μ1 , μ2 , μ3 , σ1 , σ2 , σ3 ). Hereafter we denote X 1 the latter vector derived from the trajectory associated to protocol 1 and X 2 that derived from the trajectory associated to protocol 2.

Classification procedure

We actually construct two classifiers, φ1 and φ2 , to classify subjects as hemiplegic, vestibular or normal based on either X 1 or X 2 for φ1 , and on both (X 1 , X 2 ) for φ2 .

The generic observed data structure associated to a given subject is (X 1 , X 2 , Y ), where Y ∈ {1, 2, 3} indicates the subject's group (with convention Y = 1 for hemiplegic, Y = 2 for vestibular, and Y = 3 for normal). We denote P 0 the true distribution of (X 1 , X 2 , Y ).

Let S be the set of all classifiers based on X = (X 1 , X 2 ). The misclassification risk to S ∈ S is

R (P 0 ) (S) = E P 0 [1{S(X) = Y }]. Denote R * = min S∈S R (P 0 ) (S) its minimum.
It is achieved at the Bayes classifier S * ∈ S, characterized by S * (X) = arg max y∈{1,2,3} P 0 (Y = y|X). For j = 1, 2, we also introduce the Bayes classifier S * j ∈ S which relies only on X j : S * j (X) = arg max y∈{1,2,3} P 0 (Y = y|X j ).

Our objective is to construct φ1 and φ2 which respectively estimate the better classifier among S 1 * and S 2 * (i.e., arg min S 1 * ,S 2 * R (P 0 ) (S j * )) and S * . We choose to rely on the popular methodology of random forests [START_REF] Breiman | Random forests[END_REF], which proved powerful in a variety of applications. This is made easy thanks to the R package randomForest (we used the default tuning). The construction of the estimators Ŝ1 * and Ŝ2 * of S 1 * and S 2 * , and that of φ2 is straightforward. We derive φ1 from ( Ŝ1 * , Ŝ2 * ) by V -fold cross-validation, with V = 10.

Application

This section is devoted to the application of our classification procedure to the real dataset. In Section 5.1 we apply it exactly as it is described in Section 4. In Section 5.2, we consider a slightly enhanced version, whose performances are better. We conclude on some perspectives in Section 5.3.

Performance of the classification procedure

We evaluate the performances of the classification procedure by the leave-one-out rule. We acknowledge that it may result in overly optimistic error rates. Resorting to the leaveone-out rule is notably motivated by the relatively small sample size of our dataset. The results are reported in Table 4 (second row).

With its leave-one-out performance equal to 74%, the best classifier is φ1 , which involves one protocol only. For curiosity, we also evaluate the performances of φ1 and φ2 when systematically replacing (τ 1 δ, τ2 δ) with (15, 50) (the start and end times of the perturbation phase). We report the results in Table 4 (first row). Quite satisfactorily, we note that both φ1 and φ2 are not as good as the original versions: estimating the change points proves very relevant.

Performance of an extended classification procedure

It is tempting to extend our classification procedure by simply extending the definition of the summary measure which is at its core. Following [START_REF] Chambaz | Classification in postural style[END_REF], we merely substitute (X 1 , U 1 ) and (X 2 , U 2 ) for X 1 and X 2 with U j (j = 1, 2) derived from

C 1:N as U j = ( C+ 1 -C- 1 , C- 2 -C+ 1 , C+ 2 -C- 2 ) where C- 1 = δ 5 i∈[10/δ,15/δ[ C i , C+ 1 = δ 5 i∈]15/δ,20/δ] C i , C- 2 = δ 5 i∈[45/δ,50/δ[ C i , C+ 2 = δ 5 i∈]50/δ,55/δ] C i .
We refer to [START_REF] Chambaz | Classification in postural style[END_REF] for a justification. Then we apply the extended classification procedure and report its performances in Table 4 (third row).

With its leave-one-out performance equal to 79%, the best classifier is φ2 , which both extended summary measures (X 1 , U 1 ) and (X 2 , U 2 ). This is slightly better than the best performance obtained in Section 5.1: enriching the summary measures seems to provide relevant additional information for the sake of classifying subjects as hemiplegic, vestibular or normal. performances φ1 φ2

4.2 with (τ 1 δ, τ2 δ) = (15, 50) 61% 64% procedure of Section 4.2 74% 68% 5.2 76% 79%

Table 4: Leave-one-out performances of φ1 and φ2 on the real dataset for the sake of classifying subjects as hemiplegic, vestibular or normal. First row: classification procedure of Section 4.2 when imposing (τ 1 δ, τ2 δ) = (15, 50). Second row: classification procedure of Section 4.2. Third row: extended classification procedure of Section 5.2.

Perspectives

We address the statistical challenge of classifying subjects as hemiplegic, vestibular or normal based on complex trajectories obtained through two experimental protocols which were designed to evaluate potential deficits in postural control. The classification procedure involves a dimension reduction step where the complex trajectories are summarized by finite-dimensional summary measures based on a stochastic process model for a realvalued trajectory. This allows us to retrieve from the trajectories information relative to their temporal dynamic. A leave-one-out evaluation yields a 79% performance of correct classification for a total of n = 70 subjects, with 22 hemiplegic (31%), 16 vestibular (23%) and 32 normal (46%) subjects.

In future work, we will extend the classification procedure by introducing finite-dimensional summary measures based on a stochastic process model for the original trajectories in R 4 . We will also draw advantage from our good understanding of the classification problem to tackle the closely related next statistical challenge of clustering our subjects in terms of postural style. ) (k = 1, 2, 3). For each group of subjects and each protocol, we compute over the group and for that protocol the mean and standard deviation (between parenthese) of the estimates ( θ1

, θ2 , θ3
).

Figure 1 :

 1 Figure 1: Two trajectories C 1:N respectively corresponding to protocol 1 (left) and protocol 2 (right) undergone by a single hemiplegic subject.

Figure 2 :

 2 Figure 2: Representing the estimated change points (τ 1 δ, τ2 δ) (vertical lines) as obtained based on the two (real) observed trajectories of the same single hemiplegic subject as in Figure 1 (left: protocol 1, right: protocol 2).

Figure 3 :

 3 Figure 3: In each plot (left: hemiplegic subject, right: normal subject), we represent the real (solid line) and simulated (dotted line) trajectories, where the parameters of the simulation are derived from the real trajectory by applying our two-step estimating procedure.

Table

				R 2
	respectively correspond to the left and right feet.
	protocol 1st phase (0→15s) 2nd phase (15→50s) 3rd phase (50→70s)
	1		muscular stimulation
	2	no perturbation	eyes closed	no perturbation
			muscular stimulation

Table 5 :

 5 For each subject and each protocol, we first estimate the change points (τ 1 δ, τ2 δ) then we compute the estimates ( θ1

	)
	, θ3
	, θ2

3 ) over each of the three resulting intervals [δ, τ1 δ[, [τ 1 δ, τ2 δ[ and [τ 2 δ, N δ], where θ

Acknowledgments

I wish to thank I. Bonan (Service de Médecine Physique et de réadaptation, CHU Rennes) and P-P. Vidal (CESEM, Université Paris Descartes), my PhD advisors A. Chambaz and A. Samson (MAP5, Université Paris Descartes), and my colleague V. Genon-Catalot (MAP5, Université Paris Descartes), for introducing me to this interesting problem and providing the dataset, for their attentive guidance and all the fruitful discussions.


Isabelle Corouge

Aurore Esquevin

Florence Lejeune

Jean-Christophe Ferré

Elise Bannier

Catherine Merck

Serge Belliard

Christian Barillot

Jean-Yves Gauvrit

Arterial Spin Labeling at 3T in semantic dementia: perfusion abnormalities detection and comparison with FDG-PET

Keywords: ASL, FDG-PET, MRI, semantic dementia

Arterial Spin Labeling (ASL) is a non invasive perfusion imaging technique which has shown great diagnosis potential in dementia. However, it has never been applied to semantic dementia (SD), a rare subtype of frontotemporal lobar degeneration characterized by the gradual loss of conceptual knowledge, which is actually explored by a now well established marker of SD: 18 fluorodeoxyglucose-positron emission tomography (FDG-PET) imaging. Although ASL and FDG-PET respectively measure perfusion and metabolism, they have been shown to be strongly correlated. We propose to explore the ability of ASL to detect perfusion abnormalities in SD in comparison with FDG-PET. Using patients and healthy subjects data from an ongoing clinical study, we apply our analysis framework starting with visual comparison of ASL and FDG-PET, and focusing on ASL data preprocessing and statistical analysis at the individual and group level. Preliminary results yield concordant observations between ASL and FDG-PET as well as expected hypoperfusions in SD, namely in the left temporal lobe, thus suggesting the potential of ASL to assess perfusion impairments in SD.

Introduction

Arterial spin labeling (ASL) measures brain perfusion by using blood water protons as an endogenous tracer. An inversion pulse labels the inflowing blood and after a delay called inversion time, a labeled image of the volume of interest is acquired. The subtraction of the labeled image from a control image, i.e., non labeled, reflects the quantity of spins that have perfused the imaged volume, producing a perfusion-weighted (PW) image. The signal to noise ratio of the PW map is increased by acquiring repeated pairs of label/control images. Cerebral blood flow (CBF in ml/100g/min units) is estimated by application of a quantification model.

As it is a non invasive, reproducible and quantitative technique, ASL has been used in numerous clinical applications [START_REF] Detre | Applications of Arterial Spin Labeled MRI in the brain[END_REF]. Its diagnosis potential has namely been shown in dementia, for instance to differentiate Alzheimer's disease (AD) from frontotemporal dementia or to predict cognitive decline and conversion from Mild Cognitive Impairment to AD [START_REF] Hu | Distinct cerebral perfusion patterns in FTLD and AD[END_REF][START_REF] Chao | ASL perfusion MRI predicts cognitive decline and conversion from MCI to dementia[END_REF]. Previous studies have also compared ASL with 18 fluoro-deoxyglucose positron emission tomography (FDG-PET), which remains the reference technique. Indeed, FDG-PET imaging is particularly sensitive in detecting disease-related functional changes [START_REF] Jagust | Positron emission tomography and magnetic resonance imaging in the diagnosis and prediction of dementia[END_REF] and it has been well established as a marker of degenerative dementia. Though ASL and FDG-PET image different physiological processes, perfusion on the one side and glucidic metabolism on the other side, they are generally tightly coupled at rest [START_REF] Musiek | Direct comparison of fluorodeoxyglucose positron emission tomography and arterial spin labeling magnetic resonance imaging in alzheimer's disease. Alzheimer's dementia : the journal of the Alzheimer[END_REF] and they were shown to identify similar abnormality patterns in AD [START_REF] Chen | Voxel-level comparison of arterial spin-labeled perfusion MRI and FDG-PET in Alzheimer disease[END_REF].

To our knowledge, ASL has never been applied to semantic dementia (SD). SD is a rare subtype of frontotemporal lobar degeneration [START_REF] Hodges | Semantic dementia: a unique clinicopathological syndrome[END_REF] characterized by a gradual loss of conceptual knowledge and of word meaning. The SD syndrome arises out of the temporal lobe, usually showing atrophy often bilateral but mostly left-sided on structural MRI, as well as hypometabolism, also predominantly on the left side, in FDG-PET [START_REF] Desgranges | Anatomical and functional alterations in semantic dementia: a voxel-based MRI and PET study[END_REF].

In this paper we explore the ability of ASL to detect perfusion abnormalities in SD in comparison with FDG-PET, relying on data extracted from an ongoing study. Section 2 describes the population, the imaging protocols as well the methodological framework for data processing and analysis. Section 3 presents preliminary results before discussion in Sect. 4.

Materials and Methods

Participants

Fifteen subjects have been included in the study so far: 6 patients fulfilling the SD diagnostic criteria [START_REF] Neary | Frontotemporal lobar degeneration: a consensus on clinical diagnostic criteria[END_REF] (2 males, 4 females, mean age = 63.3 ± 2.4 yo), and 9 healthy control subjects (3 males, 6 females, mean age = 60.8 ± 3.9 yo), matched for sex and age (Mann-Whitney U test with p = 0.28). All the subjects underwent an MRI scan, only patients underwent a FDG-PET examination. In addition, cognitive performance was assessed for all the subjects. Note that these neuropsychological data are not included in the analysis presented in this paper.

Data Acquisition

Anatomical and perfusion MRI Data Patients and control subjects were scanned on a 3T whole body Siemens MR scanner (Magnetom Verio, Siemens Healthcare, Erlangen, Germany) with a 32-channel head coil. Anatomical data included a high resolution 3D T1-weighted MPRAGE sequence (3D T1) with the following imaging parameters: TR/TE/TI = 1900/2.27/900 ms, 256×256 mm 2 FOV and 176 sagittal slices, 1x1x1 mm 3 resolution, parallel imaging GRAPPA2.

Perfusion data were acquired using a pulsed ASL PICORE Q2TIPS sequence with total scan time of approximately 6 minutes. The imaging parameters were: TR/TE = 3000/18 ms, flip angle 90 o , matrix size 64×64, TI 1 /TI 1 s/TI 2 = 700/1500/1800 ms. The labeling slice was 10 cm thick and was placed 3 cm below the acquisition volume. Crusher was set to achieve a cutoff velocity of 4 cm/s. Fourteen axial slices were acquired sequentially from inferior to superior in the AC-PC plane, with 3×3 mm 2 in-plane resolution, 7 mm slice thickness and 0.7 mm gap. One control volume plus 60 repetitions, i.e., label/control pairs, (121 volumes), finally composed the ASL data series.

PET Data Only patients underwent the FDG-PET protocol imaging which was conducted on a discovery ST PET/CT Imaging system, GE, Healthcare, Waukesha, WI. PET/CT imaging was initiated 30 minutes after intravenous injection of 148 MBq of 18 FDG and was conducted at rest in a dark quiet room. The in-plane PET image resolution was 2.34×2.34 mm 2 , with 3.27 mm slice thickness. Non-contrast CT acquisitions were performed (2 mm slice thickness) for anatomical co-registration, attenuation and scatter correction of the PET images.

Data Pre-Processing

Pre-processing of the image data was performed with an inhouse pipeline using MATLAB (v. 7.12.0, The MathWorks Inc., 2011) and the SPM8 toolbox (Wellcome Department of Imaging Neuroscience at University College London, UK) as follows:

Anatomical data The anatomical 3D T1 was first denoised using the NLmeans algorithm by [START_REF] Coupé | An optimized blockwise nonlocal means denoising filter for 3-d magnetic resonance images[END_REF] and corrected for intensity inhomogeneity using SPM. The images were then segmented into grey matter (GM), white matter (WM) and cerebro-spinal fluid (CSF) probability maps using the MNI template tissue probability as an a priori for brain tissue classification. Estimated by this same unified segmentation model SPM routine, spatial normalisation parameters were applied to register the 3D T1 volume to the MNI space.

ASL data The ASL data series was motion corrected by a rigid body transform minimising the sum of squared differences cost function. A two-pass procedure first realigned all the control and label volumes onto the first volume of the series, then registered the series to the mean of the images aligned in the first pass.

The motion-corrected ASL series was co-registered to the 3D T1 using a rigid transform. The latter was estimated by maximizing normalised mutual information between the mean control image, i.e., the average of all the realigned control volumes, and the 3D T1 GM map.

Letting aside the first volume of the series, pM 0 , the co-registered ASL images were pairwise subtracted (control -label images) to produce a series of perfusion weighted maps, which was subsequently averaged to produce a perfusion weighted (PW) map.

The PW map was quantified to a CBF map by applying the standard kinetic model [START_REF] Buxton | A general kinetic model for quantitative perfusion imaging with arterial spin labeling[END_REF][START_REF] Pfeuffer | Comparison of pulsed arterial spin labeling encoding schemes and absolute perfusion quantification[END_REF] :

f = λ∆M 2αT I 1 exp (-T I2+idx sl * T I sl )/T 1b M 0 ( 1
)
where f is the CBF map, ∆M is the perfusion weighted map, λ = 0.9 ml.g -1 is the blood/tissue water partition coefficient, M 0 represents the equilibrium magnetization of arterial blood and is approximated by the first volume of the series pM 0 , α = 0.95 measures the labeling efficiency, T I 2 = 1800 ms is the inversion time of the ASL sequence, idx sl is the slice index, starting from 0 for the first acquired slice, TI sl = 45 ms is the duration of acquisition of one slice, TI 1 = 700 ms is the temporal width of the bolus, T 1b = 1500 ms is the T 1 of blood [START_REF] Wang | Regional reproducibility of pulsed arterial spin labeling perfusion imaging at 3T[END_REF].

Due to its inherent low resolution, the intensity in a single voxel results from a mixture of the GM, WM, CSF measured signals and ASL suffers from severe partial volume effect. These partial volumes effects (PVE) were corrected as described in [START_REF] Petr | Denoising arterial spin labeling MRI using tissue partial volume[END_REF]. This approach considers the perfusion signal in a voxel as a weighted sum of the GM and WM perfusion signals where the weights are given by the 3D T1 brain tissues probability maps. A total variation optimization algorithm estimates the respective contributions of GM/WM perfusion signals yielding high resolution perfusion maps, henceforth called PVE CBF maps.

In order to perform inter-subject voxelwise comparison, the PVE CBF maps, already aligned to the 3D T1 image, were registered to MNI space using the spatial normalisation parameters estimated on the anatomical data.

As CBF varies greatly across individuals, the PVE CBF maps were intensity normalised to account for a global interindividual variability. The normalisation factor was set to the mean of the PVE CBF values over the voxels having a grey matter probabilty superior to 0.7.

PET Data

The attenuation corrected PET images were coregistered to the 3D T1 image using the normalised mutual information cost function. Spatial normalisation to the MNI space was then performed using the parameters estimated for the 3D T1.

Data Analysis

Data analysis is performed in 3 steps: first, the data were visually evaluated, with a focus on the concordance between perfusion and metabolism as observed in the images. Then ASL data were quantitatively analysed at an individual level to identify perfusion defects in each patient based on a template approach. Last, group analysis was done voxelwise and in regions of interest.

Note that only voxels where the brain tissues probability exceeded 0.5 were tested in the subsequent statistical analyses. For voxelwise comparisons, the data were smoothed with a 8×8×8 mm 3 Gaussian kernel. Qualitative analysis For each patient, a radiologist visually inspected the PVE CBF and PET maps seeking to detect any metabolic and perfusion abnormalities. The reading focused on the superior frontal, basifrontal, anterior and posterior temporal, parietal and occipital areas for the right and left hemispheres.

Voxelwise Individual Statistical Analysis of ASL Data Hypo-and hyperperfusions were sought at an individual level by comparing each patient ASL map to an ASL template. Assuming a Gaussian distribution in each voxel, such a template is defined as the mean and associated standard deviation map estimated from a set of input maps. Using the set of smoothed PVE corrected, spatially and intensity normalised CBF maps from the control group as inputs, the ASL template was built as described in [START_REF] Petr | Construction and evaluation of a quantitative arterial spin labeling brain perfusion template at 3T[END_REF] and is considered as a reference for brain normal perfusion for an age-and sex-matched population with respect to the patient group.

Voxelwise comparisons between the individual maps and the ASL template were run based on one-tailed t-tests1 as implemented in the SPM General Linear Model (SPM). Differences were considered significant at the cluster level at a threshold of p < 0.05 (Family-wise-error (FWE) corrected for multiple comparisons at the cluster level with a cluster defining threshold of p < 0.001). Voxelwise Group Statistical analysis The smoothed PVE corrected, spatially and intensity normalised CBF maps of both control and patient groups were fed to GLM one-tailed two-sample t-tests analysis. Between group comparisons were considered significant at the cluster level at a threshold of p < 0.05, FWE corrected, with a cluster defining threshold of p < 0.001.

ROIs Group Statistical analysis Regions of interest (ROIs) were extracted from the PVE corrected, spatially and intensity normalised CBF maps using the Automated Anatomical Labeling (AAL) template [START_REF] Tzourio-Mazoyer | Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain[END_REF]. For each hemisphere, the segmented ROIs were: the brain lobes (frontal, temporal, parietal, occipital), the cingulate gyrus, the insula and the deep grey nuclei (putamen, caudate, thalamus). Contrary to the AAL lobe definition, the hippocampus, the parahippocampus and the temporal pole were not included in the limbic but in the temporal lobe as described in [START_REF] Harnsberger | Diagnostic and Surgical Imaging Anatomy: Brain[END_REF].

For each subject, the average CBF value was computed over each ROI, left and right, excluding voxels with GM probability less than 0.7. A two-tailed twosamples t-test tested the difference between the control group and the patients group at a significant level of α = 0.05.

Results

Qualitative Analysis The visual inspection of the PVE CBF maps and FDG-PET data showed good agreement between the two modalities as illustrated in Fig. 1. For all the 6 patients, hypoperfusion and hypometabolism were observed bilaterally in the basifrontal as well as in the anterior temporal lobe though the left side was more severely affected. The posterior part of the temporal lobe was particularly impaired on the left side. For all the subjects but one, the left parietal lobe exhibited hypo-ASL and FDG-PET signals while its right counterpart presented normal imaging.

Voxelwise Individual Statistical Analysis of ASL Data Significant perfusion deviations from the ASL template were found for 3 out of the 6 patients, as illustrated in Fig. 2. Hypoperfusions could be observed in the left fusiform and lingual gyri for patient 1, in the right fusiform gyrus and in the bilateral basifrontal area for patient 2 and in the left anterior temporal area for patient 3. Hyperperfusions could also be observed in the lenticular nucleus, on the right side for patient 1 and on both sides for patient 2. As for the remaining 3 patients, no significant cluster was found.

Voxelwise Group Statistical analysis The voxelwise group statistical analysis detected hypoperfusions in the left temporo-parietal and temporo-occipital regions, and in the left and right temporal and insular regions. However, hypoperfusions globally predominate in the left anterior temporal region as shown in Fig. 3. No hyperperfusion was detected at the group level.

ROIs Group Statistical analysis

The ROIs statistical analysis found significant differences between the patients and control group in the left temporal lobe (p = 0.0004), in the left parietal lobe (p = 0.02) and in the insula also on the left side (p = 0.002).

Discussion

We have applied our framework for ASL data processing to SD, a rare neurodenerative disease. Our preliminary results are encouraging as they globally comply with the perfusion abnormalities expected in SD. First, the qualitative comparison of the PVE CBF maps and the FDG-PET data showed a good agreement between the two modalities with concordant regional hypoperfusion and hypometabolism. Second, the group analysis identified perfusion abnormalities in the frontal, insular and temporal regions, particularly pronounced on the left side as already reported in the literature. More surprising were the detected hypoperfusions in the left temporo-parietal and temporo-occipital regions. They might be related to the stage of the disease, and/or to an impairment of the fasciculi involved in SD [START_REF] Acosta-Cabronero | Atrophy, hypometabolism and white matter abnormalities in semantic dementia tell a coherent story[END_REF] (arcuate and fronto-parietal longitudinal superior tracts). An analysis of the atrophy might also shed light on this result. Similary, hyperperfusions of the lenticular globus detected at the individual level should be further investigated. Further analysis on a larger database might reveal more insight on these preliminary findings. Additional subjects are planned to be included in the next few months to increase the population of each group to at least 15. This paper focused on the quantitative analysis of ASL data. Future work will include a quantitative analysis of PET data in which hypometabolism should be automatically detected, followed by a quantitative comparison between ASL and FDG-PET. Atrophy and cortical thickness will also be analysed from the 3D T1 structural data. Eventually, correlation of the imaging data with the neuropsychological scores will be studied. Ultimately, such a study might lead to demonstrate the ability of ASL to highlight cerebral abnormalities that give rise to cognitive deficit patterns in SD. This would allow improved clinical care and functional follow-up of SD patients as ASL is a non invasive and easily repeatable technique.

Fig. 1 .

 1 Fig. 1. Visual comparison of ASL and FDG-PET illustrated for 2 patients. Left, PET data with metabolism decreasing from red to green; right, PVE CBF maps. The white arrows indicate the hypometabolic and hypoperfused areas.

Fig. 2 .

 2 Fig. 2. Voxelwise individual statistic analysis: for 3 patients, detected hypoperfusions (blue) and hyperperfusions (red) superimposed on their anatomical map.

Fig. 3 .

 3 Fig. 3. Voxelwise group statistic analysis: superposition of the detected hypoperfusions (blue) to the average anatomy.

v-μ σ√ 1/N +1 ∼ T N -1where v is the voxel value to be tested, and μ, σ and N are respectively the sample mean, variance and size.


Alain Berlinet

Rémi Servien

EMPIRICAL ESTIMATOR OF THE REGULARITY INDEX OF A PROBABILITY MEASURE

Keywords: Regularity index, Lebesgue point, Small ball probability. Classification: 62G05

The index of regularity of a measure was introduced by Beirlant, Berlinet and Biau [1] to solve practical problems in nearest neighbour density estimation such as removing bias or selecting the number of neighbours. These authors proved the weak consistency of an estimator based on the nearest neighbour density estimator. In this paper, we study an empirical version of the regularity index and give sufficient conditions for its weak and strong convergence without assuming absolute continuity or other global properties of the underlying measure.

INTRODUCTION

The subject of this paper is related to the general problem of the estimation of small ball probabilities. Beirlant et al. [START_REF] Beirlant | Higher order estimation at Lebesgue points[END_REF] introduced the notion of regularity index of a measure to specify the rate at which the ratio of ball measures converges at a Lebesgue point. Indeed, this index is the exponent appearing in the second order term of the expansion of the small ball probability. Then, they defined an estimator of this index based on the nearest neighbour density estimator and proved its weak consistency. This estimator was applied to solve practical problems in nearest neighbour density estimation such as removing bias or selecting the number of neighbours. More recently Berlinet and Servien [START_REF] Berlinet | Necessary and sufficient condition for the existence of a limit distribution of the nearest neighbour density estimator[END_REF] proved that this regularity index was the key parameter governing the limit distribution of nearest neighbour density estimators so that its estimation may be crucial in the derivation of confidence intervals. In the present paper, we study an empirical version of the regularity index and give sufficient conditions for its weak and strong convergence. Unlike Beirlant et al. [START_REF] Beirlant | Higher order estimation at Lebesgue points[END_REF] we do not assume absolute continuity of the underlying measure but only a pointwise property of small ball probabilities. Notation and former results are given in the next section. Section 3 gathers the definition of the estimator and its convergence properties. Section 4 is devoted to the proofs of the theorems and examples are given in Section 5.

NOTATION AND FORMER RESULTS

Let µ be a probability distribution and λ be the Lebesgue measure on R d equipped with the Euclidean norm ||.||. We denote by B δ (x) the open ball with center x and radius δ. To evaluate the local behaviour of µ(B δ (x)) in relation to λ(B δ (x)) one can consider the ratio of these two quantities. If, for fixed x, the following limit

(x) = lim δ→0 µ(B δ (x)) λ(B δ (x)) ( 1
)
exists and is finite, then x is called a Lebesgue point of the measure µ (see Dudley [START_REF] Dudley | Real Analysis and Probability[END_REF] and Rudin [START_REF] Rudin | Real and Complex Analysis[END_REF]). This notion of Lebesgue point is essential to state elegant results with few restrictions on the functions to be estimated. In Berlinet and Levallois [START_REF] Berlinet | Higher order analysis at Lebesgue points[END_REF], examples where the density has a bad local behaviour at Lebesgue points are examined. To evaluate rates of convergence or investigate asymptotic normality of estimators, not only the convergence of the ratio of ball measures is required but also information on its higher order behaviour. In this context, Berlinet and Levallois [START_REF] Berlinet | Higher order analysis at Lebesgue points[END_REF] define a ρ-regularity point of the measure µ as any Lebesgue point x of µ satisfying

µ(B δ (x)) λ(B δ (x)) -(x) ≤ ρ(δ), (2)
where ρ is a measurable function such that lim δ↓0 ρ(δ) = 0. To specify an exact rate of convergence of the ratio of ball measures, Beirlant et al. [START_REF] Beirlant | Higher order estimation at Lebesgue points[END_REF] assumed that a more precise relation than (2) holds at the Lebesgue point x; namely

µ(B δ (x)) λ(B δ (x)) = (x) + C x δ αx + o(δ αx ) as δ ↓ 0, (3)
where C x is a non-zero constant and α x is a positive real number. It is easy to show that Equation (3) implies ρ-regularity at the point x with ρ(δ) = D x δ αx and D x > C x . The constants C x and α x are unique (provided they exist). Examples are provided in Section 5 with an absolute continuous measure and a measure with discrete part. The index α x is a regularity index that controls the degree of smoothness of the symmetric derivative of µ with respect to λ. The larger the value of α x , the smoother the derivative of µ is at the point x. Beirlant et al. [START_REF] Beirlant | Higher order estimation at Lebesgue points[END_REF] showed the interest of estimating the regularity index to solve practical problems in nearest neighbour density estimation, such as removing bias or selecting the number of neighbours.

More recently Berlinet and Servien [START_REF] Berlinet | Necessary and sufficient condition for the existence of a limit distribution of the nearest neighbour density estimator[END_REF] analyzed the effect of the value of α x on limit distributions of nearest neighbour density estimators. They gave a necessary and sufficient condition involving α x and the number of neighbours to have a limit distribution for the estimator.

The link with the small ball probability is clear since Equation ( 3) is equivalent to the expansion

P ( X -x ≤ δ) = V d δ d ( (x) + C x δ αx + o(δ αx ))
where X has probability distribution µ and V d = π d/2 /Γ(1+d/2) denotes the volume of the unit ball in R d . In other words, the second order term in the expansion of the small ball probability of radius δ at x is equal, up to a multiplicative constant, to δ d+αx . Hence it appears that to estimate α x one needs some information on the behaviour of µ(B δ (x)) as a function of δ. This is why the following theorem, proved by Beirlant et al. [START_REF] Beirlant | Higher order estimation at Lebesgue points[END_REF] will be useful in the sequel.

Theorem 2.1. Suppose that x ∈ R d is a Lebesgue point of µ with regularity index α x . Then, for any τ > 1,

lim δ→0 ϕ τ 2 δ (x) -ϕ τ δ (x) ϕ τ δ (x) -ϕ δ (x) = τ αx
where we denote, for δ > 0,

ϕ δ (x) = µ(B δ (x)) λ(B δ (x)) .
Now let X 1 , . . . , X n denote n independent random variables with distribution µ on R d , µ being unknown. Using the k n -nearest neighbour density estimator

f kn (x) = k n nV d X (kn) (x) -x d
where X (kn) (x) is the k th n -nearest neighbour of x and V d is the volume of the unit ball in R d , Beirlant et al. [START_REF] Beirlant | Higher order estimation at Lebesgue points[END_REF] introduced an estimator ᾱn,x of the regularity index inspired by the above theorem by setting, for τ > 1,

ᾱn,x = d ln τ ln f τ 2 kn (x) -f τ kn (x) f τ kn (x) -f kn (x) , (4)
if [f τ 2 kn (x) -f τ kn (x)]/[f τ kn (x) -f kn (x)] > 1 and ᾱn,x = 0 otherwise, and proved the weak consistency of ᾱn,x .

In the paper by Beirlant et al. [START_REF] Beirlant | Higher order estimation at Lebesgue points[END_REF] most results are stated under the assumption of absolute continuity of the measure µ with respect to Lebesgue measure. This is required for instance to get a beta distribution for the random variable µ B [X(k n ) (x)-x] (x) .

Our goal in the present paper is to define an empirical estimator inspired by the same theorem. For this, we simply replace in the expression of ϕ δ (x) the unknown quantity µ(B δ (x)) by its empirical counterpart. We prove the weak and strong consistency of the resulting estimator under the sole assumption that Equation (3) holds true. The present paper stays at a theoretical level, giving conditions on the deterministic sequence (δ n ) to get consistency. This is a first step. Further work should lead to an automatic choice of this sequence from the observed data.

THE EMPIRICAL ESTIMATOR AND ITS CONVERGENCE

Let (X i ) i≥1 be a sequence of independent real d-dimensional random vectors with distribution µ. The empirical measure µ n associated with X 1 , . . . , X n is defined by

µ n (A) = 1 n n i=1 I (Xi∈A) , A ⊆ R d ,
where

I (Xi∈A) = 1 if X i ∈ A 0 otherwise
and the associated empirical estimator of ϕ δ (x) by

ϕ n,δ (x) = µ n (B δ (x)) λ(B δ (x)) .
The following theorems state the weak and strong consistency of the empirical estimator defined by

α n,x = 1 ln τ ln ϕ n,τ 2 δn (x) -ϕ n,τ δn (x) ϕ n,τ δn (x) -ϕ n,δn (x) (5)
if ϕ n,τ 2 δn (x) -ϕ n,τ δn (x) / [ϕ n,τ δn (x) -ϕ n,δn (x)] > 1 and α n,x = 0 otherwise, (δ n ) being a sequence of positive numbers which will be assumed to tend to zero. As is usually the case almost sure consistency is obtained under stronger conditions on the sequence (δ n ).

PROOFS

The weak (respectively strong) consistency of α n,x is equivalent to the weak (resp. strong) consistency, for any τ > 0, of the ratio

R n (δ n ) = ϕ n,τ 2 δn (x) -ϕ n,τ δn (x) ϕ n,τ δn (x) -ϕ n,δn (x)
to τ αx . Let us fix τ > 0 and set

S n (δ n ) = ϕ n,τ δn (x) -ϕ n,δn (x) ϕ τ δn (x) -ϕ δn (x) .
We can write

ϕ n,τ 2 δn (x) -ϕ n,τ δn (x) ϕ n,τ δn (x) -ϕ n,δn (x) = ϕ τ 2 δn (x) -ϕ τ δn (x) ϕ τ δn (x) -ϕ δn (x) ϕ n,τ 2 δn (x) -ϕ n,τ δn (x) ϕ τ 2 δn (x) -ϕ τ δn (x) × ϕ n,τ δn (x) -ϕ n,δn (x) ϕ τ δn (x) -ϕ δn (x) -1
or equivalently

R n (δ n ) = ϕ τ 2 δn (x) -ϕ τ δn (x) ϕ τ δn (x) -ϕ δn (x) S n (τ δ n ) S n (δ n ) .
Let us first look at the variance of S n (δ n ). For this let us write

S n (δ n ) = 1 + A n (τ δ n ) -A n (δ n ) ∆ n
where

A n (δ n ) = ϕ n,δn (x) -ϕ δn (x) and ∆ n = ϕ τ δn (x) -ϕ δn (x).
The following lemma gives the asymptotic variance of A n (δ n ), the asymptotic covariance of (A n (τ δ n ), A n (δ n )) and the asymptotic variance of S n (δ n ).

n δ d n E (A n (δ n )) 2 = (x) V d , lim n-→∞ n δ d n E [A n (τ δ n )A n (δ n )] = (x) τ d V d and lim n-→∞ n δ d+2αx n E (S n (δ n ) -1) 2 = (x) τ d -1 τ d V d C 2 x (τ αx -1) 2 . Proof of Lemma 4.1 First note that Equation (3) implies that ∆ n = ϕ τ δn (x) -ϕ δn (x) = C x δ αx n (τ αx -1) + o(δ αx n ) and lim n-→∞ δ 2αx n ∆ 2 n = 1 C 2
x (τ αx -1) 2 . Now, using the fact that the random variable nµ n (B δn (x)) has the binomial distribution B(n, µ(B δn (x))) we get

n δ d n E (A n (δ n )) 2 = δ d n µ(B δn (x))(1 -µ(B δn (x))) [λ(B δn (x))] 2 = [1 -µ(B δn (x))] µ(B δn (x)) λ(B δn (x)) 1 V d
which gives the asymptotic variance of A n (δ n ).

As B δn (x) ⊂ B τ δn (x) the covariance

E ([µ n (B τ δn (x)) -µ(B τ δn (x))] [µ n (B δn (x)) -µ(B δn (x))]) is equal to 1 n (1 -µ(B τ δn (x))) µ(B δn (x))
and therefore

n δ d n E [A n (τ δ n )A n (δ n )] = δ d n (1 -µ(B τ δn (x))) µ(B δn (x)) λ(B τ δn (x)) λ(B δn (x)) = [1 -µ(B τ δn (x))] µ(B δn (x)) λ(B δn (x)) 1 τ d V d which gives the asymptotic covariance of (A n (τ δ n ), A n (δ n )). Changing δ n into τ δ n as argument of A n (.) gives lim n-→∞ n (τ δ n ) d E (A n (τ δ n )) 2 = (x) V d .
Gathering the above results one gets

lim n-→∞ n δ d+2αx n E (S n (δ n ) -1) 2 = lim n-→∞ n δ d+2αx n E A n (τ δ n ) -A n (δ n ) ∆ n 2 = (x) τ d -1 τ d V d C 2 x (τ αx -1) 2 .
This ends the proof of the lemma.

Remark. Note that under the assumptions lim n→∞ δ n = 0 and (x) > 0 the condition lim

n→∞ nδ d+2αx n = ∞
is not only sufficient but also necessary for the L 2 convergence of (S n (δ n )) .

Proof of Theorem 3.1

Under the conditions of Theorem 3.1, (S n (δ n )) and (S n (τ δ n )) converge to the constant 1 in the L 2 sense and therefore also in probability. Thus, their ratio tends to 1 in probability. By Theorem 2.1 (R n (δ n )) tends to τ αx in probability. This ends the proof.

Proof of Theorem 3.2.

As already said the conclusion of Theorem 3.2, is equivalent to the following property: For any τ > 1, lim n→∞ ϕ n,τ 2 δn (x) -ϕ n,τ δn (x) ϕ n,τ δn (x) -ϕ n,δn (x) = τ αx almost surely.

From Hoeffding's inequality (see [START_REF] Devroye | Combinatorial Methods in Density Estimation[END_REF]) for a binomial distribution we have

∀t > 0, P |µ n (B δn (x)) -µ(B δn (x))| ≥ t ≤ 2 exp -2nt 2 .
Taking ε > 0 and t = ελ(B δn (x)) |∆ n | , we get,

∀ε > 0, P A n (δ n ) ∆ n ≥ ε ≤ 2 exp -2n [ελ(B δn (x))∆ n ] 2 .
By Borel-Cantelli lemma, we have the convergence

A n (δ n ) ∆ n -→0 almost completely if ∀ε > 0, ∞ n=1 exp -2n [ελ(B δn (x))∆ n ] 2 < ∞. (6)
Now, set

γ n = V 2 d ∆ 2 n δ 2αx n .
As we have from Equation (3)

∆ 2 n = δ 2αx n (C 2 x (τ αx -1) 2 + o(1))
we have

γ n = V 2 d (C 2 x (τ αx -1) 2 + o(1)
) and the summand in Condition (6) writes exp -nδ 2(d+αx)

n γ n ε 2 = exp - nδ 2(d+αx) n ln n γ n ε 2 ln n = 1 n un with (γ n ) tending to γ = V 2 d C 2
x (τ αx -1) 2 > 0 as n tends to infinity and

u n = nδ 2(d+αx) n ln n γ n ε 2 .
The condition imposed on the sequence (nδ 2(d+αx) n / ln n) implies that for any ε > 0, the sequence (u n ) tends to infinity and therefore Condition (6) is satisfied. Thus (A n (δ n )/∆ n ) converges to 0 almost completely. In the same way one proves that (A n (τ δ n )/∆ n ) converges to 0 almost completely. It follows that

S n (δ n ) = 1 + A n (τ δ n ) -A n (δ n ) ∆ n
and (S n (τ δ n )) converge to 1 almost completely. Finally, using Theorem 2.1, we get the conclusion that α n,x converges to α x almost surely.

EXAMPLES

An example with an absolutely continuous measure

First consider the measure µ, absolutely continuous with respect to the Lebesgue measure on R, with density

f (x) = 1 - √ 2 3 + |x| 1 (-1/2,1/2) (x).
The distribution function F of µ is given by

F (x) =    0 if x ≤ -1/2 (1/2) + (1 - √ 2/3) x + (2/3) x |x| if -1/2 ≤ x ≤ 1/2 1 if x ≥ 1/2
For x ∈ (0, 1/2) and δ > 0, δ small enough, one has

F (x + δ) -F (x -δ) = 2δ(1 - √ 2/3) + (2/3) (x + δ) (3/2) -(x -δ) (3/2) .
For x ∈ (-1/2, 0) and δ > 0, δ small enough, one has

F (x + δ) -F (x -δ) = 2δ(1 - √ 2/3) + (2/3) -(-x -δ) (3/2) + (-x + δ) (3/2) = 2δ(1 - √ 2/3) + (2/3) (-x + δ) (3/2) -(-x -δ) (3/2)
Now, for 0 < |u| < 1,

(1 + u)

Note that lim

x-→0

C x = -∞.

At the point x = 0 one uses the fact that for δ ∈ (0, 1) one has

Theorem 3 . 1 .

 31 (Weak consistency) Suppose that x ∈ R d is a Lebesgue point of µ with regularity index α x . Then, under the conditions lim n→∞ α n,x converges to α x in probability.

Theorem 3 . 2 .

 32 (Strong consistency) Suppose that x ∈ R d is a Lebesgue point of µ with regularity index α x . Then, under the conditions lim n→∞ δ n = 0 and lim n→∞ nδ 2(d+αx) n ln n = ∞ the empirical estimator α n,x converges to α x almost surely.

Lemma 4 . 1 .

 41 Suppose that x ∈ R d is a Lebesgue point of µ with regularity index α x . Then, under the condition lim n→∞ δ

3 + o(δ 4 )

 34 if x > 0 and δ small enough,F (x + δ) -F (x -δ) 2 + o(δ 3 ) F (x + δ) -F (x -δ) 2) δ 2 + o(δ 3 ) if x < 0 and δ small enough, F (x + δ) -F (x -δ) ) 3 δ 2 + o(δ 3 ) F (x + δ) -F (x -δ) ) (3/2) δ 2 + o(δ 3 ).Thus, for x = 0,F (x + δ) -F (x -δ) 2δ -f (x) = -

F

ACKNOWLEDGEMENT

We wish to thank the referees for their careful reading of the paper and helpful comments.

An example with a measure having a discrete part

Now, consider the discrete probability measure ν supported by the sequence

.

For any δ ∈ (0, 1) , there exists a unique positive integer k(δ) such that

and we have

From the definition of k(δ) it follows that lim

and, from (7), we have

Finally, gathering the above results, one gets

The probability measure η = (µ + ν)/2 satisfies

So, the measure η, which is clearly not absolutely continuous, satisfies Equation


Alain F Berlinet
email: berlinet@math.univ-montp2.fr

Rémi Servien

Necessary and sufficient condition for the existence of a limit distribution of the nearest neighbour density estimator

Keywords: Asymptotic normality, nearest neighbour estimator, regularity index, small ball probability AMS Subject Classification: 62G05, 60E05

In density estimation, a number of convergence results can be stated with no restriction on the function to be estimated. Unlike these universal properties, the asymptotic normality of estimators often requires hypotheses on the derivatives of the underlying density together with additional conditions on the smoothing parameter. Yet, in spite of possible bad local behaviour of the density (it is not continuous or has infinite derivative), the convergence in law of the nearest neighbour estimator can still occur and provide confidence bands for the estimated density. Therefore a natural question arises: Is it possible to give a necessary and sufficient condition for the existence of a limit distribution of the nearest neighbour estimator? We answer this question by means of the regularity index recently introduced by Beirlant, Berlinet and Biau (1). As expected the limit distribution, when existing, is gaussian. Its mean and variance are given explicitly in terms of coefficients related to the regularity index. The second order exponent in the expansion of the small ball probability is shown to be the crucial parameter. In contrast to former results on sufficiency of conditions for asymptotic normality no hypothesis is needed on the continuity of the underlying density.

Introduction

Sufficient conditions for convergence and asymptotic normality of nearest neighbour density estimators have been given a long time ago. A classical result (Loftsgaarden and Quesenberry [START_REF] Loftsgaarden | A nonparametric estimate of a multivariate density function[END_REF]) states the pointwise convergence in probability of this estimator at points where the density is positive and continuous. It can be extended to any Lebesgue point (Bosq and Lecoutre (3)). Moore and Yackel [START_REF] Moore | Large sample properties of nearest neighbor density function estimates[END_REF] proved the asymptotic normality under additional conditions on the smoothing parameter and on the derivatives of the density. Their results were extended to multivariate nearest neighbour estimators with weight function by Mack (9). It seems that no author reconsidered this question until the paper by Berlinet and Levallois (2) who proved the asymptotic normality of the nearest neighbour density estimator in cases where the density has bad local behaviour. Indeed what is important is the local behaviour of the associated measure, more precisely the rate at which the derivative of the underlying measure is approximated by ratios of ball measures. We establish in Section 2 a necessary and sufficient condition on the number of neighbours for the existence of a limit distribution of the nearest neighbour estimator. As expected, the limit distribution when existing is gaussian. Its parameters are given explicitly in terms of coefficients appearing in the expansion of small ball probabilities. In Section 3, some simulations are presented in a case where the density has no derivative.

Definitions and results

We consider a sequence (X n ) n≥1 of d-dimensional independent random vectors with probability distribution µ and density f with respect to the Lebesgue measure λ on R d equipped with the Euclidean norm . . We denote by B(x, δ) (resp. B(x, δ)) the open (resp. closed) ball with center at x and radius δ.

The problem of estimating the probability density has received considerable attention in the litterature. A lot of methods have been developped like histograms, kernel estimators and their variants. We shall consider in this paper the nearest neighbour estimator f kn of f based on (X n ) n≥1 . It is defined by

f kn (x) = k n nλ(B n (x))
where B n (x) = B(x, R n (x)) is the smallest closed ball with center x containing at least k n sample points. The integer k n plays the role of a smoothing parameter: when k n is chosen too large the data are oversmoothed; they are undersmoothed in the opposite case. The estimate f kn (x) is the ratio of the frequency of sample points falling into B kn (x) to the Lebesgue measure of B kn (x).

Berlinet and Levallois studied the asymptotic normality of this estimator and pointed out the key role of the local behaviour of the probability measure associated with the density, more precisely the role of the rate at which the derivative of the measure is approximated by the ratios of ball measures. If, for fixed x, the following limit

ℓ(x) = lim δ→0 µ(B(x, δ)) λ(B(x, δ)) (1)
exists, then x is called a Lebesgue point of the measure µ. This notion of Lebesgue point is essential to state elegant results with few restrictions on the functions to be estimated. In Berlinet and Levallois (2), examples where the density has a bad local behaviour at Lebesgue points are examined. Yet, to evaluate rates of convergence or investigate asymptotic normality of estimators, we require not only the convergence of the ratio of ball measures but also information on higher order behaviour. In this context, Berlinet and Levallois define a ρ-regularity point of the measure µ as any Lebesgue point x satisfying

µ(B(x, δ)) λ(B(x, δ)) -ℓ(x) ≤ ρ(δ), (2)
where ρ is a measurable function such that lim δ↓0 ρ(δ) = 0. As it will be seen below, Berlinet and Levallois used this notion to extend previous asymptotic normality results of Moore and Yackel.

To specify an exact rate of convergence of the ratio of ball measures, Beirlant, Berlinet and Biau (1) assumed that a more precise relation than (2) holds at the Lebesgue point x, namely

µ(B(x, δ)) λ(B(x, δ)) = f (x) + C x δ αx + o(δ αx ) as δ → 0, (3)
where C x is a non-zero constant and α x is a positive real number. It is easy to show that relation (3) implies ρ-regularity at the point x with ρ(δ) = D x δ αx and D x > C x .

The constants C x and α x are uniquely determined (provided they exist). The index α x is a regularity index which controls the degree of smoothness of the symmetric derivative of µ with respect to λ. The larger the value of α x , the smoother the derivative of µ is at the point x. Beirlant, Berlinet and Biau used this regularity index to solve practical problems in nearest neighbour density estimation, such as removing bias or selecting the number of neighbours. Equation ( 3) is clearly equivalent to the small ball probability expansion

P ( X -x ≤ δ) = V d δ d (f (x) + C x δ αx + o(δ αx ))
where X has density f and

V d = π d/2 /Γ(1 + d/2
) denotes the volume of the unit ball in R d . In other words, the second order term in the expansion of the small ball probability of radius δ at x is equal, up to a multiplicative constant, to δ d+αx .

Sufficient conditions of convergence

In discriminatory analysis, Fix and Hodges (6) introduced the classification rule based on nearest neighbour (see also Devroye, Györfi and Lugosi ( 4)). As to the nearest neighbour estimator, it was first studied by Loftsgaarden and Quesenberry (8) at points where the density is positive and continuous. These authors proved that the conditions

lim n→∞ k n = ∞ and lim n→∞ k n n = 0 (4)
ensure the convergence in probability of f kn (x) to f (x). This result can be extended to any Lebesgue point [START_REF] Bosq | Théorie de l'Estimation Fonctionnelle[END_REF]. Alain F. Berlinet and Rémi Servien Moore and Yackel proved the following asymptotic normality result.

Theorem 2.1 : Under the convergence conditions (4), if

lim n→∞ k n n 2/3 = 0,
if f is continuous at x and has bounded derivatives in a neighbourhood of x and if f (x) > 0, then

T n (x) = k n f kn (x) -f (x) f (x)
tends in distribution to the standard normal N (0, 1) as n tends to infinity.

In practice, we have no information on the existence of derivatives. It is less restrictive to suppose ρ-regularity. Assuming ρ-regularity and a suitable rate of convergence for ρ(R n (x)), Berlinet and Levallois obtained the following theorem.

Theorem 2.2 : Suppose that the convergence conditions (4) hold, that x is a ρregularity point of the measure µ and that f (x) > 0. If, as n tends to infinity,

k n ρ(R n (x)) (5)
tends to 0 in probability then

T n (x) = k n f kn (x) -f (x) f (x)
tends in distribution to N (0, 1).

When the function f satisfies a Lipschitz condition of order α > 0 at x, i.e.

|f (x) -f (y)| ≤ C x x -y α (6)
for any y in some neighbourhood of x, it is possible to give a condition on the sequence (k n ) implying condition [START_REF] Devroye | The strong uniform consistency of nearest neighbor density estimates[END_REF]. This corollary of Theorem 2.2 was given in the univariate case by Berlinet and Levallois but a quick inspection of their proof easily yields to the d-dimensional version given in Corollary 2.3 below. This proof essentially relies on the fact that condition (6) implies ρ-regularity at the point x with

ρ(δ) = D x δ α .
Hence Theorem 2.2 gives the result.

Corollary 2.3: Suppose that the convergence conditions (4) hold, that f satisfies condition [START_REF] Fix | Discriminatory analysis, nonparametric discrimination: consistency properties[END_REF] and that f (x) > 0. If

lim n→∞ k 1+d/(2α) n n = 0 then T n (x) = k n f kn (x) -f (x) f (x)
tends in distribution to N (0, 1).

Main results

Now we consider the question of existence of a limit distribution for the nearest neighbour density estimator.

Theorem 2.4 : Suppose that the convergence conditions (4) hold, that x is a Lebesgue point in R d where ( 3) is satisfied with f (x) > 0. Then, as n tends to infinity,

T n (x) = k n f kn (x) -f (x) f (x)
converges in distribution if and only if the sequence

k 1+d/(2αx) n n n≥1
has a finite limit κ. When this last condition is satisfied, the asymptotic law of T n (x) is

N C x κ αx/d V αx/d d 1 f (x) 1+αx/d , 1 .
Theorem 2.4 provides a necessary and sufficient condition for the existence of a limit distribution and gives explicitly this distribution when it does exist.

Proofs

Proof of Theorem 2.4 Theorem 2.4 is a straightforward consequence of Lemma 2.7 and 2.8 below.

The proof of the lemmas relies on the following decomposition:

T n (x) = a n (b n + c n ) where a n = k n nµ(B n (x)) , b n = k n 1 f (x) µ(B n (x)) λ(B n (x)) -f (x) and c n = n √ k n k n n -µ(B n (x)) .
Moore and Yackel [START_REF] Moore | Large sample properties of nearest neighbor density function estimates[END_REF] proved the following result.

Lemma 2.5: For the sequences of random variables (a n ) and (c n ) defined above in the decomposition of T n (x), we have, under the convergence conditions (4),

a n P → 1 and c n L → N (0, 1)
as n tends to infinity.

Let us now turn to the asymptotic behaviour of the third quantity b n introduced in the decomposition of T n (x).

Lemma 2.6: Under the convergence conditions (4), and if (3) is satisfied, we have

lim n→∞ n k n αx/d 1 √ k n b n = C x V αx/d d 1 f (x) 1+αx/d
.

Proof :

From (3) we obtain that the function

Φ(δ) = µ(B δ (x)) λ(B δ (x)) -f (x) 1 δ αx is such that lim δ→0 + Φ(δ) = C x .
Using this function Φ in the expression of b n , we have

b n = Φ(R n ) nV d R d n k n αx/d 1 V αx/d d f (x) k αx/d+1/2 n n αx/d and n k n αx/d 1 √ k n b n = Φ(R n ) nV d R d n k n αx/d 1 V αx/d d f (x) .
Then, under the convergence conditions (4), we observe that

R n P → 0 and k n nV d R d n = f kn (x) P → f (x).
The convergence in probability of (R n ) to zero implies that

Φ(R n ) P → C x .
This, together with the convergence in probability of f kn (x) to f (x), gives the result.

Lemma 2.7: Sufficient condition. Under the conditions of Theorem 2.4, if

lim n→∞ k 1+d/2αx n n = κ then T n (x) D → N C x κ αx/d V αx/d d 1 f (x) 1+αx/d , 1 . Proof : Conditions of Theorem 2.
T n (x) = k n f kn (x) -f (x) f (x)
does not converge in distribution.

Proof :

If the sequence k 1+d/2αx n n
is not bounded, then there exists a subsequence

k 1+d/2αx q(n) q(n)
which tends to infinity. So

∀M > 0, P b q(n) > M → n→∞ 1
and (b n ), and consequently (T n (x)), do not converge in distribution.

If the sequence k 1+d/2αx n n is bounded and do not converge, then we can extract two subsequences which converge to two differents limits l 1 and l 2 . Using Lemma 2.5, the Alain F. Berlinet and Rémi Servien corresponding subsequences of (T n (x)) converge in distribution respectively to N (l 1 , 1) and N (l 2 , 1). Therefore, (Z n ) does not converge in distribution.

An example

We illustrate the above considerations with the probability density f defined on (-1, 1) by

f (x) = 5 2 1 -|x| 1/4 .
The associated distribution function satisfies

F (x) =    (1/2) + (5/2) x + 2 (-x) 5/4 if -1 ≤ x ≤ 0 (1/2) + (5/2) x -2 x 5/4 if 0 ≤ x ≤ 1.
Simple calculations show that for 0 < xh < x + h < 1, the small ball probability with radius h at x is given by

F (x + h) -F (x -h) = 5 h -5 h x 1/4 + 5 32 h 3 x 7/4 + O(h 5 ) = 2 h (f (x) + C x h αx + o(h αx )) with C x = (5/64) x -7/4 α x = 2
and that for 0 < h < 1, one has

F (h) -F (-h) = 5 h -4 h 5/4 = 2 h (f (0) + C 0 h α0 ) with C 0 = -2 α 0 = 1/4.
If k n is chosen as c n β where c and β are two positive constants (β < 1), we know from Theorem 2.4 that for x = 0, (T n (x)) converges in distribution if and only if β ≤ 4/5. The asymptotic distribution is N 5 -2 2 -1/4 c 5/2

(2-2 3/4 ) 3 , 1 if β = 4/5 and N (0, 1) if 0 < β < 4/5. From the same theorem we know that (T n (0)) converges in distribution if and only if β ≤ 1/3. The asymptotic distribution is N (-4 c 3/4 5 -5/4 , 1) if β = 1/3 and N (0, 1) if 0 < β < 1/3. Note also that the optimal k n in the sense of pointwise asymptotic mean square error at any point x of (-1, 1) different from 0 is of the form c n 4/5 (see Fukunaga and Hostetler [START_REF] Fukunaga | Optimization of k-nearest-neighbor density estimates[END_REF]) and that this imply the convergence of T n (x) to a non centered gaussian distribution.

A rule of thumb for k n was early proposed by Loftsgaarden and Quesenberry [START_REF] Loftsgaarden | A nonparametric estimate of a multivariate density function[END_REF]. It consists in taking k n as an integer close to √ n. This rule ensures that the conditions of

T n (x) = k n f kn (x) -f (x) f (x)
tends to N (0, 1) as n tends to infinity whenever the conditions of Theorem 2.1 on f and its derivatives are satisfied. This is why the rule "k n close to √ n"is popular in nearest neighbour density estimation (Figure 1 described below was built with such k n ).

To illustrate the behaviour of the statistic T n (x) (x = 0 and x = 0.5) we simulated m times a set S i (1 ≤ i ≤ m) of n pseudo-random variables with density f by means of the rejection method. Each set S i of n variables was used to compute f kn (x), k n being equal to the integer part of n γ , γ ∈ 1 4 , 1 3 , 1 2 , 4 5 . Then T n (x) was computed by using the exact value of f (x). In this way we got m values of T n (x) and used them to estimate its distribution by a kernel estimate with gaussian kernel

K(t) = 1 √ 2π exp - t 2 2
and window-width

h m = 4 3m 1/5 s m
where s m is the standard deviation of the sample (see Bosq and Lecoutre (3)).

The resulting eight estimators (with n = m = 1000) are drawn in Figures 1 to 4 together with the density of N (0, 1). The estimated density of T n (1/2) (resp. T n (0)) is in dotted (resp. dashed) line and the density of N (0, 1) is in solid line. The choice of k n implies the following situations. Figure 1.

The distribution of T n (1/2) converges to N (0, 1). The distribution of T n (0) does not converge. Figure 2.

The distribution of T n (1/2) converges to N (1.043..., 1). The distribution of T n (0) does not converge. Figure 3.

The distribution of T n (1/2) converges to N (0, 1). The distribution of T n (0) converges to N (-0.534..., 1). Figure 4.

The distribution of T n (1/2) and the distribution of T n (0) converge to N (0, 1).

The lesson to draw for the practitioner is that the choice of k n should be made locally with great care (this is known but deserves reminder) and that whenever the set of data is large enough a preliminary estimation of α x , as proposed by Beirlant, Berlinet and Biau (1), may act as a safeguard against hasty conclusions.

Lemma 2 . 8 :

 28 Necessary condition. Under the conditions of Theorem 2.4, if the sequence k 1+d/(2αx) n n does not converge then the sequence of random variables

Figure 1 .Figure 2 .

 12 Figure 1. k n equal to the integer part of n 1/2 .

Figure 3 .Figure 4 .

 34 Figure 3. k n equal to the integer part of n 1/3 .