Minh-Lien Jeanne Nguyen

Nonparametric method for sparse conditional density estimation in moderately large dimensions

Keywords: conditional density, high dimension, minimax rates, kernel density estimators, greedy algorithm, sparsity, nonparametric inference

In this paper, we consider the problem of estimating a conditional density in moderately large dimensions. Much more informative than regression functions, conditional densities are of main interest in recent methods, particularly in the Bayesian framework (studying the posterior distribution, finding its modes...). Considering a recently studied family of kernel estimators, we select a pointwise multivariate bandwidth by revisiting the greedy algorithm Rodeo (Regularisation Of Derivative Expectation Operator). The method addresses several issues: being greedy and computationally efficient by an iterative procedure, avoiding the curse of high dimensionality under some suitably defined sparsity conditions by early variable selection during the procedure, converging at a quasi-optimal minimax rate.

Introduction

Motivations

In this paper, we consider the problem of the conditional density estimation. We observe a n-sample of a couple (X, Y ), in which Y is the vector of interest while X gathers auxiliary variables. We denote d the joint dimension. In particular we are interested in the inference of the d-dimensional conditional density f of Y conditionally to X. There is a growing demand for methods of conditional density estimation in a wide spectrum of applications such as Economy [START_REF] Hall | Cross-validation and the estimation of conditional probability densities[END_REF], Cosmology [START_REF] Izbicki | Nonparametric conditional density estimation in a high-dimensional regression setting[END_REF], Medicine [START_REF] Takeuchi | Nonparametric conditional density estimation using piecewise-linear solution path of kernel quantile regression[END_REF]], Actuaries [START_REF] Efromovich | Oracle inequality for conditional density estimation and an actuarial example[END_REF]], Meteorology [START_REF] Jeon | Using conditional kernel density estimation for wind power density forecasting[END_REF] among others. It can be explained by the double role of the conditional density estimation: deriving the underlying distribution of a dataset and determining the impact of the vector X of auxiliary variables on the vector of interest Y . In this aspect, the conditional density estimation is richer than both the unconditional density estimation and the regression problem. In particular, in the regression framework, only the conditional mean E [Y |X] are estimated instead of the full conditional density, which can be especially poorly informative in case of an asymmetric or multi-modal conditional density. Conversely, from the conditional density estimators, one can, e.g., derive the conditional quantiles [START_REF] Takeuchi | Nonparametric quantile estimation[END_REF]] or give accurate predictive intervals [START_REF] Fernández-Soto | Error analysis of the photometric redshift technique[END_REF]. Furthermore, since the posterior distribution in the Bayesian framework is actually a conditional density, the present paper also offers an alternative method to the ABC methodology (for Approximate Bayesian Computation) [START_REF] Beaumont | Approximate bayesian computation in population genetics[END_REF][START_REF] Marin | Approximate bayesian computation methods[END_REF][START_REF] Biau | New insights into approximate bayesian computation[END_REF] in the case of an intractable-yet-simulable model. The challenging issue in conditional density estimation is to circumvent the "curse of dimensionality". The problem is twofold: theoretical and practical. In theory, it is stigmatized by the minimax approach, stating that in a d-dimensional space the best convergence rate for the pointwise risk over a p-regular class of functions is O(n -p 2p+d ): in particular, the larger is d, the slower is the rate. In practice, the larger the dimension is, the larger the sample size is needed to control the estimation error. In order to maintain reasonable running times in moderately large dimensions, methods have to be designed especially greedy. Furthermore, one interesting question is how to retrieve the eventual relevant components in case of sparsity structure on the conditional density f . For example, if we have at disposal plenty of auxiliary variables without any indication on their dependency with our vector of interest Y , the ideal procedure will take in input the whole dataset and still achieve a running time and a minimax rate as fast as if only the relevant components were given and considered for the estimation. More precisely, two goals are simultaneously addressed : converging at rate O(n -2p 2p+r ) with r the relevant dimension, i.e. the number of components that influence the conditional density f , and detect the irrelevant components at an early stage of the procedure in order to afterwards only work on the relevant data and thus speed up the running time.

Existing methodologies

Several nonparametric methods have been proposed to estimate conditional densities: kernel density estimators [START_REF] Rosenblatt | Conditional probability density and regression estimators[END_REF][START_REF] Hyndman | Estimating and visualizing conditional densities[END_REF][START_REF] Bertin | Adaptive pointwise estimation of conditional density function[END_REF]] and various methodologies for the selection of the associated bandwidth [START_REF] Bashtannyk | Bandwidth selection for kernel conditional density estimation[END_REF][START_REF] Fan | A crossvalidation method for estimating conditional densities[END_REF][START_REF] Hall | Cross-validation and the estimation of conditional probability densities[END_REF]; local polynomial estimators [START_REF] Fan | Estimation of conditional densities and sensitivity measures in nonlinear dynamical systems[END_REF][START_REF] Hyndman | Nonparametric estimation and symmetry tests for conditional density functions[END_REF]; projection series estimators [START_REF] Efromovich | Nonparametric Curve Estimation: Methods, Theory and Applications[END_REF][START_REF] Efromovich | Conditional density estimation in a regression setting[END_REF]; piecewise constant estimator [START_REF] Györfi | Nonparametric estimation of conditional distributions[END_REF][START_REF] Sart | Estimating the conditional density by histogram type estimators and model selection[END_REF]; copula [START_REF] Faugeras | A quantile-copula approach to conditional density estimation[END_REF]]. But while most of the aforementioned works are only defined for bivariate data or at least when either X or Y is univariate, they are also computationally intractable as soon as d > 3. It is in particular the case for the kernel density methodologies [START_REF] Hall | Cross-validation and the estimation of conditional probability densities[END_REF][START_REF] Bertin | Adaptive pointwise estimation of conditional density function[END_REF]: they achieve the optimal minimax rate, and even the detection of the relevant components, thanks to an adequate choice of the bandwidth (for the two aforementioned methods by cross validation and Goldenshluger-Lepski methodology), but the computational cost of these bandwidth selections is prohibitive even for moderate sizes of n and d. To the best of our knowledge, only two kernel density methods have been proposed to handle large datasets. [START_REF] Holmes | Fast kernel conditional density estimation: A dual-tree monte carlo approach[END_REF] propose a fast method of approximated cross-validation, based on a dual-tree speed-up, but they do not establish any rate of convergence and only show the consistency of their method. For scalar Y , [START_REF] Fan | Approximating conditional density functions using dimension reduction[END_REF] proposed to perform a prior step of dimension reduction on X to bypass the curse of dimensionality, then they estimate the bivariate approximated conditional density by kernel estimators. But the proved convergence rate n -1 3 is not the optimal minimax rate n -3 8 for the estimation of a bivariate function of assumed regularity 3. Moreover, the step of dimension reduction restricts the dependency of X to a linear combination of its components, which may induce a significant loss of information. Projection series methods for scalar Y have also been proposed. [Efromovich 2010a] extends his previous work [START_REF] Efromovich | Conditional density estimation in a regression setting[END_REF]] to a multivariate X. Theoretically the method achieves an oracle inequality, thus the optimal minimax rate. Moreover it performs an automatic dimension reduction on X when there exists a smaller intrinsic dimension. To our knowledge, it is the only method which addresses datasets of dimension larger than 3 with reasonable running times and does not pay its numerical performance with non optimal minimax rates. However the computation cost is prohibitive when both n and d are large. More recently, Izbicki and Lee have proposed two methodologies using orthogonal series estimators [START_REF] Izbicki | Nonparametric conditional density estimation in a high-dimensional regression setting[END_REF][START_REF] Sart | Estimating the conditional density by histogram type estimators and model selection[END_REF]. The first method is particularly fast and can handle very large X (with more than 1000 covariates). Moreover the convergence rate adapts to an eventual smaller unknown intrinsic dimension of the support of the conditional density. The second method originally proposes to convert successful high dimensional regression methods into the conditional density estimation, interpreting the coefficients of the orthogonal series estimator as regression functions, which allows to adapt to all kind of figures (mixed data, smaller intrinsic dimension, relevant variables) in function of the regression method. However both methods converge slower than the optimal minimax rate. Moreover their optimal tunings depend in fact on the unknown intrinsic dimension. For multivariate X and Y , [START_REF] Otneim | Conditional density estimation using the local gaussian correlation[END_REF] propose a new semiparametric method, called Locally Gaussian Density Estimator: they rewrite the conditional density as a product of a function depending on the marginal distribution functions (easily estimated since univariate, then plug-in), and a term which measures the dependency between the components, which is approximated by a centred Gaussian whose covariance is parametrically estimated. Numerically, the methodology seems robust to addition of covariates of X independent of Y , but it is not proved. Moreover they only establish the asymptotic normality of their method.

Our strategy and contributions

The challenge in this paper is to handle large datasets, thus we assume at our disposal a sample of large size n and of moderately large dimension. Then our work is motivated by the following three objectives:

(i) achieving the optimal minimax rate (up to a logarithm term);

(ii) being greedy, meaning that the procedure must have reasonable running times for large n and moderately large dimensions, in particular when d > 3 ;

(iii) adapting to a potential sparsity structure of f . More precisely, in the case where f locally depends only on a number r of its d components, r can be seen as the local relevant dimension. Then the desired convergence rate has to adapt to the unknown relevant dimension r: under this sparsity assumption, the benchmark for the estimation of a p-regular function is to achieve a convergence rate of the order

O(n - 2p
2p+r ), which is the optimal minimax rate if the relevant components were given by an oracle.

Our strategy is based on kernel density estimators. The considered family has been recently introduced and studied in [START_REF] Bertin | Adaptive pointwise estimation of conditional density function[END_REF]]. This family is especially designed for conditional densities and is better adapted for the objective (iii) than the intensively studied estimator built as the ratio of a kernel estimator of the joint density over one of the marginal density of X. For example, a relevant component for the joint density and the marginal density of X may be irrelevant for the conditional density and it is the case if a component of X is independent of Y . Note though that many more cases of irrelevance exist since we define the relevance as a local property.

