<figure><img src="image_1.png" /></figure>

HAL
open science

# New estimation of Sobol' indices using kernels

Sébastien Da Veiga, Fabrice Gamboa, Agnès Lagnoux, Thierry Klein, Clémentine Prieur

► To cite this version:

Sébastien Da Veiga, Fabrice Gamboa, Agnès Lagnoux, Thierry Klein, Clémentine Prieur. New estimation of Sobol' indices using kernels. 2023. hal-04052837

HAL Id: hal-04052837

https://hal.science/hal-04052837

Preprint submitted on 30 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Sensitivity analysis from a single input/output sample

Sébastien Da Veiga, Fabrice Gamboa, Thierry Klein, Agnès Lagnoux,
Clémentine Prieur

► To cite this version:

Sébastien Da Veiga, Fabrice Gamboa, Thierry Klein, Agnès Lagnoux, Clémentine Prieur. Sensitivity analysis from a single input/output sample. 2024. hal-04052837v3

HAL Id: hal-04052837

https://hal.science/hal-04052837v3

Preprint submitted on 6 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# New estimation of Sobol' indices using kernels

Sébastien da Veiga¹, Fabrice Gamboa², Thierry Klein³, Agnès Lagnoux⁴, and Clémentine Prieur⁵

1

²Institut de Mathématiques de Toulouse and ANITI; UMR5219.

Université de Toulouse; CNRS. UT3, F-31062 Toulouse, France.

³Institut de Mathématiques de Toulouse; UMR5219. Université de Toulouse; ENAC - Ecole Nationale de l'Aviation Civile, Université de Toulouse, France

⁴Institut de Mathématiques de Toulouse; UMR5219. Université de Toulouse; CNRS. UT2J, F-31058 Toulouse, France.

⁵Laboratoire Jean Kuntzmann; CNRS, UMR 5224. Université Grenoble Alpes, F-38400 Saint Martin d'Hères, France

March 30, 2023

## Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Setting</td><td>4</td></tr><tr><td>2.1</td><td>Model and purpose</td><td>4</td></tr><tr><td>2.2</td><td>Notation and assumptions</td><td>5</td></tr><tr><td>3</td><td>Central limit theorem</td><td>6</td></tr><tr><td>4</td><td>Efficiency</td><td>9</td></tr><tr><td>5</td><td>Bandwidth selection</td><td>11</td></tr><tr><td>6</td><td>Extension to an unknown density f<sub>x</sub></td><td>12</td></tr><tr><td>7</td><td>Numerical applications</td><td>16</td></tr><tr><td>A</td><td>Proof of the results</td><td>16</td></tr><tr><td>A.1</td><td>Proof of Theorem 3.1</td><td>16</td></tr><tr><td>A.2</td><td>Proof of the remaining results</td><td>22</td></tr><tr><td>A.3</td><td>Construction of a kernel of order k in dimension d</td><td>26</td></tr></table>
# Sensitivity analysis from a single input/output sample

Sébastien Da Veiga¹, Fabrice Gamboa², Thierry Klein³, Agnès Lagnoux⁴, and Clémentine Prieur⁵

¹Université de Rennes, Ensai, CNRS, CREST - UMR 9194, F-35000 Rennes, France

²Institut de Mathématiques de Toulouse and ANITI; UMR5219. Université de Toulouse; CNRS. UT3, F-31062 Toulouse, France and ANITI Toulouse France

³Institut de Mathématiques de Toulouse; UMR5219. Université de Toulouse; ENAC - Ecole Nationale de l'Aviation Civile, Université de Toulouse, France

⁴Institut de Mathématiques de Toulouse and ANITI; UMR5219. Université de Toulouse; CNRS. UT2J, F-31058 Toulouse, France and ANITI Toulouse France

⁵Université de Grenoble Alpes, CNRS, Inria, Grenoble INP, LJK, Grenoble, France

November 5, 2025

## Abstract

The main objective of this paper is to estimate optimally Sobol' indices at any order when a unique input/output i.i.d. sample is available. Our approach stands on three main ingredients: semi-parametric estimation theory, high-order kernel estimation (inspired by the paper [15]), and mirror-type transformations as introduced in [3, 41]. We propose two different estimators. We prove that these estimators are asymptotically normal and efficient. Furthermore, we illustrate their numerical properties on standard examples.

**Keywords:** nonparametric kernel estimation; asymptotic properties; global sensitivity analysis; Sobol' indices; efficient influence function.

**AMS subject classification:** 62G05, 62G08, 62G20.

## 1 Introduction

In this paper, we consider a random variable $Y$ depending on a random vector $V = (V_1, \dots, V_p)$ through the relation $Y = G(V_1, \dots, V_p)$ where the function $G$ is defined on a compact subset of $\mathbb{R}^p$, $p \ge 2$ and is real-valued. Then we focus on the estimation of the ratio

$$S^X = \frac{\mathrm{Var}(\mathbb{E}[Y|X])}{\mathrm{Var}(Y)} = \frac{\mathbb{E}[\mathbb{E}[Y|X]^2] - \mathbb{E}[Y]^2}{\mathrm{Var}(Y)}$$

where $X$ is a $d$-dimensional vector formed by $d$ coordinates of $V$ ($d$ being a non zero integer smaller than $p$). The interest in this ratio dates back to Karl Pearson as part of analysis of variance (see, e.g., [29]), where it coincides to the square of the correlation ratio (measure of the relationship between the statistical dispersion within individual categories
# 1 Introduction

The use of complex computer models for the simulation and analysis of natural systems from physics, engineering and other fields is by now routine. These models usually depend on many input variables. Thus, it is crucial to understand which input parameter or which set of input parameters have an influence on the output. This is the aim of sensitivity analysis which has become an essential tool for system modeling and policy support (see, e.g., [39]). Global sensitivity analysis methods considers the input vector as random and propose a measure of the influence of each subset of its components on the output of interest. We refer to the seminal book by Saltelli [40] for an overview on global sensitivity analysis or to [10] for a synthesis of recent trends in the field. Among the different measures of global sensitivity analysis, variance-based measures are probably the most commonly used. The definition of the so-called Sobol' indices, introduced in [34] and later revisited in the framework of sensitivity analysis in [41] (see also [42]), is based on the Hoeffding decomposition of the variance [26]. More precisely, for $Y = f(X_1, \dots, X_p)$ where the inputs $X_i$ are assumed to be mutually independent, the Sobol' index with respect to $X$ is defined by

$$S^X = \frac{\mathrm{Var}(\mathbb{E}[Y|X])}{\mathrm{Var}(Y)} = \frac{\mathbb{E}[\mathbb{E}[Y|X]^2] - \mathbb{E}[Y]^2}{\mathrm{Var}(Y)}$$

where $X$ is a subvector of $(X_1, \dots, X_p)$, for instance $X = (X_1, \dots, X_d)$ with $1 \le d \le p$ or $X = X_i$ (first-order index) with $1 \le i \le p$. For $i = 1, \dots, p$, it is also possible to measure the total influence of $X_i$, by itself or in interaction with other inputs, with the total Sobol' index defined by $1 - S^{X-i}$ with $X_{-i} = (X_1, \dots, X_{i-1}, X_{i+1}, \dots, X_p)$ for $i = 1, \dots, p$.

In general, computing explicitly the theoretical value of $S^X$ is hopeless; therefore one of the main tasks is to provide estimators of $S^X$ that have some nice convergence properties (consistency, explicit rate of convergence, asymptotic efficiency...). Many different estimation procedures were developed in the last decades to estimate $S^X$ or $\mathbb{E}[\mathbb{E}[Y|X]^2]$ only (since one easily estimates $\mathbb{E}[Y]$ and $\mathrm{Var}(Y)$ by the empirical mean and the empirical variance respectively). Some are based on Monte Carlo, quasi Monte Carlo or nested Monte Carlo designs of experiments (see, e.g., [28, 33, 23]). Others are based on spectral approaches. Among this last class of procedures, let us cite without exhaustivity Fourier Amplitude Sensitivity Test (FAST) [8], Random Balance Design (RBD) [46], Effective Algorithm for computing global Sensitivity Indices (EASI) [36] and polynomial chaos expansions [45] (see also [10, Chapter 4] for a more complete review). In the sequel, we focus on two classes of methods that have emerged since they provide nice theoretical convergence guarantees.

The first class is based on the so-called Pick Freeze (PF) design of experiments, the theoretical properties of which (consistency, central limit theorem, concentration inequalities and Berry-Esseen bounds) have been studied in [18, 27]. Besides having very nice convergence properties, this method is very general since the only assumption needed to prove a central limit theorem is that $\mathbb{E}[Y^4] < \infty$. However, it requires the evaluation of the model on a PF design thus it cannot be computed by recycling former evaluations of the model. Also, the required sample size increases linearly with the number of indices that one wants to estimate. For instance, to estimate at rate $\sqrt{n}$ a single Sobol' index, one needs a design of experiment of size $2n$. Additionally, this size grows to $(p+1)n$ to estimate the full set of first-order Sobol' indices (still at the rate $\sqrt{n}$). By evaluating PF
and the dispersion across a whole population or a sample). It was later revisited by Sobol'
in the framework of sensitivity analysis [45, 46], under the name of closed Sobol' index.
Interestingly, it also appeared in the statistics literature under different names: measure of
explanatory power of covariates or nonparametric coefficient of determination by [15], or
the residual variance (focusing only on the numerator) which provides a lower bound for
the performance of any regression function estimator in multivariate problems, as noticed
in [13]. Finally, recent studies have emphasized the connection between machine learning
feature importance measures and this ratio (see [5] and references therein). In particular,
it is the building block of the Shapley values which are now highly popular importance
measures in machine learning for explainability [38]. It can also be used for recursive
feature elimination [2].

Here, without loss of generality we will adopt the notation and terminology from the field of sensitivity analysis, where this ratio is called the closed Sobol' index, hence the notation $S^X$ above. This research area emerged in the computer experiments community, which focuses on complex computer models to simulate and analyze natural systems in physics, engineering and other fields. These models usually depend on many input variables, and it is thus crucial to understand which input parameter or which set of input parameters have an influence on the output. This is exactly the purpose of sensitivity analysis, which has become an essential tool for systems modelling and policy support (see, e.g., [42]). Global sensitivity analysis methods consider the input vector as random and propose a measure of the influence, in terms of output fluctuations, of each subset of its components. We refer to the seminal book [43] for an overview on global sensitivity analysis or to [10] for a synthesis of recent trends in this field. Among the different measures of global sensitivity analysis, variance-based measures are probably the most commonly used. More precisely, for the output $Y$ of a computer code $G$, one of the most common measure of the sensitivity of $Y$ with respect to a vector $X$ is precisely the closed Sobol' index $S^X$ as defined at the very beginning of this section.

In recent years a myriad of different estimators have been proposed, see [10, Chapter 4] for a full review. A practical statistical method widely used to construct estimators is based on spectral analysis of the input/output functional relationship and Parseval formulae. We refer to [48] for a basic description of the method and to [10, Chapter 4] for more recent references. It should be noted that the asymptotic properties of these methods have been little studied, as they are based on the theory of non-linear (quadratic) functional estimation. Away from the spectral methods, two families of estimation methods are of particular interest because of the possibility of studying the asymptotic properties of the estimators produced.

The first family of methods is based on the so-called Pick-Freeze (PF) design of ex-
periments. The basic idea is to evaluate $G$ repeatedly for input values, fixing those of
the vector $X$, and then to calculate a Monte Carlo estimate of $S^X$ from this particular
sampling. The main advantage is that only minimal assumptions are required to derive
consistency and central limit theorems under the assumption of independent covariates.
In particular, assumptions of integrability but not regularity on $G$ are necessary (see for
example [17, 28]). As an illustration, to estimate at rate $\sqrt{n}$ a single Sobol' index, one
needs a design of experiments of size $2n$. This implies that estimating all the $p$ first-order
indices corresponding to $d=1$ involves a sample of size $(p+1)n$; a $(p+2)n$ sample is
required to estimate the full set of first-order, second-order, and total indices. This cost
estimates on replicated designs, this cost has been reduced to $2n$ to estimate the full set of first-order and second-order indices [47, 20]. However the cost for estimating total Sobol' indices with this approach remains linear in the input space dimension.

The second class of methods is based on nearest neighbors and has been studied by several authors. See e.g. [16, 29, 30, 14, 24, 15]. For instance, in [15], the authors built an estimator of $\mathbb{E}[\mathbb{E}[Y|X]^2]$ based on two independent $n$-samples. The first one allows to estimate the regression function $\mathbb{E}[Y|X = x]$ using the first nearest neighbor of $x$ among the first sample and the second one is used as a plug-in estimator. They proved a central limit theorem for their estimator. Nevertheless the recentering factor is not the unknown quantity of interest $\mathbb{E}[\mathbb{E}[Y|X]^2]$ but the expected value of their estimator. In addition, the bias is negligible if and only if the dimension $d \le 3$. In [6], the author introduced a new coefficient of correlation based on ranks for $d = 1$. It was proven that this coefficient converges to a limit that appears to be the Cramér-von Mises index with respect to $X_i$ (for $i = 1, \dots, p$) introduced in [19] as the size $n$ of the sample goes to infinity. Note that the coefficient introduced in [6] admits a representation in terms of nearest neighbor (on the right) (see [17, Section3]). A central limit theorem was proved in [6] under the assumption that $Y$ and $X_i$ are independent. This theorem can be exploited for statistical tests of independence. The central limit theorem was extended in [31] (see also references therein) to the framework where $Y$ is not a measurable function of $X_i$. The definition of the Chatterjee's coefficient was then extended in [2] including the case $d \ge 2$ with nearly optimal bounds on the rate of convergence. Chatterjee's approach was then adapted for building estimators of first-order Sobol' indices in [17], with a central limit theorem in the framework $d = 1$ and $X_i$ not necessarily independent from $Y$. Broto et al. [4] also considered a nearest neighbor approach to estimate Sobol' indices of any order. The estimators they introduced are consistent but no rate of convergence was provided. Although very powerfull it is clear thanks to the bias study in [15] that it is not possible to obtain a central limit theorem for estimating Sobol' indices of order larger than $d = 3$ using this approach. In particular, if $d \ge 5$, one can not estimate total order Sobol' indices with this kind of approaches at rate $\sqrt{n}$. It is important to note that all the procedures based on nearest neighbors require additional regularity assumptions on the model to control the bias. Close to nearest neighbor approaches are kernel methods. A natural estimator is the Nadaraya-Watson plug-in estimator or its generalization with local polynomial kernel smoothing. Based on works on the estimation of conditional covariance matrix [49, 32], kernel estimates of first-order Sobol' indices have been studied in [11, 9, 37, 43, 25]. In particular, the authors in [9] proved the asymptotic efficiency of their estimator of first-order Sobol' indices for $d = 1$. The analysis of such estimators also require regularity assumptions on the model.

Hence, in the one hand, the Pick Freeze procedure allows to estimate any Sobol' index at rate $\sqrt{n}$ with a minimal assumption on the model but requires a very specific design of experiment, the size of which growing linearly with the number of indices that one wants to estimate (except the procedure based on replicated designs to estimate first-order and second-order indices, as already mentioned). In the other hand, the nearest neighbor approach requires a single $n$-sample but needs some regularity assumptions on the model and for which the $\sqrt{n}$ parametric rate of convergence is possible only for indices up to order 3.

In the present paper, we develop a new approach mentioned in [9, Remark 1]. It consists
can be reduced to $2n$ [49, 20] if the aim is only to estimate first-order and second-order indices. Recently, the authors in [5] proposed to apply a weighting factor (called density quotient) to PF estimators to handle dependent covariates in settings in which both probability distributions of $V$ and $X$ are known. However, the main drawback of the PF estimation procedure is that it requires a specific experimental design. In particular, it cannot be used in the case where we only have a classical i.i.d. sample of $n$ input/output observations.

The second class of methods relies on local averaging and tackles this limitation. Among them are kernel estimators, which have been thoroughly studied for the case $d = 1$ [11, 9, 40, 47, 25] with central limit theorems and asymptotic efficiency as soon as $G$ satisfies regularity assumptions. Closely related are nearest neighbor approaches, which have been studied by several authors (see, e.g., [14, 33, 34, 12, 23, 13]). For instance in [13], the authors propose a plug-in estimator with statistical consistency for any $d$ and a central limit theorem with rate $\sqrt{n}$ for $d \le 3$ if, again, regularity assumptions hold. In parallel, [6] consider a variant which is consistent for any $d$ but no rate of convergence is provided. When $d = 1$, a central limit theorem for estimators based on ranks (i.e. nearest neighbors on the right) is also proved in [16].

In a nutshell, on the one hand, the class of PF estimators allows to estimate Sobol' indices at rate $\sqrt{n}$ for any $d$ with minimal assumptions on $G$ but requires a sample with highly specific structure. On the other hand, local-averaging estimators are based on a single input-output $n$-sample (standard Monte-Carlo) but need some regularity assumptions on the model and $\sqrt{n}$-parametric rate of convergence is proved only for $d \le 3$ (nearest neighbors method). In this paper, we propose and study a more general local-averaging estimator and show that the rate $\sqrt{n}$ holds for any $d$.

In our work, we build an optimal estimator. In that view, we correct a plug-in estimator using a high-order kernel and a mirror-type transformation. The correction improves the initial estimator by correcting the bias (see, e.g., [39, 37]), in line with techniques used to improve an estimator in the frame of semi-parametric inference. We refer to [51, Chapter 25] for an exhaustive overview of this theory. As a matter of fact, in our frame, the nonparametric part is the unknown regression function. It is handled by using a very specific kernel estimator. More precisely, on the one hand, we make use of high-order kernels in order to obtain $\sqrt{n}$ rate of convergence, following recent approaches of nonparametric regression (see for example [50]). Note that a high-order kernel has a non-zero negative part thus is not a probability density. Its use allows to cancel the first terms in the bias expansion (see the details of such computations in (29)). On the other hand, to remove boundary effects inherent to the kernel estimation procedure, we adapt recent mirror-type transformations [3, 41]. We propose and study two variants for the mirror transformation. In both cases, we show that the proposed estimators satisfy a central limit theorem with the optimal rate and minimal variance. Up to our knowledge, this general optimality result is new in the frame of Sobol' index estimation based on a single input-output $n$-sample. Notice that the order one case is tackled in [9] by using another estimator based on classical kernel estimation jointly with a spectral approach. Note finally that our results are closely related to the ones developed in [15] where the boundary issue is tackled using truncation leading to the estimation of a pruned version of $S^X$.

The paper is organized as follows. The setting and the notation are introduced in Section
in extending the very interesting point of view introduced in [22, 21] to estimate general nonlinear integral functionals of a density on the real line, by using empirically a kernel estimator erasing the diagonal terms. Relaxing the positiveness assumption on the kernel and choosing a kernel of order large enough, we are able to prove a central limit theorem for estimating Sobol' indices of any order (the bias is killed thanks to this signed kernel).

The paper is organized as follows. In Section 2, we introduce the setting and fix notation. Section 3 is devoted to the statement of our main results. First we introduce our kernel estimator in terms of a degenerate U-statistics. Then we state a central limit theorem. A comparison study on the efficiency is provided in Section 4. Section 5 proposes an heuristic inspired from [12] for bandwidth selection. Our methodology is extended in Section 6 to the case where the distribution of the input is unknown. Section 7 presents several numerical illustrations of our estimation procedure. Most of the proofs are postponed to Appendices A.1 and A.2. Furthermore, Appendix A.3 presents the practical construction of a kernel of fixed order.

# 2 Setting

## 2.1 Model and purpose

The quantity of interest (QoI) $Y$ is obtained from the numerical code and is regarded as a function $g$ of the vector of the random inputs $(V_i)_{i=1,...,p}$

$$Y = g(V_1, \dots, V_p),$$

where the function $g$ is defined on a compact subset of $\mathbb{R}^p$ and is real-valued. Classically, the $V_i$'s are assumed to be independent random variables and a sensitivity analysis is performed using the Hoeffding decomposition [1, 48] leading to the standard Sobol' indices [41].

Let $\mathbf{u} \subset \{1, \dots, p\}$. In the sequel, $X$ stands for a group of inputs: $X = (V_i)_{i \in \mathbf{u}}$ and $W = (V_i)_{i \notin \mathbf{u}}$ so that

$$Y = g(X, W) \tag{1}$$

(with a little abuse of notation). When $\mathbf{u}$ is reduced to a singleton, $X$ is simply one of the inputs $V_i$, $i = 1, \dots, p$ and $S^X$ is the first order Sobol' index with respect to this input. In the sequel, the cardinality of $\mathbf{u}$ is denoted by $d$: $|\mathbf{u}| = d$ with $1 \le d \le p$ and for the sake of simplicity, we assume that $\mathbf{u} = \{1, \dots, d\}$ without loss of generality.

Assume that $Y \in \mathbb{L}^2(\mathbb{R})$. We aim at estimating the Sobol' index of the output $Y$ with respect to the $d$-dimensional vector $X$; namely

$$S^X = \frac{\mathrm{Var}(\mathbb{E}[Y|X])}{\mathrm{Var}(Y)} = \frac{\mathbb{E}[\mathbb{E}[Y|X]^2] - \mathbb{E}[Y]^2}{\mathrm{Var}(Y)}. \tag{2}$$

Here, $d$ is the order of the Sobol' index while $p$ is the total number of inputs. When $\mathbf{u}$ is reduced to a singleton, $X$ is simply one of the inputs $V_i$, $i = 1, \dots, p$: $X = V_{i_0}$ and $S^X$ is the first order Sobol' index with respect to this input $V_{i_0}$.
2. Then, considering two different mirror-type transformations, we present in Section 3 two kernel-based regression estimators. Section 4 is devoted to the statement of our main results, namely central limit theorems and asymptotic efficiency. A comparison study of asymptotic estimation variance with existing methods is provided in Section 5. Section 6 presents several numerical illustrations of our estimation procedure. The proofs are postponed to Appendix A. Moreover, extended numerical studies discussing numerical stability of high-order kernels can be found as an additional file on the following link https://hal.science/hal-04052837.

# 2 Setting

## 2.1 Model and purpose

The output $Y$ is obtained from a regression model (generally computed by a numerical code or a machine learning estimated model), and is regarded as a function $G$ of random inputs $(V_i)_{i=1,...,p}$ such that $Y = G(V_1, ..., V_p)$ where the function $G$ is defined on a compact subset of $\mathbb{R}^p$, $p \ge 2$ and is real-valued. The random vector formed by the $p$ random input variables is denoted by $V$. Global sensitivity analysis is performed using Hoeffding decomposition [26] for $Y \in \mathbb{L}^2(\mathbb{R})$ which leads to standard Sobol' indices introduced below [45] assuming that the $V_i$'s are independent random variables. However, such independence assumption is only required for practical interpretation of the decomposition, but is not at all necessary for estimation purposes and asymptotic guarantees. Thus, in the rest of the paper, we do not work with such strong assumption.

For any $u \subsetneq \{1, ..., p\}$, $u \neq \emptyset$ we denote $X = (V_i)_{i \in u}$ a group of inputs with cardinality $|u| = d < p$. The closed Sobol' index of output $Y$ associated to the $d$-dimensional vector $X$ is defined as

$$S^X = \frac{\mathrm{Var}(\mathbb{E}[Y|X])}{\mathrm{Var}(Y)} = \frac{\mathbb{E}[\mathbb{E}[Y|X]^2] - \mathbb{E}[Y]^2}{\mathrm{Var}(Y)} \quad (1)$$

For the specific case $u = \{i_0\}$ with $i_0 \in \{1, ..., p\}$, that is for $X = V_{i_0}$, $S^{i_0} = S^X$ is the first-order Sobol' index associated to the input $V_{i_0}$. In addition, the total Sobol' index associated to $V_{i_0}$ writes

$$S^{\text{tot},i_0} = 1 - S^{V\sim i_0} = 1 - S^{\sim i_0} \quad (2)$$

where $V_{\sim i_0} = (V_1, ..., V_{i_0-1}, V_{i_0+1}, ..., V_p)$.

Our goal is to estimate $S^X$ from a $n$-sample $(X_j, Y_j)_{j=1,...,n}$ of the pair $(X, Y)$ with joint distribution $\mathcal{P}$. Since $\mathbb{E}[Y]$ and $\mathrm{Var}(Y)$ can naturally be estimated with classical empirical moments, we focus here on the estimation in (1) of $T = \mathbb{E}[\mathbb{E}[Y|X]^2]$.

## 2.2 Notation and assumptions

In this section, we give the general assumptions made in the paper.

(A1) The support of $(V_1, ..., V_p)$ is $[0, 1]^p$. We assume that $X$ is absolutely continuous with respect to the Lebesgue measure on $[0, 1]^d$ with density function $f_X$.

For simpler notation here we only consider $[0, 1]^p$ but results can be readily extended for support of the form $[B_1, C_1] \times \dots \times [B_p, C_p]$ where $B_i < C_i$ for all $i \in \{1, ..., p\}$.
To estimate $\mathbb{E}[Y]$ and $\mathrm{Var}(Y)$ from the $n$-sample $(X_j, Y_j)_{j=1,...,n}$ of the pair $(X, Y)$, we will naturally use the classical empirical mean and variance respectively. Thus we focus on the estimation of $\mathbb{E}[\mathbb{E}[Y|X]^2]$ in (2). In the sequel, we build a modified kernel-based procedure for the estimation of $S^X$ on a $n$-sample $(X_j, Y_j)_{j=1,...,n}$ of the pair $(X, Y)$. A natural estimator inspired from the nearest neighbor estimator as considered in [4] and the kernel-based plug-in estimator in [22] would be

