Distance and Kernel-Based Measures for Global and Local Two-Sample Conditional Distribution Testing

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

77,912 characters · 9 sections · 95 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

\affil[1]{Department of Statistics and Data Science, Cornell University} \affil[2]{Department of ISOM, Hong Kong University of Science and Technology} \affil[3]{Department of Statistics, Texas A&M University}

abstractTesting the equality of two conditional distributions is crucial in various modern applications, including transfer learning and causal inference. Despite its importance, this fundamental problem has received surprisingly little attention in the literature, with existing works focusing exclusively on global two-sample conditional distribution testing. Based on distance and kernel methods, this paper presents the first unified framework for both global and local two-sample conditional distribution testing. To this end, we introduce distance and kernel-based measures that characterize the homogeneity of two conditional distributions. Drawing from the concept of conditional U-statistics, we propose consistent estimators for these measures. Theoretically, we derive the convergence rates and the asymptotic distributions of the estimators under both the null and alternative hypotheses. Utilizing these measures, along with a local bootstrap approach, we develop global and local tests that can detect discrepancies between two conditional distributions at global and local levels, respectively. Our tests demonstrate reliable performance through simulations and real data analysis.

{\smallKeywords: Conditional distribution; Energy distance; Maximum mean discrepancy; Two-sample testing; U-statistics.}

Introduction

The canonical setting for (nonparametric) two-sample testing focuses on assessing the equality of two unconditional distributions. In many contemporary applications, however, people are instead interested in testing the equality of two conditional distributions. Consider two independent data sets $\{(Y_{i}^{(l)},X_{i}^{(l)})\}_{i=1}^{n_{l}}$ for $l=1,2$. Here $Y_{i}^{(l)}\in\mathcal{Y}$ and $X_{i}^{(l)}\in\mathbb{R}^{p}$, where $\mathcal{Y}$ is allowed to be a general metric space (see Remark (ref)). For $l=1,2$, assume that $\{(Y_{i}^{(l)},X_{i}^{(l)})\}_{i=1}^{n_{l}}$ are independent and identically distributed (i.i.d.) samples of $(Y^{(l)},X^{(l)})\sim P^{(l)}\equiv P_{Y\mid X}^{(l)}\otimes P_{X}^{(l)}$, where $P_{X}^{(l)}$ is the marginal distribution of $X^{(l)}$ and $P_{Y\mid X=x}^{(l)}$ is the conditional distribution of $Y^{(l)}$ given $X^{(l)}=x$. To ensure that the testing problem considered below is nontrivial, we assume that $P_{X}^{(1)}$ and $P_{X}^{(2)}$ are equivalent, i.e., $P_{X}^{(1)}\ll P_{X}^{(2)}$ and $P_{X}^{(2)}\ll P_{X}^{(1)}$, where $P \ll Q$ means that $P$ is absolutely continuous in reference to $Q$. We aim to test the following hypothesis, which we call the global two-sample conditional distribution testing problem,

equation[equation omitted — 176 chars of source]

Due to the equivalence between $P_{X}^{(1)}$ and $P_{X}^{(2)}$, the hypotheses in ((ref)) can be equivalently formulated by replacing $P^{(1)}_{X}$ with $P^{(2)}_{X}$.

We want to emphasize that $Y$ and $X$ denote two generic random variables and do not necessarily correspond to response and covariates, respectively. For instance, in the prior shift example below, $Y$ are covariates, and $X$ is a response in ((ref)). Moreover, the marginal distributions of $X$ from the two populations may differ (i.e., $P_{X}^{(1)}\neq P_{X}^{(2)}$), as seen in the motivating examples below. Thus, $H_{0}$ in ((ref)) is not equivalent to $P^{(1)}=P^{(2)}$, and unconditional two-sample tests for the equality of two joint distributions are generally not applicable in our context. Applying such tests would result in a failure to control the type I error when $P_{X}^{(1)}\neq P_{X}^{(2)}$. Furthermore, we mainly focus on scenarios in which both $Y$ and $X$ are continuous. When $Y$ is categorical, ((ref)) reduces to testing the equality of conditional means. Conversely, when $X$ is categorical, ((ref)) becomes equivalent to unconditional two-sample testing problems.

Hypothesis ((ref)) is central to many important problems in econometrics, machine learning, and statistics. For example, in transfer learning, the prior and covariate shift assumptions are commonly employed to tackle distributional differences between source and target populations kouw2018introduction. The prior shift assumption asserts that the conditional distribution of the covariates given the response is identical in both populations while allowing for a shift in the marginal distributions of the response. Conversely, the covariate shift assumption posits that the conditional distribution of the response given the covariates remains invariant across source and target populations, but the marginal distributions of the covariates can differ. Both assumptions are widely adopted in the literature; see, e.g., huang2024efficient,lee2024doubly for the prior shift assumption, and shimodaira2000improving,tibshirani2019conformal,liu2023augmented,ma2023optimally for the covariate shift assumption. Despite their prevalence, there is a paucity of work that formally validates these assumptions. Both of them can be framed as testing the equality of two conditional distributions as in ((ref)). Such tests are essential to the validity of methods developed under the prior or covariate shift assumptions.

Another motivating example comes from causal inference. Testing hypotheses in the context of treatment effect analysis has always been of interest imbens2009recent. Consider the standard setup based on the potential outcome framework rubin1974estimating. Suppose $\{(Y_{i},T_{i},X_{i})\}_{i=1}^{n}$ are i.i.d.\ observations of $(Y,T,X)$, where $Y$ is the observed outcome of interest, $T$ denotes a binary treatment (1: treated, 0: untreated), and $X$ are pretreatment covariates. For each subject, we define a pair of potential outcomes, $\{Y(1),Y(0)\}$, that would be observed if the subject had been given treatment, $Y(1)$, and control, $Y(0)$. One may be interested in testing the null hypothesis that the conditional distribution of $Y(1)\mid X$ is the same as that of $Y(0)\mid X$ imbens2009recent, i.e., zero conditional distributional treatment effect. Under the prevalent assumptions of consistency, $Y=TY(1)+(1-T)Y(0)$, and no unmeasured confounding, $\{Y(1),Y(0)\}\perp\!\!\!\perp T\mid X$, the conditional distributions of $Y(1)\mid X$ and $Y(0)\mid X$ can be identified as those of $Y\mid X,T=1$ and $Y\mid X,T=0$, respectively. Therefore, it can be formulated as testing hypothesis ((ref)) with $\{(Y_{i},X_{i}):T_{i}=1\}$ and $\{(Y_{i},X_{i}):T_{i}=0\}$ being the two sets of independent samples.

Global two-sample conditional distribution testing concerns the entire support of the conditioning variable. However, in many empirical applications, the primary interest lies in a subset of the population defined by a specific value of the conditioning variable bugni2025testing. For example, in regression discontinuity designs, researchers seek to compare outcome distributions conditional on the cutoff of the running variable shen2016distributional. Additional examples can be found in bugni2025testing. Formally, for a fixed value $x$ in the support of $P_{X}^{(1)}$ (or equivalently $P_{X}^{(2)}$), we are interested in testing

equation[equation omitted — 155 chars of source]

which we refer to as the local two-sample conditional distribution testing problem, in contrast to the global problem ((ref)). This problem is novel and is partially motivated by duong2013local and kim2019global in the context of unconditional two-sample testing. Notably, with a little extra effort, our proposed framework can readily handle the local testing problem ((ref)).

Recently, distance and kernel-based measures have received considerable attention in both the statistics and machine learning communities szekely2017energy,muandet2017kernel. These measures have been applied to a wide range of hypothesis testing problems, including two-sample testing szekely2004testing,baringhaus2004new,gretton2012kernel, goodness-of-fit testing szekely2005new,balasubramanian2021optimality, independence testing szekely2007measuring,gretton2007kernel,chakraborty2019distance,ke2020expected,deb2020measuring and conditional independence testing fukumizu2007kernel,wang2015conditional,sheng2023distance. This work aims to establish a distance and kernel-based framework to tackle both the global and local two-sample conditional distribution testing problems. To achieve this, we introduce the conditional energy distance ((ref)) and the conditional maximum mean discrepancy ((ref)), along with their integrated version ((ref)), which fully characterize the homogeneity of two conditional distributions. Additionally, we show the equivalence between the conditional energy distance and the conditional maximum mean discrepancy. Building on estimators of these measures, we develop global and local tests capable of detecting discrepancies between two conditional distributions at global and local levels, respectively.

Our estimation strategy employs a combination of U-statistics and kernel smoothing, initially introduced in the so-called conditional U-statistics by stute1991conditional and later applied to different problems by wang2015conditional and ke2020expected. To highlight our theoretical contributions, we summarize the distinct asymptotic behaviors of our global and local test statistics under both the null and alternative hypotheses in Table (ref). It is noteworthy that local tests based on distance and kernel measures have not been previously explored in the literature. Furthermore, while wang2015conditional and ke2020expected only provided the asymptotic distributions of their global test statistics under the null, we offer additional insights by systematically studying the properties of our statistics under both the null and alternative hypotheses. As a side note, we identify certain gaps in the derivations of the asymptotic null distribution in both wang2015conditional and ke2020expected, with details provided in Section (ref).