The main issue with kernel density estimators is the selection of the bandwidth h ∈ R d + , and in our case, we also want to complete the objective (ii), since the pre-existing methodologies of bandwidth selection does not satisfy this restriction and thus cannot handle large datasets. In this paper, it is performed by an algorithm we call CDRodeo, which is derived from the algorithm Rodeo [START_REF] Lafferty | Rodeo: Sparse, greedy nonparametric regression[END_REF][START_REF] Liu | Sparse nonparametric density estimation in high dimensions using the rodeo[END_REF]], which has respectively been applied for the regression and the unconditional density estimation. The greediness of the algorithm allows us to address datasets of large sizes while keeping a reasonable running time (see Section 3.5 for further details). We give a simulated example with a sample of size n = 10 5 and of dimension d = 5 in Section 4. Moreover, Rodeo-type algorithms ensure an early detection of irrelevant component, and thus achieve the objective (iii) while improving the objective (ii). From the theoretical point of view, if the regularity of f is known, our method achieves an optimal minimax rate (up to a logarithmic factor), which is adaptive to the unknown sparsity of f . The last property is mostly due to the Rodeo-type procedures. The improvement of our method in comparison to the paper [START_REF] Liu | Sparse nonparametric density estimation in high dimensions using the rodeo[END_REF] which estimates the unconditional density with Rodeo is twofold. First, our result is extended to any regularity p ∈ N >0 , whereas [START_REF] Liu | Sparse nonparametric density estimation in high dimensions using the rodeo[END_REF]] fixed p = 2. Secondly, our notion of relevance is both less restrictive and more natural. In [START_REF] Liu | Sparse nonparametric density estimation in high dimensions using the rodeo[END_REF], they studied the L 2 -risk of their estimator, therefore they have to consider a notion of global relevance, whereas we consider a pointwise approach, which allows us to define a local property of relevance, which can be applied to a broader class of functions. Moreover, their notion of relevance is not intrinsic to the unknown density, but in fact depends on a tuning of the method, a prior chosen baseline density which has no connexion with the density, which limits the interpretation of the relevance.

Overview

Our paper is organized as follows. We introduce the CDRodeo method in Section 2. The theoretical results are in Section 3, in which we specify the assumptions and the tunings of the procedure from which are derived the convergence rate and the complexity cost of the method. A numerical example is presented in Section 4. The proofs are in the last section.

CDRodeo method

Let W 1 , . . . , W n be a sample of a couple (X, Y ) of multivariate random vectors: for i = 1, . . . , n,

W i = (X i , Y i ), with X i valued in R d 1 and Y i in R d 2 . We denote d := d 1 + d 2 the joint dimension.
We assume that the marginal distribution of X and the conditional distribution of Y given X are absolutely continuous with respect to the Lebesgue measure, and we define f : R d → R such as for any x ∈ R d 1 , f (x, •) is the conditional density of Y conditionally to X = x. We denote f X the marginal density of X. Our method estimates f pointwisely : let us fix w = (x, y) ∈ R d the point of interest. Kernel estimators. Our method is based on kernel density estimators. More specifically, we consider the family proposed in [START_REF] Bertin | Adaptive pointwise estimation of conditional density function[END_REF], which is especially designed for the conditional density estimation. Let K : R → R be a kernel function, ie:

R K(t)dt = 1, then for any bandwidth h ∈ (R * + ) d , the estimator of f (w) is defined by:

fh (w) := 1 n n i=1 1 fX (X i ) d j=1 h -1 j K w j -W ij h j , ( 1
)
where fX is an estimator of f X , built from another sample X of X. We denote by n X the sample size of X.

The choices of K and fX are specified later (see section 3.2).

Bandwidth selection. In kernel density estimation, selecting the bandwidth is a critical choice which can be viewed as a bias-variance trade-off. In [START_REF] Bertin | Adaptive pointwise estimation of conditional density function[END_REF], it is performed by the Goldenshluger-Lepski methodology (see [START_REF] Goldenshluger | Bandwidth selection in kernel density estimation: oracle inequalities and adaptive minimax optimality[END_REF]) and requires an optimization over an exhaustive grid of couples (h, h ) of bandwidths, which leads to intractable running time when the dimension exceeds 3 (and large dataset).

That is why we focus in a method which excludes optimization over an exhaustive grid of bandwidths to rather propose a greedy algorithm derived from the algorithm Rodeo. First introduced in the regression framework [START_REF] Wasserman | Rodeo: Sparse nonparametric regression in high dimensions[END_REF][START_REF] Lafferty | Rodeo: Sparse, greedy nonparametric regression[END_REF], a variation of Rodeo was proposed in [START_REF] Liu | Sparse nonparametric density estimation in high dimensions using the rodeo[END_REF] for the density estimation. Our method we called CDRodeo (for Conditional Density Rodeo) addresses the more general problem of conditional density estimation. Like Rodeo (which means Regularisation Of Derivative Expectation Operator), the CDRodeo algorithm generates an iterative path of decreasing bandwidths, based on tests on the partial derivatives of the estimator with respect to the components of the bandwidth. Note that the greediness of the procedure leans on the selection of this path of bandwidths, which enables us to address high dimensional problems of functional inference.

Let us be more precise: we take a kernel K of class C 1 and consider the statistics Z hj for h ∈ (R * + ) d and j = 1 : d, defined by:

Z hj := ∂ ∂h j fh (w).
Z hj is easily computable, since it can be expressed by:

Z hj = -1 nh 2 j n i=1 1 fX (X i ) J w j -W ij h j d k =j h -1 k K w k -W ik h k , (2)
where J : R → R is the function defined by:

t → K(t) + tK (t). (3)
The details of the CDRodeo procedure are described in Algorithm 1 and can be summed up in one sentence: for a well-chosen threshold λ hj (specified in Section 3.3), the algorithm performs at each iteration the test |Z hj | > λ hj to determine if the component j of the current bandwidth must be shrunk or not. It can be interpreted by the following principle: the bandwidth of a kernel estimator quantifies within which distance of the point of interest w and at which degree an observation W i helps in the estimation. Heuristically, the larger the variation of f is, the smaller the bandwidth is required for an accurate estimation. The statistics Z hj = ∂ ∂h j fh (w) are used as a proxy of ∂ ∂w j f (w) to quantify the variation of f in the direction w j . Note in particular that since the partial derivatives vanish for irrelevant components, this bandwidth selection leads to an implicit variable selection, and thus to avoid the curse of dimensionality under sparsity assumptions.

Algorithm 1 CDRodeo algorithm 1. Input: the point of interest w, the data W , β ∈ (0, 1) the bandwidth decreasing factor, h 0 > 0 the bandwidth initialization value, a parameter a > 1.

Initialization:

(a) Initialize the bandwidth: for j = 1 : d, h j ← h 0 .

(b) Activate all the variables: A ← {1, . . . , d}.

While

(A = ∅) & ( d k=1 h k ≥ log n n ):
for all j ∈ A:

(a) Update Z hj and λ hj . (b) If |Z hj | ≥ λ hj : update h j ← βh j .
else: remove j from A.

4.

Output: h (and fh (w)).

Theoretical results

This section gathers the theoretical results of our method.

Assumptions

We consider K a compactly supported kernel. For any bandwidth h ∈ (R * + ) d , we define the neighbourhood U h (u) of u ∈ R d (typically, u = x or w, and d = d 1 or d) as follows:

U h (u) := u ∈ R d : ∀j = 1 : d , u j = u j -h j z j , with z ∈ (supp(K)) d .
Then we denote the CDRodeo initial bandwidth h (0) = 1 log n , . . . , 1 log n and for short, U n (u) := U h (0) (u). We also introduce the notation • ∞, U for the supremum norm over a set U.

The following first assumption ensures a certain amount of observations in the neighbourhood of our point of interest w.

Assumption 1 (f X bounded away of 0). We assume δ := inf

u∈ Un(x) f X (u) > 0.
Note that if the neighbourhood U n (x) does not contain any observation X i , the estimation of the conditional distribution of Y given the event X = x is obviously intractable.

The second assumption specifies the notions of "sparse function" and "relevant component", under which the curse of high dimensionality can be avoided.

Assumption 2 (Sparsity condition).

There exists a subset R ∈ {1, . . . , d} such that for any fixed {z j } j∈R , the function

{z k } k∈R c → f (z 1 , . . . , z d ) is constant on U n (w).
In other words, if we denote r the cardinal of R, Assumption 2 means that f locally depends on only r of its d variables. We call relevant any component in R. The notion of relevant component depends on the point where f is estimated. For example, a component w j which behaves as 1 [0,1] (w j ) in the conditional density is only relevant in the neighbourhood of 0 and 1. Note that this local property addresses a broader class of functions, which extends the application field of Theorem 2 and improves the convergence rate of the method.

Finally, the conditional density is required to be regular enough.

Assumption 3 (Regularity of f ). There exists a known integer p such that f is of class C p on U n (w) and such that ∂ p j f (w) = 0 for all j ∈ R.

Conditions on the estimator of f X

Given the definition of the estimator (1), we need an estimator fX of f X .

If f X is known. We take fX ≡ f X . This case is not completely obvious. In particular, it tackles the case of unconditional density estimation, if we set by convention d 1 = 0 and f X ≡ 1.

If f X is unknown. We need an estimator fX which satisfies the following two conditions:

(i) a positive lower bound:

δ X := inf u∈ Un(x)
fX (u) > 0

(ii) a concentration inequality in local sup norm: there exists a constant M X > 0 such that:

P sup u∈ Un(x) f X (u) -fX (u) fX (u) > M X (log n) d 2 n 1 2 ≤ exp(-(log n) 5 4
).

The following proposition proves these conditions are feasible. Furthermore, the provided estimator of f X (see the proof in Section 5.3.1) is easily implementable and does not need any optimisation.

Proposition 1. Given a sample X with same distribution as X and of size n

X = n c with c > 1, if f X is of class C p with p ≥ d 1 2(c-1)
, there exists an estimator fX which satisfies (i) and (ii).

CDRodeo parameters choice.

Kernel K. We choose the kernel function K : R → R of class C 1 , with compact support and of order p, i.e.: for = 1, . . . , p -1, R t K(t)dt = 0, and R t p K(t)dt = 0. Note that considering a compactly supported kernel is fundamental for the local approach. In particular, it relaxes the assumptions by restricting them to a neighbourhood of w. Taking a kernel of order p is usual for the control of the bias of the estimator.

Parameter β. Let β ∈ (0, 1) be the decreasing factor of the bandwidth. The larger β, the more accurate the procedure, but the longer the computational time. From the theoretical point of view, it remains of little importance, as it only affects the constant terms. In practice, we set it close to 1.

Bandwidth initialization. We recall that we set h 0 := 1 log n (and the initial bandwidth as

1 log n , . . . , 1 log n ).
Threshold λ h,j . For any bandwidth h ∈ (R * + ) d and for j = 1 : d, we set the threshold as follows:

λ hj := C λ (log n) a nh 2 j d k=1 h k , ( 4
) with C λ := 4 J 2 K d-1
2 (where J is defined in (3)) and a > 1. The expression is obtained by using concentration inequalities on Z hj . For the proof, the parameter a has to be tuned such that:

(log n) a-1 > f ∞, Un(w) δ , (5)
which is satisfied for n large enough. The influence of this parameter is discussed in the next section, once the theoretical results are stated.

Hereafter, unless otherwise specified, the parameters are chosen as described in this section.

Mains results

Let us denote ĥ the bandwidth selected by CDRodeo. In Theorem 2, we introduce a set H hp of bandwidths which contains ĥ with high probability, which leads to an upper bound of the pointwise estimation error with high probability. In Corollary 3, we deduce the convergence rate of CDRodeo from Theorem 2. More precisely, in Theorem 2, we determine lower and upper bounds (with high probability) for the stopping iteration of each bandwidth component. We set:

τ n := 1 (2p + r) log 1 β log n C τ (log n) 2p+d+a , (6)
and

T n := τ n + log C -1 T (2p + 1) log 1 β , (7)
where

C τ :=     4(p -1)!C λ min j∈R ∂ p j f (w) R t p K(t)dt     2 , C T :=   min j∈R |∂ p j f (w)| 24 max j∈R |∂ p j f (w)|   2 .
Then we define the set of bandwidths H hp by:

H hp := h ∈ R d + : h j = β θ j log n , with θ j ∈ { τ n + 1, . . . , T n } if j ∈ R, else θ j = 0 .
Theorem 2. Assume that fX satisfies Conditions (i) and (ii) of section 3.2 and Assumptions 1 to 3 are satisfied. Then, the bandwidth ĥ selected by CDRodeo belongs to H hp with high probability. More precisely, for any q > 0 and for n large enough:

P ĥ ∈ H hp ≥ 1 -n -q . ( 8
)
Moreover, with probability larger than 1 -2n -q , the CDRodeo estimator fĥ (w) verifies:

fĥ (w) -f (w) ≤ C(log n) p 2p+r (d-r+a) n -p 2p+r (9) with C := 2rC p 2p+r τ t∈R | t p p! K(t)|dt × max k∈R ∂ p k f ∞, Un(w) + 4 K d 2 f 1 2 ∞, Un(w) δ -1 2 C -r 2(2p+1) T C -r 2(2p+r) τ .
Corollary 3. Under the assumptions of Theorem 2, for any q ≥ 1:

E fĥ (w) -f (w) q 1/q ≤ C(log n) p 2p+r (d-r+a) n -p 2p+r + o n -1 .
Corollary 3 presents a generalization of the previous works on Rodeo [START_REF] Lafferty | Rodeo: Sparse, greedy nonparametric regression[END_REF]] and [START_REF] Liu | Sparse nonparametric density estimation in high dimensions using the rodeo[END_REF] whose results are restricted to the regularity p = 2 and to simpler problems, namely regression and density estimation. We compare the convergence rate of CDRodeo with the optimal minimax rate. In particular, our benchmark is the pointwise minimax rate, which is of order O n -p 2p+d , for the problem of p-regular d-dimensional density estimation, obtained by [START_REF] Donoho | Renormalization exponents and optimal pointwise rates of convergence[END_REF]. Without sparsity structure (r = d), CDRodeo achieves the optimal minimax rate, up to a logarithmic factor. The exponent of this factor depends on the parameter a. For the proofs, we need a > 1 in order to satisfy ( 5), but if an upper bound (or a pre-estimator) of

f ∞, Un(w) δ
were known, we could obtain the similar result with a = 1 and a modified constant term. Note that the logarithmic factor is a small price to pay for a computationallytractable procedure for high-dimensional functional inference, in particular see section 3.5 for the computational gain of our procedure.

Under sparsity assumptions, we avoid the curse of high dimensionality and our procedure achieves the desired rate n -p 2p+r (up to a logarithmic term), which is optimal if the relevant components were known. Note that some additional logarithmic factors could be unavoidable due to the unknown sparsity structure, which needs to be estimated. Identifying the exact order of the logarithm term in the optimal minimax rate for the sparse case remains an open challenging question.

Complexity

We now discuss the complexity of CDRodeo without taking into account the pre-computation cost of fX at the points X i , i = 1 : n (used for computing the Z hj ), but a fast procedure for fX is required, to avoid losing CDRodeo computational advantages. For CDRodeo, the main cost lies in the computation of the Z hj 's along the path of bandwidths.

The condition

d k=1 h k ≥ log n
n restricts to at most log β -1 n updates of the bandwidth across all components, leading to a worst-case complexity of order O(d.n log n). But as shown in Theorem 2, with high probability, ĥ ∈ H hp , in which only the relevant components are active after the first iteration. In first iteration, the Z h (0) j 's computation costs O(d.n) operations, while the product kernel enables us to compute the Z hj 's in following iteration with only O(r.n) operations, which leads to the complexity O(d.n + r.n log n). In order to grasp the advantage of CDRodeo greediness, we compare its complexity with optimization over an exhaustive bandwidth grid with log n values for each component of the bandwidth (which is often the case in others methods: Cross validation, Lepski methods...): for each bandwidth of (log n) d -sized grid, the computation of a statistic from the d.n-sized dataset needs at least O(d.n) operation, which leads to a complexity of order O(d.n(log β -1 n) d ). Using the parameters used in the simulated example in section 4 (n = 2.10 5 ,

d = 5, r = 3, β = 0.95),

Simulations

In this section, we test the practical performances of our method. In particular, we study CDRodeo on a 5dimensional example. The major purpose of this section is to assess if the numerical performances of our procedure. Let us describe the example. We set d 1 = 4 and d 2 = 1 and simulate an i.i.d sample {(X i , Y i )} n i=1 with the following distribution: for any i = 1, . . . , n:

-the first component X i1 of X i follows a uniform distribution on [-1, 1],
-the other components X ij , j = 2 : 4, are independent standard normal and are independent of

X i1 , -Y i is independent of X i1 , X i3 and X i4 and the conditional distribution of Y i given X i2 is exponential with survival parameter X 2 i2 .
The estimated conditional density function is then defined by:

f : (x, y) → 1 [-1,1] (x 1 ) 1 x 2 2 e -y x 2 2 .
This example enables us to test several criteria: sparsity detection, behaviour when fonctions are not continuous, bimodality estimation, robustness when f X takes small values.

In the following simulations, if not stated explicitly otherwise, Rodeo is run with sample size n = 200, 000, product Gaussian kernel, initial bandwidth value h 0 = 0.4, bandwidth decreasing factor β = 0.95 and parameter a = 1.1 and fX ≡ f X . Figure 1 illustrates CDRodeo bandwidth selection. In which, the boxplots of each selected bandwidth component are built from 200 runs of CDRodeo at the point w = (0, 1, 0, 0, 1). This figure reflects the specificity of CDRodeo to capture the relevance degree of each component, and one could compare it with variable selection (as done in [START_REF] Lafferty | Rodeo: Sparse, greedy nonparametric regression[END_REF]). The components x 3 and x 4 are irrelevant and for this point of interest, the components x 2 and y are clearly relevant while the component x 1 is barely relevant as f is constant in the direction x 1 in near neighbourhood of x 1 = 0. As expected, the irrelevant h 3 and h 4 are mostly deactivated at the first iteration, while the relevant h 2 and h 5 are systematically shrunk. The relevance degree of x 1 is also well detected as the values of h 1 are smaller than h 0 , but significantly larger than h 2 and h 5 .

Figure 1: Boxplots of each component of 200 CDRodeo selected bandwidths at the point w = (0, 1, 0, 0, 1).

Figure 2 gives CDRodeo estimation of f from one n-sample. The function f is well estimated. In particular, irrevance, jumps and bi-modality are features which are well detected by our method. As expected, main estimation errors are made on points of discontinuity for x 1 and y or at the boundaries for x 2 , x 3 and x 4 . Note that the f X values are particularly small at the boundaries of the plots in function of x, leading to lack of observations for the estimation. Note however that null value for f X does not deteriorate the estimation (cf top left plot), since the estimate of f vanishes automatically when there is no observation near the point of interest.

Running time. The simulations are implemented in R on a Macintosh laptop with a 3, 1 GHz Intel Core i7 processor. In the Figure Figure 1, the 200 runs of CDRodeo take 2952.735 seconds (around 50 minutes), or 14.8 seconds per run.

Proofs

We first give the outlines of the proofs in Section 5.1. To facilitate the lecture of the proof, we have divided the proofs of the main results (Proposition 1, Theorem 2 and Corollary 3) into intermediate results which are stated in Section 5.2 and proved in Section 5.4. The proof of the main results are in Section 5.3.

Outlines of the proofs

We first prove Proposition 1 by constructing an estimator of f X with the wanted properties. In this proof, we use some usual properties of a kernel density estimator (control of the bias, concentration inequality), which are gathered in Lemma 4. Theorem 2 states two results: the bandwidth selection (8) and the estimation error of the procedure (9). For the proof of the bandwidth selection (8), Proposition 8 mades explicit the highly probable behaviour of CDRodeo along a run, and thus the final selected bandwidth. In particular, the proof leans on an analysis of Z hj , which is made in two steps. We first consider Zhj , a simpler version of Z hj in which we substitute the estimator of f X by f X itself, and we detail its behaviour in Lemma 6. Then we control the difference Z hj -Zhj (see in 1. of Lemma 7) to ensure Zhj behaves like Z hj .

To control the estimation error of the procedure (9), we similarly analyse fh (w) in two parts: in Lemma 5, we describe the behaviour of fh (w), the simpler version of fh (w) in which we substitute the estimator of f X by f X itself, and in 2. of Lemma 7, we bound the difference f h -fh (w). Then the bandwidth selection (8) leads to the upper bound with high probability of the estimation error of fh (w) (9). Finally, we obtain the expected error of fh (w) stated in Corollary 3 by controlling the error on the residual event.

Intermediate results

For any bandwidth h X ∈ R * + , we define the kernel density estimator f K X by: for any

u ∈ R d 1 , f K X (u) := 1 n X .h d 1 X n X i=1 d 1 j=1 K X u j -X ij h X , (10)
where

K X : R → R a kernel which is compactly supported, of class C 1 , of order p X ≥ d 1 2(c-1)
, where we recall that c > 1 is defined by n X = n c . We also introduce the neighbourhood

U n (x) := {u = u -h X z : u ∈ U n (x), z ∈ supp(K X )}. ( 11
) Lemma 4 ( f K X behaviour). We assume f X is C p on U n (x) with p ≤ p X , then for any bandwidth h X ∈ R * + , 1. if we denote C bias X := K X d 1 -1 1 • p K X (•) 1 p ! d 1 max k=1:d 1 ∂ p k f X ∞, U n (x) , then E f K X -f X ∞, Un(x)
≤ C bias X h p X .

If the condition Cond

X (h X ) : h d 1 X ≥ 4 K X 2d 1 ∞ 9 K X 2d 1 2 f X ∞, U n (x) (log n) 3 2 n X is satisfied, then for λ X := 4 K X 2d 1 2 f X ∞, U n (x) h d 1 X n X (log n) 3
2 and for any u ∈ U n (x):

P f K X (u) -E f K X (u) > λ X ≤ 2 exp -(log n) 3 2
.

Lemma 5 ( fh (w) behaviour). For any bandwidth h ∈ (0, h 0 ] d , and any i = 1 : n, let us denote fhi (w) :=

K h (w-W i ) f X (X i ) .
Then, if K is chosen as in section 3.3, under Assumptions 1 to 3,

1. Let CĒ := f ∞, Un(w) K d 1 . Then E fh1 (w) ≤ E fh1 (w) ≤ CĒ.
Besides, if we denote

B h := E fh (w) -f (w) the bias of fh (w) := 1 n n i=1
fhi (w), then:

|B h | ≤ C bias k∈R h p k ,
with C bias :=

2| t∈R t p K(t)dt| p! max k∈R |∂ p k f (w)|.

Let B

h := {| fh (w) -E[ fh (w)]| ≤ σ h }, where σ h := C σ (log n) a n d k=1 h k with C σ = 2 K d 2 f 1 2 ∞, Un(w) δ 1 2 . If Cond(h): d k=1 h k ≥ 4 K 2d ∞ 9δ 2 C 2 σ (log n) a n
is satisfied, then:

P B c h ≤ 2e -(log n) a 3. Let B | f |h := {| 1 n n i=1 | fhi (w)| -E[| fh (w)|]| ≤ CĒ}. Then P B c | f |h ≤ 2e -C γ|f | n d k=1 h k , with C γ|f | := min C 2 Ē C 2 σ ; 3δCĒ 4 K d ∞ .
Lemma 6 ( Zhj behaviour). For any j ∈ {1, . . . , d} and any bandwidth h ∈ (0, h 0 ] d , we define

Zhij := 1 f X (X i ) ∂ ∂h j d k=1 h -1 k K w k -W ik h k
, and Zhj := 1 n n i=1 Zhij . If K is chosen as in section 3.3, and under Assumptions 1 to 3, 1. Under Assumptions 1 to 3, for j / ∈ R:

E Zhj = 0.
whereas, for j ∈ R, for n large enough,

1 2 C E Z,j h p-1 j ≤ E Zhj ≤ 3 2 C E Z,j h p-1 j , ( 12
)
where

C E Z,j := R t p K(t)dt (p-1)! ∂ p j f (w) . Besides, let C E| Z| := f ∞, Un(w) J 1 K d-1
1 . Then :

E | Zh1j | ≤ C E| Z| h -1 j . ( 13
) 2. Let B Z,hj := {| Zhj -E Zhj | ≤ 1 2 λ hj }.
Under Assumptions 1 to 3, if the bandwidth satisfies:

Cond Z (h):

d k=1 h k ≥ cond Z (log n) a n , with cond Z := 4 J 2 ∞ K 2(d-1) ∞ 3 2 f ∞, Un(w) J 2 2 K 2(d-1) 2 , then: P B c Z,hj ≤ 2e -γ Z,n , with γ Z,n := δ f ∞, Un(w) (log n) a . 3. Let B | Z|,hj := {| 1 n n i=1 | Zhij | -E | Zh1j | | ≤ C E| Z| h -1 j }.
Then, under Assumptions 1 to 3:

P B c | Z|,hj ≤ 2e -C γ| Z| n d k=1 h k , with C γ| Z| := min δC 2 E| Z| 4 f ∞, Un(w) J 2 2 K 2(d-1) 2 ; 3δC E| Z| 4 K d-1 ∞ J ∞ .
Lemma 7. For any h ∈ (0, h 0 ] d and any component j = 1 : d, we denote ∆ Z,hj := Z hj -Zhj and ∆ h := fh (w) -fh (w). Under Assumptions 1 to 3, if the conditions on fX are satisfied (see section 3.2), then,

1. for C M∆Z := 2C E| Z| M X C λ : 1 B | Z|,hj ∩ An |∆ Z,hj | ≤ C M∆Z (log n) a 2 λ hj 2. for C M∆ := 2CĒM X Cσ : 1 An∩B | f |h |∆ h | ≤ C M∆ (log n) a 2 σ h .
We introduce the notation h (t) , t ∈ N, the state of the bandwidth at iteration t if ĥ = h. In particular for a fixed t ∈ {0, . . . , τ n }, h (t) is identical for any h ∈ H hp . Then we consider the event:

E Z := A n ∩ j / ∈R B Z,h (0) j ∩ B | Z|,h (0) j ∩ j∈R   h∈H hp B Z,hj ∩ B | Z|,hj ∩ τn t=0 B Z,h (t) j ∩ B | Z|,h (t) j   ,
where we denote A n := sup u∈ Un(x)

f X (u)-fX (u) fX (u) ≤ M X (log n X ) b n α X .
Proposition 8 (CDRodeo behaviour). Under Assumptions 1 to 3, on E Z , ĥ ∈ H hp .

In other words, when E Z happens:

1. non relevant components are deactivated during the iteration 0;

2. at the end of the iteration τ n , the active components are exactly the relevant ones;

3. CDRodeo stops at last at the iteration T n .

Moreover, for any q > 0:

P (E c Z ) = o(n -q ).
The following lemma give a technical result to canonically obtain an upper bound of the bias of a kernel estimator. Let us denote • the multiplication terms by terms of two vectors.

(p -1), R K(t)t l dt = 0. Let f 0 : R d → R a function of class C p on U h (u) := u ∈ R d : ∀j = 1 : d , u j = u j -h j z j , with z j ∈ supp(K j ) . Then: R d   d j=1 h -1 j K u j -u j h j   f 0 (u )du -f 0 (u) R d   d j=1 K (z j )   dz = d k=1 (I k + II k ), (14)
where

I k := z∈R d   d j=1 K(z j )   ρ k dz, with the notations ρ k := ρ k (z, h, u) = (-h k z k ) p 0≤tp≤•••≤t 1 ≤1 ∂ p k f 0 (z k-1 -t p h k z k e k ) -∂ p k f 0 (z k-1 ) dt 1:p , and z k-1 := u - k-1 j=1
h j z j e j (where {e j } d j=1 is the canonical basis of R d ), and

II k := (-h k ) p t∈R t p p! K(t)dt z -k ∈R d -1 ∂ p k f 0 (z k-1 )   j =k K(z j )   dz -k .
Finally, we recall (without proof) the classical Bernstein's Inequality and Taylor's theorem with integral remainder.

Lemma 10 (Bernstein's inequality). Let U 1 , . . . , U n be independent random variables almost surely uniformly bounded by a positive constant c > 0 and such that for i = 1, . . . , n, E[U 2 i ] ≤ v . Then for any λ > 0,

P 1 n n i=1 U i -E[U i ] ≥ λ ≤ 2 exp -min nλ 2 4v , 3nλ 4c .
Note that this version is a simple consequence of Birgé and Massart (p.366 of [Birgé and Massart 1998]).

Lemma 11 (Taylor's theorem). Let g : [0, 1] → R be a function of class C q . Then we have:

g(1) -g(0) = q l=1 g (l) (0) l! + 1 t 1 =0 t 1 t 2 =0 • • • t q-1 tq=0
(g (q) (t q ) -g (q) (0))dt q dt q-1 . . . dt 1

Proofs of Proposition 1, Theorem 2 and Corollary 3

Proof of Proposition 1

We construct fX in two steps: we first construct an estimator f K X which satisfies

P f X -f K X ∞, Un(x) > M X (log n) 3 4 n 1 2 ≤ exp(-(log n) 5 4 ), (15)
then we show that if we set fX ≡ f K X ∨ (log n) -1 4 , fX satisfies Conditions (i) and (ii) for n large enough.

We take f K X as the kernel density estimator defined in(10), with a kernel K X : R → R that is compactly supported, of class C 1 , of order p X ≥ d 1 2(c-1) . and a bandwidth h

X ∈ R * + specified later.Let us control the bias E f K X -f X ∞, Un(x)
. We define p X = min(p , p X + 1). In particular, f X is of class C p X and K X has p X -1 zero moments. Therefore we can apply Lemma 4 :

E f K X -f X ∞, Un(x) ≤ C bias X h p X X ,
where C bias X :=

K X d 1 -1 1 • p X K X (•) 1 p X ! d 1 max k=1:d 1 ∂ p X k f X ∞, U n (x) . Therefore, since f K X -f X ∞, Un(x) ≤ f K X -E f K X ∞, Un(x) + E f K X -f X ∞, Un(x) ≤ f K X -E f K X ∞, Un(x) + C bias X h p X X ,
we have for any threshold λ:

P f K X -f X ∞, Un(x) ≥ λ ≤ P f K X -E f K X ∞, Un(x) ≥ λ -C bias X h p X X . (16)
Therefore, we have reduced the problem to a local concentration inequality of f K X in sup norm. In order to move from a supremum on U n (x) to a maximum on a finite set of elements of U n (x), let us construct a -net of U n (x). We denote A > 0 such that:

supp(K X ) ∪ supp(K) ⊂ -A 2 , A 2 .
We set N ( ) the smallest integer such that N ( ) ≥ A log n , i.e.:

N ( ) := A log n ,
then we introduce the notation u (l) ∈ U n (x), for a multi-index l ∈ (1 : N ( )) d 1 defined, such that the j th component of u (l) is:

u (l)j := x j - A 2 log n + (2l j -1) 2 . Then {u (l) : l ∈ (1 : N ( )) d 1 } is a -net of U n (x)
, in the meaning that for any u ∈ U n (x), there exists l ∈ {1, . . . , N ( )} d 1 such that u -u (l) ∞ := max

k=1:d 1 |u k -u (l)k | ≤ .
Therefore to obtain the desired concentration inequality, we only need to obtain the concentration inequality for each point of {u (l) : l ∈ (1 : N ( )) d 1 } and to control the following supremum

sup u∈ Un(x) min l∈(1:N ( )) d 1 f K X (u) -E f K X (u) -f K X (u (l) ) + E f K X (u (l) ) .
For this purpose, we obtain (from Taylor's Inequality): for any u, v ∈ R d 1 ,

d 1 k=1 K X (u k ) - d 1 k=1 K X (v k ) ≤ d 1 K X ∞ K X d 1 -1 ∞ u -v ∞ .
Therefore, for any u, v ∈ U n (x):

f K X (u) -f K X (v) ≤ 1 n X .h d 1 X n X i=1 d 1 k=1 K X ( u k -X ik h X ) - d 1 k=1 K X ( u k -X ik h X ) ≤ d 1 h d 1 +1 X K X ∞ K X d 1 -1 ∞ u -v ∞ . Since {u (l) : l ∈ (1 : N ( )) d 1 } is a -net of U n (x): sup u∈ Un(x) min l∈(1:N ( )) d 1 f K X (u) -f K X (u (l) ) ≤ d 1 h d 1 +1 X K X ∞ K X d 1 -1 ∞
.

Thus:

sup u∈ Un(x) min l∈(1:N ( )) d 1 E f K X (u) -E f K X (u (l) ) ≤ d 1 K X ∞ K X d 1 -1 ∞ h d 1 +1 X .
And so:

sup u∈ Un(x) min l∈(1:N ( )) d 1 f K X (u) -E f K X (u) -f K X (u (l) ) + E f K X (u (l) ) ≤ 2d 1 K X ∞ K X d 1 -1 ∞ h d 1 +1 X . We denote C diff := 2d 1 K X ∞ K X d 1 -1 ∞ . Then: f K X -E f K X ∞, Un(x) ≤ max l∈(1:N ( )) d 1 f K X (u (l) ) -E f K X (u (l) ) + sup u∈ Un(x) min l∈(1:N ( )) d 1 f K X (u) -E f K X (u) -f K X (u (l) ) + E f K X (u (l) ) ≤ max l∈(1:N ( )) d 1 f K X (u (l) ) -E f K X (u (l) ) + C diff h d 1 +1 X .
Then the inequality ( 16) becomes: for any threshold λ,

P f K X -f X ∞, Un(x) ≥ λ ≤ P f K X -E f K X ∞, Un(x) ≥ λ -C bias X h p X X ≤ P max l∈(1:N ( )) d 1 f K X (u (l) ) -E f K X (u (l) ) ≥ λ -C bias X h p X X -C diff h d 1 +1 X ≤ N ( ) d 1 max l∈(1:N ( )) d 1 P f K X (u (l) ) -E f K X (u (l) ) ≥ λ -C bias X h p X X -C diff h d 1 +1 X (17)
We want to apply 2. of Lemma 4. Therefore we fix the following settings:

• h X := n -c-1 c.d 1 X • λ := 2λ X
, where λ X is the threshold in 2. of Lemma 4;

• := h 1+ d 1 2 X n -1 2 X .
For short, we denote

C λX := 2 K X d 1 2 f X 1 2
∞, U n (x) , so:

λ X = 4 K X 2d 1 2 f X ∞, U n (x) h d 1 X n X (log n) 3 2 = C λX (log n) 3 4 h - d 1 2 X n -1 2 X = C λX (log n) 3 4 n -1 2c X .
In particular, since we take p X ≥ d 1 2(c-1) and we assume p ≥ d 1 2(c-1) , then p X = min(p , p X ) ≥ d 1 2(c-1) . Hence we obtain for n large enough:

C bias X h p X X = C bias X n - p X (c-1) c.d 1 X ≤ C bias X n -1 2c X ≤ 1 2 λ X = C λX 2 (log n) 3 4 n -1 2c X .
and also, since c > 1:

C diff h d 1 +1 X = C diff h - d 1 2 X n -1 2 X = C diff n -1 2c X ≤ 1 2 λ X = C λX 2 (log n) 3 4 n -1 2c X .
Thus:

λ -C bias X h p X X -C diff h d 1 +1 X ≥ λ X ,
and the inequality (17) becomes:

P f K X -f X ∞, Un(x) ≥ λ ≤ N ( ) d 1 max l∈(1:N ( )) d 1 P f K X (u (l) ) -E f K X (u (l) ) ≥ λ X (18)
We verify that Cond X (h X ) is satisfied for n large enough:

h d 1 X = n -c-1 c X ≥ 4 K X 2d 1 ∞ 9 K X 2d 1 2 f X ∞, U n (x) (log n) 3 2 n X .
Then we can apply 2. of Lemma 4,

P f K X (u (l) ) -E f K X (u (l) ) > λ X ≤ 2 exp -(log n) 3 2
.

Thus the inequality (18) becomes:

P f K X -f X ∞, Un(x) ≥ λ ≤ 2N ( ) d 1 exp -(log n) 3 2 . ( 19
)
Let us control 2N ( )

d 1 : 2N ( ) d 1 = 2 A log n d 1 = 2     A h 1+ d 1 2 X n -1 2 X log n     d 1 = o n d 1 +1
X Then for n large enough:

2N ( ) d 1 exp -(log n) 3 2 ≤ exp -(log n) 5 4
.

Therefore:

P f K X -f X ∞, Un(x) ≥ λ ≤ exp -(log n) 5 4 Since λ = 2C λX (log n) 3 4 n -1 2c
X , we have obtained the desired concentration inequality (15) with M X = 2C λX .

Now we consider

fX ≡ f K X ∨ (log n) -1 4
. By construction, fX satisfies Condition (i). Let us show it also satisfies Condition (ii), i.e.:

P sup u∈ Un(x) f X (u) -fX (u) fX (u) > M X (log n) d 2 n 1 2 ≤ C X exp(-(log n) 5 4 ).
We write:

P sup u∈ Un(x) f X (u) -fX (u) fX (u) > M X (log n) d 2 n 1 2 = P ∃u ∈ U n (x), f X (u) -fX (u) > fX (u)M X (log n) d 2 n 1 2 ≤ P ∃u ∈ U n (x), f X (u) -fX (u) > (log n) -1 4 M X (log n) d 2 n 1 2 ≤ P f X (u) -fX (u) ∞, Un(x) > M X (log n) d 2 -1 4 n 1 2 Since d = d 1 + d 2 ≥ 2, d 2 -1 4 ≥ 3 4
, we obtain from the previously proved concentration inequality (15):

P sup u∈ Un(x) f X (u) -fX (u) fX (u) > M X (log n) d 2 n 1 2 ≤ P f K X -f X ∞, Un(x) ≥ M X (log n) 3 4 n 1 2 ≤ exp -(log n) 5 4
.

Proof of Theorem 2

We introduce E f :=

h∈H hp B h ∩ B | f |h
and denote E := E Z ∩ E f . On E, ĥ belongs to H hp (cf Proposition 8).

Thus:

1 E fĥ (w) -f (w) = 1 E h∈H hp 1 ĥ=h fh (w) -f (w) . ( 20
)
For any h ∈ H hp , we denote ∆ h := fh (w) -fh (w) and B h := E fh (w) -f (w), and we decompose the loss as follows:

fh (w) -f (w) ≤ |∆ h | + fh (w) -E fh (w) + B h . ( 21
) Using Lemma 7, since E ⊂ A n ∩ B | f |h : 1 E |∆ h | ≤ C M∆ (log n) a 2 σ h . ( 22
)
Moreover, by Lemma 5, since E ⊂ A n ∩ B h :

fh (w) -E fh (w) ≤ σ h (23) = C σ (log n) a n d k=1 h k ≤ C σ (log n) a nh d 0 β r(Tn-τn)+rτn = C σ C -r 2(2p+1) T C -r 2(2p+r) τ (log n) p(a+d-r) 2p+r n -p 2p+r . ( 24
)
And, also:

|B h | ≤ C bias k∈R h p k ≤ rC bias β pτn h p 0 = rC bias C p 2p+r τ (log n) p(a+d-r) 2p+r n -p 2p+r . ( 25
)
To conclude, 22) and ( 23)

1 E | fĥ (w) -f (w)| ≤ 1 E h∈H hp 1 ĥ=h fh (w) -f (w) , by (20) ≤ 1 E h∈H hp 1 ĥ=h |∆ h | + fh (w) -E fh (w) + B h , by (21) ≤ 1 E h∈H hp 1 ĥ=h 1 + C M∆ (log n) a 2 σ h + B h , by (
≤ 1 E h∈H hp 1 ĥ=h C(log n) p(a+d-r) 2p+r n -p
2p+r , by ( 24) and ( 25)

= 1 E C(log n) p(a+d-r) 2p+r n -p 2p+r ,
with for n large enough (ie:

C M∆ (log n) a 2 ≤ 1), C := rC bias C p 2p+r τ + 2C σ C -r 2(2p+1) T C -r 2(2p+r) τ .
It remains to give an upper bound on P (E c ). For any q > 0:

P (E c ) ≤ P (E c Z ) + P E c f ≤ o(n -q ) + h∈H hp P B c h + P B c | f |h , using Proposition 8 ≤ o(n -q ) + h∈H hp 2e -(log n) a + 2e -C γ|f | n d k=1 h k ,
using Lemma 5, since for any h ∈ H hp , Cond(h) is satisfied. Moreover:

n d k=1 h k ≥ nβ rTn h d 0 = C r 2p+1 T C r 2p+r τ (log n) ra-2p(d-r) (2p+r) n 2p 2p+r ≥ (log n) a C γ|f | ,
for n large enough. Hence:

P (E c ) ≤ o(n -q ) + |H hp |4e -(log n) a = o(n -q ),
for n large enough, since

|H hp | = ( T n -τ n ) r == 1 (2p+1)(log( 1 β )) log( C T Cτ ) + 1 r is finite.

Proof of Corollary 3

We consider the event

E = fĥ (w) -f (w) ≤ C(log n) p 2p+r (d-r+a) n -p 2p+r
for which we proved in Theorem 2:

P (E c ) = o(n -A ),
for any A > 0. For short, we denote R h := fh (w) -f (w) for any bandwidth h ∈ (R * + ) d . Then we decompose R ĥ as follows: R ĥ = 1 E R ĥ + 1 E c R ĥ.

By definition of E, we immediately obtain:

1 E R ĥ ≤ C(log n) p 2p+r (d-r+a) n -p 2p+r (26)
For the second term, we first bound fĥ (w) a.s. In CDRodeo procedure, the loop stops when the current bandwidth becomes too small:

d k=1 h k < (log n)
n . So the final bandwidth ĥ satisfies:

d k=1 ĥk ≥ (log n) β d n . Since fĥ (w) = 1 n n i=1 1 f X (X i ) d k=1 ĥ-1 k K w j -W ij ĥj , fĥ (w) ≤ β d K d ∞ δ n log n Hence: R ĥ ≤ f (w) + β d K d ∞ δ n (log n)
Therefore, for any q > 0, using (a + b) q ≤ 2 q-1 (a q + b q ):

E 1 E c (R ĥ) q ≤ P (E c )2 q-1 f (w) q + K q ∞ n q (log n) q = o n -A , (27)
for any A > 0 (since

P (E c ) = o(n -A +q )).
We conclude by combining ( 26) and ( 27):

E fĥ (w) -f (w) q 1/q ≤ C(log n) p 2p+r (d-r+a) n -p 2p+r + o n -1 . (28)

Proof of Proposition 8 and the lemmas

Proof of Proposition 8

First, note that the final state of the bandwidth determines exactly at which iteration each component has been deactivated: for a fixed bandwidth h ∈ (R * + ) d , if ĥ = h, we denote {θ k } d k=1 such as for k = 1 : d, h k = h 0 β θ k . In particular, θ k is the iteration of deactivation of the component k. We introduce the notation h (t) , t ∈ N, the state of the bandwidth at iteration t if ĥ = h. It implies that h (t) is exactly defined by: h

(t) k = β θ k ∧t h 0 for k = 1 : d.
Notice that for a fixed t ∈ {0, . . . , τ n }, h (t) is identical for any h ∈ H hp : by definition of H hp , h

(t) j = h 0 β t if j ∈ R, else h (t) j = h 0 .
We recall the definition

E Z := A n ∩ j / ∈R B Z,h (0) j ∩ B | Z|,h (0) j ∩ j∈R   h∈H hp B Z,hj ∩ B | Z|,hj ∩ τn t=0 B Z,h (t) j ∩ B | Z|,h (t) j   .
For any component j and any bandwidth h, we decompose Z hj as follows:

1 E Z Z hj = 1 E Z Zhj + 1 E Z ∆ Z,hj = 1 E Z E Zhj + 1 E Z ( Zhj -E Zhj ) + 1 E Z ∆ Z,hj (29)
1. Let us fix j / ∈ R and h = h (0) = (h 0 , . . . , h 0 ). Using 2. of Lemma 6, E Zhj = 0. Therefore:

1 E Z |Z hj | ≤ 1 E Z Zhj -E Zhj + 1 E Z |∆ Z,hj | ≤ 1 2 λ hj + 1 E Z |∆ Z,hj | ,
using 2. of Lemma 6, since E Z ⊂ B Z,hj . Now using 1. of Lemma 7, since E Z ⊂ B | Z|,hj ∩ A n , we obtain:

1 E Z |Z hj | ≤ 1 2 λ hj + C M∆Z (log n) a 2 λ hj .
Then for n large enough (ie:

(log n) a 2 > 2C M∆Z ), 1 E Z |Z hj | < λ hj .
In other words, when E Z happens, all irrelevant components deactivate at the iteration 0.

2. Let us show that E Z implies that the relevant components remain active until iteration τ n + 1.

It suffices to prove |Z h (t) j | > λ h (t) j , for any j ∈ R and any bandwidth h (t) , t = 0 : τ n . (Indeed, by induction: (h 0 , . . . , h 0 ) = h (0) , and since the irrelevant components deactivate at the iteration 0, if the current bandwidth at the iteration t is h (t) , then the fact that all the relevant components remain active for this bandwidth implies that the bandwidth at iteration t + 1 is h (t+1) ).

Let us fix j ∈ R, t = 0 : τ n and we denote h = h (t) . Using the decomposition (29), we obtain the following lower bound:

1 E Z |Z hj | ≥ 1 E Z E Zhj -Zhj -E Zhj -|∆ Z,hj | .
Then, combining:

-

E Zhj ≥ C E Z,j 2 h p-1 j (cf 1. of Lemma 6), -Zhj -E Zhj ≤ 1 2 λ hj , since E Z ⊂ B Z,hj (cf 2. of Lemma 6), -|∆ Z,hj | ≤ C M∆Z (log n) a 2 λ hj , since E Z ⊂ B | Z|,hj ∩ A n (cf 1. of Lemma 7),
we obtain:

1 E Z |Z hj | ≥ 1 E Z C E Z,j 2 h p-1 j - 1 2 λ hj - C M∆Z (log n) a 2 λ hj . Now let us show: 1 E Z |Z hj | ≥ 1 E Z λ hj . First, if n is large enough (ie (log n) a 2 ≥ 2C M∆Z ), then C M∆Z (log n) a 2 λ hj ≤ 1 2 λ hj .
Then it suffices to prove:

C E Z,j 2 h p-1 j ≥ 2λ hj ,
i.e.:

h 2p j d k=1 h k ≥ 4 2 C 2 λ C 2 E Z,j (log n) a n .
It is ensured for t ≤ τ n , by definition of τ n in (6):

h 2p j d k=1 h k = β t(2p+r) (log n) 2p+d ≥ β τn(2p+r) (log n) 2p+d = 4 2 C 2 λ min k∈R C 2 E Z,k (log n) a n ≥ 4 2 C 2 λ C 2 E Z,j (log n) a n . (30)
Therefore, on E Z , the component j remains active until the iteration τ n .

3. Let us now prove that on E Z , each relevant component j deactivates at last at iteration T n . In particular, by definition of H hp , ĥ belongs to H hp on E Z . Assume E Z happens. We fix j ∈ R. It suffices to prove that if j is still active at iteration T n , then on E Z happens, it deactivates at the end of this iteration. We assume j is still active and we denote h the state of the bandwidth at iteration T n .

By the first point, for any k / ∈ R, h k = h 0 . Given the second point, each relevant component k was still active at the beginning of the iteration

τ n + 1, ie: for any k ∈ R, h k ≤ β τn +1 h 0 ≤ β τn h 0 .
Moreover, since j is still active, h j = β Tn h 0 . Let us prove that:

1 E Z |Z hj | < λ hj . Using the decomposi- tion (29)
:

1 E Z |Z hj | ≤ E Zhj + 1 E Z Zhj -E Zhj + 1 E Z |∆ Z,hj |
Using the points 1. and 2. of Lemma 6 and 1. Lemma 7, since

E Z ⊂ B Z,hj ∩ B | Z|,hj ∩ A n : 1 E Z |Z hj | ≤ 2C E Z,j h p-1 j + 1 2 λ hj + C M∆Z (log n) a 2 λ hj ≤ λ hj   2C E Z,j n 1 2 h p j d k=1 h 1 2 k C λ (log n) a 2 + 1 2 + C M∆Z (log n) a 2   .
Given the specific form of h :

2C E Z,j n 1 2 h p j d k=1 h 1 2 k C λ (log n) a 2 ≤ 2C E Z,j n 1 2 h 2p+d 2 0 β (2p+1) 2 (Tn-τn) β (2p+r)τn 2 C λ (log n) a 2 = 4 3 C 2 E Z,j β (2p+1)(Tn-τn) min k∈R C 2 E Z,k , by definition of τ n ≤ 1 3
, by definition of T n .

Moreover, for n large enough:

C M∆Z (log n) a 2 < 1 6 .
Therefore:

1 E Z |Z hj | < λ hj .
In other words, when E Z happens, any active component at iteration T n deactivates.

So we have proved that on E Z , ĥ ∈ H hp .

It remains to show that E Z holds with high probability.

P (E c Z ) ≤ P A c n + d k=1 P B c Z,h (0) k + P B c | Z|,h (0) k + j∈R   h∈H hp P B c Z,hj + P B c | Z|,hj + τn t=1 P B c Z,h (t) j + P B c | Z|,h (t) j  
By choice of fX :

P A c n ≤ C X e -(log n) 5 4 .
We want to apply 2. and 3. of Lemma 6 for any h ∈ H hp and any h (t) with t = 1 : τ n . These bandwidths satisfy:

d k=1 h (t) k ≥ d k=1 h k ≥ h d 0 β r Tn ≥ C r 2p+1 T C r 2p+r τ (log n) ra-2p(d-r) 2p+r n -r 2p+r ,
which ensures that for n large enough, Cond Z (h ( t)) and Cond Z (h (t) ) hold for any h ∈ H hp and any t = 0 : τ n . Note in particular that H hp ⊂ {h (t),t=0: Tn ). Therefore, for any component k = 1 : d,

P B c Z,h (0) k ≤ 2e -γ Z,n
and for any h ∈ H hp and any t = 0 : T n ,

P B c | Z|,h (t) j ≤ 2 exp -C γ| Z| n d k=1 h k ≤ 2 exp -C γ| Z| C r 2p+1 T C r 2p+r τ (log n) ra-2p(d-r) 2p+r n 2p 2p+r
≤ 2e -γ Z,n , for n large enough.

To conclude, note that

|H hp | = ( T n -τ n ) r ≤ (T n -τ n + 2) r = log(C -1 T ) (2p+1) log 1 β + 2
r is finite, so for any q > 0:

P (E c Z ) ≤ C X e -(log n) 5 4 + 2(d + r|H hp | + rτ n )2e -γ Z,n = o(n -q ), by definition γ Z,n := δ f ∞, Un(w) (log n) a .

Proof of Lemma 4

1. We control the bias

E f K X -f X ∞, Un(x)
. We write for any u ∈ U n (x):

E f K X (u) -f X (u) = 1 h d 1 X u ∈R d 1   d 1 j=1 K X u j -u j h X   f X (u )du -f X (u) R d 1   d 1 j=1 K X (z j )   dz
The kernel K X is of order p X and f X is assumed of class C p on U n (x), with in particular p -1 ≤ p X -1, then we can apply Lemma 9 with the settings u = u, d = d 1 , f 0 = f X , p = p -1, K = K X and for j = 1 : d , h k = h X . We obtain:

E f K X (u) -f X (u) = d 1 k=1 (I k + II k ). (31)
with

I k := z∈R d 1 d 1 k =1 K X (z k ) ρ k dz, ρ k := ρ k (z, h X , u) = (-h X z k ) p -1 0≤t p -1 ≤•••≤t 1 ≤1 ∂ p -1 k f X (z k-1 -t p -1 h X z k e k ) -∂ p -1 k f X (z k-1 ) dt 1:(p -1)
,

II k := (-h X ) p -1 t∈R t p -1 (p -1)! K X (t)dt z -k ∈R d 1 -1 ∂ p -1 k f X (z k-1 )   k =k K X (z k )   dz -k .
Let us control ρ k . First we write:

∂ p -1 k f X (z k-1 -t p -1 h X z k e k ) -∂ p -1 k f X (z k-1 ) = -h X z k 1 t p =0 ∂ p k f X (z k-1 -t p h X z k e k )dt p .
Therefore:

ρ k = (-h X z k ) p 0≤t p ≤•••≤t 1 ≤1 ∂ p k f X (z k-1 -t p h X z k e k )dt 1:p .
Hence:

|ρ k | ≤ |h X z k | p 0≤t p ≤•••≤t 1 ≤1 ∂ p k f X (z k-1 -t p h X z k e k ) dt 1:p = |z k | p p ∂ p k f X ∞, U n (x) h p X .
Then:

|I k | ≤ z∈R d 1 d k =1 K X (z k ) |ρ k | dz ≤ ∂ p k f X ∞, U n (x) h p X z∈R d 1 |z k | p p d k =1 K X (z k ) dz = K X d 1 -1 1 • p K X (•) 1 p ∂ p k f X ∞, U n (x) h p X (32)
Besides, K X is of order p X and p -1 < p X and so:

II k := (-h X ) p -1 (p -1)! t∈R t p -1 K X (t)dt z -k ∈R d 1 -1 ∂ p -1 k f X (z k-1 )   k =k K X (z k )   dz -k = 0.
Therefore the terms II k vanish in the equation ( 31), and with the upper bound of I k (32), we obtain:

E f K X -f X ∞, Un(x) = sup u∈ Un(x) E f K X (u) -f X (u) ≤ sup u∈ Un(x) d 1 k=1 |I k | ≤ K X d 1 -1 1 • p K X (•) 1 p ! h p X d 1 k=1 ∂ p k f X ∞, U n (x) = C bias X h p X ,
with C bias X :=

K X d 1 -1 1 • p K X (•) 1 p ! d 1 max k=1:d 1 ∂ p k f X ∞, U n (x) .
2. We apply Bernstein's inequality (see Lemma 10). We define for any u ∈ U n (x) and any i = 1 : n X :

f K Xi (u) := 1 h d 1 X d 1 j=1 K X u j -X ij h X .
Then we control f K X1 a.s.: for any u ∈ U n (x),

f K X1 (u) ≤ M h X := K X d 1 ∞ h -d 1 X .
and its variance:

Var f K X1 (u) ≤ E ( f K X1 ) 2 = h -2d 1 X u ∈R d 1   d 1 j=1 K X u j -u j h X   2 f X (u )du = h -d 1 X z∈R d 1   d 1 j=1 K X (z j )   2 f X (u -h X z)du ≤ v h X with v h X := C vX h -d 1 X and C vX := K X 2d 1 2 f X ∞, U n (x)
. Then we apply Lemma 10: for any λ > 0,

P f K X (u) -E f K X (u) > λ ≤ 2 exp -min n X λ 2 4v h X , 3n X λ 4M h X .
We set λ = λ X :=

4v h X n X (log n) 3 2 such that (log n) 3 2 = n X λ 2 4v h X
. Then we compare the rates:

n X λ 2 4v h X ≤ 3n X λ 4M h X ⇐⇒ λ 2 ≤ 3 2 C 2 vX K X 2d 1 ∞ ⇐⇒ h d 1 X ≥ 4 K X 2d 1 ∞ 9C vX (log n) 3 2 n X ,
⇐⇒ Cond X (h X ).

Proof of Lemma 5

1. We recall the notation • for the multiplication terms by terms of two vectors. Then:

|E fh1 (w) | ≤ E | fh1 (w)| = u∈R d d k=1 K(h -1 k (w k -u k )) h k f (u)du = z∈R d d k=1 K(z k ) f (w -h • z)dz ≤ f ∞, Un(w) K d 1 =: CĒ
Now let us give an upper bound on the bias of fh (w):

B h = E fh1 (w) -f (w) = u∈R d d k=1 K(h -1 k (w k -u k )) h k f (u)du -f (w) R d d k =1 K(z k )dz, since R K(t)dt = 1.
Then we apply the Lemma 9 with the settings d = d, u = w, h = h, f 0 = f , p = p and K = K. We obtain:

B h = d k=1 (I k + II k ),
where

I k := z∈R d d k =1 K(z k ) ρ k dz, ρ k := (-h k z k ) p 0≤tp≤•••≤t 1 ≤1 ∂ p k f (z k-1 -t p h k z k e k ) -∂ p k f (z k-1 ) dt 1:p , II k := (-h k ) p t∈R t p p! K(t)dt z -k ∈R d-1 ∂ p k f (z k-1 )   k =k K(z k )   dz -k .
Notice that for k / ∈ R, ∂ p k f (u) = 0 for any u ∈ U n (x), thus I k and II k vanish. Therefore:

B h = k∈R (I k + II k ),
Now let us give an equivalent of the bias. First, using Assumption 3, for any k ∈ R, we can define the modulus of continuity of ∂ p k f on U n (w) by:

Ω nk := sup z,z ∈ Un(w) ∂ p k f (z ) -∂ p k f (z)
Then we decompose II k as follows:

II k = (-h k ) p t∈R t p K(t)dt p! ∂ p k f (w) + R k , with R k := (-h k ) p t∈R t p K(t)dt p! z -k ∈R d-1 (∂ p k f (z k-1 ) -∂ p k f (w)) k =k K(z k ) dz -k such that: |R k | ≤ h p k t∈R t p p! K(t)dt Ω nk K d-1 1 (33) since ∂ p k f (z k-1 ) -∂ p k f (w) ≤ Ω nk . It remains to bound I k .
From the definition of ρ k in (50), we write:

|ρ k | ≤ |h k z k | p 0≤tp≤•••≤t 1 ≤1 ∂ p j f (z k-1 -t p h k z k e k ) -∂ p j f (z k-1 ) dt 1:p ≤ |h k z k | p Ω nk p! .
Therefore:

|I k | = z∈R d d k =1 K(z k ) ρ k dz ≤ h p k p! Ω nk z∈R d z p k d k =1 K(z k ) dz ≤ K d-1 1 t∈R t p p! K(t) dt × h p k Ω nk (34) Since U n (w) -→ n→∞ {w}, by continuity of ∂ p k f : Ω nk -→ n→∞ 0.
Therefore for n large enough, combining (33) and (34):

|I k | + |R k | ≤ | t∈R t p K(t)dt| p! max k∈R |∂ p k f (w)| × h p k
Therefore, since:

B h = k∈R (II k + I k ) = k∈R (-h k ) p t∈R t p K(t)dt ∂ p k f (w) p! + R k + I k , we obtain: B h ≤ C bias k∈R h p k , with C bias := 2| t∈R t p K(t)dt| p! max k∈R |∂ p k f (w)|.
2. We want to apply Bernstein's inequality (cf Lemma 10) to fh (w). We first obtain an almost sure upper bound:

| fh1 (w)| = 1 f X (X 1 ) d k=1 K w k -W 1k h k h k ≤ M h , (35)
where M h :=

CM d k=1 h k with CM := K d ∞ δ . Then we control the variance: Var fh1 (w) = Var   1 f X (X 1 ) d k=1 K w k -W 1k h k h k   ≤ E      1 f X (X 1 ) d k=1 K w k -W 1k h k h k   2    = u∈R d d k=1 1 h 2 k K w k -u k h k 2 f (u) f X (u 1:d 1 ) du ≤ 1 δ d k=1 h k z∈R d d k=1 K(z k ) 2 f (w -Hz)dz ≤ v h , (36)
where v h :=

C 2 σ 4 d k=1 h k
.Therefore we obtain from Bernstein's inequality (cf Lemma 10):

P B c h ≤ 2 exp -min nσ 2 h 4v h , 3nσ h 4M h
We compare the rates:

nσ 2 h 4v h ≤ 3nσ h 4M h ⇐⇒ C σ (log n) a n d k=1 h k = σ h ≤ 3v h M h = 3C 2 σ 4CM ⇐⇒ d k=1 h k ≥ 4C 2 M 9C 2 σ (log n) a n ⇐⇒ Cond(h).
Therefore, if Cond(h) is satisfied:

P B c h ≤ 2e - nσ 2 h 4v h = 2e -(log n) a .
3. We now apply Bernstein's inequality (cf Lemma 10) to 1 n n i=1

| fhi (w)|. From the upper bounds ( 35) and

(36), we obtain:

P B c | f |h ≤ 2 exp -min nC 2 Ē 4v h , 3nCĒ 4M h .
We calculate the rates: by definition of v h and M h ,

nC 2 Ē 4v h = C 2 Ē C 2 σ n d k=1 h k 3nCĒ 4M h = 3CĒ 4CM n d k=1 h k Hence: P B c | f |h ≤ 2e -C γ|f | n d k=1 h k , with C γ|f | := min C 2 Ē C 2 σ ; 3CĒ 4CM .

Proof of Lemma 6

1. First, we write Zhij more explicitly: for any bandwidth h, any observation i = 1 : n and any direction j,

Zhij = ∂ ∂h j   K( w j -W ij h j ) h j   k =j K( w k -W ik h k ) f X (X i ) k =j h k = -K( w j -W ij h j ) + w j -W ij h j K ( w j -W ij h j ) k =j K( w k -W ik h k ) f X (X i )h j d k=1 h k = -J( w j -W ij h j ) k =j K( w k -W ik h k ) f X (X i )h j d k=1 h k
where we recall J : R → R is the function t → tK (t) + K(t).

Note then that the support of J is included in the support of K, and by integration by part, we obtain for any l ∈ N: R

t l J(t)dt = R t l (tK(t)) dt = -l R t l K(t)dt (37)
In particular, since K is of order p, for l = 0 : p -1, R t l J(t)dt = 0 and R t p J(t)dt = 0. We recall the notation • for the multiplication terms by terms of two vectors. Using Assumption 2, if j / ∈ R, f (w -h • z) -f (z -j ) = 0 for any z ∈ R d . Thus we obtain:

E Zh1j = - 1 h j d k=1 h k u∈R d J( w j -u j h j )   k =j K( w k -u k h k )   f (u)du = - 1 h j z j ∈R J(z j )dz j z -j ∈R d-1   k =j K(z k )   f (w -h • z)dz -j = 0
Therefore E Zh1j = 0 for j / ∈ R. Now, we deal with the case j ∈ R. Let us fix j ∈ R. We denote z-j :

= w -(Hz) -j = w - k =j h k z k e k (with {e k } d k=1 the canonic basis of R d ).
Then we write:

E Zh1j = -1 h j d k=1 h k u -j ∈R d-1   k =j K( w k -u k h k )   u j ∈R J( w j -u j h j )du j -f (z -j ) R J(z j )dz j du -j .
Then for fixed {z k } k =j , denoting f j : z j → f (w -h • z), we apply Lemma 9 with the settings d

= 1, u = z-j , h = h j , f 0 = f j , p = p, K = J, then E Zh1j = -1 h j d k=1 h k u -j ∈R d-1   k =j K( w k -u k h k )   [I 1 + II 1 ] du -j = I j + ĨI j , (38)
where

Ĩj := (-h j ) -1 z∈R d   k =j K(z k )   J(z j )ρ j dz, (39)
with ρj := (-h j z

j ) p 0≤tp≤•••≤t 1 ≤1 ∂ p j f (z -j -t p h j z j e j ) -∂ p j f (z -j ) dt 1:p , (40)
and II j := (-h j ) p-1 t∈R

t p p! J(t)dt z -j ∈R d-1 ∂ p j f (z j-1 )   k =j K(z k )   dz -j .
Now let us determine an equivalent of E Zhj . For this purpose, let us introduce the modulus of continuity of ∂ p j f on U n (w) (which is well defined by Assumption 3):

Ω nj := sup z,z ∈ Un(w) ∂ p j f (z ) -∂ p j f (z)
Then we write:

II j = (-h j ) p-1 ∂ p j f (w) t∈R t p p! J(t)dt + R j , (41)
with

R j := (-h j ) p-1 t∈R t p p! J(t)dt z -j ∈R d-1 ∂ p j f (z -j ) -∂ p j f (w)   k =j K(z k )   dz -j .
In particular:

| R j | ≤ h p-1 j t∈R t p p! J(t)dt z -k ∈R d-1 Ω nj k =j |K(z k )|dz -j = h p-1 j Ω nj t∈R t p p! J(t)dt K d-1 1 (42)
Now let us bound Ĩj defined in (39). First, we bound ρj , defined in (40):

|ρ j | = (h j |z j |) p 0≤tp≤•••≤t 1 ≤1 ∂ p j f (z -j -t p h j z j e j ) -∂ p j f (z -j ) dt 1:p ≤ h p j |z j | p Ω nj p! ,
which lead to:

Ĩj = h -1 j z∈R d   k =j K(z k )   J(z j )ρ j dz ≤ h p-1 j Ω nj K d-1 1 z j ∈R |z j | p p! |J(z j )|dz j . (43)
Therefore using ( 41) then ( 42) and ( 43):

E Zh1j ≤ II j + Ĩj ≤ h p-1 j ∂ p j f (w) t∈R t p p! J(t)dt + R j + Ĩj ≤ C E Z,j h p-1 j + h p-1 j Ω nj t∈R t p p! J(t)dt K d-1 1 + K d-1 1 R |t| p p! |J(t)|dt
with C E Z,j := ∂ p j f (w) t∈R t p p! J(t)dt . Finally, notice that by continuity of ∂ p j f (Assumption 3), since U n (w) -→ n→∞ {w}:

Ω nj -→ n→∞ 0.
Thus for n large enough:

Ω nj t∈R t p p! J(t)dt K d-1 1 + K d-1 1 z j ∈R |z j | p p! |J(z j )|dz j ≤ 1 2 C E Z,j ,
which lead to the result (12) of Theorem 2:

1 2 C E Z,j h p-1 j ≤ E Zhj ≤ 3 2 C E Z,j h p-1 j .
To obtain the result (13) of Theorem 2, just note that:

E | Zh1j | = 1 h j d k=1 h k u∈R d J( w j -u j h j )   k =j K( w k -u k h k )   f (u)du = h -1 j z∈R d J(z j )   k =j K(z k )   f (w -Hz)dz ≤ C E| Z| h -1 j ,
with C E| Z| := f ∞, Un(w) J 1 K d-1 1 .

2. We first bound Zhij a.s. and its variance.

| Zhij | = J( w j -W ij h j ) k =j K( w k -W ik h k ) f X (X i )h j d k=1 h k ≤ J ∞ K d-1 ∞ δh j d k=1 h k = C M Z h j d k=1 h k =: M Z,hj (44)
For the variance:

Var Zhij ≤ E Z2 hij = R d J( w j -u j h j ) 2   k =j K( w k -u k h k ) 2   f XY (u) f X (u 1:d 1 ) 2 h 2 j d k=1 h 2 k du = 1 h 2 j d k=1 h k R d J(z j ) 2   k =j K(z k ) 2   f (w -Hz) f X (x -(Hz) 1:d 1 ) dz ≤ f ∞, Un(w) J 2 2 K 2(d-1) 2 δh 2 j d k=1 h k = C v Z h 2 j d k=1 h k =: v Z,hj . (45)
We apply Bernstein's inequality (cf Lemma 10) to Zhj : .

P B c Z,

Proof of Lemma 7

1. We decompose ∆ Z,hj as follows:

∆ Z,hj := Z hj -Zhj = 1 n n i=1 f X (X i ) -fX (X i ) fX (X i ) Zhij .
Using Zhij = 0 when X i / ∈ U h (x):

|∆ Z,hj | ≤ f X -fX fX ∞, Un(x) 1 n n i=1 Zhij . ( 46
)
First we deal with f X -fX fX ∞, Un(x)

. By definition of A n :

1 An f X -fX fX ∞, Un(x) ≤ M X (log n) d n 1/2 , ( 47
)
Now let us give an upper bound of 1 n n i=1

Zhij . Using Lemma 6,

1 B | Z|,hj 1 n n i=1 Zhij ≤ 1 B | Z|,hj 1 n n i=1 Zhij -E | Zh1j | + E | Zh1j | ≤ 2C E| Z| h -1 j
To conclude, combining this last result with ( 47) and ( 46):

1 B | Z|,hj ∩ An |∆ Z,hj | ≤ 2C E| Z| M X h -1 j (log n) d n 1/2 ≤ 2C E| Z| M X C λ (log n) a 2 λ hj = C M∆Z (log n) a 2 λ hj , since d k=1 h k ≤ h d 0 = 1 (log n) d .
2. We decompose ∆ h as follows:

∆ h := fh (w) -fh (w) = 1 n n i=1 f X (X i ) -fX (X i ) fX (X i ) fhi (w).
Using fhi (w) = 0 when X i / ∈ U h (x):

|∆ h | ≤ f X -fX fX ∞, Un(x)
1 n n i=1 fhi (w) .

We have proved in (47): 1 An Therefore:

1 An∩B | f |h |∆ h | ≤ 2CĒM X (log n) d n 1/2 ≤ 2CĒM X C σ (log n) -a 2 σ h , since d k=1 h k ≤ h d 0 = (log n) -d .

Proof of Lemma 9

We first denote

B := R d   d j=1 h -1 j K u j -u j h j   f 0 (u )du -f 0 (u) R d   d j=1 K (z j )   dz.
Then we obtain by integration by parts:

B := z∈R d   d j=1 K(z j )   (f 0 (u -h • z) -f 0 (u))dz (48)
For any z ∈ R d , we denote z 0 := w and for k = 1 : d , z k := u -k j=1 h j z j e j (where {e j } d j=1 is the canonical basis of R d ). Then, we write:

f 0 (u -h • z) -f 0 (u) = d k=1 f 0 (z k ) -f 0 (z k-1 ) (49)
Then we apply Taylor's theorem (cf Lemma 11) to the functions g k : t ∈ [0, 1] → f 0 (z k-1 -th k z k e k ), k ∈ (1 : d ):

f 0 (z k ) -f 0 (z k-1 ) = g k (1) -g k (0) = p l=1 (-z k h k ) l l! ∂ l k f 0 (z k-1 ) + ρ k ,
where we denote for short:

ρ k := ρ k (z, h, u) = (-h k z k ) p 0≤tp≤•••≤t 1 ≤1 ∂ p k f 0 (z k-1 -t p h k z k e k ) -∂ p k f 0 (z k-1 ) dt 1:p . ( 50
)
We introduce the notation

Figure 2 :

 2 Figure 2: CDRodeo estimator (red or blue dashed lines) VS the true density (black solid line) in function of each component, the others component being fixed following (x, y) = (0, 1, 0, 0, 1).

Lemma 9 .

 9 Let u ∈ R d and a bandwidth h ∈ R * + d . For j = 1 : d , let K : R → R be a continuous function with compact support and with at least p -1 zero moments, ie: for l = 1 :

  ) -E | fhi (w)| + E | fh1 (w)| ≤ 2CĒ.

  dz and for any z ∈ R d , we denote z -k ∈ R d -1 the vector z without its k th variable, then (48) becomes:z k )dz k dz -k  Since K has at least p -1 zero moments, the terms with l ≤ p -1 vanish, leading to:(t)dt z -k ∈R d -1 ∂ p k f 0 (z k-1 ) j =k K(z j ) dz -k .

  Let us calculate the rate: by definition of C v Z hj and M Z,hj ,

	Let us compare the rates:
				n( 4v Z,hj λ hj 2 ) 2	≤	λ hj 2 4M Z,hj 3n
		⇐⇒ C λ	(log n) a nh 2 j d k=1 h k	= λ hj ≤	6v Z,hj M Z,hj	=	6C v C M Z h j Z
		⇐⇒	d k=1	h k ≥	C 2 M 6 2 C 2 Z C 2 λ v Z (log n) a n
		⇐⇒ Cond Z (h).
	So, if Cond Z (h) is satisfied:
		P B c Z,hj ≤ 2e	-n(λ hj /2) 2 4v Z,hj	= 2e	-δ f ∞, Un(w)	(log n) a	= 2e -γ Z,n
	3. We apply Bernstein's inequality (cf Lemma 10) to 1 n	i=1 n	| Zhij | using the upper bounds (44) and (45):
		P B c | Z|,h ≤ 2 exp -min	n(C E| Z| h -1 j ) 2 4v Z,hj	,	3nC E| Z| h -1 j 4M Z,hj
					n(C E| Z| h -1 j ) 2 4v Z,hj	=	C 2 E| Z| 4C v Z n	d k=1	h k
					3nC E| Z| h -1 j 4M Z,hj	=	3C E| Z| 4C M Z n	d k=1	h k
	Hence:
					P B c | f |h ≤ 2e -C γ| Z| n d k=1 h k ,
	with C γ| Z| := min	C 2 E| Z| 4C v Z ;	3C E| Z| 4C M Z
							λ hj 2 ) 2 4v Z,hj	,	λ hj 2 4M Z,hj 3n

hj ≤ 2 exp -min n(

Acknowledgement.

The author is extremely grateful to Claire Lacour and Vincent Rivoirard for suggesting me to study this problem, for the stimulating discussions, the helpful advices and the careful proofreading.