$$ \binom{n}{2}^{-1} \sum_{1 \le j < j' \le n} \frac{Y_j Y_{j'}}{2} \left( \frac{K_{h_n}(X_{j'} - X_j)}{f_X(X_j)} + \frac{K_{h_n}(X_j - X_{j'})}{f_X(X_{j'})} \right) $$

for a bandwidth $h_n > 0$ and a kernel $K_{h_n}$. Nevertheless, as pointed in [3], boundary issues appear when the input domain is compact. To bypass this issue, we consider the following kernel-based estimator

$$ T_{n,h_n} := \binom{n}{2}^{-1} \sum_{1 \le j < j' \le n} \frac{Y_j Y_{j'}}{2} \left( \frac{K_{h_n} \circ A_{X_j}(X_{j'} - X_j)}{f_X(X_j)} + \frac{K_{h_n} \circ A_{X_{j'}}(X_j - X_{j'})}{f_X(X_{j'})} \right). \quad (3) $$

for a bandwidth $h_n > 0$, a mirror-type transformation $A_x$, and a kernel $K_{h_n}$ as defined in the next section.

## 2.2 Notation and assumptions

In this section, we give the general assumptions made in the paper.

(A1) The support of $(V_1, ..., V_p)$ is of the form $[B_1, C_1] \times \dots \times [B_p, C_p]$ where $B_i < C_i$ for all $1 \le i \le p$. We denote by $\mathcal{D}$ the support of $X$.

(A2) The distribution of the random vector $(V, Y) \in \mathbb{R}^p \times \mathbb{R}$ is absolutely continuous with respect to the Lebesgue measure. The marginal probability density function of $(X, Y)$, $V$, $X$, and $W$ are denoted by $f_{X,Y}$, $f_V$, $f_X$, and $f_W$ respectively.

The next definition is borrowed from [3, Section 2.2] and allows to circumvent the boundary issues.

(D1) For $x \in \mathcal{D}$, we define

$$ A_x: \left\{ \begin{array}{l} \mathbb{R}^d \quad \to \quad \mathbb{R} \\ u = (u_1, \dots, u_d) \quad \mapsto \quad (\sigma_1(x_1)u_1, \dots, \sigma_d(x_d)u_d) \end{array} \right. \quad (4) $$

with $\sigma_i(s) := 1 - 2\mathbb{1}_{\left(\frac{B_i+C_i}{2}, C_i\right)}(s) \in \{-1, 1\}$. Observe that $\mathcal{A} = \{A_x, x \in \mathcal{D}\}$ is a finite subset of $GL_d(\mathbb{R})$, $\mathcal{A} = \{A_1, \dots, A_\kappa\}$, with cardinal $\kappa = 2^d$. Moreover, it satisfies

(i) for any $\ell = 1, \dots, \kappa$, $|\det(A_\ell)| = 1$;

(ii) **Mirror condition:** for any $x \in \mathcal{D}$, there exists $A_\ell \in \mathcal{A}$ such that $A_x = A_\ell$ and $x + A_x^{-1}([0, 1/2]^d) \subset \mathcal{D}$.
(A2) $\exists c > 0$ such that $\inf_{x \in [0,1]^d} f_X(x) \ge c$.

(A3) $\mathbb{E}[Y^4] < \infty$ and $\sigma^2(x) = \mathrm{Var}(Y|X = x)$ is bounded on $[0, 1]^d$.

Note that under Assumption (A3) the regression function $m(\cdot) = \mathbb{E}[Y|X = \cdot]$ exists and is bounded on the support of $X$ by Cauchy-Schwartz inequality. Furthermore, we define the function $g := f_X m$.

For any integer $d$, we introduce the following multi-index notation. For any $\beta = (\beta_1, \dots, \beta_d) \in (\mathbb{R}_{>0})^d$, let $\lfloor \beta \rfloor$ be the largest integer strictly lower than $\beta$: $\lfloor \beta \rfloor = (\lfloor \beta_1 \rfloor, \dots, \lfloor \beta_d \rfloor) =: \gamma \in \mathbb{N}^d$. In addition, we introduce

$$|\gamma| = \gamma_1 + \dots + \gamma_d, \quad \gamma! = \gamma_1! \dots \gamma_d!, \quad \text{and} \quad v^\beta = v_1^{\beta_1} \dots v_d^{\beta_d} \text{ for any } v \in \mathbb{R}^d.$$

For $\alpha > 0$ let $C^\alpha([0, 1]^d)$ be the set of functions $\phi: [0, 1]^d \to \mathbb{R}$ that have derivatives up to order $\lfloor \alpha \rfloor$ and for which partial derivative of order $\lfloor \alpha \rfloor$ is $\alpha - \lfloor \alpha \rfloor$-Hölder. Namely, there exists $C_\phi > 0$ such that, for any $x$ and $x' \in [0, 1]^d$, one has

$$\left| \frac{\partial^\beta \phi}{\partial x^\beta}(x) - \frac{\partial^\beta \phi}{\partial x^\beta}(x') \right| \le C_\phi \|x - x'\|_{\infty}^{\alpha - \lfloor \alpha \rfloor} \quad (3)$$

for any $\beta \in \mathbb{N}^d$ such that $\beta = \lfloor \alpha \rfloor$ where $(\partial^\beta \cdot / \partial x^\beta)$ stands for the partial derivative of order $\beta$ and $\| \cdot \|_{\infty}$ for the supremum norm on $[0, 1]^d$.

(A4) The density $f_X$ of $X$ belongs to $C^\alpha([0, 1]^d)$ for some $\alpha > 0$.

(A5) The regression function $m$ belongs to $C^\alpha([0, 1]^d)$.

Notice that, if the model $G$ belongs to $C^\alpha$, then Assumption (A5) is satisfied.

(A6) Let $k: [0, 1] \to \mathbb{R}$ be a univariate kernel such that $\|k\|_{\infty} < \infty$ and $\int_0^1 k(v)dv = 1$. We assume that $k$ is of order $(\lfloor \alpha \rfloor + 1)$ which means that $\int_0^1 v^\ell k(v)dv = 0$ for any $\ell \in \mathbb{N}$ such that $0 < \ell \le \lfloor \alpha \rfloor$ and $\int_0^1 v^{\lfloor \alpha \rfloor+1} k(v)dv \ne 0$. Furthermore, we define the multivariate kernel $K: [0, 1]^d \to \mathbb{R}$ as: $K(v) = \prod_{k=1}^d k(v_k)$ for any $v = (v_1, \dots, v_d) \in [0, 1]^d$. Finally, we define $K_h(v) = K(v/h)/h^d = \prod_{k=1}^d k(v_k/h)/h^d$ for any $v = (v_1, \dots, v_d) \in [0, 1]^d$.

Observe that $\int_{[0,1]^d} v^\beta K(v)dv = 0$ for any $\beta \in \mathbb{N}^d$ such that $0 < |\beta| \le \lfloor \alpha \rfloor$ and $\int_{[0,1]^d} v^\beta K(v)dv \ne 0$ for some $\beta$ such that $|\beta| = \lfloor \alpha \rfloor + 1$.

(A7) The sequence $(h_n)_{n \in \mathbb{N}}$ of bandwidths is positive and such that $h_n \to 0$ as $n \to \infty$.

The previous assumptions are pretty standard. See for instance [50]. More precisely, Assumptions (A1), (A2), (A4), (A5), (A6), and (A7) are common in the literature on nonparametric kernel regression estimation. Assumption A3 naturally stems in our context where we aim to prove a central limit theorem for a quadratic functional of the regression function.

In the rest of the paper, $C$ denotes a generic constant (deterministic and finite) which may vary from line to line.
For any integer $d$, let us introduce the multi-index notation. For any $\beta = (\beta_1, \dots, \beta_d) \in \mathbb{R}_+^d$, we define the integer part of $\beta$ by

$$ \lfloor \beta \rfloor = (\lfloor \beta_1 \rfloor, \dots, \lfloor \beta_d \rfloor) =: \gamma \in \mathbb{N}^d. $$

In addition, we introduce

$$ |\gamma| = \gamma_1 + \dots + \gamma_d, \quad \gamma! = \gamma_1! \dots \gamma_d!, \quad \text{and} \quad v^\beta = v_1^{\beta_1} \dots v_d^{\beta_d} $$

for any $v \in \mathbb{R}^d$.

Let $\alpha > 0$. We define $C^\alpha(\mathcal{D})$ as the set of functions $\phi: \mathcal{D} \to \mathbb{R}$ that have derivatives up to order $\lfloor \alpha \rfloor$ and which partial derivative of order $\lfloor \alpha \rfloor$ is $\alpha - \lfloor \alpha \rfloor$-Hölder. Namely, there exists $C_\phi > 0$ such that, for any $x$ and $x' \in \mathcal{D}$, one has

$$ \left| \frac{\partial^\beta \phi}{\partial x^\beta}(x) - \frac{\partial^\beta \phi}{\partial x^\beta}(x') \right| \le C_\phi \|x - x'\|_{\infty}^{\alpha - \lfloor \alpha \rfloor} \quad (5) $$

for any $\beta \in \mathbb{N}^d$ such that $|\beta| = \lfloor \alpha \rfloor$ and where $(\partial^\beta \cdot / \partial x^\beta)$ stands for the partial derivative of order $\beta$. Here, $\| \cdot \|_\infty$ stands for the supremum norm.

(A3) Let $K: \mathbb{R}^d \to \mathbb{R}$ be a kernel with support included in $\mathcal{D}$ such that $\|K\|_\infty < +\infty$ and $\int_{\mathcal{D}} K(u)du = 1$. We assume that $K$ is of order $\lfloor \alpha \rfloor$ which means that $\int_{\mathcal{D}} u^\beta K(u)du = 0$ for any $\beta \in \mathbb{N}^d$ such that $0 < |\beta| \le \lfloor \alpha \rfloor$. Finally, we define $K_h(x) = K(x/h)/h^d$ for any $x \in \mathcal{D}$.

(A4) The sequence $(h_n)_{n \in \mathbb{N}}$ of bandwidths is positive and such that $h_n \to 0$ as $n$ goes to $\infty$.

In this article, we write $C$ for a deterministic and finite constant which value is allowed to change between occurrences.

### 3 Central limit theorem

Recall the definition of our estimator given in (3):

$$ T_{n,h_n} = \binom{n}{2}^{-1} \sum_{1 \le j < j' \le n} \frac{Y_j Y_{j'}}{2} \left( \frac{K_{h_n} \circ A_{X_j}(X_{j'} - X_j)}{f_X(X_j)} + \frac{K_{h_n} \circ A_{X_{j'}}(X_j - X_{j'})}{f_X(X_{j'})} \right). $$

The goal is to derive a central limit theorem for $T_{n,h_n}$. Since $T_{n,h_n}$ is a U-statistics, our analysis is based on the Hoeffding decomposition of U-statistics [26]. More precisely, for all $j, j' = 1, \dots, n$, we introduce the symmetric function given by

$$ R\left(\binom{x_j}{y_j}, \binom{x_{j'}}{y_{j'}}\right) = \frac{y_j y_{j'}}{2} \left( \frac{K_{h_n} \circ A_{x_j}(x_{j'} - x_j)}{f_X(x_j)} + \frac{K_{h_n} \circ A_{x_{j'}}(x_j - x_{j'})}{f_X(x_{j'})} \right). $$

Then the Hoeffding projections of $R$ are given by

$$ \begin{aligned} \pi_1 R \left( \begin{matrix} x \\ y \end{matrix} \right) &= \mathbb{E} \left[ R \left( \begin{matrix} x \\ y \end{matrix}, \begin{matrix} X_2 \\ Y_2 \end{matrix} \right) \right] - \mathbb{E} \left[ R \left( \begin{matrix} X_1 \\ Y_1 \end{matrix}, \begin{matrix} X_2 \\ Y_2 \end{matrix} \right) \right] \\ \pi_2 R \left( \begin{matrix} x_1 \\ y_1 \end{matrix}, \begin{matrix} x_2 \\ y_2 \end{matrix} \right) &= R \left( \begin{matrix} x_1 \\ y_1 \end{matrix}, \begin{matrix} x_2 \\ y_2 \end{matrix} \right) - \mathbb{E} \left[ R \left( \begin{matrix} x_1 \\ y_1 \end{matrix}, \begin{matrix} X_2 \\ Y_2 \end{matrix} \right) \right] - \mathbb{E} \left[ R \left( \begin{matrix} x_2 \\ y_2 \end{matrix}, \begin{matrix} X_2 \\ Y_2 \end{matrix} \right) \right] \\ &\quad + \mathbb{E} \left[ R \left( \begin{matrix} X_1 \\ Y_1 \end{matrix}, \begin{matrix} X_2 \\ Y_2 \end{matrix} \right) \right] \end{aligned} $$
# 3 Estimation procedure and preliminary results

In this section, we propose an estimator of $T$ based on two main ingredients: (a) estimation based on the efficient influence function of $T$ (see, e.g., [4, 15, 51]) and (b) mirror-type kernel estimators (see, e.g., [3, 41]).

Let us first exhibit the general form of the estimator of $T$ considered in this paper:

$$T_n = \frac{1}{n} \sum_{i=1}^{n} (2Y_i - m_n(X_i)) m_n(X_i). \quad (4)$$

Here, $m_n$ is an estimator of the regression function $m$ of $Y$ on $X$ (with properties to be discussed later), and $(X_i, Y_i)_{1 \le i \le n}$ are independent copies of $(X, Y)$.

Such a form can actually be justified and inferred from a heuristic based on efficient influence functions. Indeed, let $\mathcal{P}$ be the set of probability distributions on $[0, 1]^d \times \mathbb{R}$ satisfying Assumptions (A2), (A3), (A4) and (A5). Then we know from [15] that

$$(x, y) \mapsto (2y - m(x))m(x) - \mathbb{E}[\mathbb{E}[Y|X]]$$

is the efficient influence function of $T$ under the nonparametric model $\mathcal{P}$ (see [31] for detailed computations). Thus, if the probability distribution of $(X, Y)$ is in $\mathcal{P}$ and if $m$ is known, we can state, from [51, Equation (25.22)], that the estimator

$$T_n = \frac{1}{n} \sum_{i=1}^{n} (2Y_i - m(X_i)) m(X_i) \quad (5)$$

is asymptotically efficient with optimal variance equal to $\text{Var}((2Y - m(X))m(X))$. Unfortunately, the estimator in (5) is only an oracle, as $m$ is unknown in our case and needs then to be estimated, but this explains the intuition behind our focus on (4). As for the choice of the estimator of $m$, since we assume that the domain of the inputs is compact, the crucial point is to handle possible boundary effects. To do so, in [15] the authors estimate a truncated version of $T$ defined as $T^{\text{trunc},\varepsilon} = \mathbb{E}[\mathbb{E}[Y|X]^2 \mathbb{1}_{X \in (\varepsilon, 1-\varepsilon)d}]$. Even if $T^{\text{trunc},\varepsilon} \to T$ as $\varepsilon \to 0$ under mild assumptions, the practical tuning of the parameter $\varepsilon$ depends on the unknown function $G$ and its choice has a large impact, see Figure 8 in the numerical experiment section.

Here, we therefore prefer the use of plugging mirror-type kernel estimators for $m(\cdot)$ in (5). Indeed, kernel-based estimators are known to induce a mass loss near the boundary. We call a mirror-type transformation any transformation that corrects this phenomenon by pushing back the mass in the interior of the domain. Mirror-type transformations have been used, e.g., in [3, 41] in the framework of density estimation. In the following, we adapt these works to propose two nonparametric estimators of $m$, denoted as $\hat{m}$ and $\tilde{m}$. We then prove that both estimators $\tilde{T}_n$ and $\hat{T}_n$ satisfy a central limit theorem with optimal asymptotic variance $\text{Var}((2Y - m(X))m(X))$.

## 3.1 A mirror-type kernel estimator for the regression function

To estimate the regression function $m$, we take inspiration from Naradaya-Watson estimator [36, 52] and the kernel-based plug-in estimator studied in [22]. More precisely, we consider a leave-one-out kernel estimator together with a mirror-type transformation
(see [35, page 137]) leading to the Hoeffding decomposition

$$T_{n,h_n} - \mathbb{E}[T_{n,h_n}] = 2U_n^{(1)}(\pi_1 R) + U_n^{(2)}(\pi_2 R)$$

with

$$\begin{aligned}
U_n^{(1)}(\pi_1 R) &= \frac{1}{n} \sum_{j=1}^{n} (\pi_1 R) \begin{pmatrix} X_j \\ Y_j \end{pmatrix} \\
U_n^{(2)}(\pi_2 R) &= \frac{2}{n(n-1)} \sum_{1 \le j < j' \le n}^{n} (\pi_2 R) \begin{pmatrix} X_j \\ Y_j \end{pmatrix} \begin{pmatrix} X_{j'} \\ Y_{j'} \end{pmatrix}.
\end{aligned}$$

Then following [22] we write

$$\begin{aligned}
T_{n,h_n} - \mathbb{E}[T_{n,h_n}] ={}& \frac{1}{n} \sum_{j=1}^{n} 2(Y_j g_1(X_j) - \mathbb{E}[\mathbb{E}[Y|X]^2]) \\
&+ 2U_n^{(1)}(\pi_1 R) - \frac{2}{n} \sum_{j=1}^{n} (Y_j g_1(X_j) - \mathbb{E}[\mathbb{E}[Y|X]^2]) + \underbrace{U_n^{(2)}(\pi_2 R)}_{S_2}, \quad (6)
\end{aligned}$$

where $g_1(x) = \mathbb{E}[Y|X=x]$. Finally, we define $g_2(x) = \mathbb{E}[Y^2|X=x]$.

**Theorem 3.1.** Suppose that Assumptions A₁ to A₄ hold and $\mathbb{E}[Y^2] < \infty$. Let $\alpha > 0$.

I.1. Assume that $g_1 \in \mathbb{L}^1(\mathcal{D})$, $g_1 f_X \in \mathcal{C}^\alpha(\mathcal{D})$, and

$$\int_{\mathcal{D}} \left| g_1(x) \frac{\partial^\beta (g_1 f_X)}{\partial x^\beta} (x) \right| dx < \infty$$

for any $\beta$ such that $1 \le |\beta| < \lfloor \alpha \rfloor$. Then we have

$$\left| \mathbb{E}[T_{n,h_n}] - \mathbb{E}[\mathbb{E}[Y|X]^2] \right| \le Ch_n^\alpha. \quad (7)$$

I.2. Assume in addition that $g_1 \in \mathcal{C}^\alpha(\mathcal{D})$, $g_2/f_X \in \mathbb{L}^1(\mathcal{D}) \cap \mathbb{L}^2(\mathcal{D})$ and $g_2 f_X \in \mathbb{L}^2(\mathcal{D})$, and

$$\int_{\mathcal{D}} g_2(x) \left| \frac{\partial^\beta (g_1 f_X)}{\partial x^\beta} (x) \frac{\partial^{\beta'} (g_1 f_X)}{\partial x^{\beta'}} (x) \right| dx < \infty$$

for any $\beta$ and $\beta'$ such that $1 \le |\beta|, |\beta'| < \lfloor \alpha \rfloor$. Then we have

$$\mathbb{E}\left[\left(T_{n,h_n} - \mathbb{E}[T_{n,h_n}] - \frac{1}{n} \sum_{j=1}^{n} Z_j\right)^2\right] \le Ch_n^{2\alpha} + \frac{C}{h_n^d n^2}. \quad (8)$$

where, for $j = 1, \dots, n$,

$$Z_j = 2(Y_j g_1(X_j) - \mathbb{E}[\mathbb{E}[Y|X]^2]). \quad (9)$$

II. Assuming in addition that $\mathbb{E}[Y^4] < \infty$, $\alpha > d/2$, $h_n \xrightarrow{n \to \infty} 0$, $nh_n^d \xrightarrow{n \to \infty} \infty$, and $nh_n^{2\alpha} \xrightarrow{n \to \infty} 0$, we get

$$\sqrt{n} \left( T_{n,h_n} - \mathbb{E}[\mathbb{E}[Y|X]^2] \right) \xrightarrow[n \to \infty]{\mathcal{D}} \mathcal{N}(0, 4\tau^2) \quad (10)$$

with $\tau^2 = \text{Var}(Yg_1(X))$.
introduced in [3] to bypass boundary issues. More precisely, the regression function estimator is defined as follows:

$$ \hat{m}_{n,h_n,i}(X_i) = \frac{\sum_{j \neq i} Y_j K_{h_n} \circ A_{X_i}(X_j - X_i)}{\sum_{j \neq i} K_{h_n} \circ A_{X_i}(X_j - X_i)} \mathbb{1} \sum_{j \neq i} K_{h_n} \circ A_{X_i}(X_j - X_i) \neq 0 \quad (6) $$

for a bandwidth $h_n > 0$, a mirror-type transformation A, and a kernel $K_{h_n}$ satisfying Assumption (A6) presented in Section 2.2. Then, (5) is adapted as

$$ \hat{T}_{n,h_n} = \frac{1}{n} \sum_{i=1}^{n} (2Y_i - \hat{m}_{n,h_n,i}(X_i)) \hat{m}_{n,h_n,i}(X_i). \quad (7) $$

As for the mirror-transformation, for $x \in [0, 1]^d$, one may consider for instance

$$ A_x: \left\{ \begin{array}{l} \mathbb{R}^d \quad \rightarrow \quad \mathbb{R}^d \\ v = (v_1, \dots, v_d) \quad \mapsto \quad (a_1(x_1)v_1, \dots, a_d(x_d)v_d) \end{array} \right. \quad (8) $$

with $a_i(s) := 1 - 2\mathbb{1}_{\{\frac{1}{2},1\}}(s) \in \{-1, 1\}$, see Figure 1 for an illustration. Observe that $\{A_x, x \in [0, 1]^d\}$ is a finite subset of $GL_d(\mathbb{R})$ with cardinality $\kappa = 2^d$ (where $GL_d(\mathbb{R})$ is the general linear group on $\mathbb{R}$). Then we denote these $\kappa$ elements $A_1, \dots, A_\kappa$. Moreover, it satisfies

(i) for any $\ell = 1, \dots, \kappa$, $|\det(A_\ell)| = 1$;

(ii) **Mirror property:**

$$ \forall x \in [0, 1]^d, x + A_x^{-1}([0, 1/2]^d) \subset [0, 1]^d. \quad (9) $$

The $C^\alpha([0, 1]^d)$ regularity assumption can be relaxed to piecewise regularity as soon as the number $k$ of pieces is finite. The cardinality of $\{A_x, x \in [0, 1]^d\}$ is then increased to $\kappa = (2k)^d$.

Let $\hat{f}_{n,h_n,i}$ be the leave-one-out estimator of the density function $f_X$ based on the $(n-1)$-sample $(X_1, \dots, X_{i-1}, \dots, X_{i+1}, \dots, X_n)$:

$$ \hat{f}_{n,h_n,i}(x) = \frac{1}{n-1} \sum_{j \neq i} K_{h_n} \circ A_x(X_j - x). \quad (10) $$

The following lemmas establish classical controls on the bias and on the variance of $\hat{f}_{n,h_n,i}$ for all $i \in \{1, \dots, d\}$. The proofs are postponed to Appendix A.1.

**Lemma 3.1.** Under Assumptions (A1), (A4), and (A6), for all $i \in \{1, \dots, d\}$,

$$ \|\mathbb{E}[\hat{f}_{n,h_n,i}] - f_X\|_{\infty} = \|\mathbb{E}[\hat{f}_{n,h_n,1}] - f_X\|_{\infty} = O(h_n^{\alpha}). \quad (11) $$

To control the variance, one needs to define a supplementary assumption.

(A8) Assume that the family of functions

$$ \mathcal{F} = \{K_{h,x}: y \in [0, 1]^d \mapsto K_h(y-x) \in \mathbb{R}, h > 0, x \in [0, 1]^d\} \quad (12) $$
The proof of Theorem 3.1 is postponed to Appendix A.1.

Remark 3.2.

1. The assumptions of Theorem 3.1 are fulfilled as soon as $g$ and $f_X \in C^\alpha(\mathcal{D})$ and $f_X$ is bounded from below by a strictly positive constant.

2. The $C^\alpha(\mathcal{D})$ regularity assumption can be relaxed to piecewise regularity, as soon as the number of pieces $m$ is finite. The cardinality of $\mathcal{A}$ is then increased to $\kappa = (2m)^d$.

3. We still assume that $g$ and $f_X \in C^\alpha(\mathcal{D})$ but one may relax the condition of boundedness from below on $f_X$. Since the inputs $V_i$'s are independent, $f_X = \prod_{i=1}^p f_{V_i}$, it amounts to control the local behavior of each marginal $f_{V_i}$ at any point $v_0 \in \mathcal{D}$ such that $f_{V_i}(v_0) = 0$. We indeed assume that for any such $v_0$, $f_{V_i}(v) = (v - v_0)^r h_{V_i}(v)$ with $h_{V_i}(v_0) \neq 0$ and $0 < r < 1/2$. For instance, one may take all the $V_i$'s Beta-distributed:

$$f_{V_i}(v) = \frac{v^{a_i-1}(1-v)^{b_i-1}}{B(a_i, b_i)} \mathbb{I}_{[0,1]}(v) \quad (11)$$

with parameters $a_i$ and $b_i \in (1, 3/2)$ for $1 \le i \le p$ and $B(a_i, b_i)$ the normalizing constant equal to the beta function.

4. An example of convenient bandwidths $(h_n)_{n \ge 0}$ is given by $h_n = n^{-\gamma}$ with $1/(2\alpha) < \gamma < 1/d$.

Using the delta method in [48, Theorem 3.1], we are now able to get the asymptotic behavior of the estimation of $S^X$.

**Corollary 3.3.** Under all the assumptions of Theorem 3.1 (II included), one has

$$\sqrt{n} \left( \frac{T_{n,h_n} - \left( \frac{1}{n} \sum_{j=1}^{n} Y_j \right)^2}{\frac{1}{n} \sum_{j=1}^{n} Y_j^2 - \left( \frac{1}{n} \sum_{j=1}^{n} Y_j \right)^2} - S^X \right) \xrightarrow{n \to \infty} \mathcal{N}(0, \sigma^2), \quad (12)$$

where the limit variance $\sigma^2$ has an explicit expression given in (34).

Using one more time the delta method in [48, Theorem 3.1], we deduce the asymptotic behavior of the vector containing the $p$ first-order Sobol' indices. Let us denote $S^i$ the first-order Sobol index associated to $X = V_i$ and its estimator $\hat{S}^i$ given by:

$$\hat{S}^i := \frac{T_{n,h_n} - \left( \frac{1}{n} \sum_{j=1}^{n} Y_j \right)^2}{\frac{1}{n} \sum_{j=1}^{n} Y_j^2 - \left( \frac{1}{n} \sum_{j=1}^{n} Y_j \right)^2}.$$

**Corollary 3.4.** Under all the assumptions of Theorem 3.1 (II included), one has

$$\sqrt{n} \left( (\hat{S}^1, \dots, \hat{S}^p)^T - (S^1, \dots, S^p)^T \right) \xrightarrow{n \to \infty} \mathcal{N}(0, \Sigma), \quad (13)$$

where the limit variance $\Sigma$ has an explicit expression given in (35).
<figure><img src="image_2.png" /><figcaption>Figure 1: Mirror-type transformation defined in (8) with $d = 2$, for $x = (1/3, 3/4)$, and for $y = (2/3, 1/5)$.</figcaption></figure>

is a uniformly bounded Vapnik-Chervonenkis (VC)-class of functions, i.e. there exists positive numbers $A$, $B$, and $v$ such that, for all $K \in \mathcal{F}$, $\|K\|_{\infty} < \infty$ and the $\varepsilon$-covering number $\mathcal{N}(\mathcal{F}, L^2(Q), \varepsilon)$ of $\mathcal{F}$ for the $L^2(Q)$-distance satisfies

$$ \mathcal{N}(\mathcal{F}, L^2(Q), \varepsilon) \leq \left( \frac{AB}{\varepsilon} \right)^v $$

for every probability measure $Q$ on $[0, 1]^d$ and for every $\varepsilon \in (0, B)$.

Assumption (A8) is a common hypothesis, called the uniformly bounded VC class assumption, that is required to prove convergence of kernel density estimator (see, e.g., [21]). Such an assumption is satisfied for commonly used kernels as Gaussian, Epanechnikov, Uniform, as mentioned in [30].

**Lemma 3.2.** Let $\delta_n \in (0, 1)$. Under Assumptions (A1), (A6), and (A8), there exists a constant $C > 0$ such that we have, with probability $1 - \delta_n$,

$$ \| \hat{f}_{n,h_n,i} - \mathbb{E}[\hat{f}_{n,h_n,i}] \|_{\infty}^2 \leq C \frac{\log(\frac{1}{h_n}) + \log(\frac{2}{\delta_n})}{nh_n^d} \quad (13) $$

for all $i \in \{1, \dots, d\}$.

**Remark 3.3.** Choosing $h_n = n^{-1/(2\alpha+d)}$ and $\delta_n = 1/n$ in Lemmas 3.1 and 3.2, we get the following bound with probability $1 - 1/n$:

$$ \| \hat{f}_{n,h_n,i} - f_X \|_{\infty} \leq C \sqrt{\log(n)} n^{-\alpha/(2\alpha+d)} $$

which corresponds to the optimal nonparametric rate up to the log factor.

Now, we can control the mean integrated squared error of $\hat{f}_{n,h_n,i}$ together with its infimum.
Remark 3.5. The proof a Corollary 3.4 does not rely on the fact that we are dealing with order one indices, hence if $u_1, \dots, u_r$ are distinct subsets of $\{1, \dots, p\}$ we also have under the same assumptions that

$$ \sqrt{n} \left( (\hat{S}^{u_1}, \dots, \hat{S}^{u_r})^T - (S^{u_1}, \dots, S^{u_r})^T \right) \xrightarrow[n \to \infty]{D} \mathcal{N}(0, \Sigma_r), \quad (14) $$

where the limit variance $\Sigma_r$ has an explicit expression similar to (35).

# 4 Efficiency

In this section, we compare several estimators proposed in the literature to estimate $\mathbb{E}[\mathbb{E}[Y|X]^2]$ (or even $S^X$) in terms of limiting variances obtained in the central limit theorem.

## Estimating $\mathbb{E}[\mathbb{E}[Y|X]^2]$

Using a unique $n$-sample, one may compare the kernel estimator introduced in this paper and the procedure involving the estimators $\hat{T}_n$ defined in [9, page 11]. Such estimators are particularly tricky to compute and not easily tractable in practice. More precisely, the initial $n$-sample is split into two samples of sizes $n_1 = \lfloor n/\log n \rfloor$ and $n_2 = n - n_1 \approx n$. The first sample is dedicated to the estimation of the joint density of $(X, Y)$ while the second one is used to compute a Monte-Carlo estimation of the integral involved in the quantity of interest.

In the following, we give the limiting variance $\sigma_T^2$ associated to the sequence $(\hat{T}_n)_n$. By [9, Theorems 3.4 and 3.5], this sequence is asymptotically efficient to estimate $\mathbb{E}[\mathbb{E}[Y|X]^2]$; so that the variance $\sigma_T^2$ is the minimal variance expected for an estimator of $\mathbb{E}[\mathbb{E}[Y|X]^2]$ based on the previous design of experiment (a unique $n$-sample $(X_i, Y_i)_{i=1,\dots,n}$ of the input-output pair $(X, Y)$). We refer the reader to [9, 17].

**Proposition 4.1** (Theorems 3.4 and 3.5 in [9]). The minimal variance $\sigma_T^2$ is given by

$$ \sigma_T^2 = \mathrm{Var}(g_1(X)(2Y - g_1(X))) = 4\tau^2 - 3\mathrm{Var}(g_1^2(X)) $$

with $\tau^2 = \mathrm{Var}(Yg_1(X))$.

One may also compare our results to the estimation procedure proposed in [15]. Here again the initial $2n$-sample is split into two samples of sizes $n$. The first sample allows to estimate the regression function $\mathbb{E}[Y|X = x]$ using the first nearest neighbor of $x$ among the points of the first sample while the second sample is used as a plug-in estimator. They proved a central limit theorem. Nevertheless, in their result the recentering factor is not the unknown quantity of interest $\mathbb{E}[\mathbb{E}[Y|X]^2]$ as required but the expectation of the estimator. In addition, the bias is negligible if and only if the dimension $d \le 3$ so that one can not expect a central limit theorem for the difference between the estimator and the quantity of interest $\mathbb{E}[\mathbb{E}[Y|X]^2]$ as soon as $d \ge 4$. In the following proposition, we give the limiting variance obtained with this estimation procedure. The factor 2 takes into account that we have considered two samples of size $n$.
**Corollary 3.4.** Under Assumptions (A1), (A2), (A4), (A6), and (A8), one has, for all $i \in \{1, \dots, d\}$,

$$ \mathbb{E}[\int_{[0,1]^d} (\hat{f}_{n,h_n,i}(x) - f_X(x))^2 dx] = o(n^{-1/2}) \quad (14) $$

and

$$ \frac{1}{\inf_{x \in [0,1]^d} |\hat{f}_{n,h_n,i}(x)|} = O_P(1) \quad (15) $$

as soon as $nh_n^{2d} \to \infty$ and $nh_n^{4\alpha} \to 0$ as $n \to \infty$.

Assuming $h_n = n^{-\gamma}$ with $\gamma > 0$, the conditions $nh_n^{2d} \to \infty$ and $nh_n^{4\alpha} \to 0$ lead to $1/(4\alpha) < \gamma < 1/(2d)$ and $\alpha > d/2$. The proof of Corollary 3.4 is postponed to Appendix A.1.

## 3.2 An alternative regression estimator

This alternative estimator is based on the mirror-image kernel density estimator proposed in [41] (see also [35]), introduced as a generalization to dimension $d \ge 2$ of the estimator introduced in [27, 44] for the dimension $d = 1$ and generalized in [1] to the dimension $d = 2$ (see also [19]). More precisely, we consider the following transformations: for any $z \in [0, 1]$, $m^{-1}(z) = -z$, $m^0(z) = z$, and $m^1(z) = 2 - z$ and, for any $a \in \{-1, 0, 1\}^d$ and $x \in [0, 1]^d$, the $d$-dimensional vector $M^a(x) = (m^{a_1}(x_1), \dots, m^{a_d}(x_d))$. Then, let

$$ \tilde{f}_{n,h_n,i}(x) = \mathbb{1}_{[0,1]^d}(x) \frac{1}{(n-1)} \sum_{j \neq i} \sum_{a \in \{-1,0,1\}^d} \tilde{K}_{h_n}(M^a(X_j) - x), \quad (16) $$

and

$$ \tilde{g}_{n,h_n,i}(x) = \mathbb{1}_{[0,1]^d}(x) \frac{1}{(n-1)} \sum_{j \neq i} Y_j \sum_{a \in \{-1,0,1\}^d} \tilde{K}_{h_n}(M^a(X_j) - x). \quad (17) $$

Here, the bandwidth $h_n > 0$ and the multivariate kernel $\tilde{K}$ satisfy Assumption (A'6) stated below. Now we propose the following regression function estimator:

$$ \tilde{m}_{n,h_n,i}(X_i) = \frac{\tilde{g}_{n,h_n,i}(X_i)}{\tilde{f}_{n,h_n,i}(X_i)} \quad (18) $$

if $\tilde{f}_{n,h_n,i}(X_i) \neq 0$ and 0 otherwise. The associated plug-in estimator then becomes:

$$ \tilde{T}_{n,h_n} = \frac{1}{n} \sum_{i=1}^{n} (2Y_i - \tilde{m}_{n,h_n,i}(X_i)) \tilde{m}_{n,h_n,i}(X_i). \quad (19) $$

The assumptions on $f_X$ and $g$ are strengthened as follows.

(A'4) The density function $f_X \in C^\alpha([0, 1]^d)$ for some $\alpha > 0$. Its derivatives of order $\beta$, with $0 < \beta \le \lfloor \alpha \rfloor$, vanish near the boundary.
**Proposition 4.2** (Theorem 2.1 in [15]). Under the assumptions of Theorem 2.1 in [15] and for $d \le 3$, the limiting variance $\sigma_D^2$ is given by

$$ \sigma_D^2 = 2\left(\mathbb{E}[g_2(X)^2] - \mathbb{E}[g_1(X)^2]^2 + \frac{1}{2}\left(\mathbb{E}[g_2(X)g_1(X)^2] - \mathbb{E}[g_1(X)^4]\right)\right). $$

When $d = 1$, one may use the estimation procedure based on ranks introduced in [6] and studied in [17].

**Proposition 4.3** (Proof of Theorem 4.1 in [17]). Under the assumptions of Theorem 4.1 in [17], the limiting variance $\sigma_D^2$ is given by

$$ \sigma_R^2 = \mathbb{E}\left[\mathrm{Var}(YY'|X)\right] + \mathbb{E}\left[Cov(YY', YY''|X)\right] - \mathbb{E}[(Y + Y')g_x(X, W)X]^2 \\ + \mathbb{E}[(Y + Y')(\tilde{Y} + \tilde{Y}')g_x(X, W)g_x(\tilde{X}, \tilde{W})(X \wedge \tilde{X})], $$

where $Y = g(V_1, V_2, \dots, V_p) = g(X, W)$ (see (1)), $Y' = g(X, W')$, $Y'' = g(X, W'')$, $\tilde{Y} = g(\tilde{X}, \tilde{W})$, and $\tilde{Y}' = g(\tilde{X}, \tilde{W}')$ with $X$ and $\tilde{X}$ i.i.d., $W, \tilde{W}, W'$, and $W''$ i.i.d. also independent of $X$ and $\tilde{X}$. Note that $YY'$ and $Y''$ (respectively $\tilde{Y}'$) are Pick-Freeze versions of $Y$ (resp. $\tilde{Y}$). Here $g_x$ stands for the first derivative of $g$ with respect to the first coordinate.

## Estimating $S^X$

The sequence of estimators $\hat{T}_n$ being efficient to estimate $\mathbb{E}[\mathbb{E}[Y|X]^2]$, we get naturally an asymptotically efficient sequence of estimators of $S^X$. The proof of the following proposition relies on a result on efficiency in product space [48, Theorem 25.50] and another result on efficiency and delta method [48, Theorem 25.47]. The details can be found in [17].

**Proposition 4.4.** Consider the sequence of estimators $\hat{T}_n$ introduced in [9, page 11]. Assume that the joint distribution $P$ of $(X, Y)$ is absolutely continuous with respect to the product probability $P_X \otimes P_Y$, namely $P(dx, dy) = f_{X,Y}(x, y)P_X(dx)P_Y(dy)$. Then the sequence $(\hat{R}_n)_n$

$$ \hat{R}_n = \frac{\hat{T}_n - \left(\frac{1}{n} \sum_{i=1}^n Y_j\right)^2}{\frac{1}{n} \sum_{j=1}^n Y_j^2 - \left(\frac{1}{n} \sum_{j=1}^n Y_i\right)^2} $$

is asymptotically efficient in estimating $S^X$. In addition, its (minimal) variance $\sigma_{\min}^2$ is

$$ \sigma_{\min}^2 = \frac{1}{\mathrm{Var}(Y)^2} \mathrm{Var}\left(2\mathbb{E}[Y](1-S^X)Y + S^X Y^2 + \mathbb{E}[Y|X](\mathbb{E}[Y|X] - 2Y)\right). $$

The delta method can also be applied to compute the asymptotic variance for the estimation of $S^X$ with the other estimation procedures, namely the procedure in [15] and the one in [17], with the empirical mean and empirical variance to estimate respectively $\mathbb{E}[Y]$ and $\mathrm{Var}(Y)$.
<figure><img src="image_3.png" /><figcaption>Figure 2: Mirror-image transformation with $d = 2$ in red. In the left-hand side of the figure, a data-point in $[0, 1]^d$ (in green) and its 8 mirror-images (in orange). A kernel is fitted over all points of the augmented dataset (orange areas). The darker orange regions correspond to the regions where several kernels overlap. In the right-hand side of the figure, the data-point (in green) and the restriction to $[0, 1]^d$ of the augmented dataset (in green).</figcaption></figure>

Assumption (A'4) appears to be strict, and is satisfied for example by uniform and beta distributions. Nevertheless, from a practical point of view, we have performed several simulations (not reported here) using Gaussian distributions that violate this assumption, but with estimation accuracy similar to the uniform setting.

We now list the assumptions on $K$ below.

(A'6) Let $\tilde{k}: [-1, 1] \to \mathbb{R}$ be a univariate kernel such that $\|\tilde{k}\|_{\infty} < \infty$ and $\int_{-1}^{1} \tilde{k}(v)dv = 1$. We assume that $\tilde{k}$ is of order ($\lfloor \alpha \rfloor + 1$) which means that $\int_{-1}^{1} v^l \tilde{k}(v)dv = 0$ for any $l \in \mathbb{N}$ such that $0 < l \le \lfloor \alpha \rfloor$ and $\int_{-1}^{1} v^{\lfloor \alpha \rfloor+1} \tilde{k}(v)dv \ne 0$. Furthermore, we define the multivariate kernel $\tilde{K}: [-1, 1]^d \to \mathbb{R}$ as: $\tilde{K}(v) = \prod_{k=1}^{d} \tilde{k}(v_k)$ for any $v = (v_1, \dots, v_d) \in [-1, 1]^d$. Finally, we define $\tilde{K}_h(v) = \tilde{K}(v/h)/h^d = \prod_{k=1}^{d} \tilde{k}(v_k/h)/h^d$ for any $v = (v_1, \dots, v_d) \in [-1, 1]^d$.

Both numerator and denominator of the estimator defined in (18) are the sum of $3^d$ terms; one corresponds to the original data in the unit hypercube $[0, 1]^d$, and each of the remaining terms corresponds to reflected data across one of the 0-faces, 1-faces, ..., $d-1$-faces of the unit hypercube as illustrated in Figure 2 for $d = 2$. Note that Figure 2 is inspired from illustrations in [18, 41].

Note that the function $\tilde{f}_{n,h_n,i}$ is supported on $[0, 1]^d$ and $\int_{[0,1]^d} \tilde{f}_{n,h_n,i}(x)dx = 1$. Moreover, it satisfies the two following lemmas.

**Lemma 3.5.** Let $h_n \in (0, 1/2)$. Under Assumptions the (A1), (A'4), and (A'6),

$$ \forall i \in \{1, \dots, d\}, \quad \|\mathbb{E}[\tilde{f}_{n,h_n,i}] - f_X\|_{\infty} = \|\mathbb{E}[\tilde{f}_{n,h_n,1}] - f_X\|_{\infty} = O(h_n^{\alpha}). \quad (20) $$
# 5 Bandwidth selection

Recall that $X = (V_i)_{1\le i\le d}$, $W = (V_i)_{d+1\le i\le p}$ as done in Section 2 so that $Y = g(X, W)$ as in (1). In the following, we propose a bandwidth selection procedure inspired from [12] to estimate $\mathbb{E}[\mathbb{E}[Y|X]^2] = \mathbb{E}[g_1(X)^2]$ with $g_1(x) = \mathbb{E}[Y|X = x]$ by

$$T_{n,h_n} = \frac{2}{n(n-1)} \sum_{1 \le j < j' \le n} \frac{Y_j Y_{j'}}{2} \left( \frac{K_{h_n} \circ A_{X_j}(X_{j'} - X_j)}{f_X(X_j)} + \frac{K_{h_n} \circ A_{X_{j'}}(X_j - X_{j'})}{f_X(X_{j'})} \right)$$

using the sample $((V_{1,j}, \dots, V_{p,j}), Y_j)_{1 \le j \le n}$ where $X_j = (V_{i,j})_{1 \le i \le d}$.

Based on the whole sample $((V_{1,j}, \dots, V_{p,j}), Y_j)_{1 \le j \le n}$, we first build a kernel estimator of the function $g$ defined as:

$$\tilde{g}(v_1, \dots, v_p) = \frac{1}{n} \frac{1}{f_V(v_1, \dots, v_p)} \sum_{j=1}^{n} Y_j \prod_{i=1}^{p} \tilde{K}_{h_{0,i}}(V_{i,j} - v_i), \quad (15)$$

where as usual, $\tilde{K}_h(\cdot) = \tilde{K}(\cdot)/h$, where $\tilde{K}$ is a kernel defined on $\mathbb{R}$ and where the bandwidths $h_{0,i}$ for $i = 1, \dots, d$ are chosen using a classical rule of thumb. The kernel $\tilde{K}$ is different from $K$ and can be chosen arbitrarily. In the following, we will consider a Gaussian kernel. This kernel estimator $\tilde{g}$ will serve as a test function which is expected to be close to $g$, on which the bandwidths of our estimator will be calibrated with a simple procedure. First, we can get an analytical expression of the Sobol' index of $\tilde{g}$, which will be our target. Then, $\tilde{g}$ will be used to generate virtual samples that can be plugged in our estimator. Finally, we will optimize the bandwidths in order to reach the target.

In the sequel the variables $(v_1, \dots, v_p)$ will be replaced by random variables $(\tilde{V}_1, \dots, \tilde{V}_p)$ and we define $\tilde{X} := (\tilde{V}_1, \dots, \tilde{V}_d)$, hence we will have to deal with two sources of randomness one coming from our original sample $((V_{1,j}, \dots, V_{p,j}), Y_j)_{1 \le j \le n}$ and the other one from $(\tilde{V}_1, \dots, \tilde{V}_p)$. We denote by $\tilde{\mathbb{E}}$ for the conditional expectation knowing the first source of randomness (the sequence $((V_{1,j}, \dots, V_{p,j}), Y_j)_{1 \le j \le n}$ is considered as fixed and we integrate with respect to the distribution of $(\tilde{V}_1, \dots, \tilde{V}_p)$).

Then we define an approximation $\tilde{g}_1$ of $g_1$ based on $\tilde{g}$ as:

$$\tilde{g}_1(x) = \tilde{\mathbb{E}}[\tilde{g}(\tilde{V}_1, \dots, \tilde{V}_p) | \tilde{X} = x].$$

Finally, the virtual sample from $\tilde{g}$ is given by $((V_{1,j}, \dots, V_{p,j}), \tilde{Y}_j = \tilde{g}(V_{1,j}, \dots, V_{p,j}))_{1 \le j \le n}$ and is used to build our estimator

$$\tilde{T}_{n,h} = \frac{2}{n(n-1)} \sum_{1 \le j < j' \le n} \frac{\tilde{Y}_j \tilde{Y}_{j'}'}{2} \left( \frac{K_{h_n} \circ A_{X_j}(X_{j'} - X_j)}{f_X(X_j)} + \frac{K_{h_n} \circ A_{X_{j'}}(X_j - X_{j'})}{f_X(X_{j'})} \right).$$

The optimal bandwidth $h^*$ is then given by:

$$h^* = \operatorname*{argmin}_h |\tilde{T}_{n,h} - \tilde{\mathbb{E}}[\tilde{g}_1(\tilde{X})|^2| \quad (16)$$
Lemma 3.5 is an extension of [41, Proposition 5.2] to every $\alpha > 0$ in Assumption (A'4) (itself an extension to every dimension and every $\alpha \in (0, 2]$ of [35, Lemma 3.1]). Its proof is postponed to Appendix A.2.

**Lemma 3.6.** [41, Proposition 5.3] Let $h_n \in (0, 1/2)$ and $\delta_n \in (0, 1)$. Under the Assumptions (A1), (A'6), and (A8) and with the kernel $\tilde{K}$, there exists a constant $C > 0$ such that we have, with probability $1 - \delta_n$,

$$ \forall i \in \{1, \dots, d\}, \quad \|\tilde{f}_{n,h_n,i} - \mathbb{E}[\tilde{f}_{n,h_n,i}]\|_{\infty}^2 \le C \frac{\left(\log\left(\frac{1}{h_n}\right)\right)_+ + \log\left(\frac{2}{\delta_n}\right)}{nh_n^d}. \quad (21) $$

Once more, choosing $h_n = n^{-1/(2\alpha+d)}$ and $\delta_n = 1/n$ in Lemmas 3.5 and 3.6, we recover the optimal nonparametric rate up to the log factor. As before, we can deduce from Lemmas 3.5 and 3.6 the control of the mean integrated squared error of $\tilde{f}_{n,h_n,i}$ and its infimum.

**Corollary 3.7.** Assume that Assumptions (A1), (A2), (A'4), and (A'6) are satisfied, as well as Assumption (A8) with kernel $\tilde{K}$. Then one has, for all $i \in \{1, \dots, d\}$,

$$ \mathbb{E}\left[\int_{[0,1]d} (\tilde{f}_{n,h_n,i}(x) - f_X(x))^2 dx\right] = o(n^{-1/2}) \quad \text{and} \quad \frac{1}{\inf_{x \in [0,1]^d} |\tilde{f}_{n,h_n,i}(x)|} = O_P(1) \quad (22) $$

$$ nh_n^{2d} \to \infty \quad \text{and} \quad nh_n^{4\alpha} \to 0 \quad \text{as} \quad n \to \infty. $$

Assuming $h_n = n^{-\gamma}$ with $\gamma > 0$, the conditions $nh_n^{2d} \to \infty$ and $nh_n^{4\alpha} \to 0$ lead once more to $1/(4\alpha) < \gamma < 1/(2d)$ and $\alpha > d/2$.

The proof of Corollary 3.7 comes from Lemmas 3.5 and 3.6 and is similar to the proof of Corollary 3.4 from Lemmas 3.1 and 3.2. It is thus skipped.

# 4 Central limit theorems

Here, we prove a central limit theorem for both estimators defined in (7) and (19).

**Theorem 4.1** (Central limit theorem). (i) Under Assumptions (A1) to (A8), one has

$$ \sqrt{n}(\hat{T}_{n,h_n} - \mathbb{E}[\mathbb{E}[Y|X]^2]) \xrightarrow{n \to \infty} \mathcal{N}(0, \text{Var}((2Y - m(X))m(X))) \quad (23) $$

as soon as $\alpha > d/2$ and $h_n = n^{-\gamma}$ with $1/(4\alpha) < \gamma < 1/(2d)$.

(ii) Replacing Assumptions (A4) and (A6) respectively by Assumptions (A'4) and (A'6) in (i), the same conclusion holds true for $\hat{T}_{n,h_n}$.

The proof of Theorem 4.1 is postponed to Appendix A.3 while the proofs of the following results are postponed to Appendix A.4.

**Proposition 4.2** (Asymptotic efficiency for $\hat{T}_{n,h_n}$ and $\tilde{T}_{n,h_n}$). Under the assumptions of Theorems 4.1, $\hat{T}_{n,h_n}$ and $\tilde{T}_{n,h_n}$ are asymptotically efficient to estimate $\mathbb{E}[\mathbb{E}[Y|X]^2]$ from an i.i.d. sample $(X_i, Y_i)_{i=1,\dots,n}$ of the pair $(X, Y)$.
**Computation of $\tilde{\mathbb{E}}[\tilde{g}_1(\tilde{X})^2]$** First, one has

$$ \tilde{\mathbb{E}}[\tilde{g}_1(\tilde{X})^2] = \frac{1}{n^2} \sum_{1 \le j \le j' \le n} Y_j Y_{j'} \left( \prod_{i=1}^{d} \tilde{\beta}_{i,j,j'} \right) \left( \prod_{i=d+1}^{p} \tilde{\beta}_{i,j} \tilde{\beta}_{i,j'} \right) \quad (17) $$

where, for $1 \le j, j' \le n$,

$$ \tilde{\beta}_{i,j,j'} := \int_{B_i}^{C_i} \tilde{K}_{h_{0,i}}(V_{i,j} - v) \tilde{K}_{h_{0,i}}(V_{i,j'} - v) \frac{1}{f_{V_i}(v)} dv \quad \text{for } 1 \le i \le d, $$

$$ \tilde{\beta}_{i,j} := \int_{B_i}^{C_i} \tilde{K}_{h_{0,i}}(V_{i,j} - v) dv \quad \text{for } d+1 \le i \le p. $$

The details of the calculus leading to (17) are postponed to Appendix A.2.

Second, the point is to choose a kernel $\tilde{K}$ so that $\tilde{\mathbb{E}}[\tilde{g}_1(\tilde{X})^2]$ has an explicit expression. Note that if we assume that the kernel $\tilde{K}$ is Gaussian and that the inputs are uniformly distributed on $[0, 1]$, this is the case: for $1 \le j, j' \le n$, $\tilde{\beta}_{i,j}$ and $\tilde{\beta}_{i,j'}$ writes

$$ \tilde{\beta}_{i,j,j'} = \frac{1}{\sqrt{2}h_{0,i}} \varphi\left(\frac{V_{i,j} - V_{i,j'}}{\sqrt{2}h_{0,i}}\right) \left[ \Phi_{\frac{V_{i,j}+V_{i,j'}}{2}, \frac{h_{0,i}}{\sqrt{2}}} (1) - \Phi_{\frac{V_{i,j}+V_{i,j'}}{2}, \frac{h_{0,i}}{\sqrt{2}}} (0) \right] \quad \text{for } 1 \le i \le d, $$

$$ \tilde{\beta}_{i,j} = \Phi_{V_{i,j}, h_{0,i}}(1) - \Phi_{V_{i,j}, h_{0,i}}(0) \quad \text{for } d+1 \le i \le p. $$

with $\varphi$ the probability density function of the standard normal distribution and $\Phi_{\mu,\sigma}$ the cumulative distribution function of a normal distribution with mean $\mu$ and standard deviation $\sigma$.

**Algorithm 1 Bandwidth selection**

Enter the total number $n$ of calls to the computer code

Choose the kernel $\tilde{K}$

Choose the bandwidths $(h_{0,i})_{1 \le i \le p}$ using a classical rule of thumb

Simulate a $n$-sample $((V_{1,j}, \dots, V_{p,j}), Y_j)_{1 \le j \le n}$ of input/output

Compute the values of $(\tilde{\beta}_{i,j,j'})_{1 \le i \le d, 1 \le j, j' \le n}$ and $(\tilde{\beta}_{i,j})_{d+1 \le i \le p, 1 \le j \le n}$

Compute the value of $\tilde{\mathbb{E}}[\tilde{g}_1(\tilde{X})^2]$ given in (17)

Compute the $n$ sample $(\tilde{Y}_j = \tilde{g}(V_{1,j}, \dots, V_{p,j}))_{1 \le j \le n}$ where $\tilde{g}$ is defined in (15)

Solve the optimization problem given in (16)

return the optimal value $h^*$

# 6 Extension to an unknown density $f_X$

In this section, we assume that the probability density function $f_X$ is unknown and our aim is to propose a generalization of $T_{n,h_n}$ with a plug-in estimator $\hat{f}_X$:

$$ \hat{T}_{n,h} := \frac{2}{n(n-1)} \sum_{1 \le j < j' \le n} \frac{Y_j Y_{j'}}{2} \left( \frac{K_h \circ A_{X_j}(X_{j'} - X_j)}{\hat{f}_X(X_j)} + \frac{K_h \circ A_{X_{j'}}(X_j - X_{j'})}{\hat{f}_X(X_{j'})} \right). $$
Using the delta method [51, Theorem 3.1], we are now able to get the asymptotic behavior of the estimators of $S^X$ derived respectively from (7) and (19). Let

$$ \hat{S}_{n,h_n} := \frac{\hat{T}_{n,h_n} - \left(\frac{1}{n} \sum_{j=1}^{n} Y_j\right)^2}{\frac{1}{n} \sum_{j=1}^{n} Y_j^2 - \left(\frac{1}{n} \sum_{j=1}^{n} Y_j\right)^2} \quad \text{and} \quad \tilde{S}_{n,h_n} := \frac{\tilde{T}_{n,h_n} - \left(\frac{1}{n} \sum_{j=1}^{n} Y_j\right)^2}{\frac{1}{n} \sum_{j=1}^{n} Y_j^2 - \left(\frac{1}{n} \sum_{j=1}^{n} Y_j\right)^2}. $$

**Corollary 4.3** (Central limit theorem and asymptotic efficiency for $\hat{S}_{n,h_n}$ and $\tilde{S}_{n,h_n}$).

(i) Under the assumptions of Theorem 4.1 for $\hat{T}_{n,h_n}$, one has

$$ \sqrt{n} (\hat{S}_{n,h_n} - S^X) \xrightarrow[n \to \infty]{\mathcal{L}} \mathcal{N}(0, \sigma^2), \qquad (24) $$

where the limiting variance $\sigma^2$ has an explicit expression given in the proof in (38).

(ii) Under the assumptions of Theorem 4.1 for $\hat{T}_{n,h_n}$, the same result holds for $\tilde{S}_{n,h_n}$.

(iii) Moreover, $\hat{S}_{n,h_n}$ and $\tilde{S}_{n,h_n}$ are asymptotically efficient to estimate $S^X$ from an i.i.d. sample $(X_i, Y_i)_{i=1, \dots, n}$ of the pair $(X, Y)$.

Using once more the delta method, we deduce the asymptotic behavior of the vector of the $p$ first-order Sobol' indices. Let us denote by $S^i$ the first-order Sobol' index associated to $X = V_i$ and its estimators $\hat{S}^i$ and $\tilde{S}^i$ given by:

$$ \hat{S}_{n,h_n}^i := \frac{\hat{T}_{n,h_n} - \left(\frac{1}{n} \sum_{j=1}^{n} Y_j\right)^2}{\frac{1}{n} \sum_{j=1}^{n} Y_j^2 - \left(\frac{1}{n} \sum_{j=1}^{n} Y_j\right)^2}, \quad \tilde{S}_{n,h_n}^i := \frac{\tilde{T}_{n,h_n} - \left(\frac{1}{n} \sum_{j=1}^{n} Y_j\right)^2}{\frac{1}{n} \sum_{j=1}^{n} Y_j^2 - \left(\frac{1}{n} \sum_{j=1}^{n} Y_j\right)^2}. $$

Under the assumptions of Theorem 4.1, one may prove that $(\hat{S}_{n,h_n}^1, \dots, \hat{S}_{n,h_n}^p)$ and $(\tilde{S}_{n,h_n}^1, \dots, \tilde{S}_{n,h_n}^p)$ satisfy a central limit theorem with limiting variances that can be computed explicitly performing classical computations in statistics. Being useless in the rest of the paper since we do not perform statistical testing, such explicit expressions are not provided here. Moreover, $(\hat{S}_{n,h_n}^1, \dots, \hat{S}_{n,h_n}^p)$ and $(\tilde{S}_{n,h_n}^1, \dots, \tilde{S}_{n,h_n}^p)$ are asymptotically efficient to estimate $(S^1, \dots, S^p)$ from an i.i.d. sample $(\tilde{X}_i, Y_i)_{i=1, \dots, n}$ of the pair $(X, Y)$.

The proof follows from [51, Theorem 25.50] and does not rely on the fact that we are dealing with first-order Sobol' indices. Hence, if $u_1, \dots, u_r$ are distinct subsets of $\{1, \dots, p\}$, we also have under the same assumptions that $(\hat{S}_{n,h_n}^{u_1}, \dots, \hat{S}_{n,h_n}^{u_r})$ and $(\tilde{S}_{n,h_n}^{u_1}, \dots, \tilde{S}_{n,h_n}^{u_r})$ satisfy a central limit theorem with limiting variances that can be computed explicitly.

**Remark 4.4.** Assumptions (A4), (A'4) and (A5) are isotropic regularity assumptions. The results of this paper could be easily extended to anisotropic classes of regularity at the cost of more complicated notation. Then the bandwidth $h_n$ in Assumptions (A6), (A'6), (A7) and (A8) would be replaced by $d$ distinct bandwidths $h_{n,k}$, one for each dimension of $X$. In the numerical Section 6, we therefore select one bandwidth per input dimension, using leave-one-out on the regression function.

# 5 Comparison with existing methods

In this section, we compare several estimators of $\mathbb{E}[\mathbb{E}[Y|X]^2]$ proposed in the literature in terms of asymptotic variance in the central limit theorem. All of them rely on the same
**Corollary 6.1.** Assume that the assumptions of Theorem 3.1 are satisfied. Furthermore, assume that

$$ \mathbb{E}\left[\left(\frac{f_X(X_1) - \hat{f}_X(X_1)}{\hat{f}_X(X_1)}\right)^2\right] = o\left(\frac{h_n^d}{n}\right). \qquad (18) $$

Then we get

$$ \sqrt{n}\left(\hat{T}_{n,h_n} - \mathbb{E}[\mathbb{E}[Y|X]^2]\right) \xrightarrow[n \to \infty]{\mathcal{D}} \mathcal{N}(0, 4\tau^2) \qquad (19) $$

where $\tau^2$ has been defined in Theorem 3.1.

The proof of Corollary 6.1 has been postponed in Appendix A.2.

**Parametric estimation of $f_X$** Here we assume that the unknown probability density function $f_X$ belongs to a parametric family $\mathcal{F}$:

$$ \mathcal{F} := \{f_{\theta}, \theta \in \Theta \subset \mathbb{R}^m\} $$

where $f_{\theta}$ is a probability density function supported by $\mathcal{D}$ for any $\theta \in \Theta$. We denote by $\theta_0$ the true parameter: $f_X = f_{\theta_0}$. Let $\hat{\theta}$ be a consistent estimator of $\theta_0$ and estimate $f_X$ by $f_{\hat{\theta}}$.

This setting is natural. Indeed, in several industrial applications, the input variables are assumed to follow a certain distribution, the parameters of which are given by expert judgment. For instance, in [38], the authors model the fuel mass $M_{fuel}$ needed to link two fixed countries with a commercial aircraft by the Bréguet formula:

$$ M_{fuel} = \left( M_{empty} + M_{pload} \right) \left( e^{\frac{SFC \cdot g \cdot Ra}{V \cdot F}} 10^{-3} - 1 \right) \qquad (20) $$

where the cruise speed $V$ and the lift-to-drag ratio $F$ are assumed to be uniformly distributed and beta-distributed respectively with unknown parameters. See [38] for the description of the model with more details. Another example is given by the Gaussian Plume Model (GPM) that concerns a point source emitting contaminant into a unidirectional wind in an infinite domain. Such a model is also applied, for instance, to volcanic eruptions, pollen and insect dispersals (see, e.g., [5, 44]). More precisely, the contaminant concentration at location $(x, y, 0)$ (ground floor) involves several parameters and writes as:

$$ C(x, y, 0) = \frac{Q}{2\pi K x} e^{\frac{-u(y^2 + H^2)}{4Kx}}. \qquad (21) $$

Based on reality constraints and guided by the expert knowledge, the stochastic parameters $Q$, $K$, and $u$ of the model are assumed to be independent with uniform distribution on $[0, M_Q]$, $[0, M_K]$, and $[0, M_u]$ respectively.

**Examples 6.2.**
given-data design of experiment as in the setting of our work, except the Pick-Freeze estimator.

Recall also that our estimator is asymptotically efficient and its limiting variance is given in Theorem 4.1 by

$$ \mathrm{Var}((2Y - m(X))m(X)) = \mathrm{Var}(m^2(X)) + 4\mathbb{E}[m^2(X)\sigma^2(X)]. \quad (25) $$

**A seminal kernel-based method** As mentioned in Section 3, the authors of [15] consider a truncated version of $T$ defined as $T^{\mathrm{trunc},\varepsilon} = \mathbb{E}[\mathbb{E}[Y|X]^2 \mathbb{1}_{X \in (\varepsilon, 1-\varepsilon)^d}]$ to bypass boundary issues, with limiting variance

$$ \mathrm{Var}(\mathbb{1}_{X \in (\varepsilon, 1-\varepsilon)^d} m^2(X)) + 4\mathbb{E}[\mathbb{1}_{X \in (\varepsilon, 1-\varepsilon)^d} m^2(X) \sigma^2(X)] $$

As mentioned before, notice that the fine tuning of $\varepsilon$ is cumbersome. Moreover, the estimated quantity is not exactly $T$ but only its pruned version.

**An alternative kernel-based method in dimension one** An alternative procedure has been introduced in [9, page 11]. Anyway, note that the estimator $\hat{T}_n$ defined in [9, page 11] is not easily tractable in practice. More precisely, the initial $n$-sample is split into two samples of sizes $n_1 = \lfloor n/\log n \rfloor$ and $n_2 = n - n_1 \approx n$. The first sample is used to estimate the joint density of $(X, Y)$ while the second one is used in an outer loop to estimate the integral term by Monte-Carlo. The sequence of estimators $(\hat{T}_n)_{n \ge 1}$ is proved to be asymptotically efficient to estimate $\mathbb{E}[\mathbb{E}[Y|X]^2]$ and its limiting variance is naturally given by (25) as expected. We refer the reader to [9, Theorems 3.4 and 3.5] and [16].

**Nearest neighbor-based method** One may also compare our results to the estimation procedure proposed in [13], based on nearest neighbors. Here again, the initial $2n$-sample is split into two samples of equal size $n$. The first sample allows to estimate the regression function $m(\cdot) = \mathbb{E}[Y|X = \cdot]$ using the first nearest neighbor of $x$ among the points of the first sample while the second sample is used as a plug-in estimator. They proved that their estimator $S_n$ is consistent for any dimension $d$ of $X$ and that $\sqrt{n}(S_n - \mathbb{E}[S_n])$ is asymptotically Gaussian. Nevertheless, the bias term is negligible only if $d \le 3$ and in the setting where $f_X$ is Lipschitz continuous and bounded away from zero. Thus it may happen that the central limit theorem does not hold true for $\sqrt{n}(S_n - \mathbb{E}[Y|X]^2)$ if $d \ge 4$. The limiting variance obtained with this estimation procedure is

$$ 2\mathrm{Var}(m^2(X)) + 5\mathbb{E}[m^2(X)\sigma^2(X)] + 2\mathbb{E}[\sigma^2(X)], $$

the multiplicative factor 2 in front of $\mathrm{Var}(m^2(X))$ taking into account that we have considered two samples of size $n$. When $d = 1$, one may also use the estimation procedure based on ranks introduced in [7] and studied in [16] with asymptotic variance

$$ \mathrm{Var}(m^2(X)) + 4\mathbb{E}[m^2(X)\sigma^2(X)] + \mathbb{E}[\sigma^2(X)]. $$

Note that this asymptotic variance is not the efficient one, but this methodology has been proven to perform numerically particularly well.
1. Uniform distribution. For sake of clarity, we only detail the framework of first-order Sobol' indices, that is $X = V_i$ for $i \in \{1, \dots, p\}$, and we assume that $X$ is uniformly distributed on $[0, \theta]$ where $\theta$ is unknown. We estimate $\theta$ by

$$ \hat{\theta} = \max(X_1, \dots, X_n). $$

It is a well known fact that $n(\theta - \hat{\theta})$ converges in distribution towards an exponential distribution of parameter $1/\theta$. Then one has

$$ \mathbb{E}\left[\left(\frac{f_X(X_j) - \hat{f}_X(X_j)}{\hat{f}_X(X_j)}\right)^2\right] = \mathbb{E}\left[\left(\frac{\theta^{-1} - \hat{\theta}^{-1}}{\hat{\theta}^{-1}}\right)^2\right] = \frac{1}{\theta^2}\mathbb{E}[(\hat{\theta} - \theta)^2] = \frac{2}{(n+1)(n+2)}. $$

So that Condition (18) is fulfilled since $nh_n \to \infty$ ($d=1$). Note that these computations can be easily generalized to the estimation of Sobol' indices of any order with uniformly distributed inputs.

2. Beta distribution with parameters a and b in $(1, 3/2)$ unknown. Once again, we only detail the framework of first-order Sobol' index estimation. Assume thus that $X = V_i$, for some $i \in \{1, \dots, p\}$ and assume that $X$ is beta-distribution with a first unknown parameter a and a second known parameter b both in $(1, 3/2)$. Let $\hat{a} = b\bar{X}_m/(1 - \bar{X}_m)$ be the moment estimator of a, built from a m-sample $(\tilde{X}_1, \dots, \tilde{X}_m)$ of X, independent from the original n-sample. Here, $\bar{X}_m$ stands for the empirical mean of the m-sample. Since the model is not evaluated on this new sample, even for costly numerical models, the estimation cost is still the one of n model evaluations.

Now let

$$ \hat{f}_{a,b}(x) = \begin{cases} f_{\hat{a},b}(x) & \text{if } \hat{a} \in [1, 3/2] \\ \eta & \text{otherwise,} \end{cases} $$

$$ \text{where } \eta = (2 \int_0^1 \sqrt{x(1-x)} dx)^{-1}. $$

Firstly, on the event $\{\hat{a} \notin [1, 3/2]\}$, one may check that $\|f_{a,b}\|_{\infty} \le 2\eta$ so that

$$ \left| \frac{f_{a,b}(X_j) - f_{\hat{a},b}(X_j)}{f_{\hat{a},b}(X_j)} \right| = \left| \frac{f_{a,b}(X_j) - \eta}{\eta} \right| \le 1. $$

Furthermore, let us denote by $\Lambda$ the log Laplace transform of X and $\Lambda^*$ its Fenchel-Legendre transform:

$$ \Lambda(\lambda) = \log \mathbb{E}[e^{\lambda X}] \quad \text{and} \quad \Lambda^*(s) = \sup_s \{\lambda s - \Lambda(\lambda)\}. $$

See e.g. [13] for details on these functions and their properties. Now, for all $\delta > 0$, it yields, for m large enough,

$$ \begin{aligned} \mathbb{P}(\hat{a} \notin [1, 3/2]) &= \mathbb{P}(\bar{X}_m \le \frac{1}{1+b}) + \mathbb{P}(\bar{X}_m \ge \frac{3/2}{3/2+b}) \\ &= \mathbb{P}(\bar{X}_m - \mathbb{E}[X] \le \frac{1}{1+b} - \frac{a}{a+b}) + \mathbb{P}(\bar{X}_m - \mathbb{E}[X] \ge \frac{3/2}{3/2+b} - \frac{a}{a+b}) \\ &\le 2 \exp\left\{-m \min\left\{\Lambda^*\left(\frac{1}{1+b} - \frac{a}{a+b}\right), \Lambda^*\left(\frac{3/2}{3/2+b} - \frac{a}{a+b}\right)\right\} - \delta\right\} \end{aligned} \quad (22) $$
**Pick-Freeze method** The limiting variance involved in the central limit theorem of the Pick-Freeze estimation, based on the particular Pick-Freeze design, is given by

$$ \mathrm{Var}(m^2(X)) + 2\mathbb{E}[m^2(X)\sigma^2(X)] + \mathbb{E}[\sigma^2(X)]. $$

Although the Pick-Freeze estimator and ours are both asymptotically efficient, they rely on different types of design of experiments. So in some sense, the limiting variances are not comparable. Nevertheless, for a fixed sample size, one can compare numerically the numerical accuracy of these estimation methods through the variability of the estimators (see the boxplots in Section 6).

# 6 Numerical applications

In this section, we illustrate the practical performances of one of our estimators on two analytical test cases and one realistic flood case coming from the sensitivity analysis literature. Recall that both estimators are based on a high-order kernel supported on $[0, 1]^d$ for the first estimator and on $[-1, 1]^d$ for the second one. From a theoretical perspective both of them lead to equivalent convergence results, but the first one actually suffers from strong numerical instabilities as discussed in the additional file that can be found on the following link https://hal.science/hal-04052837. This explains why we focus on the second one in what follows. For all test cases:

* We focus on the estimation of first-order and total Sobol' indices. Note that the estimation of the total Sobol' index associated to input $V_i$ relies on the estimation of the closed Sobol' index of order $p-1$ associated to $X = V_{\sim i}$ (see Equation (2) in Section 2.1).

* We compute our second mirror-type estimator (19) with an Epanechnikov kernel of order 2 and 4 (see [24] for a definition), with the kernel bandwidth being optimized via leave-one-out on the regression function ("Kernel 2" and "Kernel 4").

* We also consider concurrent estimators, namely the nearest-neighbor estimator of [13] ("NN") and the asymptotically efficient version of Pick-Freeze estimator studied in [28] ("PF1") for first-order indices and for total indices, and also the replicated version of Pick-Freeze estimator proposed in [49] ("PF2"), the rank estimator of [16] ("Rank"), and the lag estimator of [32] ("Lag") for first-order indices.

* For all estimators, we generate a standard $n$-sample $(X_1, Y_1), \dots, (X_n, Y_n)$ except for the Pick-Freeze method since it relies on a structured design of experiments.

* Each experiment is repeated 100 times with a number of model evaluations fixed to $n = 500$ (then to $n = 1000$). The reference value that approximates the true value of the index is obtained from a Pick-Freeze estimation with very large sample size ($d \times 5.10^4$).

**The Bratley function** First, we consider the Bratley function defined by:

$$ g_{\text{Bratley}}(V_1, \dots, V_p) = \sum_{i=1}^{p} (-1)^i \prod_{j=1}^{i} V_j, \qquad (26) $$
since $\mathbb{E}[X] = a/(a+b)$ and using Cramér's theorem [13, Theorem 2.2.3] as $\Lambda < \infty$.
Secondly, on the event $\{\hat{a} \in [1, 3/2]\}$, we have $\hat{a} > 0$, $2a - \hat{a} > 0$, and $3a - \hat{a} > 0$ and we get

$$ \mathbb{E}\left[\left(\frac{f_{a,b}(X_j) - f_{\hat{a},b}(X_j)}{f_{\hat{a},b}(X_j)}\right)^2 | \hat{a}\right] = \frac{B(\hat{a}, b)^2}{B(a, b)^3} B(3a - 2\hat{a}, b) - 2\frac{B(\hat{a}, b)}{B(a, b)^2} B(2a - \hat{a}, b) + 1. $$

Observe that the previous expression is a function of $\bar{X}_m$ that vanishes at $\mathbb{E}[X]$.
Moreover $\{\hat{a} \in [1, 3/2]\} = \{\bar{X}_m \in [1/(1+b), 3/(3+2b)]\}$ as $\hat{a} = b\bar{X}_m/(1-\bar{X}_m)$.
Thus, the function being infinitely differentiable, by applying the mean value, inequality we obtain

$$ \begin{aligned} \mathbb{E}\left[\left(\frac{f_{a,b}(X_j) - f_{\hat{a},b}(X_j)}{f_{\hat{a},b}(X_j)}\right)^2 \mathbb{1}_{\hat{a} \in [1/2, 3/2]}\right] &= \mathbb{E}\left[\mathbb{E}\left[\left(\frac{f_{a,b}(X_j) - f_{\hat{a},b}(X_j)}{f_{\hat{a},b}(X_j)}\right)^2 | \hat{a}\right] \mathbb{1}_{\hat{a} \in [1/2, 3/2]}\right] \\ &= \mathbb{E}\left[\left(\varphi(\bar{X}_m) - \varphi(\mathbb{E}[X])\right) \mathbb{1}_{\bar{X}_m \in [1/(1+2b), 3/(3+2b)]}\right] \\ &\le C \mathbb{E}[|\bar{X}_m - \mathbb{E}[X]|] \\ &\le C \mathbb{E}[(\bar{X}_m - \mathbb{E}[X])^2]^{1/2} \\ &= \frac{C}{m^{1/2}} \mathrm{Var}(X)^{1/2}. \end{aligned} \quad (23) $$

Finally, in view of (22) and (23), Condition (18) is then fulfilled as soon as $m^{-1/2} = o(h_n n^{-1})$.

**Nonparametric estimation of $f_X$**
Here the unknown density function $f_X$ satisfies the assumptions of Theorem 3.1 and is moreover bounded from below by a constant $\eta > 0$. Then we estimate $f_X(x)$ by $\hat{f}_X(x) \vee \frac{\eta}{2}$ with $\hat{f}_X(x) = \frac{1}{m} \sum_{j=1}^m K_{\hat{h}} \circ A_x(\tilde{X}_j - m)$ the adaptive mirror-type estimator introduced in [3] built from a $m$-sample $(\tilde{X}_1, \dots, \tilde{X}_m)$ independent from the initial $n$-sample $(X_1, \dots, X_n)$. Then

$$ \begin{aligned} \left( \frac{f(X_j) - \hat{f}_X(X_j) \vee \frac{\eta}{2}}{\hat{f}_X(X_j) \vee \frac{\eta}{2}} \right)^2 &\le \frac{(f(X_j) - \hat{f}_X(X_j) \vee \frac{\eta}{2})^2}{\frac{\eta^2}{4}} \\ &\le \frac{4}{\eta^2} \left\{ \left(f(X_j) - \frac{\eta}{2}\right)^2 \mathbb{1}_{\hat{f}_X(X_j) \le \frac{\eta}{2}} + (f(X_j) - \hat{f}_X(X_j))^2 \mathbb{1}_{\hat{f}_X(X_j) > \frac{\eta}{2}} \right\} \\ &\le \frac{4}{\eta^2} \left\{ \left(\|f_X\|_{\infty} - \frac{\eta}{2}\right)^2 \mathbb{1}_{\hat{f}_X(X_j) \le \frac{\eta}{2}} + (f(X_j) - \hat{f}_X(X_j))^2 \mathbb{1}_{\hat{f}_X(X_j) > \frac{\eta}{2}} \right\} \end{aligned} $$

where the last inequality holds true as $f_X$ is continuous and thus bounded on its bounded support. Thus

$$ \mathbb{E}\left[\left(\frac{f(X_j) - \hat{f}_X(X_j) \vee \frac{\eta}{2}}{\hat{f}_X(X_j) \vee \frac{\eta}{2}}\right)^2\right] \le \frac{4}{\eta^2} \left\{ \left(\|f_X\|_{\infty} - \frac{\eta}{2}\right)^2 \mathbb{E}[\mathbb{1}_{\hat{f}_X(X_j) \le \frac{\eta}{2}}] + \mathbb{E}\left[(f(X_j) - \hat{f}_X(X_j))^2\right] \right\}. $$

In addition, we have

$$ \mathbb{1}_{\hat{f}_X(X_j) \le \frac{\eta}{2}} = \mathbb{1}_{f_X(X_j) + \hat{f}_X(X_j) - f_X(X_j) \le \frac{\eta}{2}} = \mathbb{1}_{f_X(X_j) - \frac{\eta}{2} \le f_X(X_j) - \hat{f}_X(X_j)} \le \mathbb{1}_{\frac{\eta}{2} \le f_X(X_j) - \hat{f}_X(X_j)} $$
with $V_i \sim U([0, 1])$ i.i.d. and $p = 5$. The results of the 100 experiments for all estimators of first-order indices are given as boxplots in Figure 3. We observe several trends.

* The nearest-neighbor estimator exhibits both large bias and large variance.

* The Pick-Freeze estimators and the rank one perform well except when the sensitivity index is small.

* Both the lag estimator and our mirror-type one have very small bias and variance.

<figure><img src="image_4.png" /><figcaption>Figure 3: Estimators for first-order indices of the Bratley function with $n = 500$ (left) and $n = 1000$ (right). The reference value is represented with a gray line.</figcaption></figure>

For total indices in Figure 4, recall that the only concurrent estimators are nearest-neighbor and Pick-Freeze estimators. Once again Pick-freeze estimators perform well, but the bias of the nearest-neighbor one is very large (here $d = p - 1 = 4$), this bias increasing dramatically when the sensitivity index is small. Our mirror-type estimator still has small bias and very small variance for all input variables. Experiments with lower sample size are given in Appendix B.

<figure><img src="image_5.png" /><figcaption>Figure 4: Estimators for total indices of the Bratley function with $n = 500$ (left) and $n = 1000$ (right). The reference value of the index is represented with a gray line.</figcaption></figure>

**The g-Sobol function** Then, we investigate the g-Sobol function defined by

$$g_{\text{g-Sobol}}(V_1, \dots, V_p) = \prod_{i=1}^{p} \frac{|4V_i - 2| + a_i}{1 + a_i}, \qquad (27)$$

with $V_i \sim U([0, 1])$ i.i.d., $p = 5$ and $a = (0, 1, 4.5, 9, 99)$. Notice that it is non-differentiable at any input value with a component equal to 0.5, but the impact on our estimator
yielding, using Markov's inequality,

$$
\begin{align*}
\mathbb{E}[1_{\hat{f}_X(X_j) \le \frac{\eta}{2}}] &\le \mathbb{E}\left[\mathbb{P}\left(f_X(X_j) - \hat{f}_X(X_j) \ge \frac{\eta}{2} | \widetilde{X}_1, \dots, \widetilde{X}_m\right)\right] \\
&\le \frac{4}{\eta^2} \mathbb{E}\left[\mathbb{E}\left[\left(f_X(X_j) - \hat{f}_X(X_j)\right)^2 | \widetilde{X}_1, \dots, \widetilde{X}_m\right]\right] \\
&\le \frac{4}{\eta^2} \mathbb{E}\left[\int_{\mathcal{D}} (f_X(x) - \hat{f}_X(x))^2 f_X(x) dx\right] \\
&\le \frac{4}{\eta^2} \|f_X\|_{\infty} \cdot \mathbb{E}\left[\int_{\mathcal{D}} (f_X(x) - \hat{f}_X(x))^2 dx\right].
\end{align*}
$$

In the same, way we get

$$
\mathbb{E}[(f(X_j) - \hat{f}_X(X_j))^2] \le \|f_X\|_{\infty} \mathbb{E}[\int_{\mathcal{D}} (f_X(x) - \hat{f}_X(x))^2 dx].
$$

From [3, Theorem 4] and as $f_X \in C^\alpha(\mathcal{D})$, it comes

$$
\mathbb{E} \left[ \int_{\mathcal{D}} (f_X(x) - \hat{f}_X(x))^2 dx \right] \le C m^{\frac{-2\alpha}{2\alpha+d}}.
$$

Thus, choosing the sample size $m$ such that

$$
m^{\frac{-2\alpha}{2\alpha+d}} = o\left(\frac{h_n^d}{n}\right),
$$

the central limit theorem stated in Corollary 6.1 holds true.

7 Numerical applications

A Proof of the results

In the following, $h_n$ and $L_n$ are simply denoted by $h$ and $L$ respectively. In addition,
recall that, we write $C$ for a deterministic and finite constant which value is allowed to
change between occurrences.

A.1 Proof of Theorem 3.1

*Proof of (7).* First, since $Y \in L^2(\mathbb{R})$, one has

$$
\begin{align*}
\mathbb{E}[T_{n,h}] &= \frac{1}{2} \iint_{\mathcal{D}^2} \left( \frac{K_h \circ A_{x_1}(x_2 - x_1)}{f_X(x_1)} + \frac{K_h \circ A_{x_2}(x_1 - x_2)}{f_X(x_2)} \right) f_X(x_1) f_X(x_2) g_1(x_1) g_1(x_2) dx_1 dx_2 \\
&= \iint_{\mathcal{D}^2} \frac{K_h \circ A_{x_1}(x_2 - x_1)}{f_X(x_1)} f_X(x_1) f_X(x_2) g_1(x_1) g_1(x_2) dx_1 dx_2 \\
&= \int_{\mathcal{D}} \int_{\mathcal{D}_x} K(u) g_1(x) g_1(x + h A_x^{-1}(u)) f_X(x + h A_x^{-1}(u)) du dx
\end{align*}
$$

where the last display is obtained after the variable change $(x, u) = (x_1, A_{x_1}(x_2 - x_1)/h)$,
using the fact that $|\det(A_x)| = 1$ and where $\mathcal{D} \times \mathcal{D}_x$ is the integration domain after the
performance is negligible, as can be seen in Figure 5 for first-order indices. Except for the degraded performance of the lag estimator, the conclusions are the same as for the Bratley function, even for total indices displayed in Figure 6. Experiments with lower sample size are given in Appendix B.

<figure><img src="image_6.png" /><figcaption>Figure 5: Estimators for first-order indices of the g-Sobol function with $n = 500$ (left) and $n = 1000$ (right). The reference value of the index is represented with a gray line.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Estimators for total indices of the g-Sobol function with $n = 500$ (left) and $n = 1000$ (right). The reference value of the index is represented with a gray line.</figcaption></figure>

From Theorem 4.1, we know that it is necessary to choose a kernel of order $[\alpha] + 1$ with $\alpha > d/2$ for the central limit theorem to hold. Here as $p = 5$, the dimension $d$ of $X$ for total Sobol' indices is equal to $p - 1 = 4$, thus the order of the kernel should be chosen at least equal to 3. However we do not observe on Figure 6 better practical performance for a kernel of order 4 than for a one of order 2. Although the central limit theorem requires such a strong assumption about the order of the kernel, it is important to note that this assumption is not necessary for consistency results. Also we selected the bandwidth by cross-validation, that is with a data-driven approach, which naturally adapts to the order of the kernel. Currently, we have no clear idea of when it is beneficial to use a higher-order kernel in practice. It would be interesting to propose a strategy to choose adaptively and simultaneously the bandwidth and the order of the kernel as in [3]. However this is out of the scope of the present work.
variable change. Now observe that, for h small enough, $\mathcal{D} \subset \mathcal{D}_x$ by assumption and since $K(u) = 0$ for all $u \notin \mathcal{D}$, the integration domain of the last integral is reduced to $\mathcal{D}^2$. In addition, since $\int_{\mathcal{D}} K(u)du = 1$ and $g_1(x) = \mathbb{E}[Y|X = x]$,

$$ \mathbb{E}[\mathbb{E}[Y|X]^2] = \mathbb{E}[g_1^2(X)] = \int_{\mathcal{D}} g_1(x)^2 f_X(x) dx = \iint_{\mathcal{D}^2} g_1(x)^2 f_X(x) K(u) dx du $$

leading to

$$ \begin{aligned} \mathbb{E}[T_{n,h}] - \mathbb{E}[\mathbb{E}[Y|X]^2] &= \iint_{\mathcal{D}^2} K(u)g_1(x) \left( g_1(x + hA_x^{-1}(u))f_X(x + hA_x^{-1}(u)) - g_1(x)f_X(x) \right) dxdu. \end{aligned} \quad (24) $$

The term in (24) can be handled with a Taylor expansion of $g_1 f_X$ (see, e.g., [7, Theorem 5.4]). The assumptions we use are (mainly) $g_1 f_X \in \mathcal{C}^\alpha(\mathcal{D})$ and, for any $1 \le \beta \le \lfloor \alpha \rfloor$,

$$ \int_{\mathcal{D}} \left| g_1(x) \frac{\partial^\beta (g_1 f_X)}{\partial x^\beta} \right| dx < \infty. $$

Then we get a bound in $\mathcal{O}(h^\alpha)$. More precisely, by a multivariate Taylor expansion, since $g_1 f_X$ is $\mathcal{C}^\alpha(\mathcal{D})$, we get:

$$ \begin{aligned} & g_1(x + uh)f_X(x + uh) - g_1(x)f_X(x) \\ &= \sum_{1 \le |\beta| < \lfloor \alpha \rfloor} \frac{h^{|\beta|}}{\beta!} u^\beta \frac{\partial^\beta (g_1 f_X)}{\partial x^\beta}(x) + h^{\lfloor \alpha \rfloor} \sum_{|\beta| = \lfloor \alpha \rfloor} R_\beta(x + uh)u^\beta \end{aligned} \quad (25) $$

with

$$ R_\beta(x + uh) = \frac{\lfloor \alpha \rfloor}{\beta!} \int_0^1 (1 - t)^{\lfloor \alpha \rfloor - 1} \frac{\partial^\beta (g_1 f_X)}{\partial x^\beta}(x + tuh) dt. $$

Then, one has

$$ \begin{aligned} \mathbb{E}[T_{n,h}] - \mathbb{E}[\mathbb{E}[Y|X]^2] &= \sum_{1 \le |\beta| < \lfloor \alpha \rfloor} \frac{h^{|\beta|}}{\beta!} \left( \int_{\mathcal{D}} (A_x^{-1}(u))^\beta K(u) du \right) \left( \int_{\mathcal{D}} g_1(x) \frac{\partial^\beta (g_1 f_X)}{\partial x^\beta}(x) dx \right) \\ &\quad + h^{\lfloor \alpha \rfloor} \sum_{|\beta| = \lfloor \alpha \rfloor} \int_{\mathcal{D}} g_1(x) \int_{\mathcal{D}} u^\beta K(u) R_\beta(x + hA_x^{-1}(u)) du dx \\ &= \lfloor \alpha \rfloor h^{\lfloor \alpha \rfloor} \sum_{|\beta| = \lfloor \alpha \rfloor} \frac{1}{\beta!} \int_{\mathcal{D}} g_1(x) \int_{\mathcal{D}} (A_x^{-1}(u))^\beta K(u) \int_0^1 (1-t)^{\lfloor \alpha \rfloor - 1} \frac{\partial^\beta (g_1 f_X)}{\partial x^\beta}(x + t h A_x^{-1}(u)) dt du dx \end{aligned} $$

using the fact that $(A_x^{-1}(u))^\beta$ is polynomial in $(u_1, \dots, u_d)$ of degree $\beta$ and $K$ is of order $\lfloor \alpha \rfloor$ and that

$$ \int_{\mathcal{D}} \left| g_1(x) \frac{\partial^\beta (g_1 f_X)}{\partial x^\beta} \right| dx < \infty $$

for any $\beta$ such that $1 \le |\beta| < \lfloor \alpha \rfloor$. Now, since $K$ is of order $\lfloor \alpha \rfloor$ and

$$ \int_{\mathcal{D}} u^\beta K(u) \frac{\partial^\beta (g_1 f_X)}{\partial x^\beta}(x) du = \frac{\partial^\beta (g_1 f_X)}{\partial x^\beta}(x) \int_{\mathcal{D}} u^\beta K(u) du = 0 $$
**A realistic flood model** We now consider a realistic test case from the sensitivity analysis literature. The flood model used is a simplification of the 1D Saint-Venant hydrodynamic equations, assuming constant and uniform flows and very wide rectangular cross-sections. It consists of an equation involving the characteristics of the river section upstream of the industrial site:

$$S = Z_v + H - H_d - C_b$$

where $H$ is calculated as:

$$H = \left( \frac{Q}{B K_s \sqrt{\frac{Z_m - Z_v}{L}}} \right)^{0.6}.$$

The model parameters are:

* $S$: overflow (in meters), model output;

* $H$: maximum annual water level (in meters);

* $Q$: maximum annual flow rate (in m³/s), random input following a Gumbel max distribution $Gu(1013, 558)$, truncated below at 500 and above at 3000;

* $K_s$: strickler coefficient, random input following a normal distribution $N(30, 8)$, truncated below at 15;

* $Z_v$: downstream riverbed elevation (in meters), random input following a triangular distribution $T(49, 50, 51)$;

* $Z_m$: upstream riverbed elevation (in meters), random input following a triangular distribution $T(54, 55, 56)$;

* $H_d$: dike height (in meters), random input following a uniform distribution $U(7, 9)$;

* $C_b$: bank elevation (in meters), random input following a triangular distribution $T(55, 55.5, 56)$;

* $L$: length of the river section (in meters), random input following a triangular distribution $T(4990, 5000, 5010)$;

* $B$: river width (in meters), random input following a triangular distribution $T(295, 300, 305)$.

Among the model's input variables, $H_d$ is considered random because it is studied as a design parameter. The other input variables are random due to their temporal and spatial variability, our lack of knowledge about their values, or imprecision in their estimation.

We give in Figure 7 the estimation of first-order and total Sobol' indices with a sample size $n = 500$. Similarly to the previous analytical test cases, the nearest-neighbor estimator has very large bias and variance, and the Pick-Freeze one performs well. Concerning our kernel method, we present the results obtained with a kernel of order 2 or 4, and it clearly outperforms other estimators for first-order indices. For total ones, a residual bias persists, which should decrease by increasing the sample size. Nevertheless, the ranking
for any $\beta$ such that $|\beta| = \lfloor \alpha \rfloor$, we get

$$ \mathbb{E}[T_{n,h}] - \mathbb{E}[\mathbb{E}[Y|X]^2] = \lfloor \alpha \rfloor h^{\lfloor \alpha \rfloor} \sum_{|\beta|=\lfloor \alpha \rfloor} \frac{1}{\beta!} \int_{\mathcal{D}} g_1(x) \int_{\mathcal{D}} (A_x^{-1}(u))^{\beta} K(u) \\ \int_0^1 (1-t)^{\lfloor \alpha \rfloor - 1} \left( \frac{\partial^{\beta}(g_1 f_X)}{\partial x^{\beta}} (x + t h A_x^{-1}(u)) - \frac{\partial^{\beta}(g_1 f_X)}{\partial x^{\beta}}(x) \right) dt du x. $$

Then using the fact that $g_1 f_X \in \mathcal{C}^\alpha(\mathcal{D})$ and (5), one gets

$$ \left| \frac{\partial^\beta(g_1 f_X)}{\partial x^\beta}(x + t h A_x^{-1}(u)) - \frac{\partial^\beta(g_1 f_X)}{\partial x^\beta}(x) \right| \le C_{g_1 f_X}(ht|A_x^{-1}(u)|)^{\alpha - \lfloor \alpha \rfloor} \quad (26) $$

for all $u$ and $x \in \mathcal{D}$, for all $t \in [0, 1]$, for all $h \in (0, \infty)$, and for any $\beta \in \mathbb{N}^d$ such that $|\beta| = \lfloor \alpha \rfloor$. Then,

$$ \begin{aligned} |\mathbb{E}[T_{n,h}] - \mathbb{E}[\mathbb{E}[Y|X]^2]| &\le C_{g_1 f_X} \lfloor \alpha \rfloor h^\alpha \int_{\mathcal{D}} |g_1(x)| dx \\ &\quad \sum_{|\beta|=\lfloor \alpha \rfloor} \frac{1}{\beta!} \int_{\mathcal{D}} |A_x^{-1}(u)|^\beta |A_x^{-1}(u)|^{\alpha-\lfloor \alpha \rfloor} |K(u)| du \int_0^1 (1-t)^{\lfloor \alpha \rfloor - 1} t^{\alpha-\lfloor \alpha \rfloor} dt \\ &= Ch^\alpha \int_{\mathcal{D}} |g_1(x)| dx \times \left( \sum_{|\beta|=\lfloor \alpha \rfloor} \frac{1}{\beta!} \right) \times \left( \int_{\mathcal{D}} \|u\|_{\infty}^\alpha |K(u)| du \right) \\ &\le Ch^\alpha \end{aligned} $$

since $\int_{\mathcal{D}} |g_1(x)| dx < \infty$ and

$$ |A_x^{-1}(u)|^\beta = |u|^\beta = |u_1|^{\beta_1} \dots |u_d|^{\beta_d} \le \|u\|_{\infty}^{\beta_1 + \dots + \beta_d} = \|u\|_{\infty}^{\beta} = \|u\|_{\infty}^{\lfloor \alpha \rfloor}. \quad \square $$

*Proof of (8).* Recall that, by (6), one has

$$ T_{n,h} - \mathbb{E}[T_{n,h}] = \frac{1}{n} \sum_{j=1}^{n} Z_j + S_1 + S_2 $$

where $Z_j = Y_j g_1(X_j) - \mathbb{E}[\mathbb{E}[Y|X]^2]$ has been introduced in Theorem 3.1.

(1) We consider the second moment of $S_1$. Since the pairs $(X_j, Y_j)$ are i.i.d. distributed
of the inputs in terms of influence is preserved: the most important features are $Q$, $H_d$, and $Z_v$. These results are consistent with intuition, since it is expected that the overflow is mainly driven by the maximum annual flow rate $Q$ and the dike height $H_d$. As for the previous example, the results obtained for the estimation of total Sobol' indices with the kernel of order 4 are less favourable.

<figure><img src="image_8.png" /><figcaption>Figure 7: Estimators for first-order indices (left) and total indices (right) of the flood model with $n = 500$. The reference value of the index is represented with a gray line.</figcaption></figure>
as $(X, Y)$ and the variables $Z_j = 2(Y_j g_1(X_j) - \mathbb{E}[\mathbb{E}[Y|X]^2])$ are centered, one has

$$
\begin{align*}
n \mathbb{E}[S_1^2] &= 4n \mathbb{E} \left[ \left( U_n^{(1)}(\pi_1 R) - \frac{1}{n} \sum_{j=1}^n (Y_j g_1(X_j) - \mathbb{E}[\mathbb{E}[Y|X]^2]) \right)^2 \right] \\
&= 4n \text{Var} \left( U_n^{(1)}(\pi_1 R) - \frac{1}{n} \sum_{j=1}^n (Y_j g_1(X_j) - \mathbb{E}[\mathbb{E}[Y|X]^2]) \right) \\
&= 4n \text{Var} \left( \frac{1}{n} \sum_{j=1}^n \left[ (\pi_1 R) \begin{pmatrix} X_j \\ Y_j \end{pmatrix} - (Y_j g_1(X_j) - \mathbb{E}[\mathbb{E}[Y|X]^2]) \right] \right) \\
&= 4 \text{Var} \left( \pi_1 R \begin{pmatrix} X \\ Y \end{pmatrix} - Y g_1(X) \right) \quad \text{by orthogonality of Hoeffding's decomposition} \\
&= 4 \text{Var} \left( \frac{Y}{2f_X(X)} \int_D (K_h \circ A_X)(u-X)(g_1 f_X)(u) du + \frac{Y}{2} \int_D (K_h \circ A_u)(X-u)g_1(u) du - Y g_1(X) \right) \\
&\le \mathbb{E} \left[ Y^2 \left( \frac{1}{f_X(X)} \int_D (K_h \circ A_X)(u-X)(g_1 f_X)(u) du + \int_D (K_h \circ A_u)(X-u)g_1(u) du - 2 g_1(X) \right)^2 \right] \\
&\le 2 \int_D \left( \int_D (K_h \circ A_x)(u-x)(g_1 f_X)(u) du - g_1(x) f_X(x) \right)^2 \frac{g_2(x)}{f_X(x)} dx \\
&\quad + 2 \int_D \left( \int_D (K_h \circ A_u)(x-u)g_1(u) du - g_1(x) \right)^2 g_2(x) f_X(x) dx \\
&=: 2(T_1 + T_2) \tag{27}
\end{align*}
$$

using $(a+b)^2 \le 2(a^2+b^2)$.

• We proceed following the same lines as in the proof of (7): using the variable change $(x, u) = (x_1, A_{x_1}(x_2 - x_1)/h)$, the fact that $\int_D K(u)du = 1$, and the Taylor expansion of $g_1 f_X$ in (25). Then the first term can be rewritten as

$$
\begin{align*}
T_1 &= \int_D \left( \int_D K(u)g_1(x + hA_x^{-1}(u))f_X(x + hA_x^{-1}(u))du - g_1(x)f_X(x) \right)^2 \frac{g_2(x)}{f_X(x)} dx \\
&= \int_D \left( \int_D K(u)[g_1(x + hA_x^{-1}(u))f_X(x + hA_x^{-1}(u)) - g_1(x)f_X(x)]du \right)^2 \frac{g_2(x)}{f_X(x)} dx \\
&= \int_D \left( \sum_{1 \le |\beta| < \lfloor \alpha \rfloor} \frac{h^{|\beta|}}{\beta!} \left( \int_D (A_x^{-1}(u))^{\beta} K(u)du \right) \frac{\partial^{\beta}(g_1 f_X)}{\partial x^{\beta}}(x) \right. \\
&\qquad \left. + h^{\lfloor \alpha \rfloor} \sum_{|\beta| = \lfloor \alpha \rfloor} \int_D (A_x^{-1}(u))^{\beta} K(u)R_{\beta}(x + hA_x^{-1}(u))du \right)^2 \frac{g_2(x)}{f_X(x)} dx \\
&\le 2 \int_D \left( \sum_{1 \le |\beta| < \lfloor \alpha \rfloor} \frac{h^{|\beta|}}{\beta!} \left( \int_D (A_x^{-1}(u))^{\beta} K(u)du \right) \frac{\partial^{\beta}(g_1 f_X)}{\partial x^{\beta}}(x) \right)^2 \frac{g_2(x)}{f_X(x)} dx \\
&\qquad + 2h^{2\lfloor \alpha \rfloor} \int_D \left( \sum_{|\beta| = \lfloor \alpha \rfloor} \int_D (A_x^{-1}(u))^{\beta} K(u)R_{\beta}(x + hA_x^{-1}(u))du \right)^2 \frac{g_2(x)}{f_X(x)} dx. \tag{28}
\end{align*}
$$

Since $K$ is a kernel of order $\alpha$ and $(A_x^{-1}(u))^\beta$ is a polynomial in $(u_1, \dots, u_d)$ of degree $\beta$, $\int_D (A_x^{-1}(u))^\beta K(u)du = 0$ for any $\beta$ such that $|\beta| < \alpha$ and the first term of the previous
# A Proof of the results

In the following, $h_n$ is simply denoted by $h$ and $C$ is a deterministic and finite constant, the value of which is allowed to change between occurrences. Recall also that $f_X$ is the density of $X$ with respect to the Lebesgue measure. The regression function is denoted by $m: m(x) = \mathbb{E}[Y|X = x]$ and we introduce the function $g$ defined by $g(x) = f_X(x)m(x)$. In addition, $\sigma^2(x)$ stands for $\mathrm{Var}(Y|X = x)$.

## A.1 Proof of the preliminary results of Section 3.1

Now, recall that $\hat{f}_{n,h,i}$ defined in (10) and $\hat{m}_{n,h,i}$ defined in (6) are respectively the leave-one-out estimator of $f_X$ based on $(X_1, \dots, X_{i-1}, \dots, X_{i+1}, \dots, X_n)$ and the estimator of the regression function $m$ based on the $(n-1)$-input/output sample $((X_1, Y_1), \dots, (X_{i-1}, Y_{i-1}), \dots, (X_{i+1}, Y_{i+1}), \dots, (X_n, Y_n))$. Then, the leave-one-out estimator of $g$ is naturally given by $\hat{g}_{n,h,i} = \hat{m}_{n,h,i}\hat{f}_{n,h,i}$.

To lighten notation, we write $m_i$, $\sigma_i^2$, $g_i$, and $f_i$ for $m(X_i)$, $\sigma^2(X_i)$, $g(X_i)$, and $f_X(X_i)$ respectively. Additionally, $\hat{m}_i$, $\hat{g}_i$, and $\hat{f}_i$ denote $\hat{m}_{n,h,i}(X_i)$, $\hat{g}_{n,h,i}(X_i)$, and $\hat{f}_{n,h,i}(X_i)$ respectively. By abuse of notation, we denote $\hat{f}_{n,h,i}(x)$ and $\hat{g}_{n,h,i}(x)$ by $\hat{f}_i(x)$ and $\hat{g}_i(x)$.

*Proof of Lemma 3.1.* The result can be deduced from the proof of the bound (8) in [3, Proposition 1] (see [3, Section 7.2.]) as the domain $[0, 1]^d$ is compact. More precisely,

$$ \mathbb{E}[\hat{f}_{n,h,i}(x)] = \int_{[0,1]^d} K_h \circ A_x(x' - x) f_X(x') dx' = \int_{\mathcal{D}_x} K(v) f_X(x + hA_x^{-1}(v)) dv $$

where the last display is obtained after the variable change $v = A_x(x' - x)/h$, using the fact that $|\det(A_x)| = 1$ and where $\mathcal{D}_x = A_x([0, 1]^d - x)/h$ is the integration domain after the variable change. Now observe that the mirror property in (9) ensures that, if $0 < h \le 1/2$ for all $i \in \{1, \dots, d\}$, then $[0, 1]^d \subset \mathcal{D}_x$ for all $x \in [0, 1]^d$. Thus we have

$$ \mathbb{E}[\hat{f}_{n,h,i}(x)] - f_X(x) = \int_{[0,1]^d} K(v)(f_X(x + hA_x^{-1}(v)) - f_X(x))dv \quad (28) $$

since $\int_{[0,1]^d} K(v)dv = 1$. If $\alpha \in (0, 1)$, since $f_X \in C^\alpha([0, 1]^d)$, $\|K\|_\infty < \infty$ and $\|A_x^{-1}(v)\|_\infty \le 1$, we conclude straightforwardly that

$$ |\mathbb{E}[\hat{f}_{n,h,i}(x)] - f_X(x)| \le C_{f_X} h^\alpha \int_{[0,1]^d} \|A_x^{-1}(v)\|_\infty^\alpha |K(v)| dv \le Ch^\alpha. $$

Now, for $\alpha \ge 1$, observe that, still due to the mirror property in (9), $0 < h \le 1/2$ for all $i \in \{1, \dots, d\}$ ensures that $x + hA_x^{-1}(v) \in [0, 1]^d$ for all $x \in [0, 1]^d$ and for all $v \in [0, 1]^d$. Thus the term in (28) can be handled with a Taylor expansion of $f_X$ (see, e.g., [8, Theorem 5.4]). More precisely, by a multivariate Taylor expansion, since $f_X$ is $C^\alpha([0, 1]^d)$, we get:

$$ f_X(x + vh) - f_X(x) = \sum_{1 \le |\beta| < |\alpha|} \frac{h^{|\beta|}}{\beta!} v^\beta \frac{\partial^\beta f_X}{\partial x^\beta}(x) + h^{|\alpha|} \sum_{|\beta| = |\alpha|} R_\beta(x + vh)v^\beta $$
inequality equals to

$$ \sum_{1 \le |\beta|, |\beta'| < \lfloor \alpha \rfloor} \frac{h^{|\beta|}}{\beta!} \frac{h^{|\beta'|}}{\beta'!} \left( \int_{\mathcal{D}} (A_x^{-1}(u))^\beta K(u) du \right) \left( \int_{\mathcal{D}} (A_x^{-1}(u'))^{\beta'} K(u') du' \right) \\ \times \int_{\mathcal{D}} \frac{\partial^\beta (g_1 f_X)}{\partial x^\beta}(x) \frac{\partial^{\beta'} (g_1 f_X)}{\partial x^{\beta'}}(x) \frac{g_2(x)}{f_X(x)} dx \quad (29) $$

vanishes as soon as

$$ \int_{\mathcal{D}} \frac{g_2(x)}{f_X(x)} \left| \frac{\partial^\beta (g_1 f_X)}{\partial x^\beta} \frac{\partial^{\beta'} (g_1 f_X)}{\partial x^{\beta'}} \right| dx < \infty $$

for any $\beta$ and $\beta'$ such that $1 \le |\beta|, |\beta'| < \lfloor \alpha \rfloor$. Let us turn to the remaining term in (28). Using the definition of $R_\beta(x + hA_x^{-1}(u))$, we upper bound $T_1/(2h^{2\lfloor \alpha \rfloor})$ by

$$ \int_{\mathcal{D}} \left( \sum_{|\beta|=\lfloor \alpha \rfloor} \frac{1}{\beta!} \int_{\mathcal{D}} (A_x^{-1}(u))^\beta K(u) \int_0^1 (1-t)^{\lfloor \alpha \rfloor-1} \frac{\partial^\beta (g_1 f_X)}{\partial x^\beta} (x+thA_x^{-1}(u)) dt du \right)^2 \frac{g_2(x)}{f_X(x)} dx \\ = \int_{\mathcal{D}} \left( \sum_{|\beta|=\lfloor \alpha \rfloor} \frac{1}{\beta!} \int_{\mathcal{D}} (A_x^{-1}(u))^\beta K(u) \int_0^1 (1-t)^{\lfloor \alpha \rfloor-1} \right. \\ \left. \times \left( \frac{\partial^\beta (g_1 f_X)}{\partial x^\beta} (x-tuh) - \frac{\partial^\beta (g_1 f_X)}{\partial x^\beta} (x) \right) dt du \right)^2 \frac{g_2(x)}{f_X(x)} dx. $$

It suffices to proceed as previously in (29) rewriting $(\sum_\beta z_\beta)^2$ as $\sum_{\beta,\beta'} z_\beta z_{\beta'}$. As for the first term on the right hand side of the previous inequality, we apply (26) to get

$$ \begin{aligned} |T_1| & \le Ch^{2\alpha} \left( \int_{\mathcal{D}} \frac{g_2(x)}{f_X(x)} dx \right) \times \left( \sum_{|\beta|=\lfloor \alpha \rfloor} \frac{1}{\beta!} \right)^2 \times \left( \int_{\mathcal{D}} |A_x^{-1}(u)|^\beta |A_x^{-1}(u)|^{\alpha-\lfloor \alpha \rfloor} |K(u)| du \right)^2 \\ & \times \left( \int_0^1 (1-t)^{\lfloor \alpha \rfloor-1} t^{\alpha-\lfloor \alpha \rfloor} dt \right)^2 \\ & = Ch^{2\alpha} \left( \int_{\mathcal{D}} \frac{g_2(x)}{f_X(x)} dx \right) \times \left( \sum_{|\beta|=\lfloor \alpha \rfloor} \frac{1}{\beta!} \right)^2 \times \left( \int_{\mathcal{D}} \|u\|_{\infty}^\alpha |K(u)| du \right)^2 \\ & \le Ch^{2\alpha} \end{aligned} \quad (30) $$

since $g_2/f_X \in \mathbb{L}^1(\mathcal{D})$.

• Let us turn to $T_2$. We proceed as for $T_1$ but performing a Taylor expansion of $g_1$ rather
with

$$R_{\beta}(x + vh) = \frac{\lfloor \alpha \rfloor}{\beta!} \int_{0}^{1} (1 - t)^{\lfloor \alpha \rfloor - 1} \frac{\partial^{\beta} f_{X}}{\partial x^{\beta}}(x + tvh) dt.$$

Then, recalling that $(A_x^{-1}(v))^{\beta} = (A_x^{-1}(v)_1)^{\beta_1} \dots (A_x^{-1}(v)_d)^{\beta_d}$, we have

$$ \begin{aligned} \mathbb{E}[\hat{f}_{n,h,i}(x)] - f_X(x) &= \sum_{1 \le |\beta| < \lfloor \alpha \rfloor} \frac{h^{\lfloor \beta \rfloor}}{\beta!} \left( \int_{[0,1]^d} (A_x^{-1}(v))^{\beta} K(v) dv \right) \frac{\partial^{\beta} f_X}{\partial x^{\beta}}(x) \\ &\quad + h^{\lfloor \alpha \rfloor} \sum_{|\beta| = \lfloor \alpha \rfloor} \int_{[0,1]^d} (A_x^{-1}(v))^{\beta} K(v) R_{\beta}(x + hA_x^{-1}(v)) dv \\ &= \lfloor \alpha \rfloor h^{\lfloor \alpha \rfloor} \sum_{|\beta| = \lfloor \alpha \rfloor} \frac{1}{\beta!} \int_{[0,1]^d} (A_x^{-1}(v))^{\beta} K(v) \int_0^1 (1-t)^{\lfloor \alpha \rfloor - 1} \frac{\partial^{\beta} f_X}{\partial x^{\beta}}(x + t h A_x^{-1}(v)) dt dv \quad (29) \end{aligned} $$

using the fact that $(A_x^{-1}(v))^{\beta}$ is polynomial in $(v_1, \dots, v_d)$ of degree $\beta$ and $K$ is of order $(\lfloor \alpha \rfloor + 1)$. Now, using again that $K$ is of order $(\lfloor \alpha \rfloor + 1)$, we get

$$ \begin{aligned} \mathbb{E}[\hat{f}_{n,h,i}(x)] - f_X(x) &= \lfloor \alpha \rfloor h^{\lfloor \alpha \rfloor} \sum_{|\beta|=\lfloor \alpha \rfloor} \frac{1}{\beta!} \int_{[0,1]^d} (A_x^{-1}(v))^{\beta} K(v) \\ & \qquad \int_0^1 (1-t)^{\lfloor \alpha \rfloor - 1} \left( \frac{\partial^{\beta} f_X}{\partial x^{\beta}}(x + t h A_x^{-1}(v)) - \frac{\partial^{\beta} f_X}{\partial x^{\beta}}(x) \right) dt dv. \end{aligned} $$

Then, using (3) since $f_X \in \mathcal{C}^{\alpha}([0, 1]^d)$, one gets

$$ \left| \frac{\partial^\beta f_X}{\partial x^\beta} (x + t h A_x^{-1}(v)) - \frac{\partial^\beta f_X}{\partial x^\beta} (x) \right| \le C_{f_X} (ht \|A_x^{-1}(v)\|_\infty)^{\alpha - \lfloor \alpha \rfloor} $$

for all $v \in [0, 1]^d$, $x \in [0, 1]^d$, $t \in [0, 1]$, $h \in (0, \infty)$, and $\beta \in \mathbb{N}^d$ such that $|\beta| = \lfloor \alpha \rfloor$. Then,

$$ \begin{aligned} & |\mathbb{E}[\hat{f}_{n,h,i}(x)] - f_X(x)| \\ &\le C_{f_X} \lfloor \alpha \rfloor h^{\alpha} \left( \sum_{|\beta|=\lfloor \alpha \rfloor} \frac{1}{\beta!} \right) \left( \int_{[0,1]^d} \|A_x^{-1}(v)\|_{\infty}^{\beta} \|A_x^{-1}(v)\|_{\infty}^{\alpha - \lfloor \alpha \rfloor} |K(v)| dv \right) \left( \int_0^1 (1-t)^{\lfloor \alpha \rfloor - 1} t^{\alpha - \lfloor \alpha \rfloor} dt \right) \\ & = C h^{\alpha} \left( \sum_{|\beta|=\lfloor \alpha \rfloor} \frac{1}{\beta!} \right) \left( \int_{[0,1]^d} |K(v)| dv \right) \left( \int_0^1 (1-t)^{\lfloor \alpha \rfloor - 1} t^{\alpha - \lfloor \alpha \rfloor} dt \right) \le C h^{\alpha} \end{aligned} $$

since $\|A_x^{-1}(v)\|_{\infty}^{\beta} \le 1$ and $\|K\|_{\infty} < \infty$.

<figure><img src="image_9.png" /></figure>

*Proof of Lemma 3.2.* For $u \subset \{1, \dots, d\}$, let $\sim u = \{1, \dots, d\} \setminus u$. For $x \in [0, 1]^d$,

$$ (n-1) \hat{f}_{n,h,i}(x) = \sum_{j \ne i} \prod_{\substack{1 \le k \le d \text{ s.t.} \\ x_k \in [0,1/2]}} k_h(X_{k,j} - x_k) \prod_{\substack{1 \le k \le d \text{ s.t.} \\ x_k \in (1/2,1]}} k_h(x_k - X_{k,j}) $$

$$ = \sum_{j \ne i} \sum_{u \subseteq \{1, \dots, d\}} \prod_{k \in u} \mathbb{1}_{[0,1/2]}(x_k) k_h(X_{k,j} - x_k) \prod_{k \sim u} \mathbb{1}_{(1/2,1]}(x_k) k_h(x_k - X_{k,j}) =: \sum_{j \ne i} K_{h,x}^*(X_j). $$

Arguing as in [41, Section 5.1.2] and since the class of functions $\mathcal{F}$ defined by (12) is a uniformly bounded VC-class of functions by (A8), we can say that the class $\mathcal{F}^* = \{K_{h,x}^*, h > 0, x \in [0, 1]^d\}$ is also a uniformly bounded VC-class of functions. Thus it follows from a particular case of [30, Corollary 13] that the bound in (13) holds with probability $1 - \delta_n$. $\square$

<figure><img src="image_10.png" /></figure>
than $g_1 f_X$.

$$
\begin{align*}
T_2 &= \int_{\mathcal{D}} \left( \int_{\mathcal{D}} (K_h \circ A_u)(x-u) g_1(u) du - g_1(x) \right)^2 g_2(x) f_X(x) dx \\
&= \int_{\mathcal{D}} \left( \int_{\mathcal{D}} K(u) [g_1(x+hA_x^{-1}(u)) - g_1(x)] du \right)^2 g_2(x) f_X X(x) dx \\
&= \int_{\mathcal{D}} \left( \sum_{1 \le |\beta| < \lfloor \alpha \rfloor} \frac{h^{|\beta|}}{\beta!} \left( \int_{\mathcal{D}} (A_x^{-1}(u))^\beta K(u) du \right) \frac{\partial^\beta g_1}{\partial x^\beta}(x) \right. \\
&\qquad \left. + h^{\lfloor \alpha \rfloor} \sum_{|\beta|=\lfloor \alpha \rfloor} \int_{\mathcal{D}} (A_x^{-1}(u))^\beta K(u) R_\beta(x+hA_x^{-1}(u)) du \right)^2 g_2(x) f_X(x) dx \\
&\le 2 \int_{\mathcal{D}} \left( \sum_{1 \le |\beta| < \lfloor \alpha \rfloor} \frac{h^{|\beta|}}{\beta!} \left( \int_{\mathcal{D}} (A_x^{-1}(u))^\beta K(u) du \right) \frac{\partial^\beta g_1}{\partial x^\beta}(x) \right)^2 g_2(x) f_X(x) dx \\
&\qquad + 2h^{2\lfloor \alpha \rfloor} \int_{\mathcal{D}} \left( \sum_{|\beta|=\lfloor \alpha \rfloor} \int_{\mathcal{D}} (A_x^{-1}(u))^\beta K(u) R_\beta(x+hA_x^{-1}(u)) du \right)^2 g_2(x) f_X(x) dx.
\end{align*}
$$

The first integral vanishes as soon as

$$
\int_{\mathcal{D}} g_2(x) f_X(x) \left| \frac{\partial^\beta g_1}{\partial x^\beta} \frac{\partial^{\beta'} g_1}{\partial x^{\beta'}} \right| dx < \infty
$$

for any $\beta$ and $\beta'$ such that $1 \le |\beta|, |\beta'| < \lfloor \alpha \rfloor$. As for the second, one has

$$
\begin{align*}
T_2 \le 2h^{2\lfloor\alpha\rfloor} & \int_{\mathcal{D}} \left( \sum_{|\beta|=\lfloor\alpha\rfloor} \frac{1}{\beta!} \int_{\mathcal{D}} (A_x^{-1}(u))^\beta K(u) \int_0^1 (1-t)^{\lfloor\alpha\rfloor-1} \frac{\partial^\beta g_1}{\partial x^\beta} (x+thA_x^{-1}(u)) dt du \right)^2 g_2(x) f_X(x) dx \\
& = 2h^{2\lfloor\alpha\rfloor} \int_{\mathcal{D}} \left( \sum_{|\beta|=\lfloor\alpha\rfloor} \frac{1}{\beta!} \int_{\mathcal{D}} (A_x^{-1}(u))^\beta K(u) \int_0^1 (1-t)^{\lfloor\alpha\rfloor-1} \right. \\
& \qquad \left. \times \left( \frac{\partial^\beta (g_1 f_X)}{\partial x^\beta} (x-tuh) - \frac{\partial^\beta g_1}{\partial x^\beta}(x) \right) dt du \right)^2 g_2(x) f_X(x) dx.
\end{align*}
$$

It suffices to proceed as previously in (30) to get

$$
T_2 \le Ch^{2\alpha}. \tag{31}
$$

since $\int_D g_2(x) f_X(x) dx = \mathbb{E}[g_2(X)] = \mathbb{E}[Y^2] < \infty$ by assumption on $Y$.

Finally,

$$
n \mathbb{E}[S_1^2] \le C h^{2\alpha}. \qquad (32)
$$
Proof of Corollary 3.4. First, by Lemmas 3.1 and 3.2, we get

$$
\begin{align*}
\| \hat{f}_{n,h,i} - f_X \|_{\infty}^2 &\le 2 \| \hat{f}_{n,h,i} - \mathbb{E}[\hat{f}_{n,h,i}] \|_{\infty}^2 + 2 \| \mathbb{E}[\hat{f}_{n,h,i}] - f_X \|_{\infty}^2 \\
&= O_{\mathbb{P}} \left( h^{2\alpha} + \frac{\left( \log\left(\frac{1}{h_n}\right) \right)_+ + \log\left(\frac{2}{\delta_n}\right)}{nh_n^d} \right) = o_{\mathbb{P}}(n^{-1/2})
\end{align*}
$$

as soon as $nh^{2d} \to \infty$ and $nh^{4\alpha} \to 0$ by considering, e.g., $\delta_n = 1/n$. Since the input space $[0,1]^d$ is compact by Assumption (A1), one concludes that (14) holds.

Moreover, for all $i \in \{1, \dots, d\}$, as $\hat{f}_{n,h,i}(x) = f_X(x) + \hat{f}_{n,h,i}(x) - f_X(x)$, and from Assumption (A2), we get with probability $1 - 1/\delta_n$ and for $n$ large enough:

$$
\inf_{x \in [0,1]^d} |\hat{f}_{n,h,i}(x)| \ge \delta - C \left( h^{2\alpha} + \frac{\left(\log\left(\frac{1}{h_n}\right)\right)_+ + \log\left(\frac{2}{\delta_n}\right)}{nh_n^d} \right)
$$

for some positive constant $C$. This last quantity is nonnegative as soon as $n$ is large enough when $nh^d \to \infty$ and $h^\alpha \to 0$ by considering, e.g., $\delta_n = 1/n$. $\square$

We state Lemma A.1 below as a preliminary result for the proof of (i) of Theorem 4.1.

**Lemma A.1.** Under Assumptions (A1), (A5), (A6), and (A8),

$$
\forall i \in \{1, \dots, d\}, \quad \mathbb{E} \left[ \int_{[0,1]^d} (\hat{g}_{n,h,i}(x) - g(x))^2 dx \right] = o(n^{-1/2}) \quad (30)
$$

as soon as $nh^{2d} \to \infty$, and $nh^{4\alpha} \to 0$.

Proof of Lemma A.1. Similarly as for density estimation, we first prove that $\|\mathbb{E}[\hat{g}_{n,h,1}] - g\|_{\infty} = O(h^{\alpha})$. Because $\mathbb{E}[\hat{g}_{n,h,1}(x)] = \int_{[0,1]^d} K_h \circ A_x(x' - x)g(x')dx'$, we can follow the same lines as in the proof of (11) since $g \in C^{\alpha}([0,1]^d)$. Second, we prove that there exists some constant $C > 0$ such that, with probability $1 - \delta_n$,

$$
\| \hat{g}_{n,h,1} - \mathbb{E}[\hat{g}_{n,h,1}] \|_{\infty} \le C \left( \frac{\log\left(\frac{1}{h}\right) + \log\left(\frac{2}{\delta_n}\right)}{nh^d} \right). \quad (31)
$$

One has, for all $i \in \{1, \dots, d\}$,

$$
\begin{align*}
(n-1)\hat{g}_{n,h,i}(x) &= \sum_{j \neq i} Y_j \prod_{\substack{1 \le k \le d \text{ s.t.} \\ x_k \in [0,1/2]}} k_h(X_{k,j} - x_k) \prod_{\substack{1 \le k \le d \text{ s.t.} \\ x_k \in (1/2,1]}} k_h(x_k - X_{k,j}) \\
&= \sum_{j \neq i} Y_j \sum_{u \subseteq \{1, \dots, d\}} \prod_{k \in u} \mathbb{1}_{[0,1/2]}(x_k) k_h(X_{k,j} - x_k) \prod_{k \sim u} \mathbb{1}_{(1/2,1]}(x_k) k_h(x_k - X_{k,j}) = O\left(\sum_{j \neq i} K_{h,x}^*(X_j)\right).
\end{align*}
$$

We conclude as in the proof of Lemma 3.2. Finally, we deduce (30) as soon as $nh^d \to \infty$ and $h^\alpha \to 0$ by considering, e.g., $\delta_n = 1/n$. $\square$
(2) As for the second moment of $S_2$, since $\pi_2$ is a projection, we have

$$
\begin{aligned}
\mathbb{E}[S_2^2] &\le \frac{1}{n(n-1)} \iint_{\mathcal{D}^2} \left( \frac{K_h^2 \circ A_{x_1}(x_2 - x_1)}{f_X(x_1)^2} + \frac{K_h^2 \circ A_{x_2}(x_1 - x_2)}{f_X(x_2)^2} \right) f_X(x_1) f_X(x_2) g_2(x_1) g_2(x_2) dx_1 dx_2 \\
&= \frac{2}{n(n-1)} \iint_{\mathcal{D}^2} K_h^2 \circ A_{x_1}(x_2 - x_1) \frac{f_X(x_2)}{f_X(x_1)} g_2(x_1) g_2(x_2) dx_1 dx_2 \\
&\le \frac{2\kappa}{n(n-1)} \sup_{\ell=1,\dots,\kappa} \int_{\mathcal{D}} \int_{\mathcal{D}_\ell} K_h^2 \circ A_\ell(x_2 - x_1) \frac{f_X(x_2)}{f_X(x_1)} g_2(x_1) g_2(x_2) dx_1 dx_2 \\
&\le \frac{2\kappa}{n(n-1)} \sup_{\ell=1,\dots,\kappa} \iint_{\mathcal{D}} K_h^2 \circ A_\ell(x_2 - x_1) \frac{f_X(x_2)}{f_X(x_1)} g_2(x_1) g_2(x_2) dx_1 dx_2 \\
&= \frac{2\kappa}{n(n-1)} \sup_{\ell=1,\dots,\kappa} \int_{\mathcal{D}} (K_h^2 \circ A_\ell) * \frac{g_2}{f_X}(y) g_2(y) f_X(y) dy \\
&\le \frac{2\kappa}{n(n-1)} \|g_2 f_X\|_2 \left\| \frac{g_2}{f_X} \right\|_2 \sup_{\ell=1,\dots,\kappa} \|K_h^2 \circ A_\ell\|_1 \\
&= \frac{2\kappa}{n(n-1)h^d} \|g_2 f_X\|_2 \left\| \frac{g_2}{f_X} \right\|_2 \|K\|_1^2
\end{aligned}
$$

using Young's convolution inequality again with $r = q = 2$ and $p = 1$. Since $g_2/f_X$ and $g_2f_X \in L^2(\mathcal{D})$, it yields

$$ \mathbb{E}[S_2^2] \le \frac{C}{n(n-1)h^d}. \qquad (33) $$

We conclude to (8) using (32) and (33). $\square$

*Proof of (10).* It is straightforward using (7), (8), and the standard central limit theorem for the random variables $Z_i$. $\square$

## A.2 Proof of the remaining results

*Proof of Corollary 3.3.* To prove the global CLT, we first prove a CLT for the following linear combination

$$ S_{n,h} := aT_{n,h} + b\bar{Y}_n + c\bar{Y}^2_n = aT_{n,h} + \frac{b}{n}\sum_{j=1}^{n}Y_j + \frac{c}{n}\sum_{j=1}^{n}Y_j^2 $$

for any $a, b$, and $c \in \mathbb{R}$ from which we deduce that $(T_{n,h}, \bar{Y}_n, \bar{Y}^2_n)$ is asymptotically Gaussian. Second, it suffices to conclude applying the delta method to the later vector

$$ (T_{n,h}, \bar{Y}_n, \bar{Y}^2_n) = (T_{n,h}, \frac{1}{n}\sum_{j=1}^{n}Y_j, \frac{1}{n}\sum_{j=1}^{n}Y_j^2) $$

with

$$ \Phi(x, y, z) = \frac{x - y^2}{z - y^2}. $$
## A.2 Proof of the preliminary results of Section 3.2

Now, recall that $\tilde{f}_{n,h,i}$ defined in (16) and $\tilde{g}_{n,h,i}$ defined in 17 are the leave-one-out estimators of $f_X$ and $g$ respectively based on the $(n-1)$-input/output sample $((X_1, Y_1), \dots, (X_{i-1}, Y_{i-1}), \dots, (X_{i+1}, Y_{i+1}), \dots, (X_n, Y_n))$. Then, the estimator of the regression function $m$ is naturally given by $\tilde{m}_{n,h,i} = \tilde{g}_{n,h,i}/\tilde{f}_{n,h,i}$.

**Proof of Lemma 3.5.** It follows the same lines as in the proof of Lemma 3.1. We consider the case $x \in [0, 1/2]^d$. The remaining cases can be deduced by symmetry. Let $\mathcal{A} := \{k \in \{1, \dots, d\} : x_k \le h\}$. For all $k \in \{1, \dots, d\}$ and all $t_k \in [0, 1]$, $k_h(M^1(t_k) - x_k) = k_h(2 - t_k - x_k) = 0$ as $\tilde{k}$ is supported on $[-1, 1]$ and $h < 1/2$. In addition, for all $k \in \{1, \dots, d\} \setminus \mathcal{A}$ and all $t_k \in [0, 1]$, $k_h(M^{-1}(t_k) - x_k) = \tilde{k}_h(-t_k - x_k) = 0$.

For any subset $\mathcal{B} \subset \mathcal{A}$, define $x_{\mathcal{B}}$ such that $x_{\mathcal{B},k} = x_k = M^0(x_k)$ if $k \notin \mathcal{B}$ and $x_{\mathcal{B},k} = -x_k = M^{-1}(x_k)$ if $k \in \mathcal{B}$. Then, the expected value of $\tilde{f}_{n,h,i}(x)$ can be written as:

$$ 
\begin{aligned}
\mathbb{E}[\tilde{f}_{n,h,i}(x)] &= \sum_{\mathcal{B} \subset \mathcal{A}} \int_{[0,1]^d} \prod_{k \in \mathcal{B}} \tilde{k}_h(-z_k - x_k) \prod_{k \notin \mathcal{B}} \tilde{k}_h(z_k - x_k) f_X(z) dz \\
&= \sum_{\mathcal{B} \subset \mathcal{A}} \int_{\mathcal{X}_{\mathcal{B}}} \prod_{k=1}^d \tilde{k}(u_k) f_X(z_{\mathcal{B}}) du = \sum_{\mathcal{B} \subset \mathcal{A}} \int_{\mathcal{X}_{\mathcal{B}}} \tilde{K}(v) f_X(z_{\mathcal{B}}) dv
\end{aligned}
$$

with

$$ z_{\mathcal{B},k} = \begin{cases} -x_k - v_k h, & k \in \mathcal{B} \\ x_k + v_k h, & k \notin \mathcal{B} \end{cases} \quad \text{and} \quad \mathcal{X}_{\mathcal{B},k} = \begin{cases} \{v_k \in [-1, 1) : -x_k - v_k h \in [0, 1]\}, & k \in \mathcal{B} \\ \{v_k \in [-1, 1) : x_k + v_k h \in [0, 1]\}, & k \notin \mathcal{B} \end{cases} $$

(recalling that the support of $\tilde{K}$ is $[-1, 1]^d$) and $\mathcal{X}_{\mathcal{B}} = \prod_{k=1}^d \mathcal{X}_{\mathcal{B},k}$. Note that $(\mathcal{X}_{\mathcal{B}})_{\mathcal{B} \subset \mathcal{A}}$ forms a partition of $[-1, 1]^d$. Indeed, one has

$$ \mathcal{X}_{\mathcal{B}} = \prod_{k=1}^{d} \mathcal{X}_{\mathcal{B},k} = \prod_{k \in \mathcal{B}} [-1, -\frac{x_k}{h}) \prod_{k \in \mathcal{A} \setminus \mathcal{B}} [-\frac{x_k}{h}, 1) \prod_{k \in \{1, \dots, d\} \setminus \mathcal{A}} [-1, 1). $$

Let us now define $f_X^{\text{MI}}$ on $[-1, 2]^d$ such that, for all $y \in [-1, 2]^d$,

$$ f_X^{\text{MI}}(y) = f_X(y^{\text{MI}}) \quad \text{where} \quad y_k^{\text{MI}} = \begin{cases} -y_k & \text{if } y_k \in [-1, 0] \\ y_k & \text{if } y_k \in [0, 1] \\ 2 - y_k & \text{if } y_k \in [1, 2] \end{cases}. $$

Then, for all $x \in [0, 1]^d$ and $a \in \{-1, 0, 1\}^d$, $f_X^{\text{MI}}(M^a(x)) = f_X(x)$ and similarly, for any $\mathcal{B} \subset \mathcal{A}$, $f_X^{\text{MI}}(x_{\mathcal{B}}) = f_X(x)$. Now, since $f_X^{\text{MI}} = f_X$ on $[0, 1]^d$ and for any $w \in \mathcal{X}_{\mathcal{B}}$, $f_X^{\text{MI}}(z_{\mathcal{B}}) = f_X(z_{\mathcal{B}})$, we have

$$ \mathbb{E}[\tilde{f}_{n,h,i}(x)] = \sum_{\mathcal{B} \subset \mathcal{A}} \int_{\mathcal{X}_{\mathcal{B}}} \tilde{K}(v) f_X^{\text{MI}}(z_{\mathcal{B}}) dv = \sum_{\mathcal{B} \subset \mathcal{A}} \int_{\mathcal{X}'_{\mathcal{B}}} \tilde{K}(w) f_X^{\text{MI}}(z_{\mathcal{B}}) dw \quad (32) $$

using the variable change $w = v_{\mathcal{B}}$, the symmetry of $\tilde{K}$ and with

$$ \mathcal{X}'_{\mathcal{B},k} = \begin{cases} \{w_k \in [-1, 1) : -x_k + w_k h \in [0, 1]\}, & k \in \mathcal{B} \\ \{w_k \in [-1, 1) : x_k + w_k h \in [0, 1]\}, & k \notin \mathcal{B} \end{cases}. $$
• We proceed as in the proof of Theorem 3.1: we control the bias and the variance terms as in (7) and (8) respectively. Then we conclude to the central limit theorem as in the proof of (10).

Let us first study the bias term. Let $m := a\mathbb{E}[\mathbb{E}[Y|X]^2] + b\mathbb{E}[Y] + c\mathbb{E}[Y^2]$. Then,

$$ 
\begin{aligned}
\mathbb{E}[S_{n,h}] - m &= \mathbb{E}[S_{n,h}] - (c\mathbb{E}[\mathbb{E}[Y|X]^2] + b\mathbb{E}[Y] + c\mathbb{E}[Y^2]) \\
&= a(\mathbb{E}[T_{n,h}] - \mathbb{E}[\mathbb{E}[Y|X]^2]) + b(\mathbb{E}[\bar{Y}_n] - \mathbb{E}[Y]) + c(\mathbb{E}[\bar{Y}_n^2] - \mathbb{E}[Y^2]) \\
&= a(\mathbb{E}[T_{n,h}] - \mathbb{E}[\mathbb{E}[Y|X]^2])
\end{aligned}
$$

and, by (7) in Theorem 3.1, we conclude that

$$ |\mathbb{E}[S_{n,h}] - m| \le aCh^{\alpha} $$

As for the variance, the Hoeffding decomposition given in (6) becomes

$$ 
\begin{aligned}
S_{n,h} - \mathbb{E}[S_{n,h}] &= \frac{1}{n} \sum_{j=1}^{n} \left( aZ_j + b(Y_j - \mathbb{E}[Y]) + c(Y_j^2 - \mathbb{E}[Y^2]) \right) + a(S_1 + S_2) \\
&=: \frac{1}{n} \sum_{j=1}^{n} \tilde{Z}_j + a(S_1 + S_2),
\end{aligned}
$$

where $Z_j$ has been defined in (9), $S_1$ and $S_2$ in (6). Then we deduce from (32) and (33) of Appendix A.1:

$$ \mathbb{E}\left[\left(S_{n,h} - \mathbb{E}[S_{n,h}] - \frac{1}{n}\sum_{j=1}^{n} \tilde{Z}_j\right)^2\right] = a^2 \mathbb{E}\left[\left(S_1 + S_2\right)^2\right] \le a^2 C \left(h^{2\alpha} + \frac{1}{n^2 h^d}\right). $$

Thus it remains to use a CLT for $\frac{1}{n}\sum_{j=1}^{n} \tilde{Z}_j$. By Theorem 3.1 II(b), one gets merely directly

$$ \sqrt{n}(S_{n,h} - m) \xrightarrow{n \to \infty} \mathcal{N}\left(0, \text{Var}(2aYg_1(X) + bY + cY^2)\right) $$

from which we conclude that $(T_{n,h}, \bar{Y}_n, \bar{Y}_n^2)$ is asymptotically Gaussian; namely

$$ \sqrt{n}\left(\left(\frac{T_{n,h}}{\bar{Y}_n}\right) - M\right) \xrightarrow{n \to \infty} \mathcal{N}_3(0, \Gamma), $$

where $M := (\mathbb{E}[\mathbb{E}[Y|X]^2], \mathbb{E}[Y], \mathbb{E}[Y^2])^\top$ and

$$ \Gamma := \begin{pmatrix} 4\tau^2 & 2\text{Cov}(Yg_1(X), Y) & 2\text{Cov}(Yg_1(X), Y^2) \\ 2\text{Cov}(Yg_1(X), Y) & \text{Var}(Y) & \text{Cov}(Y, Y^2) \\ 2\text{Cov}(Yg_1(X), Y^2) & \text{Cov}(Y, Y^2) & \text{Var}(Y^2) \end{pmatrix}. $$

• Now, as announced previously, we apply the delta method to $(T_{n,h}, \bar{Y}_n, \bar{Y}_n^2)$ and

$$ \Phi(x, y, z) = \frac{x - y^2}{z - y^2} $$
In addition, since $\int_{[-1,1]^d} \tilde{K}(v)dv = 1$, $(\mathcal{X}'_B)_{B\subset A}$ forms also a partition of $[-1,1)^d$, and $f_X^{\text{MI}}(x_B) = f_X(x)$, one has

$$f_X(x) = f_X(x) \int_{[-1,1]^d} \tilde{K}(v)dv = \sum_{B\subset A} \int_{\mathcal{X}'_B} \tilde{K}(v)f_X(x)dv = \sum_{B\subset A} \int_{\mathcal{X}'_B} \tilde{K}(v)f_X^{\text{MI}}(x_B)dv. \quad (33)$$

From (32) and (33), we deduce

$$\mathbb{E}[\tilde{f}_{n,h,i}(x)] - f_X(x) = \sum_{B \subset A} \int_{\mathcal{X}'_B} \tilde{K}(v) \left( f_X^{\text{MI}}(x_B + vh) - f_X^{\text{MI}}(x_B) \right) dw.$$

By construction and since Assumption (A'4) is satisfied (in particular derivatives up to order $\lfloor \alpha \rfloor$ vanish near the boundary), $f_X^{\text{MI}}$ has derivatives up to order $\lfloor \alpha \rfloor$. Then, $f_X^{\text{MI}}$ belongs to $C^\alpha([-1, 2]^d)$ with $C_{f_X^{\text{MI}}} = 3C_{f_X}$. Thus, proceeding as in the proof of Lemma 3.1 we write:

$$f_X^{\text{MI}}(x_B + vh) - f_X^{\text{MI}}(x_B) = \sum_{1 \le |\beta| < \lfloor \alpha \rfloor} \frac{h^{\lfloor \beta \rfloor}}{\beta!} v^\beta \frac{\partial^\beta f_X^{\text{MI}}}{\partial x^\beta}(x_B) + h^{\lfloor \alpha \rfloor} \sum_{|\beta| = \lfloor \alpha \rfloor} R_\beta(x_B + vh) v^\beta \quad (34)$$

with

$$R_\beta(x_B + vh) = \frac{\lfloor \alpha \rfloor}{\beta!} \int_0^1 (1-t)^{\lfloor \alpha \rfloor - 1} \frac{\partial^\beta f_X^{\text{MI}}}{\partial x^\beta}(x_B + tvh) dt. \quad (35)$$

Then, with similar arguments as in the proof of Lemma 3.1, together with Assumptions (A'4) and (A'6) and the fact that $(\mathcal{X}'_B)_{B\subset A}$ forms a partition of $[-1,1)^d$, we get

$$\begin{aligned} & |\mathbb{E}[\tilde{f}_{n,h,i}(x)] - f_X(x)| \\ & \le 3C_{f_X} \lfloor \alpha \rfloor h^\alpha \left( \sum_{|\beta|=\lfloor \alpha \rfloor} \frac{1}{\beta!} \right) \left( \int_0^1 (1-t)^{\lfloor \alpha \rfloor - 1} t^{\alpha - \lfloor \alpha \rfloor} dt \right) \left( \sum_{B \subset A} \int_{\mathcal{X}'_B} |v|^{\alpha} |\tilde{K}(v)| dv \right) \\ & \le 3C_{f_X} \lfloor \alpha \rfloor h^\alpha \left( \sum_{|\beta|=\lfloor \alpha \rfloor} \frac{1}{\beta!} \right) \left( \int_0^1 (1-t)^{\lfloor \alpha \rfloor - 1} t^{\alpha - \lfloor \alpha \rfloor} dt \right) \left( \int_{[-1,1]^d} |\tilde{K}(v)| dv \right) = Ch^\alpha. \quad \square \end{aligned}$$

We state Lemma A.2 below as a preliminary result to the proof of (ii) of Theorem 4.1.

**Lemma A.2.** Let $\alpha > 0$. Under Assumptions (A1), (A5), (A'6), and (A8),

$$\forall i \in \{1, \dots, d\}, \quad \mathbb{E}[\int (\tilde{g}_{n,h,i}(x) - g(x))^2 dx] = o(n^{-1/2}) \quad (36)$$

as soon as $nh^{2d} \to \infty$, and $nh^{4\alpha} \to 0$.

*Proof of Lemma A.2.* The proof follows similar lines as the one of Lemma A.1. Recall that for $x \in [0,1]^d$, $g(x) = f_X(x)m(x)$. We extend $m$ as $m^{\text{MI}}$ on $[-1,2]^d$ as follows. For all $y \in [-1,2]^d$, $m^{\text{MI}}(y) = m(y^{\text{MI}})$ with $y^{\text{MI}}$ defined in the proof of Lemma 3.5. Then, for all $x \in [0,1]^d$ and $a \in \{-1,0,1\}^d$, $m^{\text{MI}}(M^a(x)) = m(x)$ and similarly, for any $B \subset A$, $m^{\text{MI}}(x_B) = m(x)$. We introduce the function $g^{\text{MI}}(x) = f_X^{\text{MI}}(x)m^{\text{MI}}(x)$. We also introduce the intermediate function $\tilde{g}^{\text{MI}}(x) = f_X^{\text{MI}}(x)m(x)$. Now following the proof of Lemma A.1,
to get

$$ \sqrt{n}\left(\frac{T_{n,h}-\left(\frac{1}{n} \sum_{j=1}^{n} Y_{j}\right)^{2}}{\frac{1}{n} \sum_{j=1}^{n} Y_{j}^{2}-\left(\frac{1}{n} \sum_{j=1}^{n} Y_{j}\right)^{2}}-S^{X}\right) \xrightarrow{n \to \infty} \mathcal{N}(0, \sigma^{2}) $$

where $\sigma^2 = \ell^T \Gamma \ell$ with $\ell = \nabla \Phi(M)$. By assumption $\mathrm{Var}(Y) \neq 0$, $\Phi$ is differentiable at $M$ and we will see in the sequel that $\ell^T \Gamma \ell \neq 0$, so that the application of the delta method is justified. By differentiation, we get that, for any $x, y$, and $z$ so that $z \neq y^2$:

$$ \nabla \Phi(x, y, z) = \left( \frac{y^2}{z - y^2}, -2y \frac{z - x}{(z - y^2)^2}, -\frac{x - y^2}{(z - y^2)^2} \right)^T $$

so that

$$ \begin{aligned} \ell = \nabla \Phi(M) &= \left( \frac{1}{\mathrm{Var}(Y)}, -2\mathbb{E}[Y] \frac{\mathbb{E}[Y^2] - \mathbb{E}[\mathbb{E}[Y|X]^2]}{\mathrm{Var}(Y)^2}, -\frac{S^X}{\mathrm{Var}(Y)} \right)^T \\ &= \frac{1}{\mathrm{Var}(Y)} \left( 1, 2\mathbb{E}[Y](S^X - 1), -S^X \right)^T. \end{aligned} $$

Hence the asymptotic variance $\sigma^2$ can be computed

$$ \begin{aligned} \sigma^2 = \ell^T \Gamma \ell &= \frac{1}{\mathrm{Var}(Y)^2} \left( \begin{aligned} &4(\mathrm{Var}(Yg_1(X)) - 2\mathrm{Cov}(Yg_1(X), Y)\mathbb{E}[Y] + \mathbb{E}[Y]^2\mathrm{Var}(Y)) \\ &+ 4S^X(2\mathrm{Cov}(Yg_1(X), Y)\mathbb{E}[Y] - \mathrm{Cov}(Yg_1(X), Y^2) - 2\mathbb{E}[Y]^2\mathrm{Var}(Y) + \mathbb{E}[Y]\mathrm{Cov}(Y, Y^2)) \\ &+ (S^X)^2(4\mathbb{E}[Y]^2\mathrm{Var}(Y) - 4\mathbb{E}[Y]\mathrm{Cov}(Y, Y^2) + \mathrm{Var}(Y^2)) \end{aligned} \right). \end{aligned} \quad (34) $$

The proof of Corollary 3.3 is now complete. □

*Proof of Corollary 3.4.* For any $i = 1, \dots, p$, define $T_{n,h,i}$ as the estimator of $\mathbb{E}[\mathbb{E}[Y|V_i]^2]$. We proceed as for the proof of Corollary 3.3 to get a vectorial central limit theorem for $(T_{n,h,1}, \dots, T_{n,h,p}, \bar{Y}_n, \bar{Y}_n^2)^T$. Let $M = (S^1, \dots, S^p, \mathbb{E}[Y], \mathbb{E}[Y^2])^T$. Then, using the delta method with

$$ \Phi(x_1, \dots, x_p, y, z) = \left( \frac{x_1 - y^2}{z - y^2}, \dots, \frac{x_p - y^2}{z - y^2} \right)^T $$

we get the result with

$$ \Sigma = J_{\Phi}(M)^{T} \Gamma J_{\Phi}(M) \quad (35) $$

with $J_{\Phi}(M)$ the Jacobian of $\Phi$ at point $M$ and $\Gamma$ the covariance matrix of

$$ (Y \mathbb{E}[Y|V_1], \dots, Y \mathbb{E}[Y|V_p], Y, Y^2)^T. \quad \square $$
we get, for all $x \in [0, 1]^d$, $\mathbb{E}[\tilde{g}_{n,h,i}(x)] - g(x) = \sum_{\mathcal{B} \subset \mathcal{A}} \int_{\mathcal{X}'_B} \tilde{K}(v) (g^{\text{MI}}(x_{\mathcal{B}} + vh) - g^{\text{MI}}(x_{\mathcal{B}})) dv$.
Now,

$$
\begin{align*}
g^{\text{MI}}(x_{\mathcal{B}} + vh) - g^{\text{MI}}(x_{\mathcal{B}}) &= g^{\text{MI}}(x_{\mathcal{B}} + vh) - \tilde{g}^{\text{MI}}(x_{\mathcal{B}} + vh) + \tilde{g}^{\text{MI}}(x_{\mathcal{B}} + vh) - \tilde{g}^{\text{MI}}(x_{\mathcal{B}}) \\
&\quad + \tilde{g}^{\text{MI}}(x_{\mathcal{B}} + vh) - \tilde{g}^{\text{MI}}(x_{\mathcal{B}}) + \tilde{g}^{\text{MI}}(x_{\mathcal{B}}) - g^{\text{MI}}(x_{\mathcal{B}}) \\
&= f_{X}^{\text{MI}}(x_{\mathcal{B}} + vh) (m^{\text{MI}}(x_{\mathcal{B}} + vh) - m(x_{\mathcal{B}} + vh)) \\
&\quad + \tilde{g}^{\text{MI}}(x_{\mathcal{B}} + vh) - \tilde{g}^{\text{MI}}(x_{\mathcal{B}}) + f_{X}^{\text{MI}}(x_{\mathcal{B}}) (m(x_{\mathcal{B}}) - m^{\text{MI}}(x_{\mathcal{B}})).
\end{align*}
$$

Now by definition, we have $m^{\text{MI}}(x_{\mathcal{B}} + vh) = m(x_{\mathcal{B}} + vh)$ and $m^{\text{MI}}(x_{\mathcal{B}}) = m(x_{\mathcal{B}})$. Thus we have $g^{\text{MI}}(x_{\mathcal{B}}+vh)-g^{\text{MI}}(x_{\mathcal{B}}) = \tilde{g}^{\text{MI}}(x_{\mathcal{B}}+vh)-\tilde{g}^{\text{MI}}(x_{\mathcal{B}})$. As $f_X^{\text{MI}}$ and $m$ belong to $C^{\alpha}([-1, 2]^d)$ from Assumptions (A'4) and (A5), the function $\tilde{g}^{\text{MI}}$ also belongs to $C^{\alpha}([-1, 2]^d)$. Then, mimicking (34) and (35), we obtain for any $y \in \mathcal{Y}$,

$$
\tilde{g}^{\text{MI}}(x_{\mathcal{B}} + vh) - \tilde{g}^{\text{MI}}(x_{\mathcal{B}}) = \sum_{1 \le |\beta| < |\alpha|} \frac{h^{|\beta|}}{\beta!} v^{\beta} \frac{\partial^{\beta} \tilde{g}^{\text{MI}}}{\partial x^{\beta}}(x_{\mathcal{B}}) + h^{|\alpha|} \sum_{|\beta| = |\alpha|} R_{\beta}(x_{\mathcal{B}} + vh) v^{\beta}
$$

with

$$
R_{\beta}(x_{\mathcal{B}} + vh) = \frac{\lfloor \alpha \rfloor}{\beta!} \int_{0}^{1} (1-t)^{\lfloor \alpha \rfloor - 1} \frac{\partial^{\beta} \tilde{g}^{\text{MI}}}{\partial x^{\beta}}(x_{\mathcal{B}} + tvh) \, dt.
$$

Then, we conclude with similar arguments as in the proof of Lemma 3.1. From Assumption (A'6) and as $(\mathcal{X}'_B)_{\mathcal{B}\subset\mathcal{A}}$ forms a partition of $[-1, 1)^d$, we have

$$
\sum_{\mathcal{B} \subset \mathcal{A}} \int_{\mathcal{X}'_B} \widetilde{K}(v) \sum_{|\beta| = |\alpha|} v^\beta dv = \int_{-1}^1 \left( \sum_{|\beta| = |\alpha|} v^\beta \right) \widetilde{K}(v) dv = 0.
$$

Thus

$$
\begin{align*}
\mathbb{E}[\tilde{g}_{n,h,i}(x)] - g(x) &= \sum_{\mathcal{B} \subset \mathcal{A}} \int_{\mathcal{X}'_B} \widetilde{K}(v) \left( g^{\text{MI}}(x_{\mathcal{B}} + vh) - g^{\text{MI}}(x_{\mathcal{B}}) \right) dv \\
&= \sum_{\mathcal{B} \subset \mathcal{A}} \int_{\mathcal{X}'_B} \widetilde{K}(v) h^{|\alpha|} \sum_{|\beta|=|\alpha|} R_{\beta}(x_{\mathcal{B}} + vh) v^{\beta} dv \\
&= \sum_{\mathcal{B} \subset \mathcal{A}} \int_{\mathcal{X}'_B} \widetilde{K}(v) h^{|\alpha|} \sum_{|\beta|=|\alpha|} \frac{|\alpha|}{\beta!} \\
&\qquad \int_0^1 (1-t)^{|\alpha|} (-1)^{|\alpha|} \left( \frac{\partial^2 \tilde{g}^{\text{MI}}}{\partial x^{\beta}}(x_{\mathcal{B}} + tvh) - \frac{\partial^2 \tilde{g}^{\text{MI}}}{\partial x^{\beta}}(x_{\mathcal{B}}) \right) dt v^{\beta} dv.
\end{align*}
$$

Now, using (3) since $\tilde{g}^{\text{MI}} \in C^{\alpha}([0, 1]^d)$, we have

$$
\left| \frac{\partial^\beta \tilde{g}^{\text{MI}}}{\partial x^\beta} (x_B + tvh) - \frac{\partial^\beta \tilde{g}^{\text{MI}}}{\partial x^\beta} (x_B) \right| \le C (ht \|v\|_\infty)^{\alpha - |\alpha|} \le Ch^{\alpha - |\alpha|}
$$

for all $v \in [0, 1]^d$, $x \in [0, 1]^d$, $t \in [0, 1]$, $h \in (0, \infty)$, and $\beta \in \mathbb{N}^d$ such that $|\beta| = |\alpha|$. We hence conclude that $\mathbb{E}[\tilde{g}_{n,h,i}(x)] - g(x) \le Ch^\alpha$. $\square$

A.3 Proof of Theorem 4.1

*Proof of Theorem 4.1.* Following the same lines as in the proof of Theorem 2.1 in [15], we aim at proving that

$$
\hat{T}_{n,h} - \mathbb{E}[\mathbb{E}[Y|X]^2] = \frac{1}{n} \sum_{i=1}^{n} (2Y_i - m_i)m_i - \mathbb{E}[\mathbb{E}[Y|X]^2] + o_p(n^{-1/2}). \quad (37)
$$
Proof of (17). Let $\tilde{X} = (\tilde{V}_1, \dots, \tilde{V}_d)$ and $x = (x_1, \dots, x_d)$. First of all, we compute

$$
\begin{align*}
\tilde{g}_1(x) &= \tilde{\mathbb{E}}[\tilde{g}(\tilde{V}_1, \dots, \tilde{V}_p) | \tilde{X} = x] \\
&= \frac{1}{n} \sum_{j=1}^n \tilde{\mathbb{E}} \left[ \frac{Y_j}{f_{\tilde{V}}(\tilde{V}_1, \dots, \tilde{V}_p)} \prod_{i=1}^p \tilde{K}_{h_{0,i}}(V_{i,j} - \tilde{V}_i) \middle| \tilde{X} = x \right] \\
&= \frac{1}{n} \sum_{j=1}^n Y_j \prod_{i=1}^d \frac{\tilde{K}_{h_{0,i}}(V_{i,j} - v_i)}{f_{\tilde{V}_i}(v_i)} \prod_{i=d+1}^p \tilde{\mathbb{E}} \left[ \frac{\tilde{K}_{h_{0,i}}(V_{i,j} - \tilde{V}_i)}{f_{\tilde{V}_i}(\tilde{V}_i)} \right] \\
&= \frac{1}{n} \frac{1}{f_{\tilde{X}}(x)} \sum_{j=1}^n Y_j \prod_{i=1}^d \tilde{K}_{h_{0,i}}(V_{i,j} - v_i) \prod_{i=d+1}^p \int_{B_i}^{C_i} \tilde{K}_{h_{0,i}}(V_{i,j} - v) dv \\
&= \frac{1}{n} \frac{1}{f_{\tilde{X}}(x)} \sum_{j=1}^n Y_j \prod_{i=1}^d \tilde{K}_{h_{0,i}}(V_{i,j} - v_i) \prod_{i=d+1}^p \tilde{\beta}_{i,j}
\end{align*}
$$

recalling that the support of $V_i$ is $[B_i, C_i]$, and since the inputs are independent. Then

$$
\begin{align*}
\tilde{\mathbb{E}}[\tilde{g}_1(\tilde{X})^2] &= \tilde{\mathbb{E}}\left[\left(\frac{1}{n} \frac{1}{f_{\tilde{X}}(\tilde{X})} \sum_{j=1}^n Y_j \prod_{i=1}^d \tilde{K}_{h_{0,i}}(V_{i,j} - \tilde{V}_i) \prod_{i=d+1}^p \tilde{\beta}_{i,j}\right)^2\right] \\
&= \frac{1}{n^2} \tilde{\mathbb{E}}\left[\frac{1}{f_{\tilde{X}}(\tilde{X})^2} \sum_{1 \le j \le j' \le n} Y_j Y_{j'} \prod_{i=1}^d \tilde{K}_{h_{0,i}}(V_{i,j} - \tilde{V}_i) \tilde{K}_{h_{0,i}}(V_{i,j'} - \tilde{V}_i) \prod_{i=d+1}^p \tilde{\beta}_{i,j} \tilde{\beta}_{i,j'}\right] \\
&= \frac{1}{n^2} \sum_{1 \le j \le j' \le n} Y_j Y_{j'} \left(\prod_{i=d+1}^p \tilde{\beta}_{i,j} \tilde{\beta}_{i,j'}\right) \int_{\mathcal{D}} \frac{1}{f_{\tilde{X}}(x)^2} \prod_{i=1}^d \tilde{K}_{h_{0,i}}(V_{i,j} - v_i) \tilde{K}_{h_{0,i}}(V_{i,j'} - v_i) f_{V_i}(v_i) dv_i \\
&= \frac{1}{n^2} \sum_{1 \le j \le j' \le n} Y_j Y_{j'} \left(\prod_{i=d+1}^p \tilde{\beta}_{i,j} \tilde{\beta}_{i,j'}\right) \int_{\mathcal{D}} \frac{1}{f_{\tilde{X}}(x)} \prod_{i=1}^d \tilde{K}_{h_{0,i}}(V_{i,j} - v_i) \tilde{K}_{h_{0,i}}(V_{i,j'} - v_i) dv_i
\end{align*}
$$

recalling that $\mathcal{D} = [B_1, C_1] \times \dots \times [B_d, C_d]$. $\square$

Proof of Corollary 6.1. First, note that

$$
\sqrt{n}(\hat{T}_{n,h} - \mathbb{E}[\mathbb{E}[Y|X]^2]) = \sqrt{n}(T_{n,h} - \mathbb{E}[\mathbb{E}[Y|X]^2]) + \frac{2\sqrt{n}}{n(n-1)} \sum_{1 \le j < j' \le n} \frac{Y_j Y_{j'}}{2} \\
\times \left( K_h \circ A_{X_j}(X_{j'} - X_j) \frac{f_X(X_j) - \hat{f}_X(X_j)}{\hat{f}_X(X_j) f_X(X_j)} + K_h \circ A_{X_{j'}}(X_j - X_{j'}) \frac{f_X(X_j') - \hat{f}_X(X_j')}{\hat{f}_X(X_j') f_X(X_j')} \right).
$$

The first term in the right-hand side of the previous display is handled using Theorem 3.1. Thus it remains to prove that $\sqrt{n}\mathbb{E}[|\hat{T}_{n,h} - T_{n,h}|] \underset{n\to\infty}{\to} 0$. To do so, let us study the following expectation:

$$
\begin{align*}
& \mathbb{E}\left[\left|Y_j Y_{j'} \frac{K_h \circ A_{X_{j'}}(X_j - X_{j'})}{f_X(X_j)} \frac{f_X(X_j) - \hat{f}_X(X_j)}{\hat{f}_X(X_j)}\right|^2\right] \\
\le & \mathbb{E}\left[\left(Y_j Y_{j'} \frac{K_h \circ A_{X_{j'}}(X_j - X_{j'})}{f_X(X_j)}\right)^2\right] \mathbb{E}\left[\left(\frac{f_X(X_j) - \hat{f}_X(X_j)}{\hat{f}_X(X_j)}\right)^2\right].
\end{align*}
$$

One deduces from (33) that the first in the right-hand side of the previous display is $\mathcal{O}(1/h^d)$ and the proof is concluded using (18). $\square$
The conclusion of Theorem 4.1 will then follow directly applying the standard central limit theorem for the sum of i.i.d. random variables to the right-hand side of the previous display together with Slutsky's lemma. To establish (37), we compute

$$
\begin{align*}
\hat{T}_{n,h} - \frac{1}{n} \sum_{i=1}^{n} (2Y_i - m_i)m_i &= \frac{1}{n} \sum_{i=1}^{n} [(2Y_i - \hat{m}_i)\hat{m}_i - (2Y_i - m_i)m_i] \\
&= \frac{1}{n} \sum_{i=1}^{n} [2Y_i(\hat{m}_i - m_i) + m_i^2 - \hat{m}_i^2] \\
&= \frac{1}{n} \sum_{i=1}^{n} [2(Y_i - m_i)(\hat{m}_i - m_i) - (\hat{m}_i - m_i)^2] =: I_1 - I_2.
\end{align*}
$$

**Study of I₁** Since $\hat{m}_i - m_i = \frac{\hat{g}_i - \hat{f}_i m_i}{f_i} + \frac{(f_i - \hat{f}_i)(\hat{g}_i - \hat{f}_i m_i)}{f_i \hat{f}_i}$, $I_1$ rewrites as the sum of two terms $I_{11}$ and $I_{12}$ with

$$
\begin{align*}
I_{11} &= \frac{2}{n} \sum_{i=1}^{n} (Y_i - m_i) \frac{\hat{g}_i - \hat{f}_i m_i}{f_i} \\
&= \frac{2}{n(n-1)} \sum_{i=1}^{n} \sum_{\substack{j=1, \\ j \neq i}}^{n} (Y_i - m_i)(Y_j - m_i) \frac{K_h \circ A_{X_i}(X_j - X_i)}{f_i} \\
&= \frac{2}{n(n-1)} \sum_{i=1}^{n} \sum_{\substack{j=1, \\ j \neq i}}^{n} \varepsilon_i (\varepsilon_j + m_j - m_i) \frac{K_h \circ A_{X_i}(U_{ij})}{f_i}
\end{align*}
$$

denoting the residual $Y_i - m_i$ by $\varepsilon_i$ and the difference $X_j - X_i$ by $U_{ij}$. Conditioning by $X^{(n)} = (X_1, \dots, X_n)$ and using independence, we deduce that

$$
\begin{align*}
\mathbb{E}[I_{11}] &= 2\mathbb{E}\left[\varepsilon_1(\varepsilon_2 + m_2 - m_1) \frac{K_h \circ A_{X_1}(U_{12})}{f_i}\right] \\
&= 2\mathbb{E}\left[\mathbb{E}[\varepsilon_1|X_1](\mathbb{E}[\varepsilon_2|X_2] + m_2 - m_1) \frac{K_h \circ A_{X_1}(U_{12})}{f_i}\right]
\end{align*}
$$

that cancels since $E[\varepsilon_1|X_1] = 0$ while $\mathbb{E}[I_{11}^2]$ equals

$$
\begin{align*}
& \frac{4}{n^2(n-1)^2} \sum_{i,k=1}^{n} \sum_{\substack{j,l=1, \\ j \neq i, \\ l \neq k}}^{n} \mathbb{E} \left[ \varepsilon_i \varepsilon_k (\varepsilon_j + m_j - m_i) (\varepsilon_l + m_l - m_k) \frac{K_h \circ A_{X_i}(U_{ij})}{f_i} \frac{K_h \circ A_{X_k}(U_{kl})}{f_k} \right] \\
&= \frac{4}{n^2(n-1)^2} \sum_{i=1}^{n} \sum_{\substack{j=1, \\ j \neq i}}^{n} \mathbb{E} \left[ \varepsilon_i^2 (\varepsilon_j + m_j - m_i)^2 \frac{(K_h \circ A_{X_i}(U_{ij}))^2}{f_i^2} \right] \\
&\quad + \frac{4}{n^2(n-1)^2} \sum_{i=1}^{n} \sum_{\substack{j=1, \\ j \neq i}}^{n} \mathbb{E} \left[ \varepsilon_i \varepsilon_j (\varepsilon_j + m_j - m_i) (\varepsilon_i + m_i - m_j) \frac{K_h \circ A_{X_i}(U_{ij})}{f_i} \frac{K_h \circ A_{X_j}(U_{ji})}{f_j} \right].
\end{align*}
$$

Conditioning by $X^{(n)} = (X_1, \dots, X_n)$ once again leads to

$$
\begin{align*}
\mathbb{E}[\varepsilon_i^2 (\varepsilon_j + m_j - m_i)^2 | X^{(n)}] &= \mathbb{E}[\varepsilon_i^2 \varepsilon_j^2 | X^{(n)}] + 2\mathbb{E}[\varepsilon_i^2 \varepsilon_j (m_j - m_i) | X^{(n)}] + \mathbb{E}[\varepsilon_i^2 (m_j - m_i)^2 | X^{(n)}] \\
&= \mathbb{E}[\varepsilon_i^2 | X_i] \mathbb{E}[\varepsilon_j^2 | X_j] + 2(m_j - m_i) \mathbb{E}[\varepsilon_i^2 | X_i] \mathbb{E}[\varepsilon_j | X_j] + (m_j - m_i)^2 \mathbb{E}[\varepsilon_i^2 | X_i] \\
&= \sigma_i^2 \sigma_j^2 + (m_j - m_i)^2 \sigma_i^2
\end{align*}
$$
## A.3 Construction of a kernel of order $k$ in dimension $d$

In view of the assumptions in Section 2 on the kernel $K$ that is required to be of order $\lfloor \alpha \rfloor$, we give in this section a procedure to construct kernels of a given order. Let $f_0$ be a probability density function on $\mathbb{R}$ and let $\psi_0, \dots, \psi_k$ be the associated orthonormal polynomials: for any $m = 0, \dots, k$, $\psi_m$ is a polynomial of order $m$ and

$$ \int \psi_m(x) \psi_{m'}(x) f_0(x) dx = \delta_{mm'} $$

for any $m$ and $m' = 0, \dots, k$. For any monomial $x^m$, we define $\lambda^m = (\lambda_0^m, \dots, \lambda_k^m)$ as the coefficients of its decomposition in the $f_0$-orthonormal basis $(\psi_0, \dots, \psi_k)$ that is $x^m = \sum_{i=0}^k \lambda_i^m \psi_i(x)$. Now let $c \in \mathbb{R}^{k+1}$ solving the following linear system:

$$ \sum_{i=0}^{k} \lambda_{i}^{m} c_{i} = \delta_{m0} \quad \text{for } m = 0, \dots, k. $$

We define $K_1(x) = (\sum_{i=0}^k \psi_i(x) c_i) f_0(x)$. Then, for any $m = 0, \dots, k$, one has

$$ \begin{aligned} \int_{\mathbb{R}} x^m K_1(x) dx &= \int_{\mathbb{R}} \sum_{i=0}^{k} \lambda_i^m \psi_i(x) \sum_{i'=0}^{k} c_{i'} \psi_{i'}(x) f_0(x) dx \\ &= \sum_{i,i'=0}^{k} \lambda_i^m c_{i'} \int_{\mathbb{R}} \psi_i(x) \psi_{i'}(x) f_0(x) dx \\ &= \sum_{i=0}^{k} \lambda_i^m c_i \int_{\mathbb{R}} \psi_i(x)^2 f_0(x) dx = \sum_{i=0}^{k} \lambda_i^m c_i = \delta_{m0}. \end{aligned} $$

Thus $K_1$ is of order $k$ as required and $\int_{\mathbb{R}} K_1(u) du = 1$ (case $m = 0$). To obtain a kernel of order $k$ in dimension $d$, we proceed by tensorization defining:

$$ K(x_1, \dots, x_d) = K_1(x_1) \dots K_1(x_d). $$

The tensorization allows to transfer the properties of $K_1$ to $K$.

## References

[1] A. Antoniadis. Analysis of variance on function spaces. Statistics: A Journal of Theoretical and Applied Statistics, 15(1):59–71, 1984.

[2] M. Azadkia and S. Chatterjee. A simple measure of conditional dependence. arXiv preprint arXiv:1910.12327, 2019.

[3] K. Bertin, N. Klutchnikoff, J. R. Léon, and C. Prieur. Adaptive density estimation on bounded domains under mixing conditions. Electronic Journal of Statistics, 14(1):2198 – 2237, 2020.

[4] B. Broto, F. Bachoc, and M. Depecker. Variance reduction for estimation of shapley effects and adaptation to unknown input distribution. SIAM/ASA Journal on Uncertainty Quantification, 8(2):693–716, 2020.
and

$$
\begin{align*}
& \mathbb{E}[\varepsilon_i \varepsilon_j (\varepsilon_j + m_j - m_i)(\varepsilon_i + m_i - m_j) | X^{(n)}] \\
&= \mathbb{E}[\varepsilon_i^2 \varepsilon_j^2 | X^{(n)}] + (m_j - m_i) \mathbb{E}[\varepsilon_i^2 \varepsilon_j | X^{(n)}] + (m_i - m_j) \mathbb{E}[\varepsilon_j^2 \varepsilon_i | X^{(n)}] - (m_j - m_i)^2 \mathbb{E}[\varepsilon_i \varepsilon_j | X^{(n)}] \\
&= \sigma_i^2 \sigma_j^2.
\end{align*}
$$

Hence

$$
\begin{align*}
\mathbb{E}[I_{11}^2] &= \frac{4}{n^2(n-1)^2} \sum_{i=1}^{n} \sum_{\substack{j=1, \\ j \neq i}}^{n} \mathbb{E} \left[ (\sigma_i^2 \sigma_j^2 + (m_j - m_i)^2 \sigma_i^2) \frac{(K_h \circ A_{X_i}(U_{ij}))^2}{f_i^2} \right] \\
&\quad + \frac{4}{n^2(n-1)^2} \sum_{i=1}^{n} \sum_{\substack{j=1, \\ j \neq i}}^{n} \mathbb{E} \left[ \sigma_i^2 \sigma_j^2 \frac{K_h \circ A_{X_i}(U_{ij})}{f_i} \frac{K_h \circ A_{X_j}(U_{ji})}{f_j} \right].
\end{align*}
$$

By Assumptions (A2), (A3), (A6), and the fact that *m* is bounded from above, we
conclude that

$$
\mathbb{E}[I_{11}^2] = O\left(\frac{1}{n^2 h^{2d}}\right) = o\left(\frac{1}{n}\right)
$$

since $nh^{2d} \to \infty$ by assumption.

Let us turn to the computation of $I_{12}$ defined as follows

$$
I_{12} = \frac{2}{n} \sum_{i=1}^{n} (Y_i - m_i) \frac{(f_i - \hat{f}_i)(\hat{g}_i - \hat{f}_i m_i)}{f_i \hat{f}_i}.
$$

By (15) and since $\hat{g}_1$ does not depend on $Y_1$,

$$
\begin{align*}
|\mathbb{E}[I_{12}]| &\le \frac{C}{n} \sum_{i=1}^{n} \mathbb{E}_{Y_1, \dots, Y_n} \left[ \frac{|\varepsilon_i| |f_i - \hat{f}_i| (|g_i - \hat{g}_i| + |f_i - \hat{f}_i||m_i|)}{f_i} \right] \\
&= C \mathbb{E}_{Y_1, \dots, Y_n} \left[ |\varepsilon_1| \frac{|f_1 - \hat{f}_1| (|g_1 - \hat{g}_1| + |f_1 - \hat{f}_1||m_1|)}{f_1} \right] \\
&= C \mathbb{E}_{X_1, \dots, X_n} \left[ \frac{|f_1 - \hat{f}_1|}{f_1} \mathbb{E}_{Y_1, \dots, Y_n} [ |\varepsilon_1| (|g_1 - \hat{g}_1| + |f_1 - \hat{f}_1||m_1|) |X^{(n)}] \right] \\
&= C \mathbb{E}_{X_1, \dots, X_n} \left[ \frac{|f_1 - \hat{f}_1|}{f_1} \mathbb{E}_{Y_1} [ |\varepsilon_1| |X_1| (\mathbb{E}_{Y_2, \dots, Y_n} [ |g_1 - \hat{g}_1| |X^{(n)}| + |f_1 - \hat{f}_1||m_1|) ] \right] \\
&\le C \mathbb{E}_{X_1, \dots, X_n} \left[ \frac{|f_1 - \hat{f}_1|}{f_1} \sigma_1 (\mathbb{E}_{Y_2, \dots, Y_n} [ |g_1 - \hat{g}_1| |X^{(n)}| + |f_1 - \hat{f}_1||m_1|) \right] \\
&= C \mathbb{E}_{X_1, \dots, X_n} \left[ \frac{|f_1 - \hat{f}_1|}{f_1} \sigma_1 (|g_1 - \hat{g}_1| + |f_1 - \hat{f}_1||m_1|) \right] \\
&= C \mathbb{E}_{X_2, \dots, X_n} \left[ \int_{[0,1]^d} |f_X(x) - \hat{f}_{n,h,1}(x)| \sigma(x) (|g(x) - \hat{g}_{n,h,1}(x)| + |f_X(x) - \hat{f}_{n,h,1}(x)||m(x)|) dx \right] \\
&\le C \mathbb{E}_{X_2, \dots, X_n} \left[ \int_{[0,1]^d} |f_X(x) - \hat{f}_{n,h,1}(x)| |g(x) - \hat{g}_{n,h,1}(x)| dx \right] \\
&\quad + C \mathbb{E}_{X_2, \dots, X_n} \left[ \int_{[0,1]^d} |f_X(x) - \hat{f}_{n,h,1}(x)| |f_X(x) - \hat{f}_{n,h,1}(x)||m(x)| dx \right].
\end{align*}
$$
[5] M. Carrascal, M. Puigcerver, and P. Puig. Sensitivity of gaussian plume model to dispersion specifications. *Theoretical and Applied Climatology*, 48(2-3):147-157, 1993.

[6] S. Chatterjee. A new coefficient of correlation. *Journal of the American Statistical Association*, pages 1-26, 2020.

[7] R. Coleman. *Calculus on normed vector spaces*. Springer Science & Business Media, 2012.

[8] R. Cukier, H. Levine, and K. Shuler. Nonlinear sensitivity analysis of multiparameter model systems. *Journal of computational physics*, 26(1):1-42, 1978.

[9] S. Da Veiga and F. Gamboa. Efficient estimation of sensitivity indices. *Journal of Nonparametric Statistics*, 25(3):573-595, 2013.

[10] S. Da Veiga, F. Gamboa, B. Iooss, and C. Prieur. *Basics and Trends in Sensitivity Analysis: Theory and Practice in R*. SIAM, 2021.

[11] S. Da Veiga, F. Wahl, and F. Gamboa. Local polynomial estimation for sensitivity analysis on models with correlated inputs. *Technometrics*, 51(4):452-463, 2009.

[12] B. Delyon and F. Portier. Integral approximation by kernel smoothing. *Bernoulli*, 22(4):2177-2208, 2016.

[13] A. Dembo and O. Zeitouni. *Large deviations techniques and applications*, volume 38 of *Applications of Mathematics (New York)*. Springer-Verlag, New York, second edition, 1998.

[14] L. Devroye, P. G. Ferrario, L. Györfi, and H. Walk. Strong universal consistent estimate of the minimum mean squared error. *Empirical Inference: Festschrift in Honor of Vladimir N. Vapnik*, pages 143-160, 2013.

[15] L. Devroye, L. Györfi, G. Lugosi, and H. Walk. A nearest neighbor estimate of the residual variance. *Electronic Journal of Statistics*, 12(1):1752-1778, 2018.

[16] L. Devroye, D. Schäfer, L. Györfi, and H. Walk. The estimation problem of minimum mean squared error. *Statistics & Decisions*, 21(1):15-28, 2003.

[17] F. Gamboa, P. Gremaud, T. Klein, and A. Lagnoux. Global sensitivity analysis: A novel generation of mighty estimators based on rank statistics. *Bernoulli*, 28(4):2345-2374, 2022.

[18] F. Gamboa, A. Janon, T. Klein, A. Lagnoux, and C. Prieur. Statistical inference for Sobol Pick-Freeze Monte Carlo method. *Statistics*, 50(4):881-902, 2016.

[19] F. Gamboa, T. Klein, and A. Lagnoux. Sensitivity analysis based on Cramér von Mises distance. *SIAM/ASA Journal on Uncertainty Quantification*, 6(2):522-548, Apr. 2018.

[20] L. Gilquin, E. Arnaud, C. Prieur, and A. Janon. Making the best use of permutations to compute sensitivity indices with replicated orthogonal arrays. *Reliability Engineering & System Safety*, 187:28-39, 2019.
By applying twice Cauchy-Schwartz inequality, the first term is bounded from above by

$$
\begin{align*}
C \mathbb{E}_{\substack{X_2, \dots, X_n \\ Y_2, \dots, Y_n}} & \left[ \left( \int_{[0,1]^d} (f_X(x) - \hat{f}_{n,h,1}(x))^2 dx \right)^{1/2} \left( \int_{[0,1]^d} (g(x) - \hat{g}_{n,h,1}(x))^2 dx \right)^{1/2} \right] \\
& \le C \mathbb{E}_{\substack{X_2, \dots, X_n \\ Y_2, \dots, Y_n}} \left[ \int_{[0,1]^d} (f_X(x) - \hat{f}_{n,h,1}(x))^2 dx \right]^{1/2} \mathbb{E}_{\substack{X_2, \dots, X_n \\ Y_2, \dots, Y_n}} \left[ \int_{[0,1]^d} (g(x) - \hat{g}_{n,h,1}(x))^2 dx \right]^{1/2} \\
& = o(n^{-1/2})
\end{align*}
$$

using (14) and (30) while the second term is $o(n^{-1/2})$.

**Study of I₂** Recall that

$$
I_2 = \frac{1}{n} \sum_{i=1}^{n} (\hat{m}_i - m_i)^2 = \frac{1}{n} \sum_{i=1}^{n} \frac{(\hat{g}_i - \hat{f}_i m_i)^2}{\hat{f}_i^2} = O_P\left(\frac{1}{n} \sum_{i=1}^{n} (\hat{g}_i - \hat{f}_i m_i)^2\right)
$$

by (15). Then,

$$
\begin{align*}
\mathbb{E}[I_2] &\le C \mathbb{E}_{\substack{X_2, \dots, X_n \\ Y_2, \dots, Y_n}} \left[ \int_{[0,1]^d} (\hat{g}_{n,h,i}(x) - \hat{f}_{n,h,i}(x)m(x))^2 f_X(x)dx \right] \\
&= C \mathbb{E}_{\substack{X_2, \dots, X_n \\ Y_2, \dots, Y_n}} \left[ \int_{[0,1]^d} \left[ (\hat{g}_{n,h,i}(x) - g(x)) + (f_X(x)m(x) - \hat{f}_{n,h,i}(x)m(x)) \right]^2 f_X(x)dx \right] \\
&\le C \mathbb{E}_{\substack{X_2, \dots, X_n \\ Y_2, \dots, Y_n}} \left[ \int_{[0,1]^d} (\hat{g}_{n,h,i}(x) - g(x))^2 f_X(x)dx \right] + C \mathbb{E} \left[ \int_{[0,1]^d} (\hat{f}_{n,h,i}(x) - f(x))^2 m(x)^2 f_X(x)dx \right] \\
&= o(n^{-1/2})
\end{align*}
$$

as $f_X$ is continuous thus bounded on $[0, 1]^d$, as $m$ is also bounded on $[0, 1]^d$ (recall that it
follows from the boundedness of $\sigma^2$) and from the results stated in (14) and (30).

As for item (ii), the proof is similar except that for any $i \neq j$, $K_h \circ A_{X_i}(U_{ij})$ is replaced
by $\sum_{a \in \{-1,0,1\}^d} \prod_{k=1}^d K_h(M^{a_k}(X_{k,j}) - X_i)$. $\square$

A.4 Proof of the remaining results

*Proof of Proposition* 4.2. We denote by *P* the distribution of (*X*, *Y*) and we introduce

$$
\psi(P) = \mathbb{E}[\mathbb{E}[Y | X]^2].
$$

The influence efficient function of $\psi$ at $P$, as stated in [15], is given by $\tilde{\psi}_P(x, y) = (2y - m(x))m(x) - \mathbb{E}[Ym(X)]$ (see [31] for explicit computations). Moreover, we deduce from (37) that

$$
\hat{T}_{n,h} = \psi(P) + \frac{1}{n} \sum_{i=1}^{n} \tilde{\psi}_{P}(X_i, Y_i) + o_{\mathbb{P}}(n^{-1/2})
$$

and, by [51, Condition (25.22)], $\hat{T}_{n,h}$ is asymptotically efficient and so is $\tilde{T}_{n,h}$. $\square$

*Proof of Corollary* 4.3. In view of the definition of the asymptotic efficiency (see, e.g., [51, Lemma 25.23] or [31, Definition 2.4], the asymptotic efficiency ensures the asymptotic normality. Thus we only need to prove the asymptotic efficiency. To do so, it suffices to use the asymptotic efficiency of $\hat{T}_{n,h}$ and $\tilde{T}_{n,h}$ in Proposition 4.2, the asymptotic efficiency
[21] E. Giné and D. M. Mason. Uniform in bandwidth estimation of integral functionals of the density function. *Scandinavian Journal of Statistics*, 35(4):739–761, 2008.

[22] E. Giné, R. Nickl, et al. A simple adaptive estimator of the integrated square of a density. *Bernoulli*, 14(1):47–61, 2008.

[23] T. Goda. Computing the variance of a conditional expectation via non-nested Monte Carlo. *Operations Research Letters*, 45(1):63 – 67, 2017.

[24] L. Györfi and H. Walk. On the asymptotic normality of an estimate of a regression functional. *J. Mach. Learn. Res.*, 16:1863–1877, 2015.

[25] M. B. Heredia, C. Prieur, and N. Eckert. Nonparametric estimation of aggregated sobol'indices: application to a depth averaged snow avalanche model. *Reliability Engineering & System Safety*, 212:107422, 2021.

[26] W. Hoeffding. A class of statistics with asymptotically normal distribution. *Ann. Math. Statistics*, 19:293–325, 1948.

[27] A. Janon, T. Klein, A. Lagnoux, M. Nodet, and C. Prieur. Asymptotic normality and efficiency of two Sobol index estimators. *ESAIM: Probability and Statistics*, 18:342–364, 1 2014.

[28] S. Kucherenko and S. Song. Different numerical estimators for main effect global sensitivity indices. *Reliability Engineering & System Safety*, 165:222–238, 2017.

[29] E. Liitiäinen, F. Corona, and A. Lendasse. On nonparametric residual variance estimation. *Neural Processing Letters*, 28:155–167, 2008.

[30] E. Liitiäinen, F. Corona, and A. Lendasse. Residual variance estimation using a nearest neighbor statistic. *Journal of Multivariate Analysis*, 101(4):811–823, 2010.

[31] Z. Lin and F. Han. Limit theorems of chatterjee's rank correlation. *arXiv preprint arXiv:2204.08031*, 2022.

[32] J.-M. Loubes, C. Marteau, and M. Solís. Rates of convergence in conditional covariance matrix with nonparametric entries estimation. *Communications in Statistics-Theory and Methods*, 49(18):4536–4558, 2020.

[33] A. B. Owen. Better estimation of small sobol' sensitivity indices. *ACM Trans. Model. Comput. Simul.*, 23(2):11:1–11:17, May 2013.

[34] K. Pearson. On the partial correlation ratio. *Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character*, 91(632):492–498, 1915.

[35] V. H. Peña and E. Giné. *Decoupling: From dependence to independence*. Springer, 1999.

[36] E. Plischke. An effective algorithm for computing global sensitivity indices (easi). *Reliability Engineering & System Safety*, 95(4):354–360, 2010.
of the empirical mean $\bar{Y}_n$ and the empirical variance $\bar{Y}^2_n$ (see, e.g., [28]) together with the efficiency in product space [51, Theorem 25.50] to conclude to the joint asymptotic efficiency of $(\hat{T}_{n,h}, \bar{Y}_n, \bar{Y}^2_n)$ and $(\hat{T}_{n,h}, \bar{Y}_n, \bar{Y}^2_n)$. Finally, we follow the same lines as in the proof of [28, Proposition 2.5] using the efficiency and delta method [51, Theorem 25.47] to get the required result. It remains to perform easy computations to get the expression of the asymptotic variance

$$ \mathrm{Var}(Y)^2\sigma^2 = \mathrm{Var}((2Y - m(X))m(X)) + S^X(\mathrm{Var}(Y^2)S^X - 2\mathrm{Cov}(m(X)^2, Y^2)) \\ + 4\mathbb{E}[Y](S^X - 1)[\mathbb{E}[m(X)^2Y] - \mathbb{E}[Y]\mathbb{E}[Y^2] - \mathrm{Cov}(Y, Y^2)S^X], \quad (38) $$

from which we derive simpler expressions

* when $\mathbb{E}[Y] = 0$ or $S^X = 1$:

$$ \mathrm{Var}(Y)^2\sigma^2 = \mathrm{Var}((2Y - m(X))m(X)) + S^X(\mathrm{Var}(Y^2)S^X - 2\mathrm{Cov}(m(X)^2, Y^2)), $$

* when $S^X = 0$:

$$ \mathrm{Var}(Y)^2\sigma^2 = \mathrm{Var}((2Y - m(X))m(X)) + 4\mathbb{E}[Y](S^X - 1)[\mathbb{E}[m(X)^2Y] - \mathbb{E}[Y]\mathbb{E}[Y^2]]. $$

The proof is then complete. $\square$
[37] E. Plischke and E. Borgonovo. Fighting the curse of sparsity: Probabilistic sensitivity measures from cumulative distribution functions. Risk Analysis, 40(12):2639-2660, 2020.

[38] N. Rachdi, J.-C. Fort, and T. Klein. Stochastic inverse problem with noisy simulator-application to aeronautical model. Annales de la Faculté des Sciences de Toulouse, 6, 21:593-622, 2012.

[39] S. Razavi, A. Jakeman, A. Saltelli, C. Prieur, B. Iooss, E. Borgonovo, E. Plischke, S. Lo Piano, T. Iwanaga, W. Becker, S. Tarantola, J. H. Guillaume, J. Jakeman, H. Gupta, N. Melillo, G. Rabitti, V. Chabridon, Q. Duan, X. Sun, S. Smith, R. Sheikholeslami, N. Hosseini, M. Asadzadeh, A. Puy, S. Kucherenko, and H. Maier. The Future of Sensitivity Analysis: An essential discipline for systems modeling and policy support. Environmental Modelling and Software, 137:104954, Mar. 2021.

[40] A. Saltelli, K. Chan, and E. Scott. Sensitivity analysis. Wiley Series in Probability and Statistics. John Wiley & Sons, Ltd., Chichester, 2000.

[41] I. M. Sobol. Sensitivity estimates for nonlinear mathematical models. Math. Modeling Comput. Experiment, 1(4):407-414 (1995), 1993.

[42] I. M. Sobol. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Mathematics and Computers in Simulation, 55(1-3):271-280, 2001.

[43] M. Solís. Non-parametric estimation of the first-order sobol indices with bootstrap bandwidth. Communications in Statistics-Simulation and Computation, 50(9):2497-2512, 2021.

[44] J. M. Stockie. The mathematics of atmospheric dispersion modeling. Siam Review, 53(2):349-372, 2011.

[45] B. Sudret. Global sensitivity analysis using polynomial chaos expansions. Reliability Engineering & System Safety, 93(7):964-979, 2008.

[46] S. Tarantola, D. Gatelli, and T. Mara. Random balance designs for the estimation of first order global sensitivity indices. Reliability Engineering & System Safety, 91(6):717-727, 2006.

[47] J.-Y. Tissot and C. Prieur. A randomized orthogonal array-based procedure for the estimation of first-and second-order sobol'indices. Journal of Statistical Computation and Simulation, 85(7):1358-1381, 2015.

[48] A. W. van der Vaart. Asymptotic statistics, volume 3 of Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, Cambridge, 1998.

[49] L.-X. Zhu and K.-T. Fang. Asymptotics for kernel estimate of sliced inverse regression. The Annals of Statistics, 24(3):1053-1068, 1996.
# B Additional numerical experiments

## B.1 Influence of $\epsilon$ in [15]

We illustrate numerically that the choice of the $\epsilon$ tuning parameter of the estimator proposed in [15] is very sensitive, thus limiting its practical use as opposed to our mirror-type estimator. We consider Example 3.2 from [15] and test $\epsilon = 10^{-1}, 10^{-2}, 10^{-3}$. The comparison with our estimator with a kernel of order 2 is given in Figure 8. When $\epsilon$ is equal to $10^{-3}$, the performance of both estimators are similar. However in other cases, the bias of [15] can be very large. Since in practice such an estimation problem is unsupervised, the tuning of $\epsilon$ seems highly difficult and the non-robustness of the final estimator with respect to this parameter limits its practical use.

<figure><img src="image_11.png" /><figcaption>Figure 8: Comparison of our mirror-type estimator with the estimator of [15] for different values of $\epsilon$. The reference value of the index is represented with a gray line.</figcaption></figure>

## B.2 Analytical test functions with smaller sample size

For completeness, we include below new results on the Bratley and g-Sobol functions in dimension $p = 5$ with smaller sample sizes. As expected, the bias is larger when $n$ decreases.
<figure><img src="image_12.png" /><figcaption>Figure 9: Estimators for total indices of the Bratley function from $n = 100$ to $n = 1000$. The reference value is represented with a gray line.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 10: Estimators for total indices of the g-Sobol function from $n = 100$ to $n = 1000$. The reference value is represented with a gray line.</figcaption></figure>
## Acknowledgement

The authors would like to thank Elmar Plischke for pointing out some very relevant references related to the present work. They are also grateful to the reviewers for their valuable comments that allowed to improve this article. Support from the ANR-3IA Artificial and Natural Intelligence Toulouse Institute, the ANR GATSBII (ANR-24-CE23-6645) and the consortium in Applied Mathematics CIROQUO-2, gathering partners in technological research and academia in the development of advanced methods for Computer Experiments, are gratefully acknowledged.

## References

[1] K. Behnen, M. Hušková, and G. Neuhaus. Rank estimators of scores for testing independence. *Statistics & Risk Modeling*, 3(3-4):239–262, 1985.

[2] C. Bénard, S. Da Veiga, and E. Scornet. Mean decrease accuracy for random forests: inconsistency, and a practical solution via the Sobol-mda. *Biometrika*, 109(4):881–900, 2022.

[3] K. Bertin, N. Klutchnikoff, J. R. Léon, and C. Prieur. Adaptive density estimation on bounded domains under mixing conditions. *Electronic Journal of Statistics*, 14(1):2198 – 2237, 2020.

[4] P. J. Bickel, C. A. Klaassen, P. J. Bickel, Y. Ritov, J. Klaassen, J. A. Wellner, and Y. Ritov. *Efficient and adaptive estimation for semiparametric models*, volume 4. Springer, 1993.

[5] E. Borgonovo, E. Plischke, and C. Prieur. Total effects with constrained features. *Statistics and Computing*, 34(2):87, 2024.

[6] B. Broto, F. Bachoc, and M. Depecker. Variance reduction for estimation of shapley effects and adaptation to unknown input distribution. *SIAM/ASA Journal on Uncertainty Quantification*, 8(2):693–716, 2020.

[7] S. Chatterjee. A new coefficient of correlation. *Journal of the American Statistical Association*, pages 1–26, 2020.

[8] R. Coleman. *Calculus on normed vector spaces*. Springer Science & Business Media, 2012.

[9] S. Da Veiga and F. Gamboa. Efficient estimation of sensitivity indices. *Journal of Nonparametric Statistics*, 25(3):573–595, 2013.

[10] S. Da Veiga, F. Gamboa, B. Iooss, and C. Prieur. *Basics and Trends in Sensitivity Analysis: Theory and Practice in R*. SIAM, 2021.

[11] S. Da Veiga, F. Wahl, and F. Gamboa. Local polynomial estimation for sensitivity analysis on models with correlated inputs. *Technometrics*, 51(4):452–463, 2009.
[12] L. Devroye, P. G. Ferrario, L. Györfi, and H. Walk. Strong universal consistent estimate of the minimum mean squared error. *Empirical Inference: Festschrift in Honor of Vladimir N. Vapnik*, pages 143–160, 2013.

[13] L. Devroye, L. Györfi, G. Lugosi, and H. Walk. A nearest neighbor estimate of the residual variance. *Electronic Journal of Statistics*, 12(1):1752–1778, 2018.

[14] L. Devroye, D. Schäfer, L. Györfi, and H. Walk. The estimation problem of minimum mean squared error. *Statistics & Decisions*, 21(1):15–28, 2003.

[15] K. Doksum and A. Samarov. Nonparametric estimation of global functionals and a measure of the explanatory power of covariates in regression. *The Annals of Statistics*, pages 1443–1473, 1995.

[16] F. Gamboa, P. Gremaud, T. Klein, and A. Lagnoux. Global sensitivity analysis: A novel generation of mighty estimators based on rank statistics. *Bernoulli*, 28(4):2345–2374, 2022.

[17] F. Gamboa, A. Janon, T. Klein, A. Lagnoux, and C. Prieur. Statistical inference for Sobol' Pick-Freeze Monte Carlo method. *Statistics*, 50(4):881–902, 2016.

[18] I. Gijbels, K. Herrmann, and D. Sznajder. Flexible and dynamic modeling of dependencies via copulas. In *Modeling and Stochastic Learning for Forecasting in High Dimensions*, pages 117–146. Springer, 2015.

[19] I. Gijbels and J. Mielniczuk. Estimating the density of a copula function. *Communications in Statistics-Theory and Methods*, 19(2):445–464, 1990.

[20] L. Gilquin, E. Arnaud, C. Prieur, and A. Janon. Making the best use of permutations to compute sensitivity indices with replicated orthogonal arrays. *Reliability Engineering & System Safety*, 187:28–39, 2019.

[21] E. Giné and A. Guillou. On consistency of kernel density estimators for randomly censored data: rates holding uniformly over adaptive intervals. In *Annales de l'IHP Probabilités et statistiques*, volume 37, pages 503–522, 2001.

[22] E. Giné, R. Nickl, et al. A simple adaptive estimator of the integrated square of a density. *Bernoulli*, 14(1):47–61, 2008.

[23] L. Györfi and H. Walk. On the asymptotic normality of an estimate of a regression functional. *J. Mach. Learn. Res.*, 16:1863–1877, 2015.

[24] B. E. Hansen. Exact mean integrated squared error of higher order kernel estimators. *Econometric Theory*, 21(6):1031–1057, 2005.

[25] M. B. Heredia, C. Prieur, and N. Eckert. Nonparametric estimation of aggregated Sobol' indices: application to a depth averaged snow avalanche model. *Reliability Engineering & System Safety*, 212:107422, 2021.

[26] W. Hoeffding. A class of statistics with asymptotically normal distribution. *Ann. Math. Statistics*, 19:293–325, 1948.
[27] P. Hominal and P. Deheuvels. Estimation non paramétrique de la densité compté-tenu d'informations sur le support. *Revue de statistique appliquée*, 27(3):47-68, 1979.

[28] A. Janon, T. Klein, A. Lagnoux, M. Nodet, and C. Prieur. Asymptotic normality and efficiency of two Sobol index estimators. *ESAIM: Probability and Statistics*, 18:342-364, 1 2014.

[29] M. Kendall and A. Stuart. *The advanced theory of statistics, Vol. 2: Inference and relationship*, 2nd edition. Hafner Publishing Company, 1967.

[30] J. Kim, J. Shin, A. Rinaldo, and L. Wasserman. Uniform convergence rate of the kernel density estimator adaptive to intrinsic volume dimension. In *International Conference on Machine Learning*, pages 3398-3407. PMLR, 2019.

[31] T. Klein, A. Lagnoux, T. M. N. Nguyen, and P. Rochet. Efficient influence functions for Sobol' indices under two designs of experiments. *arXiv preprint arXiv:2407.15468*, 2024.

[32] T. Klein and P. Rochet. Efficiency of the averaged rank-based estimator for first order Sobol' index inference. *Statistics & Probability Letters*, 207:110015, 2024.

[33] E. Liitiäinen, F. Corona, and A. Lendasse. On nonparametric residual variance estimation. *Neural Processing Letters*, 28:155-167, 2008.

[34] E. Liitiäinen, F. Corona, and A. Lendasse. Residual variance estimation using a nearest neighbor statistic. *Journal of Multivariate Analysis*, 101(4):811-823, 2010.

[35] H. Liu, L. Wasserman, and J. Lafferty. Exponential concentration for mutual information estimation with application to forests. *Advances in Neural Information Processing Systems*, 25, 2012.

[36] E. A. Nadaraya. On estimating regression. *Theory of Probability & Its Applications*, 9(1):141-142, 1964.

[37] W. K. Newey, F. Hsieh, and J. Robins. Undersmoothing and bias corrected functional estimation. Work in progress. Cambridge, Mass.: Massachusetts Institute of Technology, 1998.

[38] A. B. Owen. Sobol' indices and Shapley value. *SIAM/ASA Journal on Uncertainty Quantification*, 2(1):245-251, 2014.

[39] J. Pfanzagl. Lecture notes in statistics. *Contributions to a general asymptotic statistical theory*, 13:11-15, 1982.

[40] E. Plischke and E. Borgonovo. Fighting the curse of sparsity: Probabilistic sensitivity measures from cumulative distribution functions. *Risk Analysis*, 40(12):2639-2660, 2020.

[41] L. Pujol. Nonparametric estimation of a multivariate density under Kullback-Leibler loss with ISDE. *arXiv preprint arXiv:2205.03199*, 2022.
[42] S. Razavi, A. Jakeman, A. Saltelli, C. Prieur, B. Iooss, E. Borgonovo, E. Plischke, S. Lo Piano, T. Iwanaga, W. Becker, S. Tarantola, J. H. Guillaume, J. Jakeman, H. Gupta, N. Melillo, G. Rabitti, V. Chabridon, Q. Duan, X. Sun, S. Smith, R. Sheikholeslami, N. Hosseini, M. Asadzadeh, A. Puy, S. Kucherenko, and H. Maier. The Future of Sensitivity Analysis: An essential discipline for systems modeling and policy support. *Environmental Modelling and Software*, 137:104954, Mar. 2021.

[43] A. Saltelli, K. Chan, and E. Scott. *Sensitivity analysis*. Wiley Series in Probability and Statistics. John Wiley & Sons, Ltd., Chichester, 2000.

[44] E. F. Schuster. Incorporating support constraints into nonparametric estimators of densities. *Communications in Statistics-Theory and methods*, 14(5):1123–1136, 1985.

[45] I. M. Sobol'. Sensitivity estimates for nonlinear mathematical models. *Math. Modeling Comput. Experiment*, 1(4):407–414 (1995), 1993.

[46] I. M. Sobol'. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. *Mathematics and Computers in Simulation*, 55(1-3):271–280, 2001.

[47] M. Solís. Non-parametric estimation of the first-order Sobol' indices with bootstrap bandwidth. *Communications in Statistics-Simulation and Computation*, 50(9):2497–2512, 2021.

[48] B. Sudret. Global sensitivity analysis using polynomial chaos expansions. *Reliability Engineering & System Safety*, 93(7):964–979, 2008.

[49] J.-Y. Tissot and C. Prieur. A randomized orthogonal array-based procedure for the estimation of first-and second-order Sobol' indices. *Journal of Statistical Computation and Simulation*, 85(7):1358–1381, 2015.

[50] A. B. Tsybakov. Nonparametric estimators. *Introduction to Nonparametric Estimation*, pages 1–76, 2009.

[51] A. W. van der Vaart. *Asymptotic statistics*, volume 3 of *Cambridge Series in Statistical and Probabilistic Mathematics*. Cambridge University Press, Cambridge, 1998.

[52] G. S. Watson. Smooth regression analysis. *Sankhyā: The Indian Journal of Statistics, Series A*, pages 359–372, 1964.