Our theoretical analysis reveals two key effects of applying kernel smoothing to U-statistic estimators of distance and kernel-based measures. First, kernel smoothing introduces bias into these estimators, necessitating the use of undersmoothing to control the resulting bias. Second, in contrast to the unconditional situation szekely2004testing,gretton2012kernel, kernel smoothing renders the U-statistics nondegenerate under the null hypotheses. However, when undersmoothing is used, the first-order projections in the Hoeffding decomposition of the U-statistics become asymptotically negligible under the null. Consequently, the asymptotic null distributions are determined by the second-order projections. Further details are provided in (ref)-(ref).

{

table[table omitted — 655 chars of source]

}

We now discuss related work and highlight several notable features of our framework. Although problems ((ref)) and ((ref)) are fundamental in various modern applications, surprisingly, there are very few methods available to test the equality of two conditional distributions. In the nonparametric testing literature, existing methods mainly focus on testing the equality of conditional moments of $Y$ given $X$. Specifically, most studies aim at testing the equality of conditional means, also known as the comparison of regression curves, as seen in hall1990bootstrap,kulasekera1995comparison,dette1998nonparametric,lavergne2001equality,neumeyer2003nonparametric, among others (see Section 7 in gonzalez2013updated for a detailed review). Similarly, in the causal inference literature, hypothesis testing has largely been limited to conditional average treatment effect (see, e.g., crump2008nonparametric). The methods we develop in this paper can detect general discrepancies between two conditional distributions, beyond specific moments. lee2009non and chang2015nonparametric proposed nonparametric tests for the null hypothesis of zero conditional distributional treatment effect. Their tests are built upon a Mann-Whitney statistic and cumulative distribution functions, respectively, and thus are only applicable for univariate $Y$. In our framework, $Y$ is allowed to take values in a general metric space, while $X$ can be multivariate. Very recently, utilizing techniques from conformal prediction, hu2024two proposed a test for the global testing problem ((ref)). However, their test involves (possibly unbalanced) sample splitting, and its performance relies on density ratio estimators of high quality (see Assumption 2(b) and related discussion therein). lee2024general introduced two frameworks for the hypothesis ((ref)). Their first framework converts a conditional independence test into a two-sample conditional distribution test, which effectively addresses an issue pointed out by an earlier version of this paper. Similar to hu2024two, their second framework is based on density ratio estimation and sample splitting, which are not required in our work. Notably, the local testing problem ((ref)) is new and has not been addressed in the literature. Our framework can accommodate both global and local two-sample conditional distribution testing. chen2022paired and chatterjee2024kernel considered the paired-sample problem, which is very different from our two-sample setting.

To end the introduction, we highlight the inherent challenges of two-sample conditional distribution testing in high-dimensional settings. When $Y$ is high-dimensional, even the simpler problem of unconditional two-sample testing lacks universally powerful methods. For example, distance and kernel-based approaches are known to suffer from the linearization effect as the dimension tends to infinity yan2023kernel. Conversely, when $X$ is high-dimensional, even testing the equality of conditional means becomes difficult without imposing strong structural assumptions, such as a sparse linear model xia2018two. Consequently, we do not consider high-dimensional settings in this work and instead leave them as an important direction for future research.

Preliminaries

This section provides an overview of the energy distance and maximum mean discrepancy, as well as their equivalence, which has been extensively discussed in sejdinovic2013equivalence.

For a nonempty set $\mathcal{U}$, a nonnegative function $\rho:\mathcal{U}\times\mathcal{U}\rightarrow[0,\infty)$ is called a semimetric on $\mathcal{U}$ if for any $u,u'\in\mathcal{U}$, it satisfies (i) $\rho(u,u')=0\iff u=u'$ and (ii) $\rho(u,u')=\rho(u',u)$. Then $(\mathcal{U},\rho)$ is said to be a semimetric space. The semimetric space $(\mathcal{U},\rho)$ is said to have negative type if for all $n\ge 2$, $u_{1},\ldots,u_{n}\in\mathcal{U}$ and $\alpha_{1},\ldots,\alpha_{n}\in\mathbb{R}$ with $\sum_{i=1}^{n}\alpha_{i}=0$, we have $\sum_{i,j=1}^{n}\alpha_{i}\alpha_{j}\rho(u_{i},u_{j})\le 0$. For $\theta>0$, define $\mathcal{M}_{\rho}^{\theta}(\mathcal{U})=\{\mu\in\mathcal{M}(\mathcal{U}):\int\rho^{\theta}(u,u_{0})d\mu(u)<\infty$, for some $u_{0}\in\mathcal{U}\}$, where $\mathcal{M}(\mathcal{U})$ denotes the set of all probability measures on $\mathcal{U}$. Suppose $P,Q\in\mathcal{M}_{\rho}^{1}(\mathcal{U})$, we have $\int\rho d(P-Q)^{2}\le 0$ when $(\mathcal{U},\rho)$ has negative type. We say that $(\mathcal{U},\rho)$ has strong negative type if it has negative type and the equality holds only when $P=Q$. The (generalized) energy distance between $P,Q\in\mathcal{M}_{\rho}^{1}(\mathcal{U})$ is defined as sejdinovic2013equivalence \[ \mathcal{D}_{\rho}(P,Q)=2\mathbb{E}\rho(U,V)-\mathbb{E}\rho(U,U')-\mathbb{E}\rho(V,V')=-\int\rho d(P-Q)^{2}, \] where $U,U'\overset{i.i.d.}{\sim}P$ and $V,V'\overset{i.i.d.}{\sim}Q$. If $(\mathcal{U},\rho)$ is of strong negative type, we have $\mathcal{D}_{\rho}(P,Q)\ge 0$ and the equality holds if and only if $P=Q$. Every separable Hilbert space is of strong negative type lyons2013distance. In particular, Euclidean spaces are separable.

Let $\mathcal{H}$ be a Hilbert space of real-valued functions defined on $\mathcal{U}$. A function $k:\mathcal{U}\times\mathcal{U}\rightarrow\mathbb{R}$ is a reproducing kernel of $\mathcal{H}$ if (i) $\forall u\in\mathcal{U}$, $k(\cdot,u)\in\mathcal{H}$ and (ii) $\forall u\in\mathcal{U},\forall f\in\mathcal{H}$, $\langle f,k(\cdot,u)\rangle_{\mathcal{H}}=f(u)$, where $\langle\cdot,\cdot\rangle_{\mathcal{H}}$ is the inner product associated with $\mathcal{H}$. If $\mathcal{H}$ has a reproducing kernel, it is said to be a reproducing kernel Hilbert space (RKHS). According to Moore-Aronszajn theorem berlinet2011reproducing, for every symmetric, positive definite function (henceforth kernel) $k$, there exists an associated RKHS $\mathcal{H}_{k}$ with reproducing kernel $k$. Define $\mathcal{M}_{k}^{1/2}(\mathcal{U})=\{\mu\in\mathcal{M}(\mathcal{U}):\int \sqrt{k(u,u)}d\mu(u)<\infty\}$. The kernel mean embedding of $P\in\mathcal{M}_{k}^{1/2}(\mathcal{U})$ into RKHS $\mathcal{H}_{k}$ is defined by the Bochner integral $\Pi_{k}(P)=\int k(\cdot,u)dP(u)\in\mathcal{H}_{k}$. The kernel $k$ is said to be characteristic if the mapping $\Pi_{k}$ is injective. Conditions under which kernels are characteristic have been studied by sriperumbudur2008injective,sriperumbudur2010hilbert, and examples include the Gaussian and the Laplace kernel. The maximum mean discrepancy (MMD) between $P,Q\in\mathcal{M}_{k}^{1/2}(\mathcal{U})$ is given by gretton2012kernel \[ \gamma_{k}(P,Q)=\|\Pi_{k}(P)-\Pi_{k}(Q)\|_{\mathcal{H}_k}. \] When the kernel $k$ is characteristic, we have $\gamma_{k}(P,Q)=0$ if and only if $P=Q$. Also, the following alternative representation of the squared MMD is useful \[ \gamma_{k}^{2}(P,Q)=\mathbb{E} k(U,U')+\mathbb{E} k(V,V')-2\mathbb{E} k(U,V)=\int kd(P-Q)^{2}, \] where $U,U'\overset{i.i.d.}{\sim}P$ and $V,V'\overset{i.i.d.}{\sim}Q$.

Let $\rho$ be a semimetric of negative type on $\mathcal{U}$. For any $u_{0}\in\mathcal{U}$, the function $k(u,u')=\rho(u,u_{0})+\rho(u',u_{0})-\rho(u,u')$ is positive definite, and is said to be the distance-induced kernel induced by $\rho$ and centered at $u_{0}$. Correspondingly, for a kernel $k$ on $\mathcal{U}$, the function $\rho(u,u')=\frac{1}{2}\{k(u,u)+k(u',u')\}-k(u,u')$ defines a valid semimetric $\rho$ of negative type on $\mathcal{U}$, and we will say that $k$ generates $\rho$. It is clear that every distance-induced kernel $k$ induced by $\rho$, also generates $\rho$. Theorem 22 in sejdinovic2013equivalence establishes the equivalence between the energy distance and MMD. Specifically, suppose $P,Q\in\mathcal{M}_{\rho}^{1}(\mathcal{U})$ and let $k$ be any kernel that generates $\rho$, then $\mathcal{D}_{\rho}(P,Q)=\gamma_{k}^{2}(P,Q)$.

Notation. For $l=1,2$, let $f_{l}(y,x)$, $f_{l}(x)$ and $f_{l}(y\mid x)$ be the joint probability density function of $P^{(l)}$, the marginal probability density function of $P_{X}^{(l)}$ and the conditional probability density function of $P_{Y\mid X=x}^{(l)}$, respectively (assuming the existence of these densities). For two sequences of real numbers $\{a_{n}\}_{n=1}^{\infty}$ and $\{b_{n}\}_{n=1}^{\infty}$, we write $a_{n}\asymp b_{n}$ if and only if $a_{n}=O(b_{n})$ and $b_{n}=O(a_{n})$. The symbols $\overset{p}{\rightarrow}$ and $\overset{d}{\rightarrow}$ stand for convergence in probability and in distribution, respectively.

Conditional energy distance and conditional maximum mean discrepancy

Let $\rho$ be a semimetric on $\mathcal{Y}$. We define the conditional energy distance at $x\in\mathbb{R}^{p}$ as the energy distance between $P_{Y\mid X=x}^{(1)}$ and $P_{Y\mid X=x}^{(2)}$.

definitionAssume $f_{l}(x)>0$ and $P_{Y\mid X=x}^{(l)}\in\mathcal{M}_{\rho}^{1}(\mathcal{Y})$ for $l=1,2$. The conditional energy distance (CED) between $P^{(1)}$ and $P^{(2)}$ at $x\in\mathbb{R}^{p}$ is defined as \begin{align} \begin{split} \mathcal{D}_{\rho}(x)&\coloneqq\mathcal{D}_{\rho}(P_{Y\mid X=x}^{(1)},P_{Y\mid X=x}^{(2)})=-\int_{\mathcal{Y}}\int_{\mathcal{Y}}\rho(y,y')d(P_{Y\mid X=x}^{(1)}-P_{Y\mid X=x}^{(2)})^{2}(y,y')\\ &=2\mathbb{E}\{\rho(Y^{(1)},Y^{(2)})\mid X^{(1)}=x,X^{(2)}=x\}-\mathbb{E}\{\rho(Y^{(1)},Y^{(1)\prime})\mid X^{(1)}=x,X^{(1)\prime}=x\}\\ &\quad-\mathbb{E}\{\rho(Y^{(2)},Y^{(2)\prime})\mid X^{(2)}=x,X^{(2)\prime}=x\}, \end{split} \end{align} where $(Y^{(1)\prime},X^{(1)\prime})$ and $(Y^{(2)\prime},X^{(2)\prime})$ are i.i.d.\ copies of $(Y^{(1)},X^{(1)})$ and $(Y^{(2)},X^{(2)})$, respectively.

When $\mathcal{Y}=\mathbb{R}^{q}$ and $\rho$ corresponds to the Euclidean distance, ((ref)) also has a nice interpretation in terms of conditional characteristic function. Specifically, for $l=1,2$ and $t\in\mathbb{R}^{q}$, the conditional characteristic function of $Y^{(l)}$ given $X^{(l)}=x$ is defined as $\varphi_{Y\mid X=x}^{(l)}(t)=\mathbb{E}\{\exp(\imath t^{\top}Y^{(l)})\mid X^{(l)}=x\}$, where $\imath=\sqrt{-1}$ is the imaginary unit. When $\rho(y,y')=\|y-y'\|$ with $\|\cdot\|$ being the Euclidean norm, we have \[ \mathcal{D}_{\rho}(x)=\frac{1}{c_{q}}\int_{\mathbb{R}^{q}}\frac{|\varphi_{Y\mid X=x}^{(1)}(t)-\varphi_{Y\mid X=x}^{(2)}(t)|^{2}}{\|t\|^{q+1}}dt, \] where $c_{q}=\pi^{(q+1)/2}/\Gamma((q+1)/2)$ is a constant with $\Gamma(\cdot)$ being the gamma function. The proof of this fact follows a similar approach to that of Proposition 1 in szekely2017energy, and is therefore omitted for brevity.

As the counterpart to the distance-based metric ((ref)), we now introduce a kernel-based metric to measure the discrepancy between the two conditional distributions. Let $\mathcal{H}_{k}$ be an RKHS associated with a kernel $k$ on $\mathcal{Y}$.

definitionAssume $f_{l}(x)>0$ and $P_{Y\mid X=x}^{(l)}\in\mathcal{M}_{k}^{1/2}(\mathcal{Y})$ for $l=1,2$. The conditional maximum mean discrepancy (CMMD) between $P^{(1)}$ and $P^{(2)}$ at $x\in\mathbb{R}^{p}$ is defined as the square root of \begin{align} \begin{split} \gamma_{k}^{2}(x)&\coloneqq\gamma_{k}^{2}(P_{Y\mid X=x}^{(1)},P_{Y\mid X=x}^{(2)})=\|\Pi_{k}(P_{Y\mid X=x}^{(1)})-\Pi_{k}(P_{Y\mid X=x}^{(2)})\|_{\mathcal{H}_{k}}^{2}\\ &=\int_{\mathcal{Y}}\int_{\mathcal{Y}}k(y,y')d(P_{Y\mid X=x}^{(1)}-P_{Y\mid X=x}^{(2)})^{2}(y,y')\\ &=\mathbb{E}\{k(Y^{(1)},Y^{(1)\prime})\mid X^{(1)}=x,X^{(1)\prime}=x\}+\mathbb{E}\{k(Y^{(2)},Y^{(2)\prime})\mid X^{(2)}=x,X^{(2)\prime}=x\}\\ &\quad-2\mathbb{E}\{k(Y^{(1)},Y^{(2)})\mid X^{(1)}=x,X^{(2)}=x\}, \end{split} \end{align} where $(Y^{(1)\prime},X^{(1)\prime})$ and $(Y^{(2)\prime},X^{(2)\prime})$ are i.i.d.\ copies of $(Y^{(1)},X^{(1)})$ and $(Y^{(2)},X^{(2)})$, respectively.

While conducting this research, we came across park2020measure, in which the authors present the CMMD in a slightly different form. However, their work does not address two-sample conditional distribution testing, and their estimation strategy differs entirely from ours. Specifically, the computational complexity of their estimator is $O((n_{1}+n_{2})^{3})$, whereas ours has a lower complexity of $O((n_{1}+n_{2})^{2})$ (see Remark (ref)). Moreover, park2020measure does not provide the convergence rate or the asymptotic distribution of their estimator. In contrast, we derive both the exact convergence rate and the asymptotic distribution of our statistic under both the null and alternative hypotheses. Additionally, we develop a unified framework that encompasses both distance and kernel-based measures. It is also important to note that the CMMD is a metric indexed by $x\in\mathbb{R}^{p}$, and the single measure ((ref)) introduced in Section (ref), which integrates the CMMD with a weight function, is not considered by park2020measure.

The CED and CMMD both serve to characterize the homogeneity of two conditional distributions. Besides, we can demonstrate the equivalence between the CED and CMMD, which holds for general semimetric of negative type. Similar results have been established for distance and kernel-based measures used in two-sample/independence testing sejdinovic2013equivalence and conditional independence testing sheng2023distance. Consequently, we will focus on the CED for the remainder of this paper.

theorem\begin{enumerate} • When the semimetric $\rho$ is of strong negative type, for any $x\in\mathbb{R}^{p}$ such that $f_{l}(x)>0$ and $P_{Y\mid X=x}^{(l)}\in\mathcal{M}_{\rho}^{1}(\mathcal{Y})\ (l=1,2)$, we have $\mathcal{D}_{\rho}(x)\ge 0$, and $\mathcal{D}_{\rho}(x)=0$ if and only if $H_{0}$ in ((ref)) holds. • When the kernel $k$ is characteristic, for any $x\in\mathbb{R}^{p}$ such that $f_{l}(x)>0$ and $P_{Y\mid X=x}^{(l)}\in\mathcal{M}_{k}^{1/2}(\mathcal{Y})\ (l=1,2)$, we have $\gamma_{k}(x)=0$ if and only if $H_{0}$ in ((ref)) holds. • Let $\rho$ be a semimetric of negative type on $\mathcal{Y}$. Suppose $f_{l}(x)>0$ and $P_{Y\mid X=x}^{(l)}\in\mathcal{M}_{\rho}^{1}(\mathcal{Y})$ for $l=1,2$. If $k$ is a kernel that generates $\rho$, then \[ \mathcal{D}_{\rho}(x)=\gamma_{k}^{2}(x). \] \end{enumerate}

We now present an estimator for the CED. Let $x(s)$ denote the $s$th component of $x\in\mathbb{R}^{p}$ for $s=1,\ldots,p$. For $l=1,2$, define $G_{h_{l}}:\mathbb{R}^{p}\rightarrow\mathbb{R}$ as \[ G_{h_{l}}(x)=\frac{1}{h_{l}^{p}}\prod_{s=1}^{p}g\bigg(\frac{x(s)}{h_{l}}\bigg),\quad x\in\mathbb{R}^{p}, \] where $g:\mathbb{R}\rightarrow\mathbb{R}$ is a univariate kernel function satisfying Assumption (ref) below, and $h_{1},h_{2}\in\mathbb{R}$ are the bandwidth parameters. For ease of presentation, here we use the same bandwidth $h_l$ for each component of $X^{(l)}$. In practice, one should always allow for varying bandwidths across each component of $X^{(l)}$. Motivated by the representation of CED in terms of the conditional moments given in ((ref)), we propose the following estimator for the CED:

align*[align* omitted — 892 chars of source]

which we call the sample conditional energy distance. Our estimation strategy, which combines U-statistics and kernel smoothing, has been employed in previous studies such as wang2015conditional,ke2020expected. It first appeared in the so-called conditional U-statistics stute1991conditional, which generalize the Nadaraya-Watson estimate of a regression function in the same way as Hoeffding's U-statistic is a generalization of the sample mean. Instead of V-statistics, we adopt U-statistics in this paper to avoid unnecessary bias; see the proof of Theorem (ref) for the explicit U-statistic form.

remarkThe definitions of the CED and CMMD, as well as the estimation procedures introduced, are applicable to a general metric space $\mathcal{Y}$. This implies that our framework is capable of handling cases where $Y$ is a non-Euclidean object, such as a curve, image, or even topological data encountered in various applications.

In what follows, we present the consistency of the sample CED in Theorem (ref), while Theorems (ref) and (ref) provide its asymptotic distributions under $H_{a}$ and $H_{0}$ in ((ref)), respectively. We shall make the following assumptions throughout the analysis.

assumptionThe univariate smoothing kernel $g$ is of order $\nu$ in the sense that $\int g(u)du=1$, $\int u^{r}g(u)du=0$ for $r=1,\cdots,\nu-1$, and $\int u^{\nu}g(u)du\neq 0$. Also, $g$ is bounded with $\int g^{2}(u)du<\infty$ and $\int |u|^{\nu}|g(u)|du<\infty$.

assumptionFor $l=1,2$, $h_{l}\rightarrow 0$ and $n_{l}h_{l}^{p}\rightarrow\infty$, as $n_{l}\rightarrow\infty$.

assumptionThe marginal densities $f_{1}(x)$ and $f_{2}(x)$ are $\nu$-times continuously differentiable.

assumption$\mathbb{E}\{\rho(Y^{(1)},Y^{(2)})\mid X^{(1)}=x_{1},X^{(2)}=x_{2}\}$, $\mathbb{E}\{\rho(Y^{(1)},Y^{(1)\prime})\mid X^{(1)}=x_{1},X^{(1)\prime}=x_{2}\}$, and $\mathbb{E}\{\rho(Y^{(2)},Y^{(2)\prime})\mid X^{(2)}=x_{1},X^{(2)\prime}=x_{2}\}$ are $\nu$-times continuously differentiable with respect to $(x_{1}^{\top},x_{2}^{\top})^{\top}$.

Assumptions (ref)-(ref) are standard in the literature (see, e.g., Section 1.11 in li2007nonparametric), allowing for multivariate $X$ and higher-order kernels. Assumption (ref) can be seen as the counterpart of the usual smoothness condition imposed on the regression function in classical kernel regression. A similar assumption is used in stute1991conditional.

theoremSuppose $f_{l}(x)>0$ and $P_{Y\mid X=x}^{(l)}\in\mathcal{M}_{\rho}^{2}(\mathcal{Y})$ for $l=1,2$. Under Assumptions (ref)-(ref), we have \[ \widehat{\mathcal{D}}_{\rho}(x)\overset{p}{\longrightarrow}\mathcal{D}_{\rho}(x),\quad\text{as }n_{1},n_{2}\rightarrow\infty. \]

theoremSuppose $f_{l}(x)>0$ and $P_{Y\mid X=x}^{(l)}\in\mathcal{M}_{\rho}^{2+\delta}(\mathcal{Y})$ for $l=1,2$ and some $\delta>0$. Under $H_{a}$ in ((ref)), Assumptions (ref)-(ref) and $n_{l}^{1/2}h_{l}^{p/2+\nu}\rightarrow 0$ as $n_{l}\rightarrow\infty$ for $l=1,2$, we have \[ \frac{\widehat{\mathcal{D}}_{\rho}(x)-\mathcal{D}_{\rho}(x)}{\frac{2}{\{f_{1}(x)f_{2}(x)\}^{2}}\sqrt{\frac{\xi_{1}^{2}}{n_{1}}+\frac{\xi_{2}^{2}}{n_{2}}}}\overset{d}{\longrightarrow}\mathcal{N}(0,1),\quad\text{as }n_{1},n_{2}\rightarrow\infty, \] where $\xi_{1}^{2}$ and $\xi_{2}^{2}$ are defined in equations ((ref))-((ref)) of Appendix (ref).

As shown in the proof, we have $\xi_{1}^{2}\asymp h_{1}^{-p}$ and $\xi_{2}^{2}\asymp h_{2}^{-p}$. Thus, under $H_{a}$ in ((ref)), $\widehat{\mathcal{D}}_{\rho}(x)-\mathcal{D}_{\rho}(x)=O_{p}((n_{1}h_{1}^{p})^{-1/2}+(n_{2}h_{2}^{p})^{-1/2})$, which is the same convergence rate as that of classical kernel regression ullah1999nonparametric and conditional U-statistics stute1991conditional. The condition $n_{l}^{1/2}h_{l}^{p/2+\nu}\rightarrow 0$ as $n_{l}\rightarrow\infty$ in Theorem (ref) requires undersmoothing in order to make the estimator asymptotically unbiased. Undersmoothing, coupled with the use of a higher-order kernel, is commonly employed in nonparametric and semiparametric estimation to reduce the bias of estimators involving kernel smoothing li2007nonparametric. The same condition with $p=1$ and $\nu=2$ is used in Section 3.4 of ullah1999nonparametric for the central limit theorem (CLT) of classical kernel regression and in stute1991conditional for the CLT of conditional U-statistics.

theoremDefine the double centered semimetric \[ \widetilde{\rho}(y,y')=\rho(y,y')-\mathbb{E}\{\rho(y,Y')\mid X'=x\}-\mathbb{E}\{\rho(Y,y')\mid X=x\}+\mathbb{E}\{\rho(Y,Y')\mid X=x,X'=x\}, \] where $Y\mid X=x,Y'\mid X'=x\overset{i.i.d.}{\sim}P_{Y\mid X=x}\coloneqq P_{Y\mid X=x}^{(1)}=P_{Y\mid X=x}^{(2)}$ under $H_{0}$ in ((ref)). Suppose $f_{l}(x)>0$ and $P_{Y\mid X=x}\in\mathcal{M}_{\rho}^{2+\delta}(\mathcal{Y})$ for $l=1,2$ and some $\delta>0$. Under $H_{0}$ in ((ref)), Assumptions (ref)-(ref), $n_{1}h_{1}^{p}/(n_{1}h_{1}^{p}+n_{2}h_{2}^{p})\rightarrow\tau\in(0,1)$, and $n_{l}h_{l}^{p+\nu}\rightarrow 0$ as $n_{l}\rightarrow\infty$ for $l=1,2$, we have \[ (n_{1}h_{1}^{p}+n_{2}h_{2}^{p})\widehat{\mathcal{D}}_{\rho}(x)\overset{d}{\longrightarrow}-\sum_{r=1}^{\infty}\lambda_{r}\bigg\{\bigg(\frac{1}{\sqrt{\tau}}\zeta_{r}-\frac{1}{\sqrt{1-\tau}}\eta_{r}\bigg)^{2}-\frac{(1-\tau)\sigma_{1}^{2}+\tau\sigma_{2}^{2}}{\tau(1-\tau)}\bigg\}, \] where $\{\zeta_r\}_{r=1}^{\infty}$ and $\{\eta_r\}_{r=1}^{\infty}$ are two sequences of independent normal random variables with \begin{align*} \zeta_{r}&\overset{i.i.d.}{\sim}\mathcal{N}(0,\sigma_{1}^{2}),\quad\sigma_{1}^{2}=\frac{1}{f_{1}(x)}\bigg\{\int g^{2}(u)du\bigg\}^{p},\\ \eta_{r}&\overset{i.i.d.}{\sim}\mathcal{N}(0,\sigma_{2}^{2}),\quad\sigma_{2}^{2}=\frac{1}{f_{2}(x)}\bigg\{\int g^{2}(u)du\bigg\}^{p}, \end{align*} and $\{\lambda_{r}\}_{r=1}^{\infty}$ are the eigenvalues of $\widetilde{\rho}$ with respect to the conditional distribution $P_{Y\mid X=x}$, i.e., \[ \mathbb{E}\{\widetilde{\rho}(y,Y')\phi_{r}(Y')\mid X'=x\}=\lambda_{r}\phi_{r}(y). \]

As demonstrated in Theorem (ref), the sample CED converges at a faster rate under $H_{0}$ in ((ref)): $\widehat{\mathcal{D}}_{\rho}(x)=O_{p}((n_{1}h_{1}^{p})^{-1}+(n_{2}h_{2}^{p})^{-1})$, compared to that under $H_{a}$ in ((ref)). This renders the undersmoothing condition in Theorem (ref) more stringent than that in Theorem (ref). Compared with the unconditional scenario gretton2012kernel, our results differ in several aspects. First, the convergence rates are distinct. Second, our spectral decomposition is performed with respect to the conditional distribution. Third, in contrast to i.i.d.\ standard normal random variables in gretton2012kernel, $\{\zeta_r\}_{r=1}^{\infty}$ and $\{\eta_r\}_{r=1}^{\infty}$ have different variances when $f_{1}(x)\neq f_{2}(x)$.

Integrated conditional energy distance

The CED is indexed by $x\in\mathbb{R}^{p}$, and thus is not a single number. We can obtain a single measure by integrating the CED with some weight function $\omega:\mathbb{R}^{p}\rightarrow\mathbb{R}$, i.e., \[ \mathcal{I}_{\rho}^{\omega}\coloneqq\int_{\mathbb{R}^{p}}\mathcal{D}_{\rho}(x)\omega(x)dx, \] which we call the integrated conditional energy distance (ICED). It follows directly from Theorem (ref) that $\mathcal{I}_{\rho}^{\omega}$ fully characterizes the homogeneity of two conditional distributions.

corollaryWhen $\omega$ is a nonnegative function with the same support as $P^{(1)}_{X}$ (or equivalently $P^{(2)}_{X}$) and the semimetric $\rho$ is of strong negative type, we have $\mathcal{I}_{\rho}^{\omega}=0$ if and only if $H_{0}$ in ((ref)) holds.

Motivated by su2007consistent,wang2015conditional,ke2020expected, we consider the weight function $\omega(x)=\{f_{1}(x)f_{2}(x)\}^{2}$ and define

equation[equation omitted — 128 chars of source]

This particular choice of weight function circumvents the random denominator issue. Otherwise, to deal with the density near zero and mitigate significant bias, additional stringent assumptions or trimming schemes may need to be used. For example, yin2020new assumes that the density functions $f_{1}(x)$ and $f_{2}(x)$ are bounded away from zero, which can be quite restrictive.

We employ the same estimation strategy as in Section (ref) and propose the following estimator for the ICED, which is a U-statistic by symmetrization (see the proof of Theorem (ref)):

align*[align* omitted — 816 chars of source]

We call $\widehat{\mathcal{I}}_{\rho}$ the sample integrated conditional energy distance.

remarkThe computational complexity of both the sample CED $\widehat{\mathcal{D}}_{\rho}(x)$ and sample ICED $\widehat{\mathcal{I}}_{\rho}$ is $O((n_{1}+n_{2})^{2})$, which is the same as that of the unconditional situation szekely2004testing,gretton2012kernel. The details are provided in Appendix (ref).

Theorem (ref) asserts that $\widehat{\mathcal{I}}_{\rho}$ is a consistent estimator of $\mathcal{I}_{\rho}$. By examining the Hoeffding decomposition of $\widehat{\mathcal{I}}_{\rho}$, Theorem (ref) describes its asymptotic distribution under the alternative hypothesis $H_{a}$ in ((ref)), while Theorem (ref) characterizes its asymptotic null distribution under $H_{0}$ in ((ref)). To get these results, we need to strengthen Assumptions (ref)-(ref) in the following manner.

assumptionIn addition to Assumption (ref), we require the derivatives of $f_{1}(x)$ and $f_{2}(x)$ to be bounded uniformly over $x$.

assumptionIn addition to Assumption (ref), we require the derivatives of the bivariate functions defined therein to belong to the set $\{f:\mathbb{R}^{p}\times\mathbb{R}^{p}\rightarrow\mathbb{R},\ \mathbb{E}_{P_{X}^{(1)}}|f(X,X)|<\infty,\ \mathbb{E}_{P_{X}^{(2)}}|f(X,X)|<\infty\}$.

Assumption (ref) is standard, as seen in previous studies such as lee2009non,wang2015conditional. Assumption (ref) is mild in that it does not require bounded derivatives of the conditional densities or moments. wang2015conditional assumes that the conditional densities have bounded derivatives. However, under this assumption, their proof implicitly requires the distance function $\|y_{1}-y_{2}\|$ to be integrable with respect to the Lebesgue measure (not some probability measure), which fails to hold unless the support is bounded. A careful examination of the arguments in wang2015conditional and ke2020expected suggests that conditions analogous to our Assumption (ref) are necessary.

theoremSuppose $P_{Y}^{(l)}\in\mathcal{M}_{\rho}^{2}(\mathcal{Y})$ for $l=1,2$. Under Assumptions (ref)-(ref) and (ref)-(ref), we have \[ \widehat{\mathcal{I}}_{\rho}\overset{p}{\longrightarrow}\mathcal{I}_{\rho},\quad\text{as }n_{1},n_{2}\rightarrow\infty. \]

theoremSuppose $P_{Y}^{(l)}\in\mathcal{M}_{\rho}^{2+\delta}(\mathcal{Y})$ for $l=1,2$ and some $\delta>0$. Under $H_{a}$ in ((ref)), Assumptions (ref)-(ref), (ref)-(ref) and $n_{l}^{1/2}h_{l}^{\nu}\rightarrow 0$ as $n_{l}\rightarrow\infty$ for $l=1,2$, we have \[ \frac{\widehat{\mathcal{I}}_{\rho}-\mathcal{I}_{\rho}}{2\sqrt{\frac{\delta_{1,0}^{2}}{n_{1}}+\frac{\delta_{0,1}^{2}}{n_{2}}}}\overset{d}{\longrightarrow}\mathcal{N}(0,1),\quad\text{as }n_{1},n_{2}\rightarrow\infty, \] where $\delta_{1,0}^{2}$ and $\delta_{0,1}^{2}$ are defined in equations ((ref))-((ref)) of Appendix (ref).

As shown in the proof, we have $\delta_{1,0}^{2}\asymp 1$ and $\delta_{0,1}^{2}\asymp 1$. Therefore, under $H_{a}$ in ((ref)), $\widehat{\mathcal{I}}_{\rho}-\mathcal{I}_{\rho}=O_{p}(n_{1}^{-1/2}+n_{2}^{-1/2})$. Again, undersmoothing is necessary, and the same condition has been used in lee2009non. Moreover, we emphasize that undersmoothing is not needed for establishing the consistency of $\widehat{\mathcal{D}}_{\rho}(x)$ and $\widehat{\mathcal{I}}_{\rho}$.

theoremSuppose $P_{Y}^{(l)}\in\mathcal{M}_{\rho}^{2+\delta}(\mathcal{Y})$ for $l=1,2$ and some $\delta>0$. Under $H_{0}$ in ((ref)), Assumptions (ref)-(ref), (ref)-(ref) and $n_{l}h_{l}^{p/2+\nu}\rightarrow 0$ as $n_{l}\rightarrow\infty$ for $l=1,2$, we have \[ \frac{\widehat{\mathcal{I}}_{\rho}}{\sqrt{\frac{2\delta_{2,0}^{2}}{n_{1}(n_{1}-1)}+\frac{16\delta_{1,1}^{2}}{n_{1}n_{2}}+\frac{2\delta_{0,2}^{2}}{n_{2}(n_{2}-1)}}}\overset{d}{\longrightarrow}\mathcal{N}(0,1),\quad\text{as }n_{1},n_{2}\rightarrow\infty, \] where $\delta_{2,0}^{2}$, $\delta_{1,1}^{2}$ and $\delta_{0,2}^{2}$ are defined in equations ((ref))-((ref)) of Appendix (ref).

As shown in the proof, we have $\delta_{2,0}^{2}\asymp h_{1}^{-p}$, $\delta_{1,1}^{2}\asymp h_{1}^{-p}+h_{2}^{-p}$ and $\delta_{0,2}^{2}\asymp h_{2}^{-p}$. Hence, the convergence rate of the sample ICED under $H_{0}$ in ((ref)) is $\widehat{\mathcal{I}}_{\rho}=O_{p}(n_{1}^{-1}h_{1}^{-p/2}+n_{2}^{-1}h_{2}^{-p/2})$. Because of the integration of $x$ over $\mathbb{R}^{p}$, the convergence rate of the global statistic $\widehat{\mathcal{I}}_{\rho}$ is faster than that of the local statistic $\widehat{\mathcal{D}}_{\rho}(x)$ under both the null and alternative hypotheses.

remarkInstead of a weighted sum of chi-squares, the global statistic $\widehat{\mathcal{I}}_{\rho}$ has a limiting normal distribution under the null. This distinct behavior follows from the martingale CLT (see, e.g., Corollary 3.1 of hall2014martingale). The asymptotic distributions of both $\widehat{\mathcal{D}}_{\rho}(x)$ and $\widehat{\mathcal{I}}_{\rho}$ are governed by the second-order projections in their respective Hoeffding decompositions. However, the local statistic $\widehat{\mathcal{D}}_{\rho}(x)$ does not satisfy a key condition required for the application of the martingale CLT, and thus behaves similar to its unconditional counterparts szekely2004testing,gretton2012kernel. In contrast, the global statistic $\widehat{\mathcal{I}}_{\rho}$ does meet the necessary conditions for the martingale CLT to hold. See also hall1984central for related insights.

Theorems (ref)-(ref) and (ref)-(ref) demonstrate that kernel smoothing has two key effects on the asymptotic behavior of distance and kernel-based measures:

itemize• Kernel smoothing introduces biases of order $O(h_{1}^{\nu}+h_{2}^{\nu})$. To ensure valid inference, undersmoothing or other bias reduction techniques are necessary to achieve asymptotic unbiasedness. For the local statistic $\widehat{\mathcal{D}}_{\rho}(x)$, Theorem (ref) and Assumption (ref) require the bandwidths to satisfy $h_{l}\asymp n_{l}^{\kappa_{1}}$, where $\kappa_{1}\in(-\frac{1}{p},-\frac{1}{p+\nu})$. Notably, there is no restriction on the relationship between the kernel order $\nu$ and the dimension $p$ in this case. In contrast, for the global statistic $\widehat{\mathcal{I}}_{\rho}$, Theorem (ref) and Assumption (ref) impose the condition $h_{l}\asymp n_{l}^{\kappa_{2}}$, with $\kappa_{2}\in(-\frac{1}{p},-\frac{1}{p/2+\nu})$. This requires $\nu>p/2$, indicating that higher-order kernels are necessary as the dimension $p$ increases. In particular, the use of a second-order kernel, as in wang2015conditional, becomes inadequate when $p\ge 4$ (see their Theorem 7). • As shown in the proof, due to kernel smoothing, the U-statistic estimators $\widehat{\mathcal{D}}_{\rho}(x)$ and $\widehat{\mathcal{I}}_{\rho}$ are nondegenerate under the null hypotheses. This stands in stark contrast to the unconditional setting, where the U-statistic estimators for energy distance and maximum mean discrepancy are degenerate under the null szekely2004testing,gretton2012kernel. Nonetheless, when undersmoothing is used, a more in-depth analysis reveals that the first-order projections in their Hoeffding decompositions are asymptotically negligible under the null, relative to the second-order projections. As a result, the asymptotic null distributions are governed by the second-order projections. To be specific, we have $\xi_{1}^{2}=O(h_{1}^{-p}(h_{1}^{2\nu}+h_{2}^{2\nu}))$ under $H_{0}$ in ((ref)), instead of $\xi_{1}^{2}\asymp h_{1}^{-p}$ under $H_{a}$ in ((ref)). Similarly, we have $\delta_{1,0}^{2}=O(h_{1}^{2\nu}+h_{2}^{2\nu})$ under $H_{0}$ in ((ref)), instead of $\delta_{1,0}^{2}\asymp 1$ under $H_{a}$ in ((ref)).

remarkOur work utilizes techniques similar to those used by wang2015conditional and ke2020expected, including U-statistics and kernel smoothing. However, we believe that there exists a gap in the derivations of the asymptotic null distribution in both studies. Specifically, Theorem 7 of wang2015conditional and Theorem 9 of ke2020expected both rely on Lemma B.4 of fan1996consistent, which establishes the CLT for degenerate U-statistics under certain conditions. To show that a U-statistic is degenerate, it must be demonstrated that $\mathcal{P}_{n1}(W_{1})=0$ under the null hypothesis, as per the notation used in their proofs. However, wang2015conditional merely showed that $\mathbb{E}\{\mathcal{P}_{n1}(W_{1})\}=O(h^{2})$, while ke2020expected incorrectly claimed that $\mathbb{E}\{\mathcal{P}_{n1}(W_{1})\}=0$. As a result, both studies overlook the crucial undersmoothing step necessary for valid inference.

Applications to global and local two-sample conditional distribution testing

Theorems (ref)-(ref) collectively suggest that $\widehat{\mathcal{I}}_{\rho}$ is a suitable test statistic for testing ((ref)). Nevertheless, it is impractical to use Theorem (ref) to compute the $p$-value since it is arduous to estimate $\delta_{2,0}^{2}$, $\delta_{1,1}^{2}$ and $\delta_{0,2}^{2}$. Moreover, it is widely acknowledged that a nonparametric test that relies on asymptotic normal approximation may perform poorly in finite samples su2007consistent. Consequently, we resort to the local bootstrap proposed by paparoditis2000local. This approach has been widely employed in other works involving conditional distributions; see, e.g., su2008nonparametric,huang2010testing,bouezmarni2012nonparametric,su2013nonparametric,taamouti2014nonparametric,wang2015conditional. One can follow the aforementioned references to verify the asymptotic validity of this bootstrap method in our framework. Define

equation[equation omitted — 208 chars of source]

where $\delta_{y}$ denotes a point mass at $y\in\mathcal{Y}$, and $b_1, b_2\in\mathbb{R}$ are the bandwidth parameters for local bootstrap. Essentially, $\widehat{P}_{Y \mid X = x}$ is a discrete distribution that assigns the probability $G_{b_l}(X_i^{(l)}-x) / \sum_{l=1}^{2}\sum_{i=1}^{n_l}G_{b_l}(X_i^{(l)}-x)$ to the observation $Y_i^{(l)}$. Then the following steps outline the procedure for the global two-sample conditional distribution test:

enumerate[(i)] • Calculate the global test statistic, i.e., the sample ICED, $\widehat{\mathcal{I}}_{\rho}$. • For $l=1,2$ and $i=1,\ldots,n_{l}$, draw $\widehat{Y}_{i}^{(l)}$ from $\widehat{P}_{Y \mid X = X_i^{(l)}}$. Calculate the global test statistic using the local bootstrap samples $\{(\widehat{Y}_{i}^{(l)},X_{i}^{(l)})\}_{i=1}^{n_{l}}\ (l=1,2)$, which retains the marginal distributions of $X$ but imposes the null restriction (i.e., the same conditional distribution of $Y$ given $X$). • Repeat step (ref) $B$ times, and collect $\{\widehat{\mathcal{I}}_{\rho}^{(t)}\}_{t=1}^{B}$. Then, the bootstrap-based $p$-value of the global test is given by \[ p\text{-value}=\frac{\sum_{t=1}^{B}\mathbbm{1}\{\widehat{\mathcal{I}}_{\rho}^{(t)}>\widehat{\mathcal{I}}_{\rho}\}+1}{B+1}. \]

Theorems (ref)-(ref) support the use of $\widehat{\mathcal{D}}_{\rho}(x)$ as the test statistic for testing ((ref)). Unfortunately, the asymptotic null distribution of $\widehat{\mathcal{D}}_{\rho}(x)$ in Theorem (ref) is not pivotal and involves infinite nuisance parameters. To tackle this issue, we once again utilize the local bootstrap method to calculate the $p$-value, replacing the $\widehat{\mathcal{I}}_{\rho}$ in the above steps (ref)-(ref) with $\widehat{\mathcal{D}}_{\rho}(x)$.

For the smoothing bandwidth parameters $h_l$ in the test statistics, we introduce a completely data-driven selection approach, which works well in our numerical studies. Specifically, for $l = 1, 2$, we determine $h_l$ by minimizing the cross-validation loss function \[ \operatorname{CV}(h_l(1),\ldots,h_{l}(p)) = \sum_{i=1}^{n_l} \| k(Y_i^{(l)}, \cdot) - \widehat{\Pi}_{k}^{-i}(P_{Y \mid X=X_{i}^{(l)}}^{(l)}) \|_{\mathcal{H}_k}^{2}, \] where $\widehat{\Pi}_{k}^{-i}(P_{Y \mid X=x}^{(l)}) = \{\sum_{i'\ne i} G_{h_l}(X_{i'}^{(l)} - x)\}^{-1} \sum_{i'\ne i} G_{h_l}(X_{i'}^{(l)} - x) k(Y_{i'}^{(l)}, \cdot)$ is the leave-one-out Nadaraya-Watson estimator of the conditional kernel mean embedding $\Pi_{k}(P_{Y \mid X=x}^{(l)})$. The minimization problem is solved using the limited-memory quasi-Newton method with box constraints byrd1995limited, implemented via the optim(..., method = "L-BFGS-B") function in R. The same approach is also applicable to the distance-based test statistics due to the equivalence between the CED and CMMD established in Theorem (ref). For the bandwidth parameters $b_l$ used in the local bootstrap, following the aforementioned references, we adopt the rule of thumb $b_l(s) = \min \{ \widehat{\sigma}_{X^{(l)}(s)}, \mathrm{IQR}_{X^{(l)}(s)} / 1.34 \} \cdot n_l^{-1/(p + 2\nu)} $ for $s = 1, \dots, p$, where $\widehat{\sigma}_{X^{(l)}(s)}$ and $\mathrm{IQR}_{X^{(l)}(s)}$ are the sample standard deviation and interquartile range of $X^{(l)}(s)$, respectively.

Following Remark (ref), the time complexity of both our global and local tests is $O(B(n_{1}+n_{2})^{2})$. This is the same as that of unconditional two-sample tests based on permutation approach when the number of random permutations is $B$ szekely2004testing,gretton2012kernel.

Numerical studies

We conduct numerical studies to evaluate the performance of our proposed tests in this section. The code that implements our tests and reproduces all numerical results is available at \url{https://github.com/lizhuoxi-97/TCDT}.

Monte Carlo simulations

In this subsection, we assess the finite-sample performance of our proposed methods based on the conditional energy distance (denoted as CED) and conditional maximum mean discrepancy (denoted as CMMD) through simulation examples. In Examples (ref)-(ref), we examine CED and CMMD in testing the local problem (ref), for which no competing methods currently exist. In Examples (ref)-(ref), we evaluate them in testing the global problem (ref), comparing them against the conformal prediction test of hu2024two (denoted as CONF).

In our proposed tests, we specify the semimetric $\rho$ for CED as the Euclidean distance. For the reproducing kernel $k$ in CMMD, we use the Gaussian kernel $k(y, y') = \exp\{-\|y - y'\|^2/(2\gamma^2)\}$, with the bandwidth $\gamma$ determined by the median heuristic gretton2012kernel,ke2020expected. For the smoothing kernel $g$, we adopt the Gaussian kernel. In Examples (ref)-(ref) and (ref), we use a standard second-order ($\nu = 2$) kernel $g(u) = (2\pi)^{-1/2} \exp(-u^2 / 2)$. In Example (ref), we employ a fourth-order ($\nu = 4$) kernel $g(u) = (3/2 - u^2/2) (2\pi)^{-1/2} \exp(-u^2 / 2)$ to handle the higher dimension $p = 4$, as discussed in Section (ref). We calculate the $p$-value of the proposed methods via the local bootstrap procedure with $B=299$ replications and the Gaussian kernel $g(u) = (2\pi)^{-1/2} \exp(-u^2 / 2)$ in ((ref)). The selection of bandwidth parameters, $h_l$ and $b_l$, is provided in Section (ref).

For CONF, we consider an equal data-splitting ratio and estimate the marginal and conditional density ratios by kernel logistic regression. The tuning parameter $\sigma^2$ therein is set to $200$, and $\lambda$ is selected using the out-of-sample cross entropy loss, as hu2024two suggested.

exampleTo examine the performance of the proposed local tests for problem (ref), we focus on the general regression models \begin{align} Y^{(l)} = m^{(l)}(X^{(l)}) + \epsilon^{(l)} \quad (l = 1,2), \end{align} where $m^{(l)}$ is the regression function and $\epsilon^{(l)}$ is the error satisfying $\mathbb{E}(\epsilon^{(l)}\mid X^{(l)})=0$. Let $X^{(1)} \sim \mathcal{N}(0, 1)$ and $X^{(2)} \sim t(5)$, that is, $P_X^{(1)} \ne P_X^{(2)}$. The error $\epsilon^{(l)} \sim \mathcal{N}(0, 0.5)$ is independent of $X^{(l)}$ for both $l=1,2$. For the regression functions, we consider $m^{(1)}(x) = 0.5 x$ and $m^{(2)}(x) = - 0.5 x$. It is under the null hypothesis only at $x = 0$, and the signal for the difference between the two samples gets larger when the magnitude of local point $x$ is larger. We set $n_1 = n_2 = 150$ and conducted the proposed tests at local points ranging from $-1$ to $1$, spaced by $0.1$. For each point, we obtained empirical power or size from 1000 simulations at a significance level of $\alpha=0.05$. Figure (ref) summarizes the results. At $x = 0$, the empirical sizes for CED and CMMD are $0.046$ and $0.056$, respectively, which are controlled around the nominal level. The power of both CED and CMMD increases as the local signal strengthens. \begin{figure}[htbp] \caption{ Empirical power curves under $\alpha = 0.05$, and properly scaled signals in Example (ref). } \end{figure}

exampleWe further illustrate how the local tests perform in a single data generation. Let $X^{(1)}$ follow a bivariate truncated standard normal distribution with support on $[-2, 2] \times [-2, 2]$, and let $X^{(2)} \sim \mathcal{U}(-2, 2) \times \mathcal{U}(-2, 2)$. Define $Y^{(1)} = \|X^{(1)}\|^2 + \epsilon^{(1)}$ and $Y^{(2)} = \|X^{(2)}\|^2 \mathbbm{1}\{ X^{(2)}(1) X^{(2)}(2) \ge 0 \} - \|X^{(2)}\|^2 \mathbbm{1}\{ X^{(2)}(1) X^{(2)}(2) < 0 \} + \epsilon^{(2)}$, where $\epsilon^{(l)} \sim \mathcal{N}(0, 1)$ is independent of $X^{(l)}$ for $l = 1, 2$. In this setting, the null hypothesis in (ref) holds whenever the local point $x$ lies in the first or third quadrant, and is violated when $x$ lies in the second or fourth quadrant. We set $n_1 = n_2 = 500$ and conducted the local tests on a fixed uniform grid of $21 \times 21$ points over $(x(1), x(2)) \in [-2, 2] \times [-2, 2]$ at a significance level $\alpha=0.05$. The results, based on a single data generation, are shown in Figure (ref). It demonstrates that the proposed tests correctly identify the regions of the alternative hypothesis where the signal is strong enough. The few false rejections in the first and third quadrants are due to the randomness of generating data a single time. \begin{figure}[htbp] \caption{ Significant local regions for Example (ref), obtained by the proposed tests (CED and CMMD) under $\alpha = 0.05$ based on one-time data generation. The rejection regions are filled in red, whereas the non-significant regions are filled in gray. } \end{figure}

exampleThis example assesses the performance of the proposed global tests for problem (ref). We adopt the general regression models (ref) from Example (ref). Let $X^{(1)} \sim \mathcal{N}(0, 1)$ and $X^{(2)} \sim t(5)$, resulting in $P_X^{(1)} \ne P_X^{(2)}$. In addition to CONF, we compare our methods with the regression curve testing via empirical characteristic function pardo2015nonparametric (denoted as ECF), which is applicable in the univariate case. We first consider the scenarios in which the two conditional distributions differ solely in the conditional means (i.e., $m^{(l)}$). For both $l=1,2$, let $\epsilon^{(l)} \sim \mathcal{N}(0, 0.5)$ be independent of $X^{(l)}$. Two settings of $m^{(l)}$ are considered as follows: \begin{setting} $m^{(1)}(x) = 1 + x^2$, $m^{(2)}(x) = 1 + x^2 + \frac{1}{5}c x(x+2)(x-2)$. \end{setting} \begin{setting} $m^{(1)}(x) = 1 + \exp(x)$, $m^{(2)}(x) = 1 + \exp(x) + c \sin(2\pi x)$. \end{setting} Here, $c$ is a parameter that controls the strength of the signal. The case $c = 0$ corresponds to the null hypothesis, while a larger $c$ indicates a greater difference. Besides, we consider the scenarios in which the two conditional distributions differ only in the conditional variances. Let the regression functions $m^{(1)}(x) = m^{(2)}(x) = 1 + x^2$ be identical across the two samples. The following two cases are considered: \begin{setting} $\epsilon^{(1)} \sim \mathcal{N}(0, 0.25)$, $\epsilon^{(2)} \sim \mathcal{N}(0, 0.25 + c)$. \end{setting} \begin{setting} $\epsilon^{(1)} \sim \mathcal{N}(0, 0.25^2 \exp(2X^{(1)}))$, $\epsilon^{(2)} \sim \mathcal{N}(0, (0.25 + 0.5c)^2 \exp(2X^{(2)}))$. \end{setting} We set $n_1 = n_2 = 50$ or $100$ for Settings (ref)-(ref), and set $n_1 = n_2 = 100$ or $200$ for Settings (ref)-(ref). The signal strength $c$ is varied from $0$ to $2$ in Settings (ref)-(ref), and varied from $0$ to $1$ in Settings (ref)-(ref). Table (ref) and Figure (ref) summarize the empirical sizes and powers at a significance level $\alpha = 0.05$ via 1000 simulations. \begin{table}[htbp] \caption{Empirical size ($c = 0$) under $\alpha = 0.05$ in Example (ref). } \begin{tabular}{clrrrrr} \toprule Setting & \multicolumn{1}{c}{Sample Size} & & \multicolumn{1}{c}{CED} & \multicolumn{1}{c}{CMMD} & \multicolumn{1}{c}{CONF} & \multicolumn{1}{c}{ECF} \\ \midrule \multicolumn{7}{c}{Difference in conditional means} \\ \midrule \multirow{2}[2]{*}{(ref)} & $n_1 = n_2 = 50$ & & 0.046 & 0.053 & 0.074 & 0.113 \\ & $n_1 = n_2 = 100$ & & 0.049 & 0.049 & 0.066 & 0.081 \\ \midrule \multirow{2}[2]{*}{(ref)} & $n_1 = n_2 = 50$ & & 0.043 & 0.046 & 0.043 & 0.101 \\ & $n_1 = n_2 = 100$ & & 0.043 & 0.046 & 0.037 & 0.084 \\ \midrule \multicolumn{7}{c}{Difference in conditional variances} \\ \midrule \multirow{2}[2]{*}{(ref)} & $n_1 = n_2 = 100$ & & 0.042 & 0.043 & 0.071 & 0.094 \\ & $n_1 = n_2 = 200$ & & 0.040 & 0.037 & 0.093 & 0.080 \\ \midrule \multirow{2}[2]{*}{(ref)} & $n_1 = n_2 = 100$ & & 0.027 & 0.027 & 0.066 & 0.042 \\ & $n_1 = n_2 = 200$ & & 0.023 & 0.017 & 0.082 & 0.030 \\ \bottomrule \end{tabular} \end{table} \begin{figure}[htbp] \caption{Empirical size ($c = 0$) and power ($c\neq 0$) under $\alpha = 0.05$ in Example (ref). } \end{figure} It can be observed that CONF and ECF generally yield empirical sizes higher than 0.05, whereas both CED and CMMD control the empirical size well across all scenarios. A possible reason for \texttt{CONF}'s inflated size is its reliance on high-quality density ratio estimators and a potential requirement for unbalanced sample splitting; see Assumption 2(b) and discussion in hu2024two. Although \texttt{ECF} is applicable to general distributions, its size control is often unsatisfactory when the support of $X^{(l)}$ is not the interval $[0, 1]$. In Settings (ref)-(ref) where differences in complex nonlinear regression functions are considered, the power of \texttt{CED}, \texttt{CMMD} and \texttt{ECF} are notably greater than that of \texttt{CONF}. In Settings (ref)-(ref), the proposed tests also show much higher power than \texttt{CONF}, while \texttt{ECF} shows trivial power as it is designed for comparing regression curves. The satisfactory performance of \texttt{CED} and \texttt{CMMD} in all scenarios highlights their applicability to various types of alternative hypotheses.

exampleSimilar to Example (ref), we now consider the setting where the dimension of $X$ is $p = 4$. The data generating process again follows the regression models (ref). Let $X^{(1)} \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$ and $X^{(2)} \sim \mathcal{N}(\boldsymbol{\mu}, \mathbf{I})$, where $\boldsymbol{\mu} = (1, 1, -1, 0)^{\top}$. For both $l=1,2$, the error $\epsilon^{(l)} \sim t(5)$ has a heavy tail and is independent of $X^{(l)}$. We set $m^{(1)}(x) = \{\theta(x)\}^2$ and $m^{(2)}(x) = \{\theta(x)\}^2 + c \theta(x)$, where $\theta(x) = x(1) + x(2) + x(3) + x(4) + 1$ and $c$ is the signal strength. We let $n_1 = n_2 = 100$ or $200$ and vary $c$ from $0$ to $2$. The simulations are conducted 1000 times at a significance level of $\alpha=0.05$. Table (ref) and Figure (ref) present the empirical sizes and powers. In this multivariate $X$ setting, both CED and CMMD are observed to achieve better power performance compared to CONF, while guaranteeing type I error control. \begin{table}[htbp] \caption{Empirical size ($c = 0$) under $\alpha = 0.05$ in Example (ref).} \begin{tabular}{lrrrr} \toprule \multicolumn{1}{c}{Sample Size} & & \multicolumn{1}{c}{CED} & \multicolumn{1}{c}{CMMD} & \multicolumn{1}{c}{CONF} \\ \midrule $n_1 = n_2 = 100$ & & 0.016 & 0.025 & 0.037 \\ $n_1 = n_2 = 200$ & & 0.052 & 0.049 & 0.040 \\ \bottomrule \end{tabular} \end{table} \begin{figure}[htbp] \caption{Empirical size ($c = 0$) and power ($c\neq 0$) under $\alpha = 0.05$ in Example (ref).} \end{figure}

exampleWe further consider a scenario that appears in transfer learning: testing whether the conditional distribution of covariates given the response, $P_{X\mid Y}$, is the same for two populations (i.e., prior shift assumption). To generate data, we adopt a multivariate regression model of $X$ on $Y$, conversely to (ref): \begin{align*} X^{(l)} = m^{(l)}(Y^{(l)}) + \epsilon^{(l)} \quad (l = 1, 2), \end{align*} where $m^{(l)}: \mathbb{R} \to \mathbb{R}^{p}$ is the multivariate regression function and $\epsilon^{(l)} \in \mathbb{R}^{p}$ is the multivariate error satisfying $\mathbb{E}(\epsilon^{(l)}\mid Y^{(l)})=\mathbf{0}$. Let $\epsilon^{(l)} \sim \mathcal{N}(\mathbf{0}, \pmb{\Sigma})$ be independent of $Y^{(l)}$ for both $l = 1, 2$, where $\pmb{\Sigma} = (\rho^{|i - j|})$ with $\rho = 0.5$. We set $Y^{(1)} \sim \mathcal{N}(0, 1)$ and $Y^{(2)} \sim \mathcal{N}(0, 1.5)$, that is, $P_Y^{(1)} \ne P_Y^{(2)}$. Let $\theta^{(1)}(y) = y \cdot \beta$ and $\theta^{(2)}(y) = (1 + c) y \cdot \beta$, where $\beta \in \mathbb{R}^{p}$, and $c$ controls the signal for the difference of the conditional distributions between the two samples. We then set \begin{equation*} m^{(l)}(y) = \left [\{\theta^{(l)}(y)\}^{\circ 2} + 3 \theta^{(l)}(y) + 2 \right ] \div \left [\{\theta^{(l)}(y)\}^{\circ 2} + 1 \right ]\quad (l=1,2), \end{equation*} where $^{\circ 2}$ denotes the Hadamard square and $\div$ denotes the Hadamard division. We consider $p = 5$ or $20$. For $p = 5$, we set $\beta(j) = -0.5$ for $j = 1, 2$, $\beta(j) = 0.5$ for $j = 3, 4$ and $\beta(j) = 1$ for $j = 5$. For $p = 20$, we set $\beta(j) = -0.5$ for $j = 1, 2, 3$, $\beta(j) = 0.5$ for $j = 4, 5$, $\beta(j) = 1$ for $j = 6, \dots, 10$ and $\beta(j) = 0$ for $j = 11, \dots, 20$. Let $n_1 = n_2 = 100$ or $200$, and vary the signal strength $c$ from $0$ to $1$. For each setting, we repeat 1000 simulations to obtain the empirical size ($c = 0$) and power ($c \ne 0$) of CMMD, CED and CONF at a significance level of $\alpha=0.05$. The results are presented in Table (ref) and Figure (ref). Table (ref) demonstrates that all three methods control the type I error rate. The power comparison, illustrated in Figure (ref), clearly shows an advantage of CED and CMMD over CONF. \begin{table}[htbp] \caption{Empirical size ($c = 0$) under $\alpha = 0.05$ in Example (ref). } \begin{tabular}{lrrrrrrrr} \toprule & & \multicolumn{3}{c}{$p = 5$} & & \multicolumn{3}{c}{$p = 20$} \\ \cmidrule{3-5}\cmidrule{7-9}\multicolumn{1}{c}{Sample Size} & & \multicolumn{1}{c}{CED} & \multicolumn{1}{c}{CMMD} & \multicolumn{1}{c}{CONF} & & \multicolumn{1}{c}{CED} & \multicolumn{1}{c}{CMMD} & \multicolumn{1}{c}{CONF} \\ \midrule $n_1 = n_2 = 100$ & & 0.026 & 0.026 & 0.051 & & 0.037 & 0.037 & 0.039 \\ $n_1 = n_2 = 200$ & & 0.050 & 0.054 & 0.049 & & 0.033 & 0.033 & 0.019 \\ \bottomrule \end{tabular} \end{table} \begin{figure}[htbp] \caption{Empirical size ($c = 0$) and power ($c\neq 0$) under $\alpha = 0.05$ in Example (ref). } \end{figure}

Airfoil data example

We consider the airfoil dataset from the UCI Machine Learning Repository airfoil. This dataset, collected by NASA to study airfoil sound pressure, has been previously analyzed in tibshirani2019conformal,hu2024two,huang2024efficient. It contains $N=1503$ observations, with a response $Y$ (scaled sound pressure level) and covariates $X$ with $p=5$: log-frequency, angle of attack, chord length, free-stream velocity, and suction-side log-displacement thickness.

As the dataset consists of a single sample only, we adopt a semi-synthetic approach, similar to that of tibshirani2019conformal,hu2024two,huang2024efficient, to create two-sample scenarios. We generate four different settings by partitioning the data:

enumerate• Covariate shift under the null hypothesis: The dataset is randomly partitioned into $\mathcal{D}_1$ with $n_1 = 301$ and $\mathcal{\widetilde{D}}_2$ with $\tilde{n}_2 = 1202$. The second sample, $\mathcal{D}_2$, is formed by sampling $n_2=301$ points from $\mathcal{\widetilde{D}}_2$ without replacement, with probabilities proportional to $w(X) = \exp(\alpha^\top X)$, where $\alpha = (-1, 0, 0, 0, 1)^\top$. This induces a covariate shift while preserving $P_{Y\mid X}$. • Prior shift under the null hypothesis: The procedure is identical to Setting (ref), but sampling probabilities are instead proportional to the response, $w(Y) = Y$. This induces a prior shift while preserving $P_{X\mid Y}$. • Covariate shift under the alternative hypothesis: The dataset is partitioned into two groups based on the median of the response variable. The first sample contains $n_1 = 752$ points with smaller values of $Y$, while the second sample contains the rest of $n_2 = 751$ points. To avoid singularity between two groups, we randomly select $0.05n_1$ points in each of the two samples and then flip their groups. This construction is designed to violate the equality of $P_{Y\mid X}$. • Prior shift under the alternative hypothesis: The setup is almost the same as Setting (ref), but the partition is based on the median of the “chord length" variable, creating a scenario where $P_{X\mid Y}$ differs.

For each setting, we apply our global tests at a significance level of $\alpha = 0.05$. In the covariate shift scenarios (Settings (ref) and (ref)), we test the equality of $P_{Y\mid X}$ using a Gaussian smoothing kernel of $\nu = 4$, as here $p=5$. In the prior shift scenarios (Settings (ref) and (ref)), we test the equality of $P_{X\mid Y}$ with a Gaussian smoothing kernel of $\nu = 2$. All other implementation details are the same as those in Section (ref).

The results are summarized in Table (ref). Under the null hypothesis (Settings (ref)-(ref)), we performed 500 replications and report the empirical rejection rates. The type I error of both CED and CMMD are close to the nominal 0.05 level. Under the alternative hypothesis (Settings (ref)-(ref)), which involve deterministic partitions, we report the $p$-values from a single test. In both cases, the tests correctly reject $H_0$ with very small $p$-values.

table[table omitted — 842 chars of source]

Furthermore, we apply our local tests to Setting (ref). Specifically, we test for local differences in $P_{X\mid Y=y}$ for values of the response $y$ ranging from 110 to 130. Figure (ref) plots the resulting $p$-values against $y$. The $p$-values remain consistently below the 0.05 significance level across the entire range, in agreement with the setup of Setting (ref).

figure[figure omitted — 371 chars of source]

Discussions

Several research directions warrant further exploration. First, it would be valuable to extend the framework to the $K$-sample conditional distribution testing problem for $K\geq 2$ and establish connections with distance-based multivariate analysis of variance rizzo2010disco. Beyond kernel smoothing, alternative machine learning techniques, such as random forest and neural networks, could be considered for estimating the CED and ICED. These methods are anticipated to yield improved performance for higher-dimensional $X$, but deriving the asymptotic distributions of the resulting test statistics presents significant theoretical challenges. Finally, for local testing problems, especially when inference at a boundary point is of interest, a statistic based on local polynomial fitting may be preferable to the sample CED. Such an approach could be particularly useful in the regression discontinuity design calonico2014robust.