M. Fréchet. Sur le coefficient, dit de corrélation et sur la corrélation en général. *Revue de l'Institut Interna-
tional de Statistique*, 1(4):16-23, 1934.

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.

D. Gatelli, S. Kucherenko, M. Ratto, and S. Tarantola. Calculating first-order sensitivity measures: A benchmark of some recent methodologies. *Reliability Engineering&System Safety*, 94:1212-1219, 2009.

R. Genuer, V. Michel, E. Eger, and B. Thirion. Random forests based feature selection for decoding fmri data. In *Proceedings Comstat*, volume 267, pages 1-8, 2010.

L. Gilquin, C. Prieur, and E. Arnaud. Replication procedure for grouped sobol'indices estimation in depen-
dent uncertainty spaces. *Information and Inference: A Journal of the IMA*, 4(4):354-379, 2015.

T. Goda. A simple algorithm for global sensitivity analysis with Shapley effects. *Reliability Engineering&System Safety*, 213:107702, 2021.

J. Hart and P. A. Gremaud. An approximation theoretic perspective of Sobol' indices with dependent variables. *Int. J. Uncertainty Quantification*, 8(6):483-493, 2018.

R. Helmers. On the Edgeworth expansion and the bootstrap approximation for a Studentized *U*-statistic.
*The Annals of Statistics*, 19:470–484, 1991.

W. Hoeffding. A class of statistics with asymptotically normal distribution. *The Annals of Mathematical Statistics*, 19(3):293–325, 1948.

T. Homma and A. Saltelli. Importance measures in global sensitivity analysis of nonlinear models. *Reliability Engineering&System Safety*, 52(1):1-17, 1996.

G. Hooker, L. Mentch, and S. Zhou. Unrestricted permutation forces extrapolation: variable importance requires at least one more model, or there is no free variable importance. *Statistics and Computing*, 31 (6):82:1-16, 2021.

R. L. Iman and W. J. Conover. A distribution-free approach to inducing rank correlation among input variables. *Communications in Statistics - Simulation and Computation*, 11(3):311-334, 1982.

J. Jacques, C. Lavergne, and N. Devictor. Sensitivity analysis in presence of model uncertainty and correlated inputs. *Reliability Engineering&System Safety*, 91(10-11):1126-1134, 2006.

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, 2014.

M. J. W. Jansen. Analysis of variance designs for model output. *Computer Physics Communications*, 117 (1-2):35-43, 1999.

M. J. W. Jansen, W. A. H. Rossing, and R. A. Daamen. Monte Carlo estimation of uncertainty contributions from several independent multivariate sources. *Predictability and nonlinear modelling in natural sciences and economics*, pages 334–343, 1994.

H. Joe. *Dependence Modeling with Copulas*. CRC Press, Boca Raton, 2014.

H. Knothe. Contributions to the theory of convex bodies. *Michigan Math. J.*, 4(1):39–52, 1957.

S. Kucherenko, S. Tarantola, and P. Annoni. Estimation of global sensitivity indices for models with depen-
dent variables. *Computer Physics Communications*, 183(4):937-946, 2012.

S. Kucherenko, O. V. Klymenko, and N. Shah. Sobol' indices for problems defined in non-rectangular domains. *Reliability Engineering&System Safety*, 167:218-231, 2017.

S. N. Lahiri. *Resampling Methods for Dependent Data*. Springer Science+Business Media, New York, 2003.

J. Lei, M. G'Sell, A. Rinaldo, R. J. Tibshirani, and L. Wasserman. Distribution-free predictive inference for regression. *Journal of the American Statistical Association*, 113:1094–1111, 2018.

G. Li and H. Rabitz. Relationship between sensitivity indices defined by variance- and covariance-based methods. *Reliability Engineering&System Safety*, 167:136-157, 2017.

S. Lundberg and S.-I. Lee. A unified approach to interpreting model predictions. preprint, arXiv, 2017.
https://arxiv.org/abs/1705.07874.
M. Fréchet. Sur le coefficient, dit de corrélation et sur la corrélation en général. *Revue de l'Institut Interna-
tional de Statistique*, 1(4):16-23, 1934.

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.

D. Gatelli, S. Kucherenko, M. Ratto, and S. Tarantola. Calculating first-order sensitivity measures: A benchmark of some recent methodologies. *Reliability Engineering&System Safety*, 94:1212-1219, 2009.

R. Genuer, V. Michel, E. Eger, and B. Thirion. Random forests based feature selection for decoding fmri data. In *Proceedings Comstat*, volume 267, pages 1-8, 2010.

L. Gilquin, C. Prieur, and E. Arnaud. Replication procedure for grouped sobol'indices estimation in depen-
dent uncertainty spaces. *Information and Inference: A Journal of the IMA*, 4(4):354-379, 2015.

T. Goda. A simple algorithm for global sensitivity analysis with Shapley effects. *Reliability Engineering&System Safety*, 213:107702, 2021.

J. Hart and P. A. Gremaud. An approximation theoretic perspective of Sobol' indices with dependent variables. *Int. J. Uncertainty Quantification*, 8(6):483-493, 2018.

R. Helmers. On the Edgeworth expansion and the bootstrap approximation for a Studentized *U*-statistic.
*The Annals of Statistics*, 19:470–484, 1991.

W. Hoeffding. A class of statistics with asymptotically normal distribution. *The Annals of Mathematical Statistics*, 19(3):293–325, 1948.

T. Homma and A. Saltelli. Importance measures in global sensitivity analysis of nonlinear models. *Reliability Engineering&System Safety*, 52(1):1-17, 1996.

G. Hooker, L. Mentch, and S. Zhou. Unrestricted permutation forces extrapolation: variable importance requires at least one more model, or there is no free variable importance. *Statistics and Computing*, 31 (6):82:1–16, 2021.

R. L. Iman and W. J. Conover. A distribution-free approach to inducing rank correlation among input variables. *Communications in Statistics - Simulation and Computation*, 11(3):311-334, 1982.

J. Jacques, C. Lavergne, and N. Devictor. Sensitivity analysis in presence of model uncertainty and correlated inputs. *Reliability Engineering&System Safety*, 91(10-11):1126-1134, 2006.

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, 2014.

M. J. W. Jansen. Analysis of variance designs for model output. *Computer Physics Communications*, 117 (1-2):35-43, 1999.

M. J. W. Jansen, W. A. H. Rossing, and R. A. Daamen. Monte Carlo estimation of uncertainty contributions from several independent multivariate sources. *Predictability and nonlinear modelling in natural sciences and economics*, pages 334–343, 1994.

H. Joe. *Dependence Modeling with Copulas*. CRC Press, Boca Raton, 2014.

H. Knothe. Contributions to the theory of convex bodies. *Michigan Math. J.*, 4(1):39–52, 1957.

S. Kucherenko, S. Tarantola, and P. Annoni. Estimation of global sensitivity indices for models with depen-
dent variables. *Computer Physics Communications*, 183(4):937-946, 2012.

S. Kucherenko, O. V. Klymenko, and N. Shah. Sobol' indices for problems defined in non-rectangular domains. *Reliability Engineering&System Safety*, 167:218-231, 2017.

S. N. Lahiri. *Resampling Methods for Dependent Data*. Springer Science+Business Media, New York, 2003.

J. Lei, M. G'Sell, A. Rinaldo, R. J. Tibshirani, and L. Wasserman. Distribution-free predictive inference for regression. *Journal of the American Statistical Association*, 113:1094–1111, 2018.

G. Li and H. Rabitz. Relationship between sensitivity indices defined by variance- and covariance-based methods. *Reliability Engineering&System Safety*, 167:136-157, 2017.

S. Lundberg and S.-I. Lee. A unified approach to interpreting model predictions. preprint, arXiv, 2017.
https://arxiv.org/abs/1705.07874.
G. Mainik. Risk aggregation with empirical margins: Latin hypercubes, empirical copulas, and convergence of sum distributions. *Journal of Multivariate Analysis*, 141:197-216, 2015.

T. A. Mara and S. Tarantola. Variance-based sensitivity indices for models with dependent inputs. *Reliability Engineering & System Safety*, 107:115-121, 2012.

T. A. Mara, S. Tarantola, and P. Annoni. Non-parametric methods for global sensitivity analysis of model output with dependent inputs. *Environmental Modelling & Software*, 72:173-183, 2015.

J. Matoušek. On the $L^2$-discrepancy for anchored boxes. *Journal of Complexity*, 14(4):527-556, 1998.

W. J. Murdoch, C. Singh, K. Kumbier, R. Abbasi-Asl, and B. Yu. Definitions, methods and applications in interpretable machine learning. *Proceedings of the National Academy of Sciences*, 116(44):22071-22080, 2019.

M. H. Neumann. A central limit theorem for triangular arrays of weakly dependent random variables, with applications in statistics. *ESAIM: Probability and Statistics*, 17:120-134, 2013.

J. E. Oakley and A. O'Hagan. Probabilistic sensitivity analysis of complex models: A Bayesian approach. *Journal of the Royal Statistical Society, Series B*, 66(3):751-769, 2004.

A. B. Owen. Scrambled net variance for integrals of smooth functions. *The Annals of Statistics*, 25(4): 1541-1562, 1997.

A. B. Owen and C. R. Hoyt. Efficient estimation of the ANOVA mean dimension, with an application to neural net classification. *SIAM/ASA Journal on Uncertainty Quantification*, 9(2), 2021.

A. B. Owen and C. Prieur. On Shapley value for measuring importance of dependent inputs. *SIAM/ASA Journal on Uncertainty Quantification*, 5(1):986-1002, 2017.

K. Pearson. Notes on regression and inheritance in the case of two parents. *Proc. Royal Soc. London*, 58: 240-242, 1895.

E. Plischke. How to compute variance-based sensitivity indicators with your spreadsheet software. *Environmental Modelling & Software*, 35:188-191, 2012.

E. Plischke, G. Rabitti, and E. Borgonovo. Has the spell been broken? Estimating global sensitivity measures via nearest neighbors. 2022. In Preparation.

C. Prieur and S. Tarantola. Variance-based sensitivity analysis: Theory and estimation algorithms. In *Handbook of Uncertainty Quantification*, pages 1217-1239. Springer Verlag, Cham, 2017.

H. Rabitz and Ö. F. Alış. General foundations of high-dimensional model representations. *J. Math. Chem.*, 25(2-3):197-233, 1999.

M. Rosenblatt. Remarks on a multivariate transformation. *Ann. Math. Statist.*, 23(3):470-472, 1952.

A. Saltelli and S. Tarantola. On the relative importance of input factors in mathematical models: Safety assessment for nuclear waste disposal. *Journal of the American Statistical Association*, 97(459):702-709, 2002.

A. Saltelli, K. Chan, and E. M. Scott. *Sensitivity Analysis*. John Wiley&Sons, Chichester, 2000.

I. M. Sobol'. Sensitivity estimates for nonlinear mathematical models. *Mathematical Modelling & Computational Experiments*, 1:407-414, 1993.

I. M. Sobol', S. Tarantola, D. Gatelli, S. S. Kucherenko, and W. Mauntz. Estimating the approximation error when fixing unessential factors in global sensitivity analysis. *Reliability Engineering & System Safety*, 92:957-960, 2007.

D. M. Sparkman, J. E. Garza, H. R. Millwater, Jr., and B. P. Smarslok. Importance sampling-based post-processing method for global sensitivity analysis. In *18th AIAA Non-Deterministic Approaches Conference. 4-8 January 2016, San Diego, California, USA*. AIAA SciTech, 2016. Paper #AIAA 2016-1444.

X. Sun, W. Zhong, and P. Ma. An asymptotic and empirical smoothing parameters selection method for smoothing spline anova models in large samples. *Biometrika*, 108(1):149-166, 2021.

A. W. van der Vaart. *Asymptotic Statistics*, volume 3 of *Cambridge Series in Statistical and Probabilistic Mathematics*. Cambridge University Press, 1998.
G. Mainik. Risk aggregation with empirical margins: Latin hypercubes, empirical copulas, and convergence of sum distributions. *Journal of Multivariate Analysis*, 141:197-216, 2015.

T. A. Mara and S. Tarantola. Variance-based sensitivity indices for models with dependent inputs. *Reliability Engineering&System Safety*, 107:115-121, 2012.

T. A. Mara, S. Tarantola, and P. Annoni. Non-parametric methods for global sensitivity analysis of model output with dependent inputs. *Environmental Modelling&Software*, 72:173-183, 2015.

J. Matoušek. On the $L^2$-discrepancy for anchored boxes. *Journal of Complexity*, 14(4):527-556, 1998.

W. J. Murdoch, C. Singh, K. Kumbier, R. Abbasi-Asl, and B. Yu. Definitions, methods and applications in interpretable machine learning. *Proceedings of the National Academy of Sciences*, 116(44):22071-22080, 2019.

M. H. Neumann. A central limit theorem for triangular arrays of weakly dependent random variables, with applications in statistics. *ESAIM: Probability and Statistics*, 17:120-134, 2013.

J. E. Oakley and A. O'Hagan. Probabilistic sensitivity analysis of complex models: A Bayesian approach. *Journal of the Royal Statistical Society, Series B*, 66(3):751-769, 2004.

A. B. Owen. Scrambled net variance for integrals of smooth functions. *The Annals of Statistics*, 25(4): 1541-1562, 1997.

A. B. Owen and C. R. Hoyt. Efficient estimation of the ANOVA mean dimension, with an application to neural net classification. *SIAM/ASA Journal on Uncertainty Quantification*, 9(2), 2021.

A. B. Owen and C. Prieur. On Shapley value for measuring importance of dependent inputs. *SIAM/ASA Journal on Uncertainty Quantification*, 5(1):986-1002, 2017.

K. Pearson. Notes on regression and inheritance in the case of two parents. *Proc. Royal Soc. London*, 58: 240-242, 1895.

E. Plischke. How to compute variance-based sensitivity indicators with your spreadsheet software. *Environmental Modelling&Software*, 35:188-191, 2012.

E. Plischke, G. Rabitti, and E. Borgonovo. Has the spell been broken? Estimating global sensitivity measures via nearest neighbors. 2022. In Preparation.

C. Prieur and S. Tarantola. Variance-based sensitivity analysis: Theory and estimation algorithms. In *Handbook of Uncertainty Quantification*, pages 1217-1239. Springer Verlag, Cham, 2017.

H. Rabitz and Ö. F. Alış. General foundations of high-dimensional model representations. *J. Math. Chem.*, 25(2-3):197-233, 1999.

M. Rosenblatt. Remarks on a multivariate transformation. *Ann. Math. Statist.*, 23(3):470-472, 1952.

A. Saltelli and S. Tarantola. On the relative importance of input factors in mathematical models: Safety assessment for nuclear waste disposal. *Journal of the American Statistical Association*, 97(459):702-709, 2002.

A. Saltelli, K. Chan, and E. M. Scott. *Sensitivity Analysis*. John Wiley&Sons, Chichester, 2000.

I. M. Sobol'. Sensitivity estimates for nonlinear mathematical models. *Mathematical Modelling & Computational Experiments*, 1:407-414, 1993.

I. M. Sobol', S. Tarantola, D. Gatelli, S. S. Kucherenko, and W. Mauntz. Estimating the approximation error when fixing unessential factors in global sensitivity analysis. *Reliability Engineering&System Safety*, 92:957-960, 2007.

D. M. Sparkman, J. E. Garza, H. R. Millwater, Jr., and B. P. Smarslok. Importance sampling-based post-processing method for global sensitivity analysis. In *18th AIAA Non-Deterministic Approaches Conference. 4-8 January 2016, San Diego, California, USA*. AIAA SciTech, 2016. Paper #AIAA 2016-1444.

X. Sun, W. Zhong, and P. Ma. An asymptotic and empirical smoothing parameters selection method for smoothing spline anova models in large samples. *Biometrika*, 108(1):149-166, 2021.

A. W. van der Vaart. *Asymptotic Statistics*, volume 3 of *Cambridge Series in Statistical and Probabilistic Mathematics*. Cambridge University Press, 1998.
B. D. Williamson, P. B. Gilbert, M. Carone, and N. Simon. Nonparametric variable importance assessment using machine learning techniques. *Biometrics*, 77:9–22, 2021. With discussion.

B. D. Williamson, P. B. Gilbert, N. R. Simon, and M. Carone. A general framework for inference on algorithm-agnostic variable importance. *Journal of the American Statistical Association*, 118(543):1645–1658, 2023.

# A Proofs

## A.1 More on conditional independence

In this appendix we detail more the role of conditional independence for the computation of total Sobol' indices. In particular, we provide the proof of Lemma 1, Proposition 2 stated in Section 2.

*Proof of Lemma 1.* From Fréchet [1934] we know that:

$$V[E[Y|X_{-u}]] = E[E[Y|X_{-u}]^2] - E[Y]^2 = E[E[Y \cdot E[Y|X_{-u}] - E[Y]^2|X_{-u}]] = \text{cov}(Y, E[Y|X_{-u}]).$$

Then, we deduce the following covariance representation of $\tau_u$:

$$\tau_u = E[V[Y|X_{-u}]] = V[Y] - \text{cov}(Y, E[Y|X_{-u}]). \qquad (45)$$

Now, let $Y'$ be an independent replicate of $Y$ conditionally on $X_{-u}$. As a consequence of conditional independence, $E[Y \cdot Y'|X_{-u}] = E[Y|X_{-u}]E[Y'|X_{-u}]$. Hence

$$\text{cov}(Y, Y') = E[E[Y \cdot Y' - E[Y]^2|X_{-u}]] = E[E[Y|X_{-u}] \cdot E[Y'|X_{-u}]] - E[Y]^2 = V[E[Y|X_{-u}]]$$

so that

$$\tau_u = V[Y] - \text{cov}(Y, Y') = \frac{1}{2}(V[Y] + V[Y']) - \text{cov}(Y, Y') = \frac{1}{2}E[(Y - Y')^2]$$

where $V[Y] = V[Y']$ as both are identically distributed. □

Let us shed some more light on conditional independence. For this, let $Z$ be a general $m$-dimensional random vector that takes its values on a Cartesian product space. Again for notation simplicity, we assume that the joint probability distribution of $Z$ has a density function $h(z)$ with respect to Lebesgue measure. For $u, v$ two disjoint subsets of $[m] = \{1, 2, \dots, m\}$, we denote the disjoint union of $u$ and $v$ by $u + v$.

**Definition 15.** Let $u, v, w$ be pairwise disjoint index sets in $[m]$. Then $Z_u$ and $Z_v$ are conditionally independent given $Z_w$ ($u \perp v|w$, for short) if for all $z_u \in \mathcal{Z}_u, z_v \in \mathcal{Z}_v, z_w \in \mathcal{Z}_w$, the density $h$ satisfies

$$h_{u+v|w}(z_{u+v}|z_w) = h_{u|w}(z_u|z_w) \cdot h_{v|w}(z_v|z_w). \qquad (46)$$

Multiplying (46) by $h_w(z_w)$ we find

$$h_{u+v+w}(z_{u+v+w}) = h_{u|w}(z_u|z_w) \cdot h_{v|w}(z_v|z_w)h_w(z_w) = \begin{cases} h_{v|w}(z_v|z_w) \cdot h_{u+w}(z_{u+w}), \\ h_{u|w}(z_u|z_w) \cdot h_{v+w}(z_{v+w}). \end{cases} \qquad (47)$$

*Proof of Proposition 2.* The proof of Proposition 2 follows straightforwardly from these observations. □

## A.2 Proof of the results stated in Section 3

*Proof of Proposition 5.* Consider an independent copy $X'$ of $X$. Projections onto index subsets $u$ and $-u$ keep independence intact, i.e., $X'_u = P_u(X')$ and $X_{-u} = P_{-u}(X)$ are independent. The random vector obtained by glueing these two vectors together therefore has a density $f_u \cdot f_{-u}$, breaking the inter-block dependence. Hence generally for a measurable function $h: \mathbb{R}^d \to \mathbb{R}$, we have

$$E[h(X'_u : X_{-u})] = \iint h(x'_u : x_{-u}) f_{-u}(x_{-u}) f_u(x'_u) dx_{-u} dx'_u.$$

Now, in order to compare the expectation of $h(X)$ with $X$ possibly being dependent, we split the argument $X$ into two arguments via projections onto subdimensions indexed by $u$ and $-u$, so that we may write the joint density as product of marginal and conditional density,

$$E[h(X_u : X_{-u})] = \iint h(x_u : x_{-u}) f_{u|-u}(x_u|x_{-u}) f_{-u}(x_{-u}) dx_u dx_{-u}.$$
B. D. Williamson, P. B. Gilbert, M. Carone, and N. Simon. Nonparametric variable importance assessment using machine learning techniques. Biometrics, 77:9-22, 2021. With discussion.

B. D. Williamson, P. B. Gilbert, N. R. Simon, and M. Carone. A general framework for inference on algorithm-agnostic variable importance. Journal of the American Statistical Association, 118(543):1645-1658, 2023.

# A Proofs

## A.1 More on conditional independence

In this appendix we detail more the role of conditional independence for the computation of total Sobol' indices. In particular, we provide the proof of Lemma 1, Proposition 2 stated in Section 2.

*Proof of Lemma 1.* From Fréchet [1934] we know that:

$$V[E[Y|X_{-u}]] = E[E[Y|X_{-u}]^2] - E[Y]^2 = E[E[Y \cdot E[Y|X_{-u}] - E[Y]^2|X_{-u}]] = \text{cov}(Y, E[Y|X_{-u}]).$$

Then, we deduce the following covariance representation of $\tau_u$:

$$\tau_u = E[V[Y|X_{-u}]] = V[Y] - \text{cov}(Y, E[Y|X_{-u}]). \qquad (45)$$

Now, let $Y'$ be an independent replicate of $Y$ conditionally on $X_{-u}$. As a consequence of conditional independence, $E[Y \cdot Y'|X_{-u}] = E[Y|X_{-u}]E[Y'|X_{-u}]$. Hence

$$\text{cov}(Y, Y') = E[E[Y \cdot Y' - E[Y]^2|X_{-u}]] = E[E[Y|X_{-u}] \cdot E[Y'|X_{-u}]] - E[Y]^2 = V[E[Y|X_{-u}]]$$

so that

$$\tau_u = V[Y] - \text{cov}(Y, Y') = \frac{1}{2} (V[Y] + V[Y']) - \text{cov}(Y, Y') = \frac{1}{2} E[(Y - Y')^2]$$

where $V[Y] = V[Y']$ as both are identically distributed. □

Let us shed some more light on conditional independence. For this, let $Z$ be a general $m$-dimensional random vector that takes its values on a Cartesian product space. Again for notation simplicity, we assume that the joint probability distribution of $Z$ has a density function $h(z)$ with respect to Lebesgue measure. For $u, v$ two disjoint subsets of $[m] = \{1, 2, \dots, m\}$, we denote the disjoint union of $u$ and $v$ by $u + v$.

**Definition 15.** Let $u, v, w$ be pairwise disjoint index sets in $[m]$. Then $Z_u$ and $Z_v$ are conditionally independent given $Z_w$ ($u \perp v|w$, for short) if for all $z_u \in \mathcal{Z}_u$, $z_v \in \mathcal{Z}_v$, $z_w \in \mathcal{Z}_w$, the density $h$ satisfies

$$h_{u+v|w}(z_{u+v}|z_w) = h_{u|w}(z_u|z_w) \cdot h_{v|w}(z_v|z_w). \qquad (46)$$

Multiplying (46) by $h_w(z_w)$ we find

$$h_{u+v+w}(z_{u+v+w}) = h_{u|w}(z_u|z_w) \cdot h_{v|w}(z_v|z_w)h_w(z_w) = \begin{cases} h_{v|w}(z_v|z_w) \cdot h_{u+w}(z_{u+w}), \\ h_{u|w}(z_u|z_w) \cdot h_{v+w}(z_{v+w}). \end{cases} \qquad (47)$$

*Proof of Proposition 2.* The proof of Proposition 2 follows straightforwardly from these observations. □

## A.2 Proof of the results stated in Section 3

*Proof of Proposition 5.* Consider an independent copy $X'$ of $X$. Projections onto index subsets $u$ and $-u$ keep independence intact, i.e., $X'_u = P_u(X')$ and $X_{-u} = P_{-u}(X)$ are independent. The random vector obtained by glueing these two vectors together therefore has a density $f_u \cdot f_{-u}$, breaking the inter-block dependence. Hence generally for a measurable function $h: \mathbb{R}^d \to \mathbb{R}$, we have

$$E[h(X'_u : X_{-u})] = \iint h(x'_u : x_{-u}) f_{-u}(x_{-u}) f_u(x'_u) dx_{-u} dx'_u.$$

Now, in order to compare the expectation of $h(X)$ with $X$ possibly being dependent, we split the argument $X$ into two arguments via projections onto subdimensions indexed by $u$ and $-u$, so that we may write the joint density as product of marginal and conditional density,

$$E[h(X_u : X_{-u})] = \iint h(x_u : x_{-u}) f_{u|-u}(x_u|x_{-u}) f_{-u}(x_{-u}) dx_u dx_{-u}.$$
Considering all three terms in a function $h_2: \mathbb{R}^{|u|} \times \mathbb{R}^{|d-|u|} \times \mathbb{R}^{|u|} \to \mathbb{R}$ and taking its expectation then leads to

$$ \mathbb{E}[h_2(X_u, X_{-u}, X'_u)] = \iiint h_2(x_u, x_{-u}, x'_u) f_{u|-u}(x_u|x_{-u}) f_{-u}(x_{-u}) f_u(x'_u) dx_u dx_{-u} dx'_u. $$

Let us now consider the weighted Jansen's estimator,

$$ h_2(X_u, X_{-u}, X'_u) = \frac{1}{2} \frac{f_{u,-u}(X'_u : X_{-u})}{f_u(X'_u) f_{-u}(X_{-u})} (g(X_u : X_{-u}) - g(X'_u : X_{-u}))^2. $$

Then by (9),

$$ \begin{aligned} \mathbb{E}[h_2(X_u, X_{-u}, X'_u)] &= \frac{1}{2} \iiint \frac{f_{u,-u}(x'_u : x_{-u})}{f_u(x'_u) f_{-u}(x_{-u})} (g(x_u : x_{-u}) - g(x'_u : x_{-u}))^2 \cdot \\ & \qquad f_{u|-u}(x_u|x_{-u}) f_{-u}(x_{-u}) f_u(x'_u) d(x_u, x_{-u}, x'_u) \\ &= \frac{1}{2} \iiint f_{u|-u}(x'_u|x_{-u}) (g(x_u : x_{-u}) - g(x'_u : x_{-u}))^2 \cdot \\ & \qquad f_{u|-u}(x_u|x_{-u}) f_{-u}(x_{-u}) d(x_u, x_{-u}, x'_u) = \tau_u. \end{aligned} \quad (48) $$

The last equality follows from (9). By definition, a sample consists of realizing $n$ independent copies of $X$. Hence two different copies of $X$, $X^i$ and $X^j$, $i, j = 1, \dots, n$, $i \neq j$, are independent. Then, an estimator of $\tau_u$ is

$$ \hat{\tau}_{u,n}^{U} = \frac{1}{n(n-1)} \sum_{i=1}^{n} \sum_{j \neq i} h_2(X_u^i, X_{-u}^i, X_u^j), \quad (49) $$

which yields (18).

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

**Proof of Lemma 6.** Let us define the random vector $W = (W_1, W_2) = (X_u, X_{-u})$ (the random vector $X$ split according to the index set $u$). Let $W^i$, $i = 1, \dots, n$, be identical copies of $W$. We then write the estimator as follows:

$$ \hat{\tau}_{u,n}^{U} = \frac{1}{n(n-1)} \sum_{i=1}^{n} \sum_{j \neq i} \Phi(W^i, W^j) = \left(\frac{n}{2}\right)^{-1} \sum_{i=1}^{n} \sum_{j>i} \phi^s(W^i, W^j) \quad (50) $$

with

$$ \Phi^s(W^i, W^j) = \frac{1}{2} (\Phi(W^i, W^j) + \Phi(W^j, W^i)), \quad (51) $$

$$ \Phi(W^i, W^j) = \frac{1}{2} i_u (W_1^i, W_2^j) \left( g(W_1^j : W_2^j) - g(W_1^i : W_2^j) \right)^2. \quad (52) $$

$$ \text{Hence, } \hat{\tau}_{u,n}^{U} \text{ defines a U-statistic of order 2 for } \tau_u = \mathbb{E}[\Phi(W^i, W^j)] = \mathbb{E}[\hat{\tau}_{u,n}^{U}]. \quad \square $$

**Proof of Lemma 7.** The proof of (19) follows directly from (50) using the computations in [Bose and Chatterjee, 2018, Section 1.2.1]. The asymptotic normality is a consequence of [Bose and Chatterjee, 2018, Theorem 1.1].

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

## A.3 Proof of the results stated in Section 4

**Proof of Theorem 8.** First, the $i$th realization and the $s_n(i)$th one in the input sample are independent. We thus deduce from (48) that for any $1 \le i \le n$,

$$ \frac{1}{2} \mathbb{E} \left[ \iota_j (X^{s_n(i)}, X^i) \left( g(X_j^i : X_{-j}^i) - g(X_j^{s_n(i)} : X_{-j}^i) \right)^2 \right] = \tau_j \quad (53) $$

thus the estimator $\hat{\tau}_{j,n}$ is unbiased. To prove the central limit theorem, we first decompose $\hat{\tau}_{j,n}^S$ as:

$$ \hat{\tau}_{j,n}^{S} = \frac{n-1}{n} \hat{\tau}_{j,n-1} + \frac{1}{2n} \frac{f_{j,-j}(X_j^1 : X_{-j}^n)}{f_j(X_j^1) f_{-j}(X_{-j}^n)} (g(X_j^n : X_{-j}^n) - g(X_j^1 : X_{-j}^n))^2 $$

where $\hat{\tau}_{j,n-1}$ stands for the estimator built from Formula (21) on $(X_1, \dots, X_{n-1})$ and by replacing the cyclic shift-by-one $s_n$ by the acyclic one $s_{n-1}^a$. Then, $\hat{\tau}_{j,n-1} = \frac{1}{n-1} \sum_{i=1}^{n-1} V_j^i$ and the sequence $(V_j^i)_{i \ge n}$ is stationary and 1-dependent. Thus, the limit $\sqrt{n}(\hat{\tau}_{j,n-1} - \tau_j) \to N(0, \sigma_j^2)$ follows from [Diananda, 1954, Theorem 5]. Finally, noting that

$$ \mathbb{E} \left[ \left[ \sqrt{n} \frac{1}{2n} \frac{f_{j,-j}(X_j^1 : X_{-j}^n)}{f_j(X_j^1) f_{-j}(X_{-j}^n)} \left( g(X_j^n : X_{-j}^n) - g(X_j^1 : X_{-j}^n) \right)^2 \right] \right] = \frac{\tau_j}{\sqrt{n}} \xrightarrow{n \to +\infty} 0 \quad (54) $$

and applying Slutsky's Theorem we get (22).

<figure><img src="image_15.png" /></figure>
Considering all three terms in a function $h_2: \mathbb{R}^{|u|} \times \mathbb{R}^{d-|u|} \times \mathbb{R}^{|u|} \to \mathbb{R}$ and taking its expectation then leads to

$$ \mathbb{E}[h_2(X_u, X_{-u}, X'_u)] = \iiint h_2(x_u, x_{-u}, x'_u) f_{u|-u}(x_u|x_{-u}) f_{-u}(x_{-u}) f_u(x'_u) dx_u dx_{-u} dx'_u. $$

Let us now consider the weighted Jansen's estimator,

$$ h_2(X_u, X_{-u}, X'_u) = \frac{1}{2} \frac{f_{u,-u}(X'_u : X_{-u})}{f_u(X'_u) f_{-u}(X_{-u})} (g(X_u : X_{-u}) - g(X'_u : X_{-u}))^2. $$

Then by (9),

$$ \begin{aligned} \mathbb{E}[h_2(X_u, X_{-u}, X'_u)] &= \frac{1}{2} \iiint \frac{f_{u,-u}(x'_u : x_{-u})}{f_u(x'_u) f_{-u}(x_{-u})} (g(x_u : x_{-u}) - g(x'_u : x_{-u}))^2 \cdot \\ & \qquad f_{u|-u}(x_u|x_{-u}) f_{-u}(x_{-u}) f_u(x'_u) d(x_u, x_{-u}, x'_u) \\ &= \frac{1}{2} \iiint f_{u|-u}(x'_u|x_{-u}) (g(x_u : x_{-u}) - g(x'_u : x_{-u}))^2 \cdot \\ & \qquad f_{u|-u}(x_u|x_{-u}) f_{-u}(x_{-u}) d(x_u, x_{-u}, x'_u) = \tau_u. \end{aligned} \quad (48) $$

The last equality follows from (9). By definition, a sample consists of realizing $n$ independent copies of $X$. Hence two different copies of $X$, $X^i$ and $X^j$, $i, j = 1, \dots, n$, $i \neq j$, are independent. Then, an estimator of $\tau_u$ is

$$ \hat{\tau}_{u,n}^{U} = \frac{1}{n(n-1)} \sum_{i=1}^{n} \sum_{j \neq i} h_2(X_u^i, X_{-u}^i, X_u^j), \quad (49) $$

which yields (18). □

**Proof of Lemma 6.** Let us define the random vector $W = (W_1, W_2) = (X_u, X_{-u})$ (the random vector $X$ split according to the index set $u$). Let $W^i$, $i = 1, \dots, n$, be identical copies of $W$. We then write the estimator as follows:

$$ \hat{\tau}_{u,n}^{U} = \frac{1}{n(n-1)} \sum_{i=1}^{n} \sum_{j \neq i} \Phi(W^i, W^j) = \left(\frac{n}{2}\right)^{-1} \sum_{i=1}^{n} \sum_{j>i} \phi^s(W^i, W^j) \quad (50) $$

with

$$ \Phi^s(W^i, W^j) = \frac{1}{2} (\Phi(W^i, W^j) + \Phi(W^j, W^i)), \quad (51) $$

$$ \Phi(W^i, W^j) = \frac{1}{2} i_u (W_1^i, W_2^j) \left( g(W_1^j : W_2^j) - g(W_1^i : W_2^j) \right)^2. \quad (52) $$

$$ \text{Hence, } \hat{\tau}_{u,n}^{U} \text{ defines a U-statistic of order 2 for } \tau_u = \mathbb{E}[\Phi(W^i, W^j)] = \mathbb{E}[\hat{\tau}_{u,n}^{U}]. \quad \square $$

**Proof of Lemma 7.** The proof of (19) follows directly from (50) using the computations in [Bose and Chatterjee, 2018, Section 1.2.1]. The asymptotic normality is a consequence of [Bose and Chatterjee, 2018, Theorem 1.1]. □

## A.3 Proof of the results stated in Section 4

**Proof of Theorem 8.** First, the $i$th realization and the $s_n(i)$th one in the input sample are independent. We thus deduce from (48) that for any $1 \le i \le n$,

$$ \frac{1}{2} \mathbb{E} \left[ \iota_j (X^{s_n(i)}, X^i) \left( g(X_j^i : X_{-j}^i) - g(X_j^{s_n(i)} : X_{-j}^i) \right)^2 \right] = \tau_j \quad (53) $$

thus the estimator $\hat{\tau}_{j,n}$ is unbiased. To prove the central limit theorem, we first decompose $\hat{\tau}_{j,n}^S$ as:

$$ \hat{\tau}_{j,n}^S = \frac{n-1}{n} \hat{\tau}_{j,n-1} + \frac{1}{2n} \frac{f_{j,-j}(X_j^1 : X_{-j}^n)}{f_j(X_j^1) f_{-j}(X_{-j}^n)} (g(X_j^n : X_{-j}^n) - g(X_j^1 : X_{-j}^n))^2 $$

where $\hat{\tau}_{j,n-1}$ stands for the estimator built from Formula (21) on $(X_1, \dots, X_{n-1})$ and by replacing the cyclic shift-by-one $s_n$ by the acyclic one $s_{n-1}^a$. Then, $\hat{\tau}_{j,n-1} = \frac{1}{n-1} \sum_{i=1}^{n-1} V_j^i$ and the sequence $(V_j^i)_{i \ge n}$ is stationary and 1-dependent. Thus, the limit $\sqrt{n}(\hat{\tau}_{j,n-1} - \tau_j) \to N(0, \sigma_j^2)$ follows from [Diananda, 1954, Theorem 5]. Finally, noting that

$$ \mathbb{E} \left[ \left[ \sqrt{n} \frac{1}{2n} \frac{f_{j,-j}(X_j^1 : X_{-j}^n)}{f_j(X_j^1) f_{-j}(X_{-j}^n)} \left( g(X_j^n : X_{-j}^n) - g(X_j^1 : X_{-j}^n) \right)^2 \right] \right] = \frac{\tau_j}{\sqrt{n}} \xrightarrow{n \to +\infty} 0 \quad (54) $$

and applying Slutsky's Theorem we get (22). □
*Proof of Corollary 9.* The result follows from Theorem 8 and by applying the Delta method (see, e.g., van der Vaart [1998]). First, by mimicking the proof of Theorem 8, it is possible to prove that for any $\alpha, \beta, \gamma$, a central limit theorem holds true for $\alpha\hat{\tau}_{j,n}^S + \beta\frac{1}{n}\sum_{i=1}^n g(X^i) + \gamma\frac{1}{n}\sum_{i=1}^n g^2(X^i)$. It yields a central limit theorem for $(U_{1,n}, U_{2,n}, U_{3,n}) = (\hat{\tau}_{j,n}^S, \frac{1}{n}\sum_{i=1}^n g(X^i), \frac{1}{n}\sum_{i=1}^n g^2(X^i))$, namely

$$ \sqrt{n} \left( (U_{1,n}, U_{2,n}, U_{3,n})^T - \theta_j \right) \xrightarrow{n \to +\infty} \mathcal{N}(0, \Sigma_j) $$

with $\theta_j = (\tau_j, \mathbb{E}[Y], \mathbb{E}[Y^2])$ and

$$ \Sigma_j = \begin{pmatrix} \mathrm{V}[V_j^1] & \mathrm{cov}(V_j^1, Y^1) & \mathrm{cov}(V_j^1, (Y^1)^2) \\ * & \mathrm{V}(Y^1) & \mathrm{cov}(Y^1, (Y^1)^2) \\ * & * & \mathrm{V}[(Y^1)^2] \end{pmatrix}. $$

Then we can prove the central limit theorem for $T_j$, using the Delta method on $\psi(x, y, z) = \frac{x}{z-y^2}$ and $\theta_j = (\tau_j, \mathbb{E}[Y], \mathbb{E}[Y^2])$. More precisely, let $\rho_j$ denote the gradient of $\psi$ at $\theta_j$. We have $\rho_j = \nabla\psi(\theta_j) = \frac{1}{\mathrm{V}[Y]} [1, 2T_j \mathbb{E}[Y], -T_j]^T$. Thus

$$ \sqrt{n} (\hat{\tau}_{j,n}^S - T_j) \xrightarrow{n \to +\infty} \mathcal{N}(0, \rho_j^T \Sigma_j \rho_j). $$

It concludes the proof of Corollary 9. $\square$

*Proof of Theorem 10.* To prove that $\hat{\tau}_{j,n}^D$ is unbiased, we use the same arguments as the ones used to prove that $\hat{\tau}_{j,n}^S$ defined by (21) in Theorem 8 is unbiased, additionally noting that $\pi_n(i) \neq i$ for all $i = 1, \dots, n$. To prove the central limit theorem, we first decompose, for each $n$, the permutation $\pi_n$ in cycles $C_{1,n}, \dots, C_{m_n,n}$. Let us arbitrarily fix the first element in each cycle. We then form $p_n$ blocks, with $p_n = \max_{1 \le k \le m_n} \ell_{k,n}$ and $\ell_{k,n}$ the length of cycle $C_{k,n}$. For each $n$, we re-order the $X^i$'s so that the first $b_{1,n} = m_n$ re-ordered variables are the first element of each cycle, the next $b_{2,n}$ re-ordered variables are the second element of each cycle with length at least two and so on until the $n - \sum_{v=1}^{p_n-1} b_{v,n}$ last $b_{p_n,n}$ re-ordered variables which are the last element in each cycle of length $p_n$. Here, $1 \le b_{p_n,n} \le \dots \le b_{1,n} = m_n$. For each $n$, we denote the re-ordered sequence of $X^i$'s by $X^{i,n}$, $1 \le i \le n$, $n \ge 1$. We then define $S_{v,n} = \sum_{i=k_{v-1,n}+1}^{k_{v,n}} \tilde{V}_j^{i,n}$, with $k_0 = 0$, $k_{v,n} = \sum_{w=1}^v b_{w,n}$ and $\tilde{V}_j^{i,n}$ defined as $V_j^i - \tau_j/\sqrt{n}$ but with the $X^{i,n}$s in place of the $X^i$s. More precisely, we have to use the trick in the proof of Theorem 8 by replacing first the last variable $X^i$ in cycle $C_{1,n}$ by $X^{n+1}, \dots$, the last variable $X^i$ in cycle $C_{m_n,n}$ by $X^{n+m_n}$. As $\lim_{n \to +\infty} m_n/\sqrt{n} = 0$, we prove that the remaining term decreases fast enough not to perturb the result of the central limit theorem (see the proof of Theorem 8 for more details). Now, as $\lim_{n \to +\infty} m_n/\sqrt{n} = 0$, then $p_n \ge n/m_n \ge \sqrt{n}/m_n \to +\infty$. For each $n \ge 1$, the sequence $(S_{v,n})_{1 \le v \le p_n}$ is a sequence of 1-dependent variables. Then applying [Neumann, 2013, Theorem 2.1], we get the central limit theorem, as soon as there exists $\delta > 0$ such that $\mathbb{E}[|V_j^1|^{2+\delta}|] < +\infty$ for some $\delta > 0$. Indeed, as variables inside each block are centered, independent and identically distributed,

$$ \sum_{v=1}^{p_n} \mathbb{E}[S_{v,n}^2] = \sum_{v=1}^{p_n} \frac{b_{v,n}}{n} \mathrm{V}[V_j^1] = \mathrm{V}[V_j^1] < +\infty. $$

Then, Assumption (2.1) in [Neumann, 2013, Theorem 2.1] is true due to the stationarity of $V_j^i V_j^{i+1}$. Indeed, due to 1-dependence,

$$ \begin{aligned} \mathrm{V} \left[ \sum_{v=1}^{p_n} S_{v,n} \right] &= \sum_{v=1}^{p_n} \mathrm{V}[S_{v,n}] + 2 \sum_{1 \le v < w \le p_n} \mathrm{cov}(S_{v,n}, S_{w,n}) = \sum_{v=1}^{p_n} b_{v,n} \frac{\mathrm{V}[V_j^1]}{n} + 2 \sum_{v=1}^{p_n-1} \mathrm{cov}(S_{v,n}, S_{v+1,n}) \\ &= \frac{\mathrm{V}[V_j^1]}{n} \sum_{v=1}^{p_n} b_{v,n} + \frac{2}{n} \sum_{v=1}^{p_n-1} b_{v+1,n} \mathrm{cov}(V_j^1, V_j^2) = \mathrm{V}[V_j^1] + 2 \frac{n-m_n}{n} \mathrm{cov}(V_j^1, V_j^2) \xrightarrow{n \to +\infty} \sigma_j^2 < +\infty. \end{aligned} $$
*Proof of Corollary 9.* The result follows from Theorem 8 and by applying the Delta method (see, e.g., van der Vaart [1998]). First, by mimicking the proof of Theorem 8, it is possible to prove that for any $\alpha, \beta, \gamma$, a central limit theorem holds true for $\alpha\hat{\tau}_{j,n}^S + \beta\frac{1}{n}\sum_{i=1}^n g(X^i) + \gamma\frac{1}{n}\sum_{i=1}^n g^2(X^i)$. It yields a central limit theorem for $(U_{1,n}, U_{2,n}, U_{3,n}) = (\hat{\tau}_{j,n}^S, \frac{1}{n}\sum_{i=1}^n g(X^i), \frac{1}{n}\sum_{i=1}^n g^2(X^i))$, namely

$$ \sqrt{n} \left( (U_{1,n}, U_{2,n}, U_{3,n})^T - \theta_j \right) \xrightarrow{n \to +\infty} \mathcal{N}(0, \Sigma_j) $$

with $\theta_j = (\tau_j, \mathbb{E}[Y], \mathbb{E}[Y^2])$ and

$$ \Sigma_j = \begin{pmatrix} \mathrm{V}[V_j^1] & \mathrm{cov}(V_j^1, Y^1) & \mathrm{cov}(V_j^1, (Y^1)^2) \\ * & \mathrm{V}(Y^1) & \mathrm{cov}(Y^1, (Y^1)^2) \\ * & * & \mathrm{V}[(Y^1)^2] \end{pmatrix}. $$

Then we can prove the central limit theorem for $T_j$, using the Delta method on $\psi(x, y, z) = \frac{x}{z-y^2}$ and $\theta_j = (\tau_j, \mathbb{E}[Y], \mathbb{E}[Y^2])$. More precisely, let $\rho_j$ denote the gradient of $\psi$ at $\theta_j$. We have $\rho_j = \nabla\psi(\theta_j) = \frac{1}{\mathrm{V}[Y]} [1, 2T_j \mathbb{E}[Y], -T_j]^T$. Thus

$$ \sqrt{n} (\hat{\tau}_{j,n}^S - T_j) \xrightarrow{n \to +\infty} \mathcal{N}(0, \rho_j^T \Sigma_j \rho_j). $$

It concludes the proof of Corollary 9. $\square$

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

*Proof of Theorem 10.* To prove that $\hat{\tau}_{j,n}^D$ is unbiased, we use the same arguments as the ones used to prove that $\hat{\tau}_{j,n}^S$ defined by (21) in Theorem 8 is unbiased, additionally noting that $\pi_n(i) \neq i$ for all $i = 1, \dots, n$. To prove the central limit theorem, we first decompose, for each $n$, the permutation $\pi_n$ in cycles $C_{1,n}, \dots, C_{m_n,n}$. Let us arbitrarily fix the first element in each cycle. We then form $p_n$ blocks, with $p_n = \max_{1 \le k \le m_n} \ell_{k,n}$ and $\ell_{k,n}$ the length of cycle $C_{k,n}$. For each $n$, we re-order the $X^i$'s so that the first $b_{1,n} = m_n$ re-ordered variables are the first element of each cycle, the next $b_{2,n}$ re-ordered variables are the second element of each cycle with length at least two and so on until the $n - \sum_{v=1}^{p_n-1} b_{v,n}$ last $b_{p_n,n}$ re-ordered variables which are the last element in each cycle of length $p_n$. Here, $1 \le b_{p_n,n} \le \dots \le b_{1,n} = m_n$. For each $n$, we denote the re-ordered sequence of $X^i$'s by $X^{i,n}$, $1 \le i \le n$, $n \ge 1$. We then define $S_{v,n} = \sum_{i=k_{v-1,n}+1}^{k_{v,n}} \tilde{V}_j^{i,n}$, with $k_0 = 0$, $k_{v,n} = \sum_{w=1}^v b_{w,n}$ and $\tilde{V}_j^{i,n}$ defined as $V_j^i - \tau_j/\sqrt{n}$ but with the $X^{i,n}$s in place of the $X^i$s. More precisely, we have to use the trick in the proof of Theorem 8 by replacing first the last variable $X^i$ in cycle $C_{1,n}$ by $X^{n+1}, \dots$, the last variable $X^i$ in cycle $C_{m_n,n}$ by $X^{n+m_n}$. As $\lim_{n \to +\infty} m_n/\sqrt{n} = 0$, we prove that the remaining term decreases fast enough not to perturb the result of the central limit theorem (see the proof of Theorem 8 for more details). Now, as $\lim_{n \to +\infty} m_n/\sqrt{n} = 0$, then $p_n \ge n/m_n \ge \sqrt{n}/m_n \to +\infty$. For each $n \ge 1$, the sequence $(S_{v,n})_{1 \le v \le p_n}$ is a sequence of 1-dependent variables. Then applying [Neumann, 2013, Theorem 2.1], we get the central limit theorem, as soon as there exists $\delta > 0$ such that $\mathbb{E}[|V_j^1|^{2+\delta}|] < +\infty$ for some $\delta > 0$. Indeed, as variables inside each block are centered, independent and identically distributed,

$$ \sum_{v=1}^{p_n} \mathbb{E}[S_{v,n}^2] = \sum_{v=1}^{p_n} \frac{b_{v,n}}{n} \mathrm{V}[V_j^1] = \mathrm{V}[V_j^1] < +\infty. $$

Then, Assumption (2.1) in [Neumann, 2013, Theorem 2.1] is true due to the stationarity of $V_j^i V_j^{i+1}$. Indeed, due to 1-dependence,

$$ \begin{aligned} \mathrm{V} \left[ \sum_{v=1}^{p_n} S_{v,n} \right] &= \sum_{v=1}^{p_n} \mathrm{V}[S_{v,n}] + 2 \sum_{1 \le v < w \le p_n} \mathrm{cov}(S_{v,n}, S_{w,n}) = \sum_{v=1}^{p_n} b_{v,n} \frac{\mathrm{V}[V_j^1]}{n} + 2 \sum_{v=1}^{p_n-1} \mathrm{cov}(S_{v,n}, S_{v+1,n}) \\ &= \frac{\mathrm{V}[V_j^1]}{n} \sum_{v=1}^{p_n} b_{v,n} + \frac{2}{n} \sum_{v=1}^{p_n-1} b_{v+1,n} \mathrm{cov}(V_j^1, V_j^2) = \mathrm{V}[V_j^1] + 2 \frac{n-m_n}{n} \mathrm{cov}(V_j^1, V_j^2) \xrightarrow{n \to +\infty} \sigma_j^2 < +\infty. \end{aligned} $$
Let us now prove that Assumption (2.2) of [Neumann, 2013, Theorem 2.1] holds. For any $\varepsilon > 0$ we have

$$
\begin{align*}
\sum_{v=1}^{p_n} \mathbb{E} [S_{v,n}^2 \mathbb{I}_{|S_{v,n}| > \varepsilon}] &\le \sum_{v=1}^{p_n} \mathbb{E} [|S_{v,n}|^{2+\delta}]^{\frac{2}{2+\delta}} (\Pr(|S_{v,n}| > \varepsilon))^{\frac{\delta}{2+\delta}} \quad \text{(using Hölder Inequality)} \\
&\le \varepsilon^{-\delta} \sum_{v=1}^{p_n} \mathbb{E} [|S_{v,n}|^{2+\delta}] \quad \text{(using Markov Inequality)} \\
&\le \varepsilon^{-\delta} \sum_{v=1}^{p_n} \mathbb{E} \left[ \left( \sum_{i=k_{v-1,n}+1}^{k_{v,n}} n^{-1/2} |V_j^i - \tau_j| \right)^{2+\delta} \right] \\
&= \varepsilon^{-\delta} \sum_{v=1}^{p_n} \frac{b_{v,n}^{2+\delta}}{n^{1+\frac{\delta}{2}}} \mathbb{E} \left[ \left( \sum_{i=k_{v-1,n}+1}^{k_{v,n}} b_{v,n}^{-1} |V_j^i - \tau_j| \right)^{2+\delta} \right] \\
&\le \varepsilon^{-\delta} \sum_{v=1}^{p_n} \frac{b_{v,n}^{2+\delta}}{n^{1+\frac{\delta}{2}}} \mathbb{E} \left[ \sum_{i=k_{v-1,n}+1}^{k_{v,n}} b_{v,n}^{-1} |V_j^i - \tau_j|^{2+\delta} \right] \quad \text{(using Jensen Inequality)} \\
&\le \frac{1}{\varepsilon^{\delta}} \frac{1}{n^{1+\frac{\delta}{2}}} \sum_{v=1}^{p_n} b_{v,n}^{1+\delta} \sum_{i=k_{v-1,n}+1}^{k_{v,n}} \mathbb{E} [|V_j^i - \tau_j|^{2+\delta}] \\
&\le \frac{1}{\varepsilon^{\delta}} \frac{1}{n^{1+\frac{\delta}{2}}} m_n^{1+\delta} \sum_{v=1}^{p_n} b_{v,n} \mathbb{E} [|V_j^1 - \tau_j|^{2+\delta}] \quad \text{(by stationarity of } (|V_j^i - \tau_j|)_{i \ge 1} \text{)} \\
&\le \frac{1}{\varepsilon^{\delta}} \frac{m_n^{1+\delta}}{n^{\frac{\delta}{2}}} \mathbb{E} [|V_j^1 - \tau_j|^{2+\delta}] \underset{n \to +\infty}{\longrightarrow} 0
\end{align*}
$$

The limit holds as $\mathbb{E} [|V_j^1|^{2+\delta}] < +\infty$ and $\lim_{n \to +\infty} \frac{m_n^{1+\delta}}{\sqrt{n}^{\delta}} \to 0$. This concludes the proof of Theorem 10. $\square$

## A.4 Link with Breiman's Permutation Importance

*Proof of Proposition 13.* Let us write the squared loss function as $\mathcal{L}(Y, g(X; \theta)) = (Y - g(X; \theta))^2$. Then, (31) becomes

$$ \mathrm{MDA}_{j}^{\mathrm{FR}} = \mathbb{E}[\iota_{j}(X', X)(Y - g(X' : X_{-j}; \theta))^{2}] - \mathbb{E}[(Y - g(X; \theta))^{2}]. \quad (55) $$

Using the assumption that $g(X; \theta)$ is a perfect predictor, we have $Y = g(X; \theta)$ for all values of $X$, so that $\mathbb{E}[(Y - g(X; \theta))^2] = \mathbb{E}[0] = 0$. Then, (55) becomes

$$ \mathrm{MDA}_j = \mathbb{E}[\iota_j(X', X)(g(X, \theta) - g(X'_j : X_{-j}; \theta))^2] = \tau_j. $$

Hence total effects and MDA with a weighted squared loss are the same. $\square$

## A.5 Proof of Theorem 14

*Proof of Theorem 14.* Main effect: We consider the Gaussian multivariate distribution of $(X, Y)$. Setting $u = \{d+1\}, v = \emptyset, w = \{j\}$ in (32), then the variance of $Y$ conditionally on $X_j$ is

$$ \Sigma_{Y,Y} - \Sigma_{Y,j} \Sigma_{j,j}^{-1} \Sigma_{j,Y} = \beta^T \Sigma \beta - \beta^T (\Sigma_{[d],j} \Sigma_{j,j}^{-1} \Sigma_{j,[d]}) \beta, \quad (56) $$

using here the special structure of the augmented matrix. This variance is constant (as it does not depend on the value of $X_j$), so that $\mathbb{E}[V[Y|X_j]] = V[Y|X_j]$. However, for main effects we are interested in the variance of the conditional expectation, so we have to subtract the value in (56) from the total variance which yields (33).

**Total effect:** The Rosenblatt transform in the Gaussian case uses the Cholesky decomposition of the covariance matrix. The Cholesky matrix $C = \mathrm{chol}(\Sigma)$ is an upper triangular matrix such that $C^T C = \Sigma$. If $Z \sim \mathcal{N}(0, I)$ then $X = \mu + C^T Z \sim \mathcal{N}(\mu, \Sigma)$. One can define the Cholesky decomposition recursively,

$$ \mathrm{chol}(\Sigma) = \begin{pmatrix} \Sigma_{1,1}^{1/2} & \Sigma_{1,1}^{-1/2} \cdot \Sigma_{1,-1} \\ 0 & \mathrm{chol}(\Sigma_{-1,-1} - \Sigma_{1,-1}^T \Sigma_{1,1}^{-1} \Sigma_{1,-1}) \end{pmatrix} $$

where the index $-1$ denotes all coordinates but the first (if $\Sigma$ is a scalar then $\mathrm{chol}(\Sigma) = \Sigma^{1/2}$). By reordering the input factors, we may assume without loss of generality that $j = d$. Then define $Y = \beta^0 + \beta^T(\mu + C^T Z)$ and $Y' = \beta^0 + \beta^T(\mu + C^T Z')$ with $Z, Z' \sim \mathcal{N}(0, I)$, differing only in their last coordinate independent from each other. Then $Y$ and $Y'$ are identically (but not independently) distributed, and we know from Section 2.2 that $\mathrm{cov}(Y, Y') = \frac{1}{2} \mathbb{E}[(Y - Y')^2] = \frac{1}{2} \beta_d^2 C_{d,d}^2 \mathbb{E}[(Z_d - Z_d')^2]$. But $Z_d$ and $Z_d'$ are iid. standard normal, i.e., $\frac{1}{2} \mathbb{E}[(Z_d - Z_d')^2] = 1$. We are left with the identification of the last diagonal entry of the Cholesky matrix. Because of its hierarchical triangular structure, it keeps subdeterminants intact, and $C_{d,d} = \sqrt{\frac{\det \Sigma}{\det \Sigma_{-d,-d}}}$.
Let us now prove that Assumption (2.2) of [Neumann, 2013, Theorem 2.1] holds. For any $\varepsilon > 0$ we have

$$
\begin{align*}
\sum_{v=1}^{p_n} \mathbb{E} [S_{v,n}^2 \mathbb{I}_{|S_{v,n}| > \varepsilon}] &\le \sum_{v=1}^{p_n} \mathbb{E} [|S_{v,n}|^{2+\delta}]^{\frac{2}{2+\delta}} (\Pr(|S_{v,n}| > \varepsilon))^{\frac{\delta}{2+\delta}} \quad \text{(using Hölder Inequality)} \\
&\le \varepsilon^{-\delta} \sum_{v=1}^{p_n} \mathbb{E} [|S_{v,n}|^{2+\delta}] \quad \text{(using Markov Inequality)} \\
&\le \varepsilon^{-\delta} \sum_{v=1}^{p_n} \mathbb{E} \left[ \left( \sum_{i=k_{v-1,n}+1}^{k_{v,n}} n^{-1/2} |V_j^i - \tau_j| \right)^{2+\delta} \right] \\
&= \varepsilon^{-\delta} \sum_{v=1}^{p_n} \frac{b_{v,n}^{2+\delta}}{n^{1+\frac{\delta}{2}}} \mathbb{E} \left[ \left( \sum_{i=k_{v-1,n}+1}^{k_{v,n}} b_{v,n}^{-1} |V_j^i - \tau_j| \right)^{2+\delta} \right] \\
&\le \varepsilon^{-\delta} \sum_{v=1}^{p_n} \frac{b_{v,n}^{2+\delta}}{n^{1+\frac{\delta}{2}}} \mathbb{E} \left[ \sum_{i=k_{v-1,n}+1}^{k_{v,n}} b_{v,n}^{-1} |V_j^i - \tau_j|^{2+\delta} \right] \quad \text{(using Jensen Inequality)} \\
&\le \frac{1}{\varepsilon^{\delta}} \frac{1}{n^{1+\frac{\delta}{2}}} \sum_{v=1}^{p_n} b_{v,n}^{1+\delta} \sum_{i=k_{v-1,n}+1}^{k_{v,n}} \mathbb{E} [|V_j^i - \tau_j|^{2+\delta}] \\
&\le \frac{1}{\varepsilon^{\delta}} \frac{1}{n^{1+\frac{\delta}{2}}} m_n^{1+\delta} \sum_{v=1}^{p_n} b_{v,n} \mathbb{E} [|V_j^1 - \tau_j|^{2+\delta}] \quad \text{(by stationarity of } (|V_j^i - \tau_j|)_{i \ge 1} \text{)} \\
&\le \frac{1}{\varepsilon^{\delta}} \frac{m_n^{1+\delta}}{n^{\frac{\delta}{2}}} \mathbb{E} [|V_j^1 - \tau_j|^{2+\delta}] \underset{n \to +\infty}{\longrightarrow} 0
\end{align*}
$$

The limit holds as $\mathbb{E} [|V_j^1|^{2+\delta}] < +\infty$ and $\lim_{n \to +\infty} \frac{m_n^{1+\delta}}{\sqrt{n}^{\delta}} \to 0$. This concludes the proof of Theorem 10. $\square$

## A.4 Link with Breiman's Permutation Importance

*Proof of Proposition 13.* Let us write the squared loss function as $\mathcal{L}(Y, g(X; \theta)) = (Y - g(X; \theta))^2$. Then, (31) becomes

$$ \text{MDA}_{j}^{\text{FR}} = \mathbb{E}[\iota_{j}(X', X)(Y - g(X'_{j} : X_{-j}; \theta))^{2}] - \mathbb{E}[(Y - g(X; \theta))^{2}]. \quad (55) $$

Using the assumption that $g(X; \theta)$ is a perfect predictor, we have $Y = g(X; \theta)$ for all values of $X$, so that $\mathbb{E}[(Y - g(X; \theta))^2] = \mathbb{E}[0] = 0$. Then, (55) becomes

$$ \text{MDA}_j = \mathbb{E}[\iota_j(X', X)(g(X, \theta) - g(X'_j : X_{-j}; \theta))^2] = \tau_j. $$

Hence total effects and MDA with a weighted squared loss are the same. $\square$

## A.5 Proof of Theorem 14

*Proof of Theorem 14.* Main effect: We consider the Gaussian multivariate distribution of $(X, Y)$. Setting $u = \{d+1\}, v = \emptyset, w = \{j\}$ in (32), then the variance of $Y$ conditionally on $X_j$ is

$$ \Sigma_{Y,Y} - \Sigma_{Y,j} \Sigma_{j,j}^{-1} \Sigma_{j,Y} = \beta^T \Sigma \beta - \beta^T (\Sigma_{[d],j} \Sigma_{j,j}^{-1} \Sigma_{j,[d]}) \beta, \quad (56) $$

using here the special structure of the augmented matrix. This variance is constant (as it does not depend on the value of $X_j$), so that $\mathbb{E}[V[Y|X_j]] = V[Y|X_j]$. However, for main effects we are interested in the variance of the conditional expectation, so we have to subtract the value in (56) from the total variance which yields (33).

**Total effect:** The Rosenblatt transform in the Gaussian case uses the Cholesky decomposition of the covariance matrix. The Cholesky matrix $C = \text{chol}(\Sigma)$ is an upper triangular matrix such that $C^T C = \Sigma$. If $Z \sim \mathcal{N}(0, I)$ then $X = \mu + C^T Z \sim \mathcal{N}(\mu, \Sigma)$. One can define the Cholesky decomposition recursively,

$$ \text{chol}(\Sigma) = \begin{pmatrix} \Sigma_{1,1}^{1/2} & \Sigma_{1,1}^{-1/2} \cdot \Sigma_{1,-1} \\ 0 & \text{chol}(\Sigma_{-1,-1} - \Sigma_{1,-1}^T \Sigma_{1,1}^{-1} \Sigma_{1,-1}) \end{pmatrix} $$

where the index $-1$ denotes all coordinates but the first (if $\Sigma$ is a scalar then $\text{chol}(\Sigma) = \Sigma^{1/2}$). By reordering the input factors, we may assume without loss of generality that $j = d$. Then define $Y = \beta^0 + \beta^T(\mu + C^T Z)$ and $Y' = \beta^0 + \beta^T(\mu + C^T Z')$ with $Z, Z' \sim \mathcal{N}(0, I)$, differing only in their last coordinate independent from each other. Then $Y$ and $Y'$ are identically (but not independently) distributed, and we know from Section 2.2 that $\text{cov}(Y, Y') = \frac{1}{2} \mathbb{E}[(Y - Y')^2] = \frac{1}{2} \beta_d^2 C_{d,d}^2 \mathbb{E}[(Z_d - Z_d')^2]$. But $Z_d$ and $Z_d'$ are iid. standard normal, i.e., $\frac{1}{2} \mathbb{E}[(Z_d - Z_d')^2] = 1$. We are left with the identification of the last diagonal entry of the Cholesky matrix. Because of its hierarchical triangular structure, it keeps subdeterminants intact, and $C_{d,d} = \sqrt{\frac{\det \Sigma}{\det \Sigma_{-d,-d}}}$.
B Software Availability

The repository

https://gitlab.gwdg.de/elmar.plischke/global-sensitivity-analysis-collection/
contains implementations of the generalized winding stairs total estimator (windsi.m, windsi.R), of the mix-
and-reweight U-statistics estimator (totalsdep.m), and of the derange-and-reweight and shift-and-reweight
pick-and-freeze estimators (totalsdep_pnf.m, totalsdep_pnfshift.m) for use with MATLAB/OCTAVE or R.
B Software Availability

The repository

https://gitlab.gwdg.de/elmar.plischke/global-sensitivity-analysis-collection/
contains implementations of the generalized winding stairs total estimator (windsi.m, windsi.R), of the mix-
and-reweight U-statistics estimator (totalsdep.m), and of the derange-and-reweight and shift-and-reweight
pick-and-freeze estimators (totalsdep_pnf.m, totalsdep_pnfshift.m) for use with MATLAB/OCTAVE or R.