In this paper we introduce a kernel-based measure for detecting differences between two conditional distributions. Using the `kernel trick' and nearest-neighbor graphs, we propose a consistent estimate of this measure which can be computed in nearly linear time (for a fixed number of nearest neighbors). Moreover, when the two conditional distributions are the same, the estimate has a Gaussian limit and its asymptotic variance has a simple form that can be easily estimated from the data. The resulting test attains precise asymptotic level and is universally consistent for  detecting differences between two conditional distributions. We also provide a resampling based test using our estimate that applies to the conditional goodness-of-fit problem, which controls Type I error in finite samples and is asymptotically consistent with only a finite number of resamples. A method to de-randomize the resampling test is also presented.   The proposed methods can be readily applied to a broad range of problems, ranging from classical nonparametric statistics to modern machine learning. Specifically, we explore three applications: testing model calibration, regression curve evaluation, and validation of emulator models in simulation-based inference. We illustrate the superior performance of our method for these tasks, both in simulations as well as on real data. In particular, we apply our method to (1) assess the calibration of neural network models trained on the CIFAR-10 dataset, (2) compare regression functions for wind power generation across two different turbines, and (3) validate emulator models on benchmark examples with intractable posteriors and for generating synthetic `redshift' associated with galaxy images.

§ INTRODUCTION

The conditional 2-sample problem is to test whether the conditional distributions of two response variables $X$ and $Y$ are the same, given a set of covariates $Z$. More formally, given independent samples $(X_1,Y_1,Z_1), (X_2,Y_2,Z_2),\ldots,(X_n,Y_n,Z_n)$ from a joint distribution $ P_{XYZ}$, we want to test the following null hypothesis:

\begin{align}\label{eq:H0H1}
\bH_0: P_{X|Z} =  P_{Y|Z} \text{ almost everywhere } P_Z ,
\end{align}

where $P_{X|Z}$, $P_{Y|Z}$ denote the conditional distributions of $X|Z$ and $Y|Z$, respectively, and $P_Z$ is the marginal distribution of $Z$. This problem appears in different guises in the literature, both classical and contemporary. The following are three examples:

(Calibration tests for predictive models)
Many problems in supervised learning involve estimating a conditional probability distribution $P_{Y|Z}$ of a response $Y$ given a feature vector $Z$. For models which output such predictive distributions, statistical guarantees  beyond accuracy are necessary for quantifying predictive uncertainties. One such guarantee is calibration, which, loosely stated, ensures that almost every prediction matches the conditional distribution of the response given this prediction. Calibration has been used to provide predictive guarantees in meteorological and statistical studies for many years [DeGroot and Fienberg, 1983, Murphy and Winkler, 1977].
This notion has seen a resurgence in modern machine learning, following the breakthrough work of Guo et al., 2017, which showed that
common neural network architectures trained on image and text data are often miscalibrated. This prompted a slew of work on different notions of calibration for classification (where $Y$ is categorical) [DeGroot and Fienberg, 1983, Murphy and Winkler, 1977, Kull and Flach, 2015, Kumar et al., 2018, Guo et al., 2017, Wenger et al., 2020, Widmann et al., 2019, Lee et al., 2023, Vaicenavicius et al., 2019] and regression problems (where $Y$ is continuous)
[Ho and Lee, 2005, Rueda et al., 2007, Gneiting et al., 2007, Taillardat et al., 2016, Song et al., 2019]. Specifically, for a $r$-class classification problem one can define calibration as follows [Widmann et al., 2019]: Consider a pair of random of variables $(Y, Z)$ with joint distribution $P_{YZ}$ over a space $\cY \times \cZ$, with $\cY = \{1,2,\ldots, r\}$ and a predicted model $f:\cZ\ra \cP_{\cY}$, where

\begin{align*}
\cP_{\cY} = \left\{\bp = (p_{1},p_{2},\ldots, p_{r}):0\leq p_{s}\leq 1, \sum_{s=1}^{r}p_{s}=1\right\} ,
\end{align*}

is the set of distributions over $\cY$. Then $f$ is said to be calibrated  if and only if, for $1 \leq s \leq r$,

\begin{align}\label{eq:classcalibration}
\P\left[Y = s \middle| f(Z) = \bp\right] = p_{s} \quad \text{ for all } \bp\in f(\cZ)\subseteq\cP_{\cY} ,
\end{align}

where $f(Z)$ is the vector of predicted probabilities. Calibration testing can be formulated as a conditional 2-sample problem by noting that if $X\sim \mathrm{Multinomial}(r, f(Z))$, then by definition of the multinomial distribution: $\P\left[X = s \middle| f(Z) = \bp \right] = p_{s}$, for $1 \leq s \leq r$.
Hence, from (<ref>), we say the predictive model $f$ is calibrated if and only if,

\begin{align}\label{eq:classification}
P_{Y|f(Z)} =  P_{X|f(Z)} \text{ almost everywhere } P_Z.
\end{align}

Note that $f$ is a known trained model, hence given data $(Y_1, Z_1), (Y_2, Z_2), \ldots, (Y_n, Z_n)$ from $P_{YZ}$, it is usually easy to draws samples $X_i \sim \mathrm{Multinomial}(r, f(Z_i))$, for $1 \leq i \leq n$. Hence, (<ref>) can be implemented as a conditional 2-sample problem based on the samples $(X_1,Y_1,f(Z_1)),\ldots,(X_n,Y_n,f(Z_n))$. For details on the implementation, comparison with existing methods, and extension to regression problems, see Section <ref>.

(Comparing regression curves) A prototypical example of conditional 2-sample testing is the problem of comparing two regression curves. The curves usually correspond to the mean functions of a control and a treatment outcome given a collection of covariates $Z$. For instance, suppose one has $n$ pairs of independent observations,

\begin{align}\label{eq:compH0}
(X_i, Y_i) = (f(Z_i) + \vep_i, g(Z_i) + \delta_i),
\end{align}

for $1\leq i\leq n$, where  $\vep_{1},\ldots,\vep_n$ and $\delta_1,\ldots,\delta_n$ are i.i.d. error variables. Then the hypothesis of no treatment effect, that is, testing whether or not $f=g$ is equivalent to (<ref>). Note that when $f$ and $g$ are linear and the errors are normally distributed, this is the classical analysis of covariance problem [Hocking, 2013]. This problem is also well-studied the nonparametric setting (<ref>), beginning with the works of Hardle and Marron, 1990 and King et al., 1991, where tests based on the kernel density estimates of the regression functions were proposed. Many variations have been considered over the years, focusing primarily on testing equality of conditional moments  (see [Dette and Munk, 1998, Neumeyer and Dette, 2003, Chaudhuri and Marron, 1999, Hardle and Marron, 1990, Kulasekera, 1995, Pardo-Fernández et al., 2015, Fan and Lin, 1998] and the references therein). The framework considered in this paper extends beyond comparing the mean and variance functions, to the entire conditional distribution of the responses given the covariates. For instance, suppose
$X = f(Z,\vep)$ and $Y = g(Z,\delta)$,
where $\vep,\delta$ are independent and identically distributed error variables. Then testing the hypothesis $f=g$ is equivalent to $\bH_0$ in (<ref>) (see Section <ref> for a specific example).

(Validation of emulators models in Simulation-Based  Inference)
Likelihood Free Inference (LFI) and Simulation-Based  Inference (SBI) broadly refer to the collection of methods that use simulations to infer the posterior in situations where the likelihood function is intractable. This approach has found widespread success in many scientific domains (see Cranmer et al., 2020 and the references therein). In a typical setup, the simulator takes as input a vector of parameters $\theta$ generated from a prior $p(\theta)$, samples a set of latent variable $\eta\sim p(\eta|\theta)$ and generates a data vector $X\sim p(X|\eta,\theta)$. In this case, the likelihood function $p(X|\theta)$ is implicitly defined as

\begin{align*}
p(X|\theta) = \int p(X|\eta,\theta)p(\eta)\d\eta ,
\end{align*}

which is intractable for most real-life simulators. To overcome this issue, several methods like Approximate Bayesian Computation [Marin et al., 2012] and Synthetic Likelihood [Wood, 2010] have been proposed. These methods make repeated calls to the simulator and use the simulated data to provide an estimate to the posterior distribution of $\theta$.
However, the accuracy of the estimate depends on the number of simulation calls, leading to increased computational cost, especially for expensive simulators. Additionally, inference chain depends on the choice of hyperparameters and low dimensional summary statistics, which can potentially reduce the quality of inference. Moreover, directly using the simulator presents a lack of amortization, that is, the inference chain has to be restarted every time new data is available. To mitigate these issues, one often trains a faster surrogate or emulator $q(\cdot|\theta)$ for the computationally expensive simulator. Examples of emulator models include Gaussian mixture density networks [Lueckmann et al., 2019], density ratio estimators [Hermans et al., 2020, Durkan et al., 2020, Dinev and Gutmann, 2018], and, more recently, neural conditional density estimators, such as normalizing flows [Trippe and Turner, 2018, Glöckler et al., 2022, Dirmeier et al., 2023], autoregressive models [Hansen, 1994, Papamakarios et al., 2017, Bruinsma et al., 2023, Papamakarios et al., 2019],  conditional diffusion models [Batzolis et al., 2021, Shi et al., 2022], and conditional deep generative networks [Zhou et al., 2022, Liu et al., 2021, Song et al., 2023, Alfonso et al., 2023], among others. To evaluate the validity of the emulator one needs to understand the extent to which it can imitate the simulator.  The basic diagnostic check towards this is to test the following hypothesis:

\begin{align}\label{eq:SBItest}
p(\cdot|\theta) = q(\cdot|\theta)\text{ almost everywhere }p(\theta) ,
\end{align}

This can be implemented as an instance of the conditional 2-sample hypothesis (<ref>) as follows: Given independent samples $\{\theta_{1}, \theta_{2}, \ldots, \theta_{n}\}$ from the prior $p(\theta)$, generate a sample $X_{i}\sim p(\cdot|\theta_{i})$ from the simulator and a sample $Y_{i}\sim q(\cdot|\theta_{i})$ from the emulator, for $1 \leq i \leq n$. Then we can test the hypothesis (<ref>) using the samples $(X_{1}, Y_{1}, {\theta}_{1}), (X_{2}, Y_{2}, {\theta}_{2}), \ldots, (X_{n}, Y_{n}, {\theta}_{n})$ (see, for example, [Dalmasso et al., 2020, Dalmasso et al., 2020, Dalmasso et al., 2021, Lemos et al., 2023]).

In addition to emulating the simulator to provide posterior inferences about the parameter $\theta$, there have been approaches to directly emulate the true posterior from observed samples without additional calls to the simulator. This method offers an amortized way to estimate the posterior distribution. It allows for faster analysis in situations where time is limited or a large amount of data needs to be processed [Gonçalves et al., 2020, Dax et al., 2021]. Moreover, it facilitates rapid application of diagnostic techniques that rely on obtaining posterior samples for numerous observations [Cook et al., 2006, Talts et al., 2018]. Popular approaches towards approximating the posterior include mixture density network [Papamakarios and Murray, 2016, Lueckmann et al., 2017], masked autoregressive flows [Greenberg et al., 2019, Deistler et al., 2022], normalising flows [Papamakarios et al., 2021, Rodrigues et al., 2021, Wiqvist et al., 2021], generative adversarial networks [Ramesh et al., 2022], and diffusion models [Gloeckler et al., 2024, Linhart et al., 2024], among others. Once again the basic diagnostic check towards understanding the validity of approximate posterior is to test,

\begin{align}\label{eq:sbiposterior}
p\left(\cdot|X\right) = q\left(\cdot|X\right)\text{ almost everywhere } p(X),
\end{align}

where $p(\cdot|X)$ is the true posterior distribution and $q(\cdot|X)$ is the proposed approximation. Given independent samples $(X_1, \theta_1), \ldots, (X_n, \theta_n)$ from the joint distribution $p(X, \theta)$ the test from (<ref>) can be considered as an instance of (<ref>) as follows: Generate a sample $\theta_i^\prime\sim q(\cdot|X_i)$ from the emulator for $1\leq i\leq n$. Now, we can test the hypothesis (<ref>) in a similar manner as (<ref>) using the samples $(\theta_1^\prime, \theta_1,X_1),\ldots(\theta_n^\prime, \theta_n, X_n)$ (see, for example [Zhao et al., 2021, Linhart et al., 2024]).

§.§ Summary of Results

In this paper we propose a measure of discrepancy between 2 conditional distributions, by adapting the well-known kernel Maximum Mean Discrepancy (MMD) [Gretton et al., 2012] to the conditional setting. Towards this, we first embed the conditional distributions $P_{X|Z}$ and $P_{Y|Z} $ in a reproducing kernel Hilbert space (RKHS), through their conditional kernel mean embeddings [Park and Muandet, 2020]. Then we quantify the discrepancy between $P_{X|Z}$ and $P_{Y|Z} $ in terms of the norm difference (in the RKHS) between the conditional mean embeddings averaged over the marginal distribution of $Z$. We refer to this measure as the Expected Conditional Mean Embedding (ECMMD). This measure characterizes the equality of two conditional distributions, that is, the $\EMMD$ measure is zero if and only if  $P_{X|Z} = P_{Y|Z} $ almost surely (Proposition <ref>). Moreover,  leveraging the reproducing property of the Hilbert space (the well-known `kernel-trick') we can express the $\EMMD$ in terms of the kernel dissimilarities averaged over the respective conditional distributions (Proposition <ref>). This allows us to estimate the $\EMMD$ efficiently based on the observed data. Specifically, we propose a nearest-neighbor graph based estimate of the $\EMMD$ that has the following properties:

* The $\EMMD$ estimate has a simple, interpretable form, which does not require any estimation of density or distribution functions. Moreover, it encompasses both categorial and continuous responses and, consequently, can be easily applied to a range of data types.

* The estimate can be computed in near-linear time (with a fixed number of nearest neighbors) irrespective of the dimension of the data (see Remark <ref>).

* The estimate is consistent for the population $\EMMD$ measure under mild moment conditions on the kernel (Theorem <ref>).

* Under $\bH_0$ as in (<ref>) the estimate is asymptotically Gaussian (Theorem <ref>) and its  variance has a simple tractable form (Proposition <ref>). In particular, we can consistently estimate the variance under $\bH_0$ (Proposition <ref>), using which we obtain a universally consistent test for (<ref>) that requires no nuisance parameter estimation (see Remark <ref> for details). Moreover, both the consistency and the null distribution hold for any fixed kernel bandwidth, which is in contrast to several results on conditional inference, where density estimation with dimension dependent smoothing bandwidths are required for consistent estimation/testing.

Next, motivated by the model-$X$ framework for conditional independence testing [Candes et al., 2018], in Section <ref> we propose a resampling based test for the hypothesis (<ref>) when there is sample access from one of the conditional distributions. This also provides a test for the conditional goodness-of-fit problem that controls Type I error in finite samples and is asymptotically consistent with only a finite number of resamples (see Proposition <ref>). To further improve the stability of the procedure, we propose a de-randomized test which attains precise asymptotic level and is universally consistent (see Section <ref>). Our method outperforms competing methods for conditional goodness-of-fit testing, such as the test based on the Kernel Conditional Stein Discrepancy (KCSD) [Jitkrittum et al., 2020], in simulations, even in simple Gaussian models (see Section <ref>). Finally, we return to the examples mentioned before and illustrate how our proposed method performs both in simulations and real-data. The following is a summary of our findings:

* In Section <ref> and Section <ref> we apply the $\EMMD$ test for assessing calibration in classification and regression models, respectively. We compare our method with the test based on squared kernel calibration error (SKCE) proposed recently in [Widmann et al., 2019, Widmann et al., 2022].   In terms of power, the $\EMMD$ and SKCE are comparable for classification model but with the proposed $\EMMD$ we obtain a significantly more power for regression model. Moreover, under the null hypothesis the SKCE statistic has a non-Gaussian asymptotic distribution (an infinite weighted sum of $\chi^2$ distributions) with no closed form expressions for the quantiles, hence the rejection threshold has to be chosen based on bootstrap/permutation methods. In contrast, the null distribution of the $\EMMD$ statistic is Gaussian, hence, the rejection threshold can be obtained readily.

* We also apply our method to test whether convolution neural network models for classification are calibrated using the CIFAR 10 dataset. The $\EMMD$ test is able to successfully predict the insufficient calibration and detect the significant change in calibration performance when recalibration technique is applied (see Section <ref>).

* For comparing regression curves, we apply the $\EMMD$ method on the wind energy dataset [Ding, 2019, Prakash et al., 2022, Hwangbo et al., 2017], to test if the effect wind speed on wind power generation remain the same across different turbines. In contrast to existing methods for comparing regression curves, which  are primarily focused on comparing the mean functions, the $\EMMD$ based test can be used to detect arbitrary differences between the two models (see Section <ref> for details).

* In Section <ref> we use the $\EMMD$ to test the validity of posterior approximations on benchmark examples arising in simulation-based inference (SBI). Specifically, we apply the $\EMMD$ to test how well methods based on Mixture Density Networks (MDNs) and Neural Spline Flow (NSF) approximate the posterior distributions in the well-known Two Moons and Simple Likelihood Complex Posterior (SLCP) tasks.

* In Section <ref> we apply the $\EMMD$ test to validate an emulator for the conditional density of the synthetic `redshift' associated with photometric galaxy images. Using a Gaussian convolutional mixture density network (ConvMDN) as an emulator for the synthetic `redshift' (as in [Zhao et al., 2021, Dey et al., 2022, LeRoy and Zhao, 2021, Schmidt et al., 2020]), we illustrate the efficacy of the $\EMMD$ test in accurately detecting similarities and differences between the emulator and the simulator.

§.§ Related Work

Although, as discussed above, specific instance of the conditional two-sample problem are abundant, very few rigorous statistical methods have been proposed in the generality considered in this paper.  Very recently, Chen et al., 2022 proposed a test for the hypothesis (<ref>) using a conditional version of the celebrated energy distance test [Székely, 2003, Wang et al., 2015, Székely and Rizzo, 2004, Baringhaus and Franz, 2004]. Their method uses kernel density estimation techniques, which, as mentioned before, require a careful choice of the bandwidth for the asymptotic properties to hold, that
can be difficult to control beyond low dimensions. The $\EMMD$ measure has also appeared in the recent paper by Huang et al., 2022, where it is used as a metric for empirically quantifying the discrepancy between two conditional generative models. However, their estimation method and final aim are very different (see Remark <ref> are further details). In the context of simulation-based inference, Linhart et al., 2024 recently proposed a test for comparing 2 conditional distributions locally at a given observation using the classifier 2-sample test [Friedman, 2003, Lopez-Paz and Oquab, 2022].

In another variant of the conditional two-sample problem one assumes that the responses $X$ and $Y$ are conditioned on (potentially) different covariates and the response covariate pairs are independently generated from their respective joint distributions. In this setting, Hu and Lei, 2023 proposed a test using techniques from conformal prediction and a classifier based estimation of the marginal density ratio of the covariates.
However, this requires sample-splitting and the performance of the test depends on the accuracy of the classifier. A refinement based on de-biased two-sample $U$-statistics has been proposed very recently by  Chen and Lei, 2024. Another test has been proposed by Yan and Zhang, 2022 based on the integrated conditional energy distance. This requires estimation of the marginal density of the covariates based on kernels with smoothing bandwidths. Consequently, their test statistic has a non-negligible bias for dimensions greater than 4 and the asymptotic variance is also intractable, hence, the rejection threshold is chosen using a bootstrap resampling approach.

A similarity measure between the mean functions of 2 conditional distributions
sharing the same set of covariates, based on the minimum mean square error (mMSE) gap, also appears in <cit.>.

§ CONDITIONAL TWO SAMPLE TEST AND EXPECTED CONDITIONAL MMD

In this section we introduce the Expected Conditional MMD (ECMMD) measure for the conditional two-sample sample problem
(<ref>). We begin by introducing the necessary formalism. To this end, suppose $\mathcal X$ and $\mathcal Z$ are Polish spaces, that is, complete and separable metric spaces, and $\mathscr B(\mathcal X)$ and $\mathscr B(\mathcal Z)$ be the $\sigma$-algebras generated by the open sets of $\mathcal X$ and $\mathcal Z$, respectively. Denote by $\mathcal P(\mathcal X \times \mathcal X \times \mathcal Z)$ the collection of all probability distributions on $(\mathcal X \times \mathcal X \times \mathcal Z, \mathcal B(\mathcal X) \times \mathcal B(\mathcal X) \times \mathcal B(\mathcal Z))$. Suppose $P_{XYZ}  \in \mathcal P(\mathcal X \times \mathcal X \times \mathcal Z)$ and $(X,Y,Z)\sim P_{XYZ}$ be a random variable with distribution $P_{XYZ}$. Denote by $P_X, P_Y, P_Z, P_{XY}, P_{YZ}, P_{XZ}$ the marginal distributions of $X, Y, Z, (X, Y), (Y, Z), (X, Z)$, respectively. Also, denote by $P_{X|Z}$, $P_{Y|Z}$, and $P_{XY|Z}$, the regular conditional distributions of $X|Z$, $Y|Z$, and $(X, Y)|Z$, respectively, which exist by <cit.>.
Now, the conditional two-sample sample hypothesis (<ref>) can be stated more formally as follows:

\begin{align}\label{eq:H0H1hypothesis}
\bH_0: P_{Z}\left[P_{X|Z} = P_{Y|Z}\right] = 1 \text{ versus } \bH_1: P_{Z}\left[P_{X|Z} = P_{Y|Z}\right]<1 ,
\end{align}

where $P_{Z}$ is the marginal distribution of $Z$. Our goal is to test the above hypothesis based on i.i.d. samples $\{(X_1,Y_1,Z_1),\ldots,(X_n,Y_n,Z_n)\}$ from the joint distribution $ P_{XYZ}$.

Note that when $Z$ is almost surely a constant, then (<ref>) reduces to the familiar unconditional two-sample testing problem between the distributions $P_X$ and $P_Y$ which are the marginal distributions of $X$ and $Y$, respectively.

To introduce the $\EMMD$ measure we first recall the fundamentals of the kernel MMD from Gretton et al., 2012.

§.§ Kernel Maximum Mean Discrepancy

Denote by $\mathcal P(\mathcal X)$ the collection of all probability measures on $(\mathcal X, \mathscr B(\mathcal X))$. The maximum mean discrepancy (MMD) between two probability measures $P_X, P_Y \in \mathcal P(\mathcal X)$ is defined as,

\begin{align}\label{eq:FPQ}
\MMD\left[\mathcal{F}, P_X, P_Y \right] = \sup_{f\in\mathcal{F}}\left\{ \mathbb{E}_{X \sim P_X } [ f( X ) ] -\mathbb{E}_{Y \sim P_Y } [f( Y ) ] \right\} ,
\end{align}

where $\mathcal{F}$ is the unit ball of a reproducing kernel Hilbert space (RKHS) $\mathcal H$ defined on $\mathcal X$ [Aronszajn, 1950]. Since $\mathcal H$ is an RKHS, by the Riesz representation theorem <cit.> there exists a positive definite kernel $\sfK : \mathcal X \times \mathcal X \rightarrow \mathbb R$ such that the feature map $\psi_x\in\cH$ for all $x\in \cX$ satisfies $\sfK(x,\cdot) = \psi_{x}(\cdot)$ and $\sfK (x,y) = \langle \psi_x,\psi_y\rangle_{\mathcal{H}}$. The notion of feature map can now be extended to define the kernel mean embedding $\mu_{P}$ for any distribution $P \in \cP(\cX)$ as,

\begin{align}\label{eq:Hf}
\langle f,\mu_{P}\rangle_{\mathcal{H}} = \mathbb{E}_{X\sim  P}[f(X)] ,
\end{align}

for all $f\in \mathcal{H}$. By the canonical form of the feature map it follows that

\begin{align}\label{eq:muP}
\mu_{P}(t) = \E_{X \sim  P}  [\sfK(X,t)] ,
\end{align}

for all $t\in \cX$.

Throughout we will assume the following, which ensures that a RKHS is rich enough to distinguish the two distributions from their corresponding mean embeddings.

The kernel $\sfK : \mathcal X \times \mathcal X \rightarrow \mathbb R$ is positive definite and satisfies the following:

$(1)$  $\mathbb{E}_{X \sim P_X}[\sfK (X,X)^{\frac{1}{2}}]<\infty$ and $\mathbb{E}_{Y \sim P_Y}[\sfK (Y,Y)^{\frac{1}{2}}]<\infty$.

$(2)$ The kernel $\sfK$ is characteristic, that is, the mean embedding $\mu: \mathcal{P}(\mathcal{X})\rightarrow\mathcal{H}$ is a one-to-one (injective) function.Moreover, the RKHS $\cH$ generated by $\sfK$ is separable.

Assumption <ref> ensures that $\mu_{P_X}, \mu_{P_Y}\in\mathcal{H}$ and the MMD can then be expressed as the distance in the RKHS between the corresponding kernel mean embeddings (see <cit.>):

\begin{align}\label{eq:KPQ}
\MMD^2\left[\mathcal{F}, P_X, P_Y\right] = \left\|\mu_{ P_X}-\mu_{ P_Y}\right\|_{\mathcal{H}}^2 ,
\end{align}

where $\|\cdot\|_{\mathcal{H}}$ is the norm corresponding to the inner product $\langle\cdot,\cdot\rangle_{\mathcal{H}}$.
This, in particular, implies that $ \mathrm{MMD}^2\left[\mathcal{F}, P_X, P_Y\right]  = 0$ if and only if $ P_X= P_Y$.

§.§ Conditional Kernel Mean Embedding and Expected Conditional MMD

Using the representation of the kernel mean embedding from (<ref>), one can define the conditional kernel mean embeddings of $X|Z$ and $Y|Z$ as follows (see <cit.>):

\begin{align}\label{eq:conditionalmean}
\mu_{ P_{X|Z}}(t) = \E_{P_{X|Z}}\left[\sfK(X,t)\middle|Z\right] \text{ and } \mu_{ P_{Y|Z}}(t) = \E_{P_{Y|Z}}\left[\sfK(Y,t)\middle|Z\right]  ,
\end{align}

for all $t\in \cX$, where $ P_{X|Z}$ and $P_{Y|Z}$ are the regular conditional distributions of $X|Z$ and $Y|Z$, respectively. By Assumption <ref>, the conditional mean embeddings are well-defined and for all $f \in \cH$, similar to (<ref>) (see <cit.>),

$$\langle f,\mu_{P_{X|Z}}\rangle_{\mathcal{H}} = \mathbb{E}_{ P_{X|Z}}[f(X)|Z] \text{ and }  \langle f, \mu_{P_{Y|Z}}\rangle_{\mathcal{H}} = \mathbb{E}_{P_{Y|Z}}[f(Y)|Z] ,$$

almost surely. Also, note that, unlike the (unconditional) kernel mean embeddings which are fixed elements in $\cH$, the conditional kernel mean embedding $\mu_{ P_{X|Z}}$ and $\mu_{ P_{Y|Z}}$ are random variables taking values in $\cH$. Consequently, $\|\mu_{ P_{X|Z}} - \mu_{ P_{Y|Z}}\|_{\cH}$ is random metric that measures the difference between the distributions $ P_{X|Z}$ and $P_{Y|Z}$, given a specific value of $Z$. Averaging over $Z$ leads to the following definition:

The Expected Conditional $\MMD$ ($\EMMD$) between the conditional distributions $P_{X|Z}$ and $P_{Y|Z}$ is defined as:

\begin{align}\label{eq:defEMMD}
\EMMD^2\left[\cF,  P_{X|Z}, P_{Y|Z}\right] := \E_{Z\sim P_Z}\left[\left\|\mu_{ P_{X|Z}} - \mu_{ P_{Y|Z}}\right\|_{\cH}^2\right] ,
\end{align}

where the expectation taken over the marginal distribution $P_Z$.

The following result shows that $\EMMD$ characterizes the equality of two conditional distributions (see Appendix <ref> for the proof).

Suppose the kernel $\sfK$ satisfies Assumption <ref>. Then

\begin{align*}
\EMMD^2\left[\cF,  P_{X|Z}, P_{Y|Z}\right] = 0\text{ if and only if }\P\left[ P_{X|Z} =  P_{Y|Z}\right] = 1.
\end{align*}

The above result shows that the conditional 2-sample hypothesis (<ref>) can be equivalently reformulated as:

$$\bH_0: \EMMD^2\left[\cF,  P_{X|Z}, P_{Y|Z}\right]= 0 \quad \text{ versus } \quad \bH_1: \EMMD^2\left[\cF,  P_{X|Z}, P_{Y|Z}\right] > 0.$$

Hence, to obtain a consistent test for the conditional 2-sample problem it suffices to consistently estimate $\EMMD^2\left[\cF,  P_{X|Z}, P_{Y|Z}\right]$ based on the data. While at first glance this might seem difficult, because the  MMD involves a supremum over functions in the unit ball of a RKHS, the well-known `kernel'-trick [Gretton et al., 2012] allows us to express the $\EMMD$ in the following more tractable form:

Suppose the kernel $\sfK$ satisfies Assumption <ref>. Then,

\begin{align}\label{eq:EMMDK}
\EMMD^2\left[\cF,  P_{X|Z}, P_{Y|Z}\right] = \E\left[\sfK(X,X') + \sfK(Y,Y') - \sfK(X,Y') - \sfK(X',Y)\right] ,
\end{align}

where $(X,X',Y,Y',Z)$ has the following distribution: Generate $Z\sim P_Z$ and then sample $(X,Y)$ and $(X',Y')$ independently from the conditional distribution $ P_{XY|Z}$.

The proof of Proposition <ref> is given in Appendix <ref>. In the next section, we use the representation in (<ref>) to estimate $\EMMD$ based on nearest-neighbors.

As mentioned in the Introduction, the $\EMMD$ measure has been used in  [Huang et al., 2022] to quantify the discrepancy between two conditional generative models (see also earlier related work [Ren et al., 2016] where the conditional kernel mean embeddings was used to measure the pointwise difference between 2 conditional distributions). In [Huang et al., 2022] authors use multiple rounds of resamples from the generative models to estimate the $\EMMD$ and demonstrate numerically the performance of the estimate in measuring discrepancies between conditional distributions. However, the distributional properties of the estimate was not explored. On the other hand, we use the $\EMMD$ measure to develop a test for the conditional 2-sample problem with asymptotic guarantees. Our estimation is based on a simple nearest-neighbor technique that requires no resampling. We establish the consistency and asymptotic normality of our estimate and, as a consequence, obtain a universally consistent test for the conditional 2-sample problem (see Section <ref>).

§ ESTIMATING $\EMMD$ USING NEAREST-NEIGHBORS

In this section we propose an estimate of the $\EMMD$ measure (<ref>) given i.i.d. samples $(X_1,Y_1,Z_1),\ldots,$
$(X_n,Y_n,Z_n)$ from the joint distribution $ P_{XYZ}$. Towards this, suppose $\sfK$ is a kernel satisfying Assumption <ref> and for $\bw = (x,y), \bw' = (x',y')\in \cX\times\cX$ define,

\begin{align}\label{eq:defH}
\sfH(\bw,\bw') = \sfK(x,x') + \sfK(y,y') - \sfK(x,y') - \sfK(x',y) .
\end{align}

Then using Proposition <ref> the $\EMMD$ can be expressed as:

\begin{align}\label{eq:HWW}
\EMMD^2\left[\cF,  P_{X|Z}, P_{Y|Z}\right] = \E\left[\sfH(\bW,\bW')\right] = \E\left[\E\left[\sfH(\bW,\bW')\middle|Z\right]\right] ,
\end{align}

where $Z\sim  P_Z$ and $\bW = (X,Y)$ and $\bW' = (X',Y')$ are generated independently from $ P_{XY|Z}$. To estimate the RHS of (<ref>), we first fix $Z=Z_u$, for $1 \leq u \leq n$, and consider the condition expectation $\E\left[\sfH(\bW,\bW')\middle|Z= Z_u\right]$.  To estimate this conditional expectation the idea is to average the discrepancies of the centered kernel $\sfH$ over indices which are `close'  to $Z_u$. One natural way to capture such proximity is through neighbors-neighbor graphs. This motivates the following construction:

* Fix $K=K_n \geq 1$ and construct the directed $K$-nearest neighbor ($K$-NN) graph $G(\sZ_n)$ of the data points $\sZ_n := \{Z_1,Z_2,\ldots,Z_n\}$. Unless otherwise specified, we will use the notation $(u, v)$ to denote the directed edge $Z_u\ra Z_v$, for $Z_u, Z_v \in \sZ_n$.
Also, we will denote the edge set of $G(\sZ_n)$ by
$$E(G(\sZ_n)) = \{(u, v) \in  [n]^2: \text{such that } Z_u\ra Z_v \text{ is a directed  edge in } G(\sZ_n) \},$$

where $[n]  := \{1, 2, \ldots, n\}$.

* Then the $K$-NN based estimate of (<ref>) is given by:

\begin{align}\label{eq:estEMMD}
\EMMD^2[\sfK, \bcW_n, \sZ_n] = \frac{1}{n}\sum_{u=1}^{n}\frac{1}{K}\sum_{v \in N_{G(\sZ_n)}(u)}\sfH(\bW_{u},\bW_{v}) ,
\end{align}

where $\bcW_n:=\left\{\bW_i = (X_i,Y_i):1\leq i\leq n\right\}$ and $N_{G(\sZ_n)}(u) = \{v:
(u,v) \in E(G(\sZ_n))\}$.

Note that the estimate (<ref>) can be computed easily in $O(K n \log n)$ time, in any dimensions. This is because $K$-NN graph can be computed in $O(K n \log n)$ time (see, for example, [Friedman et al., 1977]) and, given the graph, the sum in (<ref>) can be computed in $O(K n)$ time, since the $K$-NN graph has $O(Kn)$ edges.

To establish the consistency of (<ref>) we assume the following on the conditioning variable $Z$.

The random variable $Z$ takes values in $\cZ = \R^d$, for some $d \in \N$, and $\|Z - Z'\|_2$ has a continuous distribution, where $Z, Z'$ are i.i.d. samples from $P_Z$.

This assumption ensures that the $K$-NN graph constructed using $\|\cdot\|_2$ norm is well-defined and the degrees of its vertices scales proportional to $K$ (see [Jaffe et al., 2020, Deb et al., 2020]). One can easily relax this to include any finite dimensional inner product space over $\R$, by isometrically isometric embedding such spaces into the Euclidean space. It is worth noting here that we do not require the space $\cX$ (on which the random variables $X, Y$ are defined) to be Euclidean or even finite-dimensional. (For example, while calibration testing in classification models the variables $X, Y$ are categorical.) The following theorem establishes the consistency of

Suppose Assumption <ref> and Assumption <ref> hold. Moreover, suppose the kernel $\sfK$ satisfies $\int \sfK(x,x)^{2+\delta}\mathrm{d} P_X(x)<\infty$, and  $\int \sfK(x,x)^{2+\delta}\mathrm{d} P_Y(x)<\infty$, for some $\delta>0$. Then with $K = o(n/\log n)$,

\begin{align*}
\EMMD^2[\sfK,\bcW_n,\sZ_n]\pto \EMMD^2\left[\cF,  P_{X|Z}, P_{Y|Z}\right] ,
\end{align*}

where $\EMMD^2[\sfK,\bcW_n,\sZ_n]$ and $\EMMD^2\left[\cF,  P_{X|Z}, P_{Y|Z}\right]$ are defined in (<ref>) and (<ref>),  respectively.

The proof of Theorem <ref> is given in Appendix <ref>. The proof involves the following steps:

* First we show that the expectation of $\EMMD^2[\sfK,\bcW_n,\sZ_n]$ converges to the population $\EMMD^2$ (Lemma <ref>). The main idea here is that averaging the centered kernel $\sfH$ (recall (<ref>)) around the nearest neighbors of a point $Z=z$, provides an asymptotically unbiased estimate of the $\MMD^2$ distance between conditional mean embeddings $P_{X|Z=z}$ and $P_{Y|Z=z}$.

* Next, using the Efron-Stein inequality [Efron and Stein, 1981], we show that the variance of (<ref>) converges to zero (Lemma <ref>). This leverages the local dependence of the $K$-NN graph, which is controlled by the condition $K = o(n/\log n)$.

§ ASYMPTOTIC TEST BASED ON $\EMMD$

To use the estimate (<ref>) for testing the conditional two-sample hypothesis, we need to derive its asymptotic distribution under $\bH_{0}: P_{X|Z} = P_{Y|Z}$. To this end, consider the scaled version of the $\emmdstat$ statistic,

\begin{align}\label{eq:defetan}
\eta_{n}:=\sqrt{nK}\emmdstatsq = \frac{1}{\sqrt{nK}}\sum_{u=1}^{n}\sum_{v \in N_{G(\sZ_n)}(u)} \sfH(\bW_{u},\bW_{v}).
\end{align}

Also, denote by $\cF(\sZ_n)$ the $\sigma$-algebra generated by $Z_1,Z_2,\ldots,Z_n$.  To begin with, note that

\begin{align}\label{eq:EH0}
\E_{\bH_{0}}\left[\eta_{n}|\cF(\sZ_n)\right] = \frac{1}{\sqrt{nK}}\sum_{1 \leq u,v \leq n} \E_{\bH_{0}}\left[\sfH(\bW_{u},\bW_{v})|\cF(\sZ_n)\right]\one\left\{(u,v)\in E(G(\sZ_n))\right\} = 0 ,
\end{align}

since, for $1 \leq u\neq v \leq n$, recalling (<ref>),

\begin{align}
\E_{\bH_{0}}\left[\sfH(\bW_{u},\bW_{v})|\cF_{n}(\sZ_n)\right]
= 0 , \label{eq:condexpij0}
\end{align}

almost surely.
Next, we compute the conditional variance of $\eta_{n}$ under $\bH_{0}$.

Denote by $\sigma^2_n := \Var_{\bH_0}\left[\eta_n\middle|\cF(\sZ_n)\right]$, conditional variance of $\eta_n$ given $\cF(\sZ_n)$  under  $\bH_{0}$. Then

\begin{align}\label{eq:condeta2exp}
\sigma^2_n = \frac{1}{nK}\sum_{ 1 \leq u, v \leq n} f(Z_u,Z_v) \left( \one\left\{(u,v)\in E(G(\sZ_n))\right\} + \one\left\{(u,v),(v,u)\in E(G(\sZ_n))\right\} \right) ,
\end{align}

where $f(Z,Z'):= \E_{\bH_0}\left[\sfH^2((X,Y),(X',Y')) |Z,Z'\right]$ with $(X,Y,Z),(X',Y',Z')$ i.i.d. samples from $ P_{XYZ}$.

The proof of Proposition <ref> is given in Appendix <ref>. The proof relies on the observation that for $1 \leq u \ne u'  \ne v \ne v' \leq n$,

$$\mathrm{Cov}_{\bH_0}[ \sfH(\bW_{u},\bW_{v}),  \sfH(\bW_{u'},\bW_{v'}) | \cF(\sZ_n) ] = 0, $$
unless $\{u, v\} = \{ u', v'\}$.
In other words, under the null hypothesis, the summands in (<ref>) are pairwise conditionally uncorrelated, hence, the only terms that contribute to the conditional variance are when (1) $u=u'$ and $v = v'$ which corresponds to the first term in (<ref>) or (2) $u=v'$ and $v = u'$ which corresponds to the second term in (<ref>).  Consequently, $\sigma^2_n$ has a simple form as the average of the conditional expectation of $\sfH^2$ over the edges of the $K$-NN graph $G(\sZ_n)$ (counted twice when edges are present in both directions).
This representation is particularly convenient because just by replacing the (unknown) conditional expectation $f(Z_u, Z_v)$ with kernel values $\sfH^2(\bW_{u},\bW_{v})$, we get following natural estimate of $\sigma_{n}^2$:

\begin{align}\label{eq:defSn}
\hat{\sigma}_{n}^2:= \frac{1}{nK}\sum_{u=1}^{n}\sum_{v=1}^{n}\sfH^{2}\left(\bW_{u},\bW_{v}\right)\left(\one\left\{(u,v)\in E(G(\sZ_n))\right\}+\one\left\{(u,v),(v,u)\in E(G(\sZ_n))\right\}\right).
\end{align}

In the following we establish the consistency of $\hat{\sigma}_{n}^2$ (see Appendix <ref> for the proof) and subsequently use $\hat{\sigma}_n$ to define a studentized version of $\eta_n$ for constructing the asymptotic test.

Suppose Assumption <ref> and Assumption <ref> holds. Furthermore, for $(X, Y, Z)\sim P_{XYZ}$ assume that $P_{XY}(X\neq Y)>0$ and the kernel $\sfK$ satisfies $\int\sfK(x,x)^{4+\delta}\rmd P_{X}(x)<\infty$ and $\int\sfK(y,y)^{4+\delta}\rmd P_{Y}(y)<\infty$, for some $\delta>0$. Then under $\bH_0$ with $K = o(n/\log n)$, as $n \rightarrow \infty$

\begin{align}\label{eq:Snapproxvar}
\left|\frac{\hat{\sigma}_n^2}{\sigma_n^2}-1\right| = o_{p}\left(n^{-\frac{\delta}{32+4\delta}}\right).
\end{align}

where $\sigma_n^2$ and $\hat{\sigma}_n^2$ are defined in (<ref>) and (<ref>), respectively.

We are now ready to state the result about the asymptotic null distribution of $\eta_{n}$. Specifically, in the following theorem we show that $\eta_n$ scaled by $\hat{\sigma}_n$ converges to $N(0, 1)$, under $\bH_0$, in the Kolmogorov distance.

Suppose Assumption <ref> and Assumption <ref> holds. Furthermore, for $(X, Y, Z)\sim P_{XYZ}$ assume that $P_{XY}(X\neq Y)>0$ and the kernel $\sfK$ satisfies $\int\sfK(x,x)^{4+\delta}\rmd P_{X}(x)<\infty$ and $\int\sfK(y,y)^{4+\delta}\rmd P_{Y}(y)<\infty$, for some $\delta>0$. Then under $\bH_0$, with $K = o(n^{1/44})$, as $n \rightarrow \infty$,

\begin{align}\label{eq:distributionH0}
\sup_{z\in\R}\left|\P_{\bH_0}\left[\frac{\eta_{n}}{\hat{\sigma}_n}\leq z\right] - \Phi(z)\right|\ra 0 ,
\end{align}

where $\Phi(\cdot)$ is the CDF of the standard Gaussian distribution.

The proof of Theorem <ref> is given in Appendix <ref>. The proof entails showing that

\begin{align}\label{eq:distributionH0pf}
\sup_{z\in\R} \left|\P_{\bH_0}\left[\frac{\eta_{n}}{\sqrt{\Var_{\bH_0}[\eta_n|\cF(\sZ_n)]}}\leq z\right] - \Phi(z)\right|\ra 0 .
\end{align}

The result in (<ref>) follows by combining the above with Proposition <ref>. To show (<ref>) we use the Stein's method based on dependency graphs [Chen and Shao, 2004], which allows us to control the
Kolmogorov distance between $\eta_{n}/\hat{\sigma}_n$ and $N(0, 1)$ in terms
of maximum degree of the $K$-NN graph.

To choose the rejection threshold for $\hat \eta_n$ based on Theorem <ref>, fix $\alpha \in (0, 1)$ and consider the test function:

\begin{align}\label{eq:defphin}
\phi_{n} := \bm{1}\left\{\left|\eta_n/\hat{\sigma}_n\right|>z_{\alpha/2}\right\} ,
\end{align}

where $\eta_n$ and $\hat{\sigma}_n$ are as defined in (<ref>) and (<ref>), respectively. Theorem <ref> directly implies that $\phi_{n}$ is asymptotically level $\alpha$ for $\bH_0$ as in (<ref>), that is,

$$\lim_{n\ra\infty} \P_{\bH_0}\left[\phi_{n} = 1\right] = \alpha.$$

Since $\EMMD$ characterizes the equality of the conditional distributions and (<ref>) is a consistent estimate of $\EMMD$ (recall Theorem <ref>), the test $\phi_n$ is consistent for fixed alternatives. This is summarized in the following result (see Appendix <ref> for the proof).

Suppose the assumptions of Theorem <ref> hold. Then for any $P_{XYZ} \in \bH_1$,
$$\lim_{n\ra\infty} \P_{\bH_1}\left[\phi_n = 1\right] = 1.$$

The ECMMD test has several interesting features which are in contrast to other kernel and nearest-neighbor based nonparametric tests:

* The limiting distribution of the statistic $\eta_{n}$ is normal under $\bH_0$ and the asymptotic variance has a simple form, which can be easily estimated from the data. Consequently, the rescaled test statistic $\eta_{n}/\hat{\sigma}_n$ converges to $N(0, 1)$ under $\bH_0$ and we can obtain the rejection threshold as in (<ref>), without having to estimate any nuisance parameter. In other words, $\eta_{n}/\hat{\sigma}_n$ is asymptotically distribution-free, that is, its limiting distribution under  $\bH_0$ does not depend on the unknown distribution of the data. In contrast, the familiar kernel MMD statistic for the (unconditional) two-sample problem has a non-Gaussian (specifically, an infinite mixture of chi-squares) limiting distribution under the null [Gretton et al., 2012]. Closed form estimates for the quantiles of such distributions are not available, in general,  which necessitates the use permutation/bootstrap resampling techniques or conservative approximations based on concentration inequalities for determining the rejection thresholds [Schrab et al., 2023, Gretton et al., 2009, Chatterjee and Bhattacharya, 2023]. We circumvent this issue through the use of nearest neighbor graphs (on the space of the covariate variable), which mitigates the dependence among the summands in (<ref>) in such a way that, although the kernel $\sfH$ is degenerate under $\bH_0$, the asymptotic distribution of $\eta_n$ is normal.

* Another important property of the $\EMMD$ test is that the statistic $\eta_n$ is unbiased (has mean zero) under $\bH_0$ (recall (<ref>)). This is different from the recent work on conditional independence testing based on nearest neighbors [Azadkia and Chatterjee, 2021, Shi et al., 2021, Huang et al., 2022], where the test statistics have a non-zero bias under the null, and, hence, cannot be directly used for inference, without additional de-biasing. The bias issue also appears in density-estimation based methods, both for the conditional independence [Wang et al., 2015] and the conditional 2-sample problem [Yan and Zhang, 2022], which is usually handled by choosing dimension dependent smoothing bandwidths. We, on the other hand, are able to cancel the bias because of the paired nature of the samples (recall that for each $Z_i$ we have paired samples $(X_i, Y_i)$) and through the use of nearest-neighbors.

§ A RESAMPLING BASED CONDITIONAL GOODNESS-OF-FIT TEST USING $\EMMD$

In this section, drawing parallel from the model-$X$ framework for conditional independence testing [Candes et al., 2018], we design a resampling based test the hypothesis  (<ref>) that controls Type I error in finite samples when it is possible to efficiently sample from one of the conditional distributions $P_{X|Z}$ or $P_{Y|Z}$. This principle applies more broadly to the conditional goodness-of-fit problem which entails testing
the hypothesis in (<ref>), when one of the conditional distributions is specified. Specifically, suppose we are given on i.i.d. samples $(Y_1, Z_1),\ldots,(Y_n, Z_n)$ from the joint distribution $ P_{YZ}$ and we wish to test the hypothesis

\begin{align}\label{eq:H0H1finite}
\bH_0: P_{Z}\left[P_{Y|Z}=P_{X|Z}\right] = 1 \quad \text{ versus } \quad \bH_1: P_{Z}\left[P_{Y|Z} \ne P_{X|Z}\right]>0,
\end{align}

where $P_{X|Z}$ is a specified conditional distribution.

A finite sample conditional goodness-of-fit test

* For each $1\leq u \leq n$, generate i.i.d. samples $\left(X_{u}^{(1)},\ldots,X_{u}^{(M+1)}\right)$ from the distribution $ P_{X|Z=Z_u}$.

* Denote by $\bm W_{u}^{(m)}= (X_u^{(M+1)}, X_u^{(m)})$, for $1 \leq m \leq M$, and $\bm W_{u}^{(M+1)}= (X_u^{(M+1)}, Y_u)$, for $1\leq u \leq n$. Define

\begin{align}\label{eq:defetanm}
\eta_{n}^{(m)}=\frac{1}{\sqrt{nK}}\sum_{u=1}^{n}\sum_{v \in N_{G(\sZ_n)}(u)}\sfH(\bW_{u}^{(m)},\bW_{v}^{(m)}),
\end{align}

for $1 \leq m \leq M+1$, where $\sfH$ is defined in (<ref>) and $G(\sZ_n)$ is the $K$-NN graph of the data points $\sZ_n = \{Z_1,\ldots,Z_n\}$.

* Report the $p$-value

\begin{align}\label{eq:Mfinite}
p_{M}:=\frac{1}{M+1}\left[1+\sum_{m=1}^M \mathbf{1}\left\{\eta_{n}^{(m)} \geq \eta_{n}^{(M+1) }\right\}\right] .
\end{align}

In Algorithm <ref> we develop a resampling based test for (<ref>). Note that (<ref>) in Algorithm <ref> is a valid $p$-value because the collection $\{\eta_{n}^{(1)},\eta_{n}^{(2)},\ldots,\eta_{n}^{(M+1)}\}$ is exchangeable conditional on $\sZ_n$ when $ P_{X|Z} =  P_{Y|Z}$ almost surely.

Consequently, the resulting test controls Type I error in finite samples.
This is formalized in the following result which also establishes the asymptotic consistency of test with a finite number of resamples (see Appendix <ref> for the proof):

Fix $\alpha \in (0, 1)$  and consider the test function $\tilde \phi_{n, M} = \bm{1}\left\{p_{M}\leq \alpha\right\}$, with $p_M$ as in (<ref>). Then the following hold:

$(1)$ $\P_{\bH_0}[\tilde \phi_{n, M} = 1]\leq \alpha$.

$(2)$ For any $P_{XYZ} \in \bH_1$ (that is, $P_{Y|Z}\ne P_{X|Z}$), $\lim_{n \rightarrow \infty} \P_{\bH_1}[\tilde \phi_{n, M} = 1] = 1$, whenever $M>\frac{1}{\alpha}-1$.

In the following we summarize the current state-of-the-art in nonparametric conditional goodness-of-fit testing (Section <ref>), discuss how the resampling based test in Algorithm <ref> fits into this literature and its relevance in modern machine learning problems (Section <ref>), and propose a de-randomized version of the resampling test and study its asymptotic properties (Section <ref>).

§.§ Prior Work on Conditional Goodness-of-Fit Testing

The conditional goodness-of-fit problem has its roots in the econometrics literature, beginning with the work of Andrews [Andrews, 1997], which extended the classic Kolmogorov test to the conditional case. Thereafter, other methods for the conditional goodness-of-fit problem have been proposed, however, these tests either involve density estimation [Zheng, 2000], which require decaying smoothing bandwidths that can be difficult to control, or are designed for specific families of conditional models, such as structural equation models [Moreira, 2003] or generalized linear models [Stute and Zhu, 2002]. Recently, Jitkrittum et al., 2020 proposed a general nonparametric test for the conditional goodness-of-fit problem that does not require any density estimation or knowledge of the normalizing constant of the conditioning distribution. Specifically, the method extends the well-known Kernel Stein Discrepancy (KSD) [Liu et al., 2016] to the conditional setting, referred to as the Kernel Conditional Stein Discrepancy (KCSD), and only requires knowledge of the score-function of the conditional distribution $P_{X|Z}$. While the KCSD method circumvents several of the limitations of previous density estimation based methods, it still remains inapplicable in situations where $P_{X|Z}$ is implicitly defined and one only has access to samples from the conditional distribution. In contrast, the resampling based method described above can be readily applied in such situations. We elaborate on this in the next section.

§.§ When is Resampling Useful?

Any goodness-of-fit problem can be transformed into a 2-sample problem by repeatedly sampling from the known null distribution. Therefore, it is no surprise that the $\EMMD$ statistic, which is a measure of difference between two conditional distributions, can be calibrated for the conditional goodness-of-fit problem through resampling. Friedman <cit.> summarizes this principle succinctly as follows: `Using this additional
information has the potential for increased power at the expense of having to generate many Monte Carlo samples, instead of just one.' We are not advocating that one should always resort to resampling for conditional goodness-of-fit testing, but there are important cases where it might be reasonable to apply the finite sample test in Proposition <ref>, over  existing conditional goodness-of-fit methods. This is indeed the case for two of the examples considered in this paper: (1) calibration testing and (2) validation of emulator models in SBI, as explained below:

(1) Calibration testing for classification is an example of a conditional goodness-of-fit problem that can be implemented as a conditional 2-sample problem by sampling from the known distribution $P_{X|f(Z)}$  (recall (<ref>)). We illustrate the numerical performance of  this method in Section <ref>. We also compare our method with the test in [Widmann et al., 2019], which is based on directly estimating kernel versions of the expected calibration error (ECE) (see Appendix <ref> for the definition of ECE). Beyond the classification setting, that is, for continuous response (see Section <ref> for more details), kernel based ECE measures for testing calibration require taking expectations against a specified (usually intractable) probabilistic model and involve kernels on the space of distributions [Widmann et al., 2022]. Although these issues have been mitigated in the recent work based on KCSD [Glaser et al., 2023], it still requires the score function of the generative model to be available in closed form. On the other hand, our formulation of the calibration problem as a conditional 2-sample hypothesis only requires sample access from the conditional distribution $P_{X|f(Z)}$, which is usually readily available given the trained model $f(Z)$. Consequently, our proposed method can be applied easily to categorical and continuous responses. In fact, from the numerical experiments in Section <ref> we will see that even the asymptotic test, which only requires sampling a single $X_i$ for each observation $(Y_i, f(Z_i))$, is powerful in a variety of examples.

(2) Validation tests of emulator models or approximate posteriors is another instance of a conditional goodness-of-fit problem that can be naturally operationalized as a conditional 2-sample problem. Here, $P_{X|Z}$ in (<ref>) corresponds to $p(\cdot| \theta)$ (the density/score function of simulator) or $p(\cdot|X)$ (the true posterior distribution), and
$P_{Y|Z}$ corresponds to either the emulator distribution $q(\cdot| \theta)$ as in (<ref>) or the approximate posterior $q(\cdot|X)$ as in (<ref>), respectively. In either case, there is usually no tractable form $P_{X|Z}$. Moreover, it is usually easier to sample from $q(\cdot| \theta)$ or $q(\cdot|X)$ than to deal with their actual functional forms. In this situation, following Friedman's aphorism, we can validate the performance of the emulator (test the hypothesis (<ref>)) or the approximate posterior (test the hypothesis (<ref>)) using the finite-sample test (<ref>), by repeatedly sampling from the emulator, or the asymptotic test (<ref>), with a single set of samples (recall the discussion after (<ref>) and (<ref>)). We discuss examples in Section <ref>.

§.§ A Derandomized  Asymptotic Test Based on ECMMD

One issue with the test in Algorithm <ref> is that it is a randomized procedure, that is, different runs of the algorithm produce different $p$-values, which can lead to inconsistent conclusions. In this section we propose a de-randomization method that, instead of calculating the $\EMMD$  statistic for each run of the algorithm as in (<ref>), computes a single test statistic by averaging the kernel discrepancies over $M = M_n$ resamples for each given $Z= Z_u$, for $1 \leq u \leq n$, as described below:

(1) For each $1\leq u\leq n$ generate i.i.d. samples $(X_{u}^{(1)},\ldots, X_{u}^{(M_n)})$ from the distribution $P_{X|Z = Z_u}$ independent of $Y_u$ given $Z_u$.

(2) Denoting $\bW_{u}^{(m)} = (X_u^{(m)}, Y_{u})$, for $1\leq m\leq M_n$ and $1\leq u\leq n$ define the de-randomized test statistic as follows:

\begin{align}\label{eq:defDn}
D_n = \frac{1}{nK}\sum_{u=1}^{n}\sum_{v\in N_{G(\sZ_n)}(u)}\frac{1}{M_n}\sum_{m=1}^{M_n}\sfH(\bW_{u}^{(m)},\bW_{v}^{(m)}) .
\end{align}

Note that for $M_n=1$ the statistic $D_n$ equals the estimate $\emmdstatsq$ defined in (<ref>). The averaging step in (<ref>) is meant to mitigate the sensitivity to the resampling uncertainty. To choose the rejection threshold for $D_n$, we now investigate its asymptotic properties. First, we show that $D_n$ consistently estimates the population $\EMMD$, as $n \rightarrow \infty$, irrespective of the choice of $M_n$ (see Appendix <ref> for the proof):

Suppose Assumption <ref> and Assumption <ref> hold. Moreover, suppose the kernel $\sfK$ satisfies $\int \sfK(x,x)^{2+\delta}\mathrm{d} P_X(x)<\infty$, and  $\int \sfK(x,x)^{2+\delta}\mathrm{d} P_Y(x)<\infty$, for some $\delta>0$. Then with $K = o(n/\log n)$,

\begin{align*}
D_n\pto \EMMD^2\left[\cF,P_{X|Z},P_{Y|Z}\right],
\end{align*}

where $\EMMD^2\left[\cF,P_{X|Z},P_{Y|Z}\right]$ is defined in  (<ref>).

Next, we show that a studentized version of $D_n$ converges to $N(0, 1)$ under the null hypothesis. Towards this, define the studentization factor:

\begin{align}\label{eq:estimateM}
\hat{\tau}_n^2 := \frac{1}{nK}\sum_{\substack{1\leq u\leq n\\v\in N_{G(\sZ_n)}(u)}}\left(\frac{1}{M_n}\sum_{m=1}^{M_n}\sfH\left(\bW_{u}^{(m)},\bW_{v}^{(m)}\right)\right)^2 \left(\one\{\cE_{u, v}\} +  \one\{\cE_{u, v}^+\} \right)  ,
\end{align}

where $\cE_{u, v} := \{(u,v)\in E(G(\sZ_n)) \}$ and  $\cE_{u, v}^+ := \{(u,v),(v,u)\in E(G(\sZ_n)) \}$.

Suppose Assumption <ref> and Assumption <ref> hold. Also, assume that the kernel $\sfK$ satisfies $\int\sfK(x,x)^{4+\delta}\rmd P_{X}(x)<\infty$ and $\int\sfK(y,y)^{4+\delta}\rmd P_{Y}(y)<\infty$, for some $\delta>0$. Then, for $K = o(n^{1/44})$ and $M_n\ra\infty$, as $n\ra\infty$, the following holds under $\bH_0$,

\begin{align}\label{eq:distributionH0Dn}
\sup_{z\in\R}\left|\P_{\bH_0}\left[\frac{\sqrt{nK}D_{n}}{\hat{\tau}_n}\leq z\right] - \Phi(z)\right|\ra 0 ,
\end{align}

where $\Phi(\cdot)$ is the CDF of standard Gaussian distribution.

The proof of Theorem <ref> is given in Appendix <ref>. Theorem <ref> shows that the test function

\begin{align}\label{eq:defphinM}
\tilde{\phi}_{n}^* := \bm{1}\left\{\left| \sqrt{nK}D_{n} \right|>z_{\alpha/2} \hat{\tau}_n \right\} ,
\end{align}

is asymptotically level $\alpha$, that is, $\lim_{n \rightarrow \infty} \P_{\bH_0}[\tilde{\phi}_{n}^* =1] = \alpha$. Moreover, since $D_n$ is a consistent estimate of $\EMMD$ (recall Theorem <ref>), then following the proof of Corollary <ref>, the test $\tilde{\phi}_{n}^*$ is also consistent for fixed alternatives, that is, $\lim_{n \rightarrow \infty} \P_{\bH_1}[\tilde{\phi}_{n}^* =1] = 1$. This shows that $\tilde{\phi}_{n}^*$ is a derandomized test for the conditional goodness-of-fit problem (<ref>) which attains precise asymptotic level and is universally consistent.

§ APPLICATIONS

In this section we apply the $\EMMD$ test to the examples discussed in the Introduction. The section is organized as follows: In Section <ref> we use the $\EMMD$ method for testing calibration in both classification and regression models.  We investigate the finite-sample performance of the resampling based test in Section <ref>. We  compare regression functions in the wind energy dataset in Section <ref>. In Section <ref> we apply the ECMMD measure to validate emulators in benchmark SBI examples and for simulating redshifts associated with galaxy images.

§.§ Calibration Tests

In Section <ref> we discuss $\EMMD$ based calibration tests for classification and in Section <ref> we discuss $\EMMD$ based calibration tests for regression. Calibration of convolutional neural network models is tested on the CIFAR-10 dataset using the $\EMMD$ measure in Section <ref>

§.§.§ Calibration Tests for Classification

Calibration tests for classification: (a) Type I error for $n=75$, (b) empirical power for $n=75$, (c) Type I error for $n=100$, and (d) empirical power for $n=100$, as a function of $\rho$.

For testing calibration in classification we consider the following data
generating mechanism (as in [Widmann et al., 2019]),

\begin{align*}
\bm f(Z) = (f_1(Z), 1-f_1(Z))^\top \sim \textnormal{Dir}(\rho, 1-\rho) \text{ and } X\sim \mathrm{Ber}(f_1(Z)) .
\end{align*}

To examine the Type-I error rate and the power we consider the following setups:

* Null hypothesis:  $ Y\sim \textnormal{Bern}(f_1(Z))$;

* Alternative hypothesis: $ Y\sim \textnormal{Bern}(f_1( Z) - f_1(Z)^2)$.

We implement the $\EMMD$ test with the linear kernel, that is, $\sfK(x,y)=x\cdot y$ by varying $\rho\in\{0.3, 0.4, 0.5, 0.6, 0.7\}$, the number of nearest-neighbor $K \in \{5, 10, 15, 20, 25\}$, and the sample size $n \in \{75, 100\}$. Note that, under the alternative hypothesis, as $\rho$ (the signal strength) grows, $f_1^2(Z)$ will tend to be larger, which makes distribution of $Y|Z$ further apart from that of $X|Z$. For comparison we also implement the test based on squared kernel calibration error (SKCE) [Widmann et al., 2019]. The asymptotic null distribution of the test based on SKCE is an infinite weighted sum of $\chi^2$ distributions [Widmann et al., 2019], whose quantiles are intractable.  Hence, to chose the cut-off under the null hypothesis we use a parametric bootstrap procedure as in [Widmann et al., 2019]. Throughout, the nominal level is to be $0.05$. The results are shown in Figure <ref>. The following observations are revealed:

* The $\EMMD$ test is well-calibrated under the null hypothesis for all choices of $K$. The SKCE based test also produces reasonably well-calibrated results (see Figure <ref>(a) and <ref>(c)).

* Overall the $\EMMD$ test has comparable power with the SKCE test, even though SKCE slightly outperforms the $\EMMD$ test when the sample size $n=75$ (see Figure <ref>(b) and <ref>(d)).

This shows that the ECMMD which can be readily implemented without any bootstrap resampling, hence is computationally much more efficient that the SKCE, has comparable statistical performance with the SKCE. In particular, the results in Figure <ref> shows that the computational advantage of the ECMMD comes with almost no sacrifice in statistical performance in comparison to the SKCE.

§.§.§ Calibration Tests for Regression

In the regression framework calibration is often defined in terms of the quantiles of response distribution [Liu et al., 2023]. Specifically, suppose $Q(\rho, z)$ is a pre-trained quantile prediction model, which gives the prediction of the $\rho$-th quantile of the conditional distribution of $Y$ given $Z = z$. The quantile prediction model $Q$ is said to be calibrated if and only if,

\begin{align}\label{eq:regcalibration}
\P\left(Y\leq Q(\rho,Z)|Z = z\right) = \rho ,
\end{align}

for all $0\leq \rho\leq 1$ and $z\in \text{Supp}( P_{Z})$. Notice that for $U\sim \text{Unif}[0,1]$ independent of all previous data, $Q(\cdot, Z)$ is a quantile function for $Q(U,Z)|Z$ and hence, from (<ref>), $Q$ is calibrated if and only if,

\begin{align*}
P_{Y|Z} =  P_{Q(U,Z)|Z}\text{ almost surely } P_Z,
\end{align*}

reducing the test for calibration of $Q$ to the form (<ref>). For more general predictive models, (<ref>) can be further extended as follows (see [Widmann et al., 2022]):  A model $F_Z$ of the conditional distribution $ P_{Y|Z}$ is called calibrated if and only if,

\begin{align}\label{eq:calibration}
P_{Y|F_Z} = F_Z\text{ almost surely } P_{F_Z}.
\end{align}

Note that when $F_Z$ is a classification model, (<ref>) matches with (<ref>). Similar to the classification setting, testing for calibration of a general predictive model $F_Z$ can be framed in terms of  (<ref>) as follows: For $X\sim F_Z$ notice that $ P_{X|F_Z} = F_Z$ almost surely $ P_{F_Z}$. Then by (<ref>), the predictive model $F_Z$ is calibrated if and only if,

\begin{align*}
P_{Y|F_Z} =  P_{X|F_Z}\text{ almost surely } P_{F_Z}.
\end{align*}

Thus, testing for calibration of $F_Z$ is now equivalent to the hypothesis test in (<ref>) with samples $(X_1,Y_1, F_{Z_1}),\ldots, (X_n,Y_n, F_{Z_n})$, where $X_i\sim F_{Z_i}$, for all $1\leq i\leq n$.

Empirical power of calibration tests for the regression model (<ref>) as a function of $\rho$, with (a) $n=50$ and (b) $n=75$.

A special case of interest is when $F_Z$ is a Gaussian linear model. In this case, denote the conditional mean and conditional variance as $\E[X|F_Z]$ and $\mathrm{Var}[X|F_Z]$, respectively. Then testing $X|F_Z\overset{D}{=}Y|F_Z$ is the same as:

since in the mean and the variance determines a Gaussian distribution. In particular, if one uses the homoscedastic linear model, then it suffices to condition only on the conditional mean.  We will use such a model in the following simulation.

To evaluate the performance of the $\EMMD$ test for regression calibration we consider the following model inspired by Widmann et al., 2022:

\begin{align}\label{eq:YZ}
Y=\rho\times\sin(\pi Z)+|1+Z| \varepsilon, \text{ where } \varepsilon\sim N(0,0.15^2) \text{ and } Z\sim \text{Unif}([-1,1]) ,
\end{align}

for $\rho > 0$. Considering a training set of size $n_{\mathrm{train}} = 200$ we fit an ordinary least squares (OLS) to $\{(Y_i, Z_i)\}_{1 \leq i \leq n_{\mathrm{train}}}$ generated i.i.d from (<ref>), by varying $\rho \in \{0.2, 0.4, 0.6, 0.8, 1\}$. We denote the fitted regression coefficient as $\hat\beta$.  For testing this model we generate data as follows: given $Z_i$ we use (<ref>) to generate $Y_i$ and the conditional normal model to generate $X_i \sim N(Z_i\hat\beta,\hat\sigma^2)$, where $\hat\sigma^2=\sum_{i=1}^{n_{\mathrm{train}}} (Y_i-Z_i \hat\beta)^2/(n_{\mathrm{train}}-1)$, for $1 \leq i \leq n_{\mathrm{test}}$. We will use the samples $\{(X_i, Y_i, Z_i\hat \beta)\}_{1 \leq i \leq n_{\mathrm{test}}}$ to perform the calibration test, since the training model uses the homoscedastic error (recall Remark <ref>).  Specifically, we implement the $\EMMD$ method with the Gaussian kernel with the bandwidth chosen as the median of absolute differences $\{|X_i-Y_i|:1\leq i\leq n_{\mathrm{test}}\}$. The number of nearest neighbors varies as $K \in \{2, 4, 6, 8, 10\}$. We also implement the SKCE test from [Widmann et al., 2022] for comparison.  The results are shown in Figure <ref> for $n_{\mathrm{test}} \in \{50, 75\}$. Note that due to model misspecification, the OLS model is clearly not calibrated in this case (see <cit.> for further discussion). Moreover, the true model is similar to a heteroscedastic linear model for small values of $\rho$ (since, by a Taylor expansion, $\rho\sin(\pi Z)\approx \rho \pi Z$, when $\rho$ is small), but becomes more non-linear as $\rho$ increases. Hence, the power of the tests are expected to increase as the signal strength $\rho$ increases. This aligns with the results in Figure <ref>. The plots also reveal the following:

* The $\EMMD$ tests have better power than the SKCE test and the power of the $\EMMD$ tests increase as $K$ increases.

* Even when $\rho$ is small, the $\EMMD$ test has non-trivial power, whereas the SKCE almost has no power.  This shows that the ECMMD is more sensitive to detecting the heteroscedasticity in the true model (which is approximately linear for small values of $\rho$) than the SKCE.

§.§.§ Calibration on Real Data

In this section, we apply the $\EMMD$ for testing calibration of convolutional neural networks using the CIFAR-10 dataset (<https://www.cs.toronto.edu/ kriz/cifar.html>). The dataset
consists of 10 classes of objects: airplane, automobile, bird, cat, deer, dog, frog, horse, ship, and truck. Each class has 6000 images and the total number of images is 60000. For our experiment, we use the following pairs for binary classification

\begin{align}\label{eq:imagepairs}
\{\text{bird, cat}\}, \{\text{cat, dog}\}, \{\text{cat, deer}\}, \{\text{cat, frog}\}, \{\text{cat, horse}\}.
\end{align}

The ratio of training data and test data is chosen to be $3:2$. The training data is used to learn a convolutional neural network (CNN) classifier, and the test data will used for assessing calibration based on the $\EMMD$ measure.

§.§.§ Model Setup

To train the classifier, we use three convolutional layers with 32, 64, and 64 filters, respectively, interspersed with $2\times 2$ max-pooling layers to reduce spatial dimensions. Post convolution, a flattening layer
transforms the 2D feature maps into a 1D vector, followed by a dense layer with 64 neurons. The output  layer uses a softmax activation function for binary classification.  The model employs the Adam optimizer, categorical cross-entropy loss, and tracks accuracy as its metric.

§.§.§ Calibration Test

We implement the ECMMD test with the linear kernel, which is a characteristic kernel under the Bernoulli distribution class. To test for calibration we implement the asymptotic test (<ref>) as well as the finite-sample test proposed in Algorithm <ref> with the number of resamples $M=200$. Given the relatively large sample size (around $1000$), we vary the number of nearest neighbors as $K \in \{40, 60, 80, 100\}$.

§.§.§ Potential Miscalibration and Recalibration

Following the influential paper of Guo et al., 2017, it is now common knowledge that deep-learning models tend to overfit the data which can lead potential miscalibration. To check this we will use the $\EMMD$ measure to test if the convolutional neural network trained as described above is calibrated using the pairs in (<ref>). We then split the test data according to $2:1$ ratio, use the first part of the data to recalibrate the prediction probabilities using isotonic regression (see Appendix <ref> for details), and the second part to again test using the ECMMD if  the recalibrated probabilities are indeed calibrated.

Results for bird-cat classification: (a) reliability plot before recalibration and after recalibration, and (b) $p$-values of the $\EMMD$ test for different values of $K$ and ECE before and after recalibration.

Results for cat-dog classification: (a) reliability plot before recalibration and after recalibration, and (b) $p$-values of the $\EMMD$ test for different values of $K$ and ECE before and after recalibration.

Results for cat-deer classification: (a) reliability plot before recalibration and after recalibration, and (b) $p$-values of the $\EMMD$ test for different values of $K$ and ECE before and after recalibration.

Results for cat-frog classification: (a) reliability plot before recalibration and after recalibration, and (b) $p$-values of the $\EMMD$ test for different values of $K$ and ECE before and after recalibration.

Results for cat-horse classification: (a) reliability plot before recalibration and after recalibration, and (b) $p$-values of the $\EMMD$ test for different values of $K$ and ECE before and after recalibration.

§.§.§ Results and Interpretation

To illustrate our results and the effect of recalibration, we present the reliability plots before after recalibration for the 5 pairs in (<ref>) are shown in Figure <ref>(a)-<ref>(a).

A curve aligned with the diagonal line in the reliability diagram indicates that the model is well-calibrated (see Appendix <ref> for details). The  $p$-values of the ECMMD tests, both before and after recalibration, are reported in the tables in Figure <ref>(b)-<ref>(b). We use $\tt{FS}$ and $\tt{Asymp}$ to denote the finite sample and the asymptotic test respectively. The tables also show the estimated values of the expected calibration error (ECE) computed with the number of bins set to be $100$ (see Appendix <ref> for details). The following is a summary of our findings:

* Before recalibration, the reliability diagrams and the $p$-values indicate that  CNN classifier is mis-calibrated for all the 5 pairs in (<ref>). In fact, the $p$-values of the asymptotic tests in all the 5 cases are smaller than $10^{-4}$.

* After recalibration, the ECMMD tests accept the null hypothesis (that is, the model is calibrated) at significance level  $0.05$, except for $\{\text{cat},~\text{deer}\}$ case with $K=80, 100$. Even in this case, the $p$-values are significantly higher after recalibration. The improved calibration can also be seen from the significant reduction of the ECE values.

§.§ Conditional Goodness-of-Fit Test

In this section we illustrate the performance of the resampling based conditional goodness-of-fit test from Section <ref>. We compare our method with the KCSD test, by considering the following example from [Jitkrittum et al., 2020]: Suppose

$$P_{X|Z} = N( Z^\top\one_d, g( Z)), \text{ where } g( Z)=1 + \rho\exp\left(-\frac{\| Z-1.5I_d\|_2^2 }{2\times 0.8^2}\right) \text{ and }  Z\sim N_d(\bm 0, I_d) ,$$

for $\rho > 0$. We wish to test the hypothesis (<ref>) with $P_{X|Z}$ as above, based on samples $(Y_1, Z_1)$, $(Y_2, Z_2), \ldots,$ $(Y_n, Z_n)$ from $P_{Y|Z} = N( Z^\top\one_d, 1)$ with $Z\sim N_d(\bm 0, I_d)$.[Here, $\bm 1_d$ denotes the vector all 1 in dimension $d$. Also, $N_d( \mu, \Sigma)$ denotes the $d$-dimensional normal distribution with mean vector $\mu$ and covariance matrix $\Sigma$. In dimension 1, the subscript $d$ is dropped.] Note that $\rho = 0$ corresponds to null hypothesis. In the simulations we set $d=3$ and, as in [Jitkrittum et al., 2020], we consider $\rho = 0$ and $\rho =10$.

Figure <ref> shows the empirical Type I error (with $\rho =0)$ and power (with $\rho = 10)$, averaged over 200 repetitions, of the KCSD test and the finite sample ECMMD test (from Proposition <ref>) as the sample size varies over $n \in \{100,150, 200, 250, 300, 350\}$. For the ECMMD test we use the Gaussian kernel with the bandwidth chosen as the median of absolute differences $\{ |X_i-Y_i|:1\leq i\leq n \}$. Also, the number of nearest-neighbors vary over $K \in \{5, 10, 15, 20, 25\}$. For implementing the finite sample test we use $M=300$ many resamples. The results show that the proposed test is more powerful than the KCSD, while maintaining similar Type I control. Also, as expected, power improves for larger values of $K$.

Conditional goodness-of-fit test: (a) Type-I error and (b) power of the test in (<ref>) and the KCSD test as a function of the sample size, for different values of $K$.

§.§ Comparing Regression Curves: Wind Energy Data Application

To illustrate the application of the $\EMMD$ based test in comparing two regression curves, we consider the wind energy dataset [Ding, 2019, Prakash et al., 2022, Hwangbo et al., 2017].
The dataset is publicly available at <https://aml.engr.tamu.edu/book-dswe/dswe-datasets/> under  folder.
There are four datasets in the folder, each corresponding to one wind turbine. Here, we consider the datasets WT3 and WT4 which correspond to offshore wind turbines. Both these datasets consist of 4 years (2007-2010) of data and the main variable of interest (the outcome) is the wind power generated by different turbines across different time periods. The data also has the following five covariates:  wind speed (V), wind direction (D), air density ($\rho$), turbulence intensity (I), and humidity (S).

The question of interest to researchers in this field is to understand how the wind speed (V) affects the power generation in the presence of other covariates in order to further improve the efficiency of wind turbines. Previous literature considered the problem of testing if the effect of wind speed on power generation is the same across different years. For example, [Prakash et al., 2022] adopted a Bayesian nonparameteric approach
and provide a confidence band on the regression functions compared across different years. Here, we investigate if the effect of wind speed on wind power generation remain the same across different turbines within these two periods: 2007-2008 and 2009-2010?
To alleviate the temporal effect, we preprocess the data by averaging the observations within each day. This results in a sample of size $475$ and $626$ for periods $2007$-$2008$ and $2009$-$2010$, respectively. Since the data are recorded over time and the two turbines may be located in adjacent areas, we expect the covariate information to be similar in each recording point. Denote the outcome (averaged within each day as mentioned above) for WT3 and WT4 by $X$ and $Y$, respectively. We then concatenate the covariate vector $Z_1$ for WT3 and the covariate vector $Z_2$ for WT4 as a joint vector $Z\equiv (Z_1,Z_2)$. The combined covariate vector is of size 10 (5 from each turbine). The final preprocessed data can be represented as $\{X_i, Y_i, Z_i\}_{1 \leq i \leq n}$. Based on this data, we now test the hypothesis (<ref>) using the $\EMMD$ statistic (<ref>) with the Gaussian kernel. We consider the bandwidth to be the median of absolute difference: $\mathrm{Median}\{|X_i-Y_i|:1\leq i\leq n\}$. We also test the equality of conditional expectations, that is, $\E[Y|Z]=\E[X|Z]$ almost surely, by choosing the linear kernel $\sfK(x, y) = x\cdot y$ in (<ref>). We implement the asymptotic test (<ref>) for this application. Notice whenever the null hypothesis $X|Z\overset{D}{=}Y|Z$ is accepted, we can conclude that the effect of wind speed on wind power generation remain the same across turbines WT3 and WT4, given the other covariates are fixed.

Frequency plots of wind power in WT3 ($X$) and WT4 ($Y$) during the 2007-2008 period (left) and 2009-2010 period (right).

§.§.§ Results and Interpretation

We first check the marginal distributions of $X$ and $Y$ by comparing their frequency plots based on the samples $\{X_i\}_{1 \leq i \leq n}$ and $\{Y_i\}_{1\leq i \leq n}$ (see Figure <ref>). The frequency plots show that the marginal distributions of $X$ and $Y$ are more different in the 2007-2008 period than the 2009-2010 period (although the overall trend in both cases are very similar). This suggests that there is potentially more heterogeneity between the conditional distributions $X|Z$ and $Y|Z$ in 2007-2008 period compared to the 2009-2010 period.

$K$     $\E[X|Z]=\E[Y|Z]$     $\P_{X|Z}=\P_{Y|Z}$

20     0.265     0.115

25     0.187     0.080

30     0.350      0.124

35     0.318     0.110

40     0.545     0.037

45     0.419      0.046

$K$     $\E[X|Z]=\E[Y|Z]$     $\P_{X|Z}=\P_{Y|Z}$

20     0.113     0.273

25     0.083     0.297

30     0.200      0.147

35     0.251      0.206

40     0.141     0.233

45     0.056     0.265

The $p$-values for testing equality of conditional means and  equality of conditional distributions for WT3 versus WT4 in the (a) 2007-2008 and (b) 2009-2010 period, for varying $K$.

To confirm this, in Table <ref> we report the $p$-values of the test (<ref>) with the linear kernel (for testing conditional mean equality) and the Gaussian kernel (for testing conditional distribution equality). We see that during the 2007-2008 period (see Table <ref>(a)), at significance level $0.05$, the hypothesis $\P_{X|Z} = \P_{Y|Z}$ is rejected when $K$ (the number of nearest-neighbors) is large, but the conditional mean hypothesis $\E[X|Z]=\E[Y|Z]$ is not rejected. This indicates there is higher order heterogeneity in the conditional distribution across turbines beyond the conditional mean. However, this heterogeneity is not severe since it takes large number (for example, $K=40$ and $45$) of nearest neighbors to detect the difference. For the time period 2009-2010 (see Table <ref>(b)), both the hypotheses $\P_{X|Z} = \P_{Y|Z}$ and  $\E[X|Z]=\E[Y|Z]$ are accepted, suggesting that during this period both turbines have similar conditional power generation given the wind speed and the other covariates. Similar comparison has been considered in [Prakash et al., 2022]. They provide the confidence band for the conditional mean difference of the wind power of WT3 and WT4 over wind speed for each year between 2007-2010, after adjusting the other covariates. Our result on conditional mean testing aligns with the conclusion in <cit.> that there is no significant statistical difference in conditional mean between the wind power of WT3 and WT4 between the periods 2007-2008 and 2009-2010.

§.§ Validation of Emulators in SBI

In this section we illustrate the applicability of $\EMMD$ measure in testing validity of emulators and approximate posteriors arising from for the true posterior distribution in simulation based inference (SBI).

§.§.§ Approximate Posteriors in Benchmark Examples

Examples of posterior distribution in the Two Moons task. Samples are generated using the reference posterior from the  [Lueckmann et al., 2021] package with (a) $x = (0.0030109584, -0.6457662)$ and (b) $x = (0.19458583, 1.0400153)$.

We consider the following 2 examples from the literature:

* Two Moons: This is a two-dimensional example with a posterior that exhibits both global (bimodality) and local
(crescent shape) structure (see [Greenberg et al., 2019, Lueckmann et al., 2021, Papamakarios et al., 2019, Linhart et al., 2024, Wildberger et al., 2024, Glaser et al., 2022]). Here, the data generating process has the following description:

* Generate $\theta = (\theta_1,\theta_2)^\top$ from the uniform distribution on the unit square $[-1,1]\times [-1,1]$.

* Given $\theta$ generate $X$ from the simulator:

\begin{align*}
X|\theta =
\begin{pmatrix}
r\cos(\alpha) + 0.25\\
\end{pmatrix} +
\begin{pmatrix}
\left(-\theta_1+\theta_2\right)/\sqrt{2}
\end{pmatrix} ,
\end{align*}

where $\alpha\sim \text{Unif}(-\pi/2,\pi/2)$ and $r\sim N(0.1, 0.01^2)$.

In Figure <ref> we show the heatmap of the posterior distribution $\theta | X=x$ for 2 different values of $x$ illustrating the bimodality and the cresent shaped structure.

* Simple Likelihood Complex Posterior (SLCP): This is an example with a simple likelihood that exhibits complex posterior phenomena. Here, we consider the following data generating process, which is a
simplification of the SLCP task considered in [Papamakarios et al., 2019, Greenberg et al., 2019, Hermans et al., 2020, Durkan et al., 2020, Lueckmann et al., 2021].

* Generate $\theta = (\theta_1,\theta_2,\theta_3,\theta_4,\theta_5)^\top$ from Uniform distribution on the box $[-3,3]^5$.

* Given $\theta$ generate  $ X|\theta\sim N\left({m}_{\theta}, {S}_{\theta}\right)$ where,

\begin{align*}
{m}_{\theta} =
\begin{pmatrix}
\theta_1\\
\theta_2
\end{pmatrix}\text{ and }
{S}_{\theta} =
\begin{pmatrix}
\theta_3^4 & \theta_3^2\theta_4^2\tanh(\theta_5)\\
\theta_3^2\theta_4^2\tanh(\theta_5) & \theta_4^4
\end{pmatrix} .
\end{align*}

For both the above examples the posterior distributions are intractable. To approximate the posteriors we consider the following emulators:

* Mixture Density Networks (MDNs) [Bishop, 1994, Papamakarios and Murray, 2016]: In this case, the emulator $q_\phi(\theta|\bm X)$ takes the form of a mixture of $\ell$ Gaussian components where the mixing coefficients, the means, and the covariance matrices are parametrized by a feed-forward neural network (which we denote by $\phi$) taking $X$ as an input. For the  example we consider MDNs with $1, 3, 5, 7$ mixtures and for  we consider MDNs with $5, 10, 15, 25$ mixtures.

* Conditional density estimates using Neural Spline Flow (NSF) [Durkan et al., 2019], as in [Papamakarios and Murray, 2016, Lueckmann et al., 2021]. In our experiments we implement a NSF with $5$ transforms.

Empirical power of the ECMMD test for validating the emulators of the posterior distributions in (a)  and (b)  examples as a function of training budget.

The MDN and NSF based posterior approximations are trained using their implementation in the  [Tejero-Cantero et al., 2020] package. We use a training budget (sample size) of $N\in \{10^2, 5\times 10^2, 10^3, 5\times 10^3, 10^4, 5\times 10^4, 10^5\}$ for the  example and $N\in \{10^4, 2.5\times 10^4, 5\times 10^4, 7.5\times 10^4, 10^5\}$ for the  example.  For testing the validity of the emulators we consider a test budget of $n=1000$ with $K = 50$ nearest neighbors. For each $1\leq i\leq 1000$, we generate $(\theta_i, X_i)$, where $\theta_i$ is generated from the prior and $X_i$ is generated from the simulator given $\theta_i$. Then with $X_i$ as input to the posterior emulator we generate a sample $\theta_i^\prime$, for $1\leq i\leq 1000$. Based on samples $(\theta_i^\prime, \theta_i, X_i)$, for $1\leq i\leq 1000$, we implement the asymptotic ECMMD based test from (<ref>). We use a nominal level of $\alpha = 0.05$ and the Gaussian kernel with bandwidth $\lambda = \text{Median}\left\{\left\|\theta_i-\theta_i^\prime\right\|:1\leq i\leq 1000\right\}$. The test is repeated $100$ times (with the same trained model) to estimate the proportion of rejections as a function of the training budget, for different models.

The results are shown in Figure <ref>.
In the  example, MDNs with $5$ and $7$ mixtures are able to emulate the posterior distribution (with moderate training budget). On the other hand, in the more complex  example, even with $25$ mixtures the MDN is unable to fully capture the complexities of the posterior distribution. In contrast, the NSF model with $5$ transforms is able to emulate the posterior distribution in both examples with moderate training budget.
This showcases the power of the NSF in emulating complex posterior distribution compared to MDNs.

§.§.§ Conditional Density for Photometric Redshift

In this section we illustrate the applicability of the $\EMMD$ measure in testing  emulators for simulating redshifts associated with galaxy images. Specifically, our goal is to validate conditional density estimates of the synthetic `redshift' $X$ associated with photometric or `photo-z' galaxy images $Z$. This example is motivated by the diagnostic studies on conditional density of photometric redshifts in [Zhao et al., 2021, Dey et al., 2022, LeRoy and Zhao, 2021] and presents the $\EMMD$ test as a powerful method for validating photo-z probability densities.

In our experiments, we use , the open source toolkit for simulating realistic images of astronomical objects introduced in [Rowe et al., 2015],  to generate $20\times 20$-pixel images of galaxies. The parameters are set as follows:

* The galaxy's rotational angle $\alpha$ with respect to the $x$-axis is chosen from $\text{Unif}\left[-\pi/2,\pi/2\right]$.

* The parameter $\lambda$, which is the ratio between the minor and major axes of the projection of the elliptical galaxy, is set to be $0.75$ with a small uniform noise generated from $\text{Unif}\left[-0.1,0.1\right]$
added to it.

These random parameters $\alpha$ and $\lambda$ are used as input to  to generate a galaxy image $\bm Z \in [0, 1]^{20 \times 20}$. Example of galaxy images with $(\alpha, \lambda) = (-\pi/4,  0.75)$ and $(\alpha, \lambda) = (\pi/4,  0.75)$ are shown in Figure <ref>.

Galaxy images generated by  with $\lambda = 0.75$ and (a) $\alpha = -\pi/4$, (b) $\alpha = \pi/4$.

To simulate the `redshift' $X$ associated with $\bm Z$ we consider the following multimodal setup:

\begin{align}\label{eq:Xgivenalpha}
\begin{cases}
\frac{1}{2} N \left(2,1\right) + \frac{1}{2} N \left(-2,1\right) & \text{ if }\alpha\geq \theta\\
\frac{1}{4} N\left(3,0.0625\right) + \frac{1}{4}N\left(\mu,0.0625\right) + \frac{1}{4}N\left(-3,0.0625\right) + \frac{1}{4}N\left(-\mu,0.0625\right) & \text{otherwise,}
\end{cases}
\end{align}

for fixed parameters $\theta\in \left\{-\pi/3, -\pi/4, -\pi/6, -\pi/12, 0, \pi/12, \pi/6, \pi/4, \pi/3\right\}$ and $\mu\in\{0, 0.5, 1, 1.5\}$. The marginal distribution of $X$ for $\theta \in  \{-\pi/3, 0, \pi/3\}$ and different values of $\mu$ are shown in Figure <ref>.

Marginal distribution of `redshift' $X$ generated according to (<ref>) for (a) $\theta = -\pi/3$, (b) $\theta = 0$, and (c) $\theta = \pi/3$. The density and histogram are generated using $10000$ samples.

In our experiment we use the `redshift' value $X$ and the galaxy image $\bm Z$ as observations, and the parameter $\theta$ remains implicit and hence,  unobserved. Based on the marginal distributions and following the experiments from [Zhao et al., 2021, LeRoy and Zhao, 2021] we model the conditional distribution of $X|\bm Z$ as a mixture of Gaussians. In particular,  we consider a Gaussian convolutional mixture density network (ConvMDN) [D’Isanto and Polsterer, 2018] with $2$ mixtures. This will play the role of the emulator. To fit the ConvMDN we use a convolutional neural network with two convolutional layers interspersed with ReLU layers and having filters of size $3$ with $6$ and $12$ channels,  respectively. Post convolution, after passing through a $2\times 2$ max-pooling layer, a flattening layer transforms the 2D feature maps into a 1D vector, followed by two fully connected layers with $128$ and $10$ neurons and an output layer with $2$ neurons. The means, standard deviation, and weights of the mixtures model are determined using the output layer and the mechanism from <cit.>. To train the models we use the negative log-likelihood loss and optimize through the Adam optimizer  [Kingma and Ba, 2014] with learning rate $10^{-3}$ and parameters $\beta_1 = 0.9$ and $\beta_2 = 0.999$.

For training the ConvMDN, we use $N = 20000$ training samples generated according to the above schematics and to test the accuracy of the trained model we generate $n = 2500$ testing samples. For each $1\leq i\leq 2500$, we
generate $(X_i, \bm Z_i)$, where $\bm Z_i$ is a galaxy image generated from  and $X_i$ follows (<ref>), with parameters chosen as above. Also, with the galaxy image $\bm Z_i$ as input to the trained ConvMDN model, we generate a sample $Y_i$, for $1\leq i\leq 2500$. To evaluate the ConvMDN model in emulating redshifts, we now test whether or not $X|\bm Z$ (the simulator) has the same distribution as $Y |\bm Z$ (the emulator) using the asymptotic $\EMMD$ test from (<ref>), based on the samples $(X_i, Y_i, \bm Z_i)$, for $1 \leq i \leq n$. We use a nominal level of $\alpha=0.05$, the number of nearest-neighbors $K= 100$, and the Gaussian kernel $\sfK_{\lambda}(x,y) = \exp\left(-(x-y)^2/\lambda^2\right)$ with bandwidth $\lambda = \text{Median}\{|X_i-Y_i|: 1\leq i\leq 2500\}$. The test is repeated $100$ times (with the same trained model) to estimate the proportion of rejections as a function of $\theta$, for different values of $\mu$.

Empirical power of the $\EMMD$ test for validating the  Gaussian ConvMDN emulator for generating redshifts as a function of $\theta$:  (a) for $n=2500$, $K=100$, and different values of $\mu$, and (b) for $\mu=0.5$ and different values of $n$ and $K$.

The results are shown in Figure <ref>. We observe the following:

$1.$ the power of the $\EMMD$ test is more for larger values of $\theta$; and

$2.$ the power decreases as $\mu$ increases.

This can be explained by looking at the marginal distribution of $X$.
Note that the marginal distribution of $X$ has 2 modes for small values of $\theta$ (see Figure <ref>(a)), hence, the ConvMDN model provides a good fit in this case. Therefore, the $\EMMD$ test accepts the null hypothesis $\P_{X|\bm Z} = \P_{Y |\bm Z}$ more often when $\theta$ is small. As $\theta$ increases, the number of modes increase, which leads to more rejections, hence, the increasing trends in the power curves in Figure <ref>(a). Further, even though a $2$ component Gaussian mixture seems inadequate for larger values of $\theta$, when $\mu$ is also large, the modes tend to flatten out (see the histograms for $\mu=1.5$ in Figure <ref>), which makes it harder to distinguish the true distribution from a $2$-mixture of Gaussian. Hence, the rates of increase of the power curves in Figure <ref>(a) are slower for larger values of $\mu$.
In Figure <ref>(b) we fix $\mu=0.5$ and show the proportion of rejections as a function of $\theta$, for different values of $n$ and $K$. As expected, the power improves as $n$ and $K$ increases and the $\EMMD$ test accurately detects the differences between the 2 distributions.
Hence, the $\EMMD$ test successfully validates whether the emulator model is a reasonable surrogate for simulating photometric redshift. This also illustrates how the proposed method can be more generally applicable as a model diagnostic tool within the realm of simulation-based inference.

§.§ Acknowledgements
B. B. Bhattacharya was supported by NSF CAREER grant DMS 2046393, NSF grant DMS 2113771, and a Sloan Research Fellowship. The authors also thank Lucas Janson for helpful discussions.

[Alfonso et al., 2023]
J. Alfonso, R. Baptista, A. Bhakta, N. Gal, A. Hou, I. Lyubimova, D. Pocklington, J. Sajonz, G. Trigila, and R. Tsai.
A generative flow for conditional sampling via optimal transport.
arXiv:2307.04102, 2023.

[Andrews, 1997]
D. W. Andrews.
A conditional kolmogorov test.
Econometrica: Journal of the Econometric Society, pages 1097–1128, 1997.

[Aronszajn, 1950]
N. Aronszajn.
Theory of reproducing kernels.
Transactions of the American mathematical society, 680 (3):0 337–404, 1950.

[Arrieta-Ibarra et al., 2022]
I. Arrieta-Ibarra, P. Gujral, J. Tannen, M. Tygert, and C. Xu.
Metrics of calibration for probabilistic predictions.
Journal of Machine Learning Research, 230 (351):0 1–54, 2022.

[Azadkia and Chatterjee, 2021]
M. Azadkia and S. Chatterjee.
A simple measure of conditional dependence.
The Annals of Statistics, 490 (6):0 3070–3102, 2021.

[Baringhaus and Franz, 2004]
L. Baringhaus and C. Franz.
On a new multivariate two-sample test.
J. Multivariate Anal., 880 (1):0 190–206, 2004.

[Batzolis et al., 2021]
G. Batzolis, J. Stanczuk, C.-B. Schönlieb, and C. Etmann.
Conditional image generation with score-based diffusion models.
arXiv:2111.13606, 2021.

[Berlinet and Thomas-Agnan, 2011]
A. Berlinet and C. Thomas-Agnan.
Reproducing kernel Hilbert spaces in probability and statistics.
Springer Science & Business Media, 2011.

[Berta et al., 2024]
E. Berta, F. Bach, and M. Jordan.
Classifier calibration with roc-regularized isotonic regression.
In International Conference on Artificial Intelligence and Statistics, pages 1972–1980. PMLR, 2024.

[Bishop, 1994]
C. M. Bishop.
Mixture density networks.

[Boucheron et al., 2013]
S. Boucheron, G. Lugosi, and P. Massart.
Concentration Inequalities: A Nonasymptotic Theory of Independence.
Oxford University Press, 2013.

[Bröcker and Smith, 2007]
J. Bröcker and L. A. Smith.
Increasing the reliability of reliability diagrams.
Weather and forecasting, 220 (3):0 651–661, 2007.

[Bruinsma et al., 2023]
W. P. Bruinsma, S. Markou, J. Requiema, A. Y. Foong, T. R. Andersson, A. Vaughan, A. Buonomo, J. S. Hosking, and R. E. Turner.
Autoregressive conditional neural processes.
arXiv:2303.14468, 2023.

[Candes et al., 2018]
E. Candes, Y. Fan, L. Janson, and J. Lv.
Panning for gold: `Model-X' knockoffs for high dimensional controlled variable selection.
Journal of the Royal Statistical Society Series B: Statistical Methodology, 800 (3):0 551–577, 2018.

[Chatterjee and Bhattacharya, 2023]
A. Chatterjee and B. B. Bhattacharya.
Boosting the power of kernel two-sample tests.
arXiv:2302.10687, 2023.

[Chaudhuri and Marron, 1999]
P. Chaudhuri and J. S. Marron.
Sizer for exploration of structures in curves.
Journal of the American Statistical Association, 940 (447):0 807–823, 1999.

[Chen and Shao, 2004]
L. H. Chen and Q.-M. Shao.
Normal approximation under local dependence.
The Annals of Probability, 320 (3A):0 1985–2028, 2004.

[Chen et al., 2022]
M. Chen, T. Tian, J. Zhu, W. Pan, and X. Wang.
Paired-sample tests for homogeneity with/without confounding variables.
Statistics and Its Interface, 150 (3):0 335–348, 2022.

[Chen and Lei, 2024]
Y. Chen and J. Lei.
De-biased two-sample $u$-statistics with application to conditional distribution testing.
arXiv:2402.00164, 2024.

[Cook et al., 2006]
S. R. Cook, A. Gelman, and D. B. Rubin.
Validation of software for bayesian models using posterior quantiles.
Journal of Computational and Graphical Statistics, 150 (3):0 675–692, 2006.

[Cranmer et al., 2020]
K. Cranmer, J. Brehmer, and G. Louppe.
The frontier of simulation-based inference.
Proceedings of the National Academy of Sciences, 1170 (48):0 30055–30062, 2020.

[Dalmasso et al., 2020]
N. Dalmasso, R. Izbicki, and A. Lee.
Confidence sets and hypothesis testing in a likelihood-free inference setting.
In International Conference on Machine Learning, pages 2323–2334. PMLR, 2020a.

[Dalmasso et al., 2020]
N. Dalmasso, A. Lee, R. Izbicki, T. Pospisil, I. Kim, and C.-A. Lin.
Validation of approximate likelihood and emulator models for computationally intensive simulations.
In International Conference on Artificial Intelligence and Statistics, pages 3349–3361. PMLR, 2020b.

[Dalmasso et al., 2021]
N. Dalmasso, L. Masserano, D. Zhao, R. Izbicki, and A. B. Lee.
Likelihood-free frequentist inference: Confidence sets with correct conditional coverage.
arXiv:2107.03920, 2021.

[Dax et al., 2021]
M. Dax, S. R. Green, J. Gair, J. H. Macke, A. Buonanno, and B. Schölkopf.
Real-time gravitational wave science with neural posterior estimation.
Physical review letters, 1270 (24):0 241103, 2021.

[Deb et al., 2020]
N. Deb, P. Ghosal, and B. Sen.
Measuring association on topological spaces using kernels and geometric graphs.
arXiv:2010.01768, 2020.

[DeGroot and Fienberg, 1983]
M. H. DeGroot and S. E. Fienberg.
The comparison and evaluation of forecasters.
Journal of the Royal Statistical Society: Series D (The Statistician), 320 (1-2):0 12–22, 1983.

[Deistler et al., 2022]
M. Deistler, P. J. Goncalves, and J. H. Macke.
Truncated proposals for scalable and hassle-free simulation-based inference.
Advances in Neural Information Processing Systems, 35:0 23135–23149, 2022.

[Dette and Munk, 1998]
H. Dette and A. Munk.
Nonparametric comparison of several regression functions: exact and asymptotic theory.
The Annals of Statistics, 260 (6):0 2339–2368, 1998.

[Dey et al., 2022]
B. Dey, D. Zhao, B. Andrews, J. Newman, R. Izbicki, and A. Lee.
Calibrated predictive distributions for photometric redshifts.
In Machine Learning for Astrophysics, proceedings of the Thirty-ninth International Conference on Machine Learning (ICML 2022), 2022.

[Dinev and Gutmann, 2018]
T. Dinev and M. U. Gutmann.
Dynamic likelihood-free inference via ratio estimation (dire).
arXiv preprint arXiv:1810.09899, 2018.

[Ding, 2019]
Y. Ding.
Data science for wind energy.
CRC Press, 2019.

[Dirmeier et al., 2023]
S. Dirmeier, C. Albert, and F. Perez-Cruz.
Simulation-based inference using surjective sequential neural likelihood estimation.
arXiv preprint arXiv:2308.01054, 2023.

[Durkan et al., 2019]
C. Durkan, A. Bekasov, I. Murray, and G. Papamakarios.
Neural spline flows.
Advances in neural information processing systems, 32, 2019.

[Durkan et al., 2020]
C. Durkan, I. Murray, and G. Papamakarios.
On contrastive learning for likelihood-free inference.
In International conference on machine learning, pages 2771–2781. PMLR, 2020.

[D’Isanto and Polsterer, 2018]
A. D’Isanto and K. L. Polsterer.
Photometric redshift estimation via deep learning-generalized and pre-classification-less, image based, fully probabilistic redshifts.
Astronomy & Astrophysics, 609:0 A111, 2018.

[Efron and Stein, 1981]
B. Efron and C. Stein.
The jackknife estimate of variance.
The Annals of Statistics, pages 586–596, 1981.

[Fan and Lin, 1998]
J. Fan and S.-K. Lin.
Test of significance when data are curves.
Journal of the American Statistical Association, 930 (443):0 1007–1021, 1998.

[Frank, 2015]
E. H. Frank.
Regression modeling strategies with applications to linear models, logistic and ordinal regression, and survival analysis, 2015.

[Friedman, 2003]
J. H. Friedman.
On multivariate goodness–of–fit and two–sample testing.
Statistical Problems in Particle Physics, Astrophysics, and Cosmology, 1:0 311, 2003.

[Friedman et al., 1977]
J. H. Friedman, J. L. Bentley, and R. A. Finkel.
An algorithm for finding best matches in logarithmic expected time.
ACM Transactions on Mathematical Software (TOMS), 30 (3):0 209–226, 1977.

[Glaser et al., 2022]
P. Glaser, M. Arbel, A. Doucet, and A. Gretton.
Maximum likelihood learning of energy-based models for simulation-based inference.

[Glaser et al., 2023]
P. Glaser, D. Widmann, F. Lindsten, and A. Gretton.
Fast and scalable score-based kernel calibration tests.
In Uncertainty in Artificial Intelligence, pages 691–700. PMLR, 2023.

[Glöckler et al., 2022]
M. Glöckler, M. Deistler, and J. H. Macke.
Variational methods for simulation-based inference.
arXiv preprint arXiv:2203.04176, 2022.

[Gloeckler et al., 2024]
M. Gloeckler, M. Deistler, C. Weilbach, F. Wood, and J. H. Macke.
All-in-one simulation-based inference.
arXiv preprint arXiv:2404.09636, 2024.

[Gneiting et al., 2007]
T. Gneiting, F. Balabdaoui, and A. E. Raftery.
Probabilistic forecasts, calibration and sharpness.
Journal of the Royal Statistical Society Series B: Statistical Methodology, 690 (2):0 243–268, 2007.

[Gonçalves et al., 2020]
P. J. Gonçalves, J.-M. Lueckmann, M. Deistler, M. Nonnenmacher, K. Öcal, G. Bassetto, C. Chintaluri, W. F. Podlaski, S. A. Haddad, T. P. Vogels, et al.
Training deep neural density estimators to identify mechanistic models of neural dynamics.
Elife, 9:0 e56261, 2020.

[Greenberg et al., 2019]
D. Greenberg, M. Nonnenmacher, and J. Macke.
Automatic posterior transformation for likelihood-free inference.
In International Conference on Machine Learning, pages 2404–2414. PMLR, 2019.

[Gretton et al., 2009]
A. Gretton, K. Fukumizu, Z. Harchaoui, and B. K. Sriperumbudur.
A fast, consistent kernel two-sample test.
Advances in neural information processing systems, 22, 2009.

[Gretton et al., 2012]
A. Gretton, K. M. Borgwardt, M. J. Rasch, B. Schölkopf, and A. Smola.
A kernel two-sample test.
The Journal of Machine Learning Research, 130 (1):0 723–773, 2012.

[Guo et al., 2017]
C. Guo, G. Pleiss, Y. Sun, and K. Q. Weinberger.
On calibration of modern neural networks.
In International conference on machine learning, pages 1321–1330. PMLR, 2017.

[Gweon and Yu, 2019]
H. Gweon and H. Yu.
How reliable is your reliability diagram?
Pattern Recognition Letters, 125:0 687–693, 2019.

[Hansen, 1994]
B. E. Hansen.
Autoregressive conditional density estimation.
International Economic Review, pages 705–730, 1994.

[Hardle and Marron, 1990]
W. Hardle and J. S. Marron.
Semiparametric comparison of regression curves.
The Annals of Statistics, pages 63–89, 1990.

[Hermans et al., 2020]
J. Hermans, V. Begy, and G. Louppe.
Likelihood-free mcmc with amortized approximate ratio estimators.
In International conference on machine learning, pages 4239–4248. PMLR, 2020.

[Ho and Lee, 2005]
Y. H. Ho and S. M. Lee.
Calibrated interpolated confidence intervals for population quantiles.
Biometrika, 920 (1):0 234–241, 2005.

[Hocking, 2013]
R. R. Hocking.
Methods and applications of linear models: regression and the analysis of variance.
John Wiley & Sons, 2013.

[Hu and Lei, 2023]
X. Hu and J. Lei.
A two-sample conditional distribution test using conformal prediction and weighted rank sum.
Journal of the American Statistical Association, pages 1–19, 2023.

[Huang et al., 2022]
Z. Huang, N. Deb, and B. Sen.
Kernel partial correlation coefficient?a measure of conditional dependence.
The Journal of Machine Learning Research, 230 (1):0 9699–9756, 2022a.

[Huang et al., 2022]
Z. Huang, H. Lam, and H. Zhang.
Evaluating aleatoric uncertainty via conditional generative models.
arXiv:2206.04287, 2022b.

[Hwangbo et al., 2017]
H. Hwangbo, A. Johnson, and Y. Ding.
A production economics analysis for quantifying the efficiency of wind turbines.
Wind Energy, 200 (9):0 1501–1513, 2017.

[Jaffe et al., 2020]
A. Jaffe, Y. Kluger, G. C. Linderman, G. Mishne, and S. Steinerberger.
Randomized near-neighbor graphs, giant components and applications in data science.
Journal of applied probability, 570 (2):0 458–476, 2020.

[Jitkrittum et al., 2020]
W. Jitkrittum, H. Kanagawa, and B. Schölkopf.
Testing goodness of fit of conditional density models with kernels.
In Conference on Uncertainty in Artificial Intelligence, pages 221–230. PMLR, 2020.

[Kallenberg, 1997]
O. Kallenberg.
Foundations of modern probability, volume 2.
Springer, 1997.

[King et al., 1991]
E. King, J. D. Hart, and T. E. Wehrly.
Testing the equality of two regression curves using linear smoothers.
Statistics & Probability Letters, 120 (3):0 239–247, 1991.

[Kingma and Ba, 2014]
D. P. Kingma and J. Ba.
Adam: A method for stochastic optimization.
arXiv:1412.6980, 2014.

[Klenke, 2013]
A. Klenke.
Probability theory: a comprehensive course.
Springer Science & Business Media, 2013.

[Kulasekera, 1995]
K. Kulasekera.
Comparison of regression curves using quasi-residuals.
Journal of the American Statistical Association, 900 (431):0 1085–1093, 1995.

[Kull and Flach, 2015]
M. Kull and P. Flach.
Novel decompositions of proper scoring rules for classification: Score adjustment as precursor to calibration.
In Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2015, Porto, Portugal, September 7-11, 2015, Proceedings, Part I 15, pages 68–85. Springer, 2015.

[Kumar et al., 2018]
A. Kumar, S. Sarawagi, and U. Jain.
Trainable calibration measures for neural networks from kernel mean embeddings.
In International Conference on Machine Learning, pages 2805–2814. PMLR, 2018.

[Lee et al., 2023]
D. Lee, X. Huang, H. Hassani, and E. Dobriban.
T-cal: An optimal test for the calibration of predictive models.
Journal of Machine Learning Research, 240 (335):0 1–72, 2023.

[Lemos et al., 2023]
P. Lemos, A. Coogan, Y. Hezaveh, and L. Perreault-Levasseur.
Sampling-based accuracy testing of posterior estimators for general inference.
arXiv:2302.03026, 2023.

[LeRoy and Zhao, 2021]
B. LeRoy and D. Zhao.
Md-split+: Practical local conformal inference in high dimensions.
arXiv:2107.03280, 2021.

[Linhart et al., 2024]
J. Linhart, G. V. Cardoso, A. Gramfort, S. L. Corff, and P. L. Rodrigues.
Diffusion posterior sampling for simulation-based inference in tall data settings.
arXiv preprint arXiv:2404.07593, 2024a.

[Linhart et al., 2024]
J. Linhart, A. Gramfort, and P. Rodrigues.
L-c2st: Local diagnostics for posterior approximations in simulation-based inference.
Advances in Neural Information Processing Systems, 36, 2024b.

[Liu et al., 2016]
Q. Liu, J. Lee, and M. Jordan.
A kernelized stein discrepancy for goodness-of-fit tests.
In International conference on machine learning, pages 276–284. PMLR, 2016.

[Liu et al., 2021]
S. Liu, X. Zhou, Y. Jiao, and J. Huang.
Wasserstein generative learning of conditional distribution.
arXiv:2112.10039, 2021.

[Liu et al., 2023]
S. Liu, Z. Cai, and X. Li.
Distribution-free model-agnostic regression calibration via nonparametric methods.
arXiv:2305.12283, 2023.

[Lopez-Paz and Oquab, 2022]
D. Lopez-Paz and M. Oquab.
Revisiting classifier two-sample tests.
In International Conference on Learning Representations, 2022.

[Lueckmann et al., 2017]
J.-M. Lueckmann, P. J. Goncalves, G. Bassetto, K. Öcal, M. Nonnenmacher, and J. H. Macke.
Flexible statistical inference for mechanistic models of neural dynamics.
Advances in neural information processing systems, 30, 2017.

[Lueckmann et al., 2019]
J.-M. Lueckmann, G. Bassetto, T. Karaletsos, and J. H. Macke.
Likelihood-free inference with emulator networks.
In Symposium on Advances in Approximate Bayesian Inference, pages 32–53. PMLR, 2019.

[Lueckmann et al., 2021]
J.-M. Lueckmann, J. Boelts, D. Greenberg, P. Goncalves, and J. Macke.
Benchmarking simulation-based inference.
In International conference on artificial intelligence and statistics, pages 343–351. PMLR, 2021.

[Marin et al., 2012]
J.-M. Marin, P. Pudlo, C. P. Robert, and R. J. Ryder.
Approximate bayesian computational methods.
Statistics and computing, 220 (6):0 1167–1180, 2012.

[Moreira, 2003]
M. J. Moreira.
A conditional likelihood ratio test for structural models.
Econometrica, 710 (4):0 1027–1048, 2003.

[Murphy and Winkler, 1977]
A. H. Murphy and R. L. Winkler.
Reliability of subjective probability forecasts of precipitation and temperature.
Journal of the Royal Statistical Society Series C: Applied Statistics, 260 (1):0 41–47, 1977.

[Naeini et al., 2015]
M. P. Naeini, G. Cooper, and M. Hauskrecht.
Obtaining well calibrated probabilities using bayesian binning.
In Proceedings of the AAAI conference on artificial intelligence, volume 29, 2015.

[Neumeyer and Dette, 2003]
N. Neumeyer and H. Dette.
Nonparametric comparison of regression curves: an empirical process approach.
The Annals of Statistics, 310 (3):0 880–920, 2003.

[Niculescu-Mizil and Caruana, 2005]
A. Niculescu-Mizil and R. Caruana.
Predicting good probabilities with supervised learning.
In Proceedings of the 22nd international conference on Machine learning, pages 625–632, 2005.

[Papamakarios and Murray, 2016]
G. Papamakarios and I. Murray.
Fast $\varepsilon$-free inference of simulation models with bayesian conditional density estimation.
Advances in neural information processing systems, 29, 2016.

[Papamakarios et al., 2017]
G. Papamakarios, T. Pavlakou, and I. Murray.
Masked autoregressive flow for density estimation.
Advances in neural information processing systems, 30, 2017.

[Papamakarios et al., 2019]
G. Papamakarios, D. Sterratt, and I. Murray.
Sequential neural likelihood: Fast likelihood-free inference with autoregressive flows.
In The 22nd International Conference on Artificial Intelligence and Statistics, pages 837–848. PMLR, 2019.

[Papamakarios et al., 2021]
G. Papamakarios, E. Nalisnick, D. J. Rezende, S. Mohamed, and B. Lakshminarayanan.
Normalizing flows for probabilistic modeling and inference.
Journal of Machine Learning Research, 220 (57):0 1–64, 2021.

[Pardo-Fernández et al., 2015]
J. C. Pardo-Fernández, M. D. Jiménez-Gamero, and A. El Ghouch.
Tests for the equality of conditional variance functions in nonparametric regression.
Electronic Journal of Statistics, 9:0 1826–1851, 2015.

[Park and Muandet, 2020]
J. Park and K. Muandet.
A measure-theoretic approach to kernel conditional mean embeddings.
Advances in neural information processing systems, 33:0 21247–21259, 2020.

[Prakash et al., 2022]
A. Prakash, R. Tuo, and Y. Ding.
Gaussian process-aided function comparison using noisy scattered data.
Technometrics, 640 (1):0 92–102, 2022.

[Ramesh et al., 2022]
P. Ramesh, J.-M. Lueckmann, J. Boelts, Á. Tejero-Cantero, D. S. Greenberg, P. J. Gonçalves, and J. H. Macke.
Gatsbi: Generative adversarial training for simulation-based inference.
arXiv preprint arXiv:2203.06481, 2022.

[Reed and Simon, 1980]
M. Reed and B. Simon.
Methods of Modern Mathematical Physics. vol. 1. Functional Analysis.
Academic New York, 1980.

[Ren et al., 2016]
Y. Ren, J. Zhu, J. Li, and Y. Luo.
Conditional generative moment-matching networks.
Advances in Neural Information Processing Systems, 29, 2016.

[Rodrigues et al., 2021]
P. Rodrigues, T. Moreau, G. Louppe, and A. Gramfort.
Hnpe: Leveraging global parameters for neural posterior estimation.
Advances in Neural Information Processing Systems, 34:0 13432–13443, 2021.

[Rowe et al., 2015]
B. T. Rowe, M. Jarvis, R. Mandelbaum, G. M. Bernstein, J. Bosch, M. Simet, J. E. Meyers, T. Kacprzak, R. Nakajima, J. Zuntz, et al.
Galsim: The modular galaxy image simulation toolkit.
Astronomy and Computing, 10:0 121–150, 2015.

[Rueda et al., 2007]
M. Rueda, S. Martínez-Puertas, H. Martínez-Puertas, and A. Arcos.
Calibration methods for estimating quantiles.
Metrika, 66:0 355–371, 2007.

[Schmidt et al., 2020]
S. Schmidt, A. Malz, J. Soo, I. Almosallam, M. Brescia, S. Cavuoti, J. Cohen-Tanugi, A. Connolly, J. DeRose, P. Freeman, et al.
Evaluation of probabilistic photometric redshift estimation approaches for the rubin observatory legacy survey of space and time (lsst).
Monthly Notices of the Royal Astronomical Society, 4990 (2):0 1587–1606, 2020.

[Schrab et al., 2023]
A. Schrab, I. Kim, M. Albert, B. Laurent, B. Guedj, and A. Gretton.
Mmd aggregated two-sample test.
Journal of Machine Learning Research (JMLR), 24, 2023.

[Sejdinovic et al., 2013]
D. Sejdinovic, B. Sriperumbudur, A. Gretton, and K. Fukumizu.
Equivalence of distance-based and rkhs-based statistics in hypothesis testing.
The annals of statistics, pages 2263–2291, 2013.

[Shi et al., 2021]
H. Shi, M. Drton, and F. Han.
On azadkia-chatterjee's conditional dependence coefficient.
arXiv:2108.06827, 2021.

[Shi et al., 2022]
Y. Shi, V. De Bortoli, G. Deligiannidis, and A. Doucet.
Conditional simulation using diffusion schrödinger bridges.
In Uncertainty in Artificial Intelligence, pages 1792–1802. PMLR, 2022.

[Song et al., 2019]
H. Song, T. Diethe, M. Kull, and P. Flach.
Distribution calibration for regression.
In International Conference on Machine Learning, pages 5897–5906. PMLR, 2019.

[Song et al., 2023]
S. Song, T. Wang, G. Shen, Y. Lin, and J. Huang.
Wasserstein generative regression.
arXiv:2306.15163, 2023.

[Sriperumbudur et al., 2008]
B. K. Sriperumbudur, A. Gretton, K. Fukumizu, G. Lanckriet, and B. Schölkopf.
Injective hilbert space embeddings of probability measures.
In 21st annual conference on learning theory (COLT 2008), pages 111–122. Omnipress, 2008.

[Stute and Zhu, 2002]
W. Stute and L.-X. Zhu.
Model checks for generalized linear models.
Scandinavian Journal of Statistics, 290 (3):0 535–545, 2002.

[Székely, 2003]
G. J. Székely.
E-statistics: The energy of statistical samples.
Bowling Green State University, Department of Mathematics and Statistics Technical Report, 30 (05):0 1–18, 2003.

[Székely and Rizzo, 2004]
G. J. Székely and M. L. Rizzo.
Testing for equal distributions in high dimension.
InterStat, 50 (16.10):0 1249–1272, 2004.

[Taillardat et al., 2016]
M. Taillardat, O. Mestre, M. Zamo, and P. Naveau.
Calibrated ensemble forecasts using quantile regression forests and ensemble model output statistics.
Monthly Weather Review, 1440 (6):0 2375–2393, 2016.

[Talts et al., 2018]
S. Talts, M. Betancourt, D. Simpson, A. Vehtari, and A. Gelman.
Validating bayesian inference algorithms with simulation-based calibration.
arXiv preprint arXiv:1804.06788, 2018.

[Tejero-Cantero et al., 2020]
A. Tejero-Cantero, J. Boelts, M. Deistler, J.-M. Lueckmann, C. Durkan, P. J. Gonçalves, D. S. Greenberg, and J. H. Macke.
Sbi–a toolkit for simulation-based inference.
arXiv preprint arXiv:2007.09114, 2020.

[To and Yip, 1975]
T. O. To and K. W. Yip.
A generalized jensen's inequality.
Pacific J. Math., 610 (2):0 255–259, 1975.

[Trippe and Turner, 2018]
B. L. Trippe and R. E. Turner.
Conditional density estimation with bayesian normalising flows.
arXiv:1802.04908, 2018.

[Vaicenavicius et al., 2019]
J. Vaicenavicius, D. Widmann, C. Andersson, F. Lindsten, J. Roll, and T. Schön.
Evaluating model calibration in classification.
In The 22nd International Conference on Artificial Intelligence and Statistics, pages 3459–3467. PMLR, 2019.

[Wang et al., 2015]
X. Wang, W. Pan, W. Hu, Y. Tian, and H. Zhang.
Conditional distance correlation.
Journal of the American Statistical Association, 1100 (512):0 1726–1734, 2015.

[Wenger et al., 2020]
J. Wenger, H. Kjellström, and R. Triebel.
Non-parametric calibration for classification.
In International Conference on Artificial Intelligence and Statistics, pages 178–190. PMLR, 2020.

[Widmann et al., 2019]
D. Widmann, F. Lindsten, and D. Zachariah.
Calibration tests in multi-class classification: A unifying framework.
Advances in neural information processing systems, 32, 2019.

[Widmann et al., 2022]
D. Widmann, F. Lindsten, and D. Zachariah.
Calibration tests beyond classification.
arXiv:2210.13355, 2022.

[Wildberger et al., 2024]
J. Wildberger, M. Dax, S. Buchholz, S. Green, J. H. Macke, and B. Schölkopf.
Flow matching for scalable simulation-based inference.
Advances in Neural Information Processing Systems, 36, 2024.

[Wiqvist et al., 2021]
S. Wiqvist, J. Frellsen, and U. Picchini.
Sequential neural posterior and likelihood approximation.
arXiv preprint arXiv:2102.06522, 2021.

[Wood, 2010]
S. N. Wood.
Statistical inference for noisy nonlinear ecological dynamic systems.
Nature, 4660 (7310):0 1102–1104, 2010.

[Yan and Zhang, 2022]
J. Yan and X. Zhang.
A nonparametric two-sample conditional distribution test.
arXiv:2210.08149, 2022.

[Zadrozny and Elkan, 2002]
B. Zadrozny and C. Elkan.
Transforming classifier scores into accurate multiclass probability estimates.
In Proceedings of the eighth ACM SIGKDD international conference on Knowledge discovery and data mining, pages 694–699, 2002.

[Zhang, 2022]
L. Zhang.
Inference on nonparametric targets and discrete structures.
PhD thesis, Harvard University, 2022.

[Zhao et al., 2021]
D. Zhao, N. Dalmasso, R. Izbicki, and A. B. Lee.
Diagnostics for conditional density models and bayesian inference algorithms.
In Uncertainty in Artificial Intelligence, pages 1830–1840. PMLR, 2021.

[Zheng, 2000]
J. X. Zheng.
A consistent test of conditional parametric distributions.
Econometric Theory, 160 (5):0 667–691, 2000.

[Zhong and Kwok, 2013]
W. Zhong and J. T. Kwok.
Accurate probability calibration for multiple classifiers.
In Twenty-third international joint conference on artificial intelligence, 2013.

[Zhou et al., 2022]
X. Zhou, Y. Jiao, J. Liu, and J. Huang.
A deep generative approach to conditional sampling.
Journal of the American Statistical Association, pages 1–12, 2022.

§ PROOF OF PROPOSITION <REF> AND PROPOSITION <REF>

§.§ Proof of Proposition <ref>

To begin with suppose $\P_{Z \sim P_Z}\left[ P_{X|Z} =  P_{Y|Z}\right] = 1$. , by recalling (<ref>), the conditional mean embeddings $\mu_{ P_{X|Z}}$ and $\mu_{ P_{X|Z}}$ are equal on a $P_Z$-almost sure set. Hence, (<ref>) implies, $\EMMD^2\left[\cF,  P_{X|Z}, P_{Y|Z}\right] = 0$.

Now, suppose $\EMMD^2\left[\cF,  P_{X|Z}, P_{Y|Z}\right] = 0$. Then by (<ref>),

\begin{align}\label{eq:condmeanembed0}
\left\|\mu_{ P_{X|Z}} - \mu_{ P_{Y|Z}}\right\|_{\cH} = 0\text{ almost surely } P_Z .
\end{align}

Since $\sfK$ is characteristic by Assumption <ref>, recalling the definition of characteristic kernels from Sriperumbudur et al., 2008 the identity from (<ref>) implies, $P_{X|Z} =  P_{Y|Z}$, on a $P_Z$-almost sure set, which completes the proof of Proposition <ref>.

§.§ Proof of Proposition <ref>

From the definition of the conditional mean embeddings in (<ref>),

\begin{align}\label{eq:EMMDequivKXY}
\E\left[\left\|\mu_{ P_{X|Z}}(\cdot) - \mu_{ P_{Y|Z}}(\cdot)\right\|_{\cH}^2\right] = \E\left[\left\|\E\left[\sfK(X,\cdot) - \sfK(Y,\cdot)\middle|Z\right]\right\|_{\cH}^2\right].
\end{align}

Now, the construction of $(X,X',Y,Y')$ gives,

\begin{align*}
\E\left[\left\|\E\left[\sfK(X,\cdot) - \sfK(Y,\cdot)\middle|Z\right]\right\|_{\cH}^2\right] = \E\left[\left\langle\E\left[\sfK(X,\cdot) - \sfK(Y,\cdot)\middle|Z\right],\E\left[\sfK(X',\cdot) - \sfK(Y',\cdot)\middle|Z\right]\right\rangle_{\cH}\right]
\end{align*}

Interchanging the inner product and conditional expectation gives,

\begin{align}\label{eq:gequiv}
\E\left[\left\|\E\left[\sfK(X,\cdot) - \sfK(Y,\cdot)\middle|Z\right]\right\|_{\cH}^2\right]
& = \E\left[\left\langle\sfK(X,\cdot) - \sfK(Y,\cdot), \sfK(X',\cdot) - \sfK(Y',\cdot)\right\rangle_{\cH}\right]\nonumber\\
& = \E\left[\sfK(X,X') + \sfK(Y,Y') - \sfK(X,Y') - \sfK(X',Y)\right].
\end{align}

This together with the equality in (<ref>) completes the proof of Proposition <ref>.

§ PROOF OF THEOREM <REF>

To begin with, for notational convenience, define

\begin{align}\label{eq:defTn}
T_n = T_n(\bcW_n,\sZ_n) := \EMMD^2[\sfK, \bcW_n, \sZ_n] = \frac{1}{n}\sum_{u=1}^{n}\frac{1}{K}\sum_{v \in N_{G(\sZ_n)}(u)}\sfH(\bW_{u},\bW_{v}) .
\end{align}

To show Theorem <ref> we will prove that $T_n$ converges to $\EMMD\left[\cF, P_{X|Z}, P_{Y|Z}\right]$ in $L^2$. The proof is based on the following 2 lemmas:

Under the assumptions of Theorem <ref>,

\begin{align*}
\E[T_n] \ra\E\left[\left\|g_{Z_{1}}\right\|_{\cH}^2\right] ,
\end{align*}

where $g_{Z_1}(\cdot) := \E\left[\sfK(\cdot, X_{1}) - \sfK(\cdot, Y_{1})| Z_{1}\right]$.

Under the assumptions of Theorem <ref>,

\begin{align}\label{eq:Tnvariance}
\mathrm{Var}[T_n] = o\left(n^{-\frac{\delta}{4 + \delta}}\right).
\end{align}

The proofs of Lemma <ref> and Lemma <ref> are given in Appendix <ref> and Appendix <ref>, respectively. Using these results the proof of Theorem <ref> can be completed as follows: Note that Lemma <ref> and Lemma <ref> together implies,

\begin{align}\label{eq:TgZ}
T_n\pto \E\left[\left\|g_{Z_{1}}\right\|_{\cH}^2\right].
\end{align}

Now, consider $(X,X',Y,Y',Z)$ such that $Z\sim P_Z$ and $(X,Y), (X',Y')$ are generated i.i.d. from $ P_{XY|Z}$. Then recalling (<ref>) gives,

\begin{align}\label{eq:EgEMMD}
\E\left[\left\|g_{Z_{1}}\right\|_{\cH}^2\right] = \E\left[\sfK(X,X') + \sfK(Y,Y') - \sfK(X,Y') - \sfK(X',Y)\right].
\end{align}

Collecting (<ref>), (<ref>), (<ref>), and Proposition <ref>, the result in Theorem <ref> follows.

§.§ Proof of Lemma <ref>

To begin with, note that for all $g_{Z_u} := \E\left[\sfK(\cdot, X_{u}) - \sfK(\cdot, Y_{u})| Z_{u}\right] \in \cH$ is well defined for all $1\leq u\leq n$ by <cit.>. Now, recalling the definitions of $T_n$ and $\sfH$ from (<ref>) and
(<ref>), respectively, and by exchangeability,

\begin{align}\label{eq:TnXY}
\E\left[T_{n}\right]
= & \E\left[\frac{1}{K}\sum_{u=1}^{n}\left\langle\sfK(\cdot, X_{1}) - \sfK(\cdot, Y_{1}), \sfK(\cdot, X_{u}) - \sfK(\cdot, Y_{u})\right\rangle_{\cH}\I\left\{(1,u)\in E(G(\sZ_n))\right\}\right]  \nonumber \\
= & \E\left[\frac{1}{K}\sum_{u=1}^{n} \E\left[ \left\langle\sfK(\cdot, X_{1}) - \sfK(\cdot, Y_{1}), \sfK(\cdot, X_{u}) - \sfK(\cdot, Y_{u}) \right\rangle_{\cH} | \cF(\sZ_n) \right] \I\left\{(1,u)\in E(G(\sZ_n))\right\}\right] \nonumber \\
= &\  \E\left[\frac{1}{K}\sum_{u=1}^{n} \langle g_{Z_{1}}, g_{Z_{u}}\rangle_{\cH} \I\left\{(1,u)\in E(G(\sZ_n))\right\}\right] ,
\end{align}

where the last step follows by interchanging the inner product and expectation. Now, suppose $N(1)$ is an index selected from the neighbors of vertex $1$ in $G(\sZ_n)$ uniformly at random. Then from (<ref>),

\begin{align}
\left|\E\left[T_{n}\right] - \E\left[\left\|g_{Z_{1}} \right\|_{\cH}^2\right]\right|
&\ \leq \E\left[\frac{1}{K}\sum_{u=1}^{n}\left|\langle g_{Z_{1}}, g_{Z_{u}}\rangle_{\cH} - \|g_{Z_{1}}\|_{\cH}^2\right|\I\left\{(1,u)\in E(G(\sZ_n))\right\}\right] \label{eq:gZH}\\
&\ =\E\left[\frac{1}{K}\sum_{u=1}^{n}\left|\langle g_{Z_{1}}, g_{Z_{u}} - g_{Z_{1}}\rangle_{\cH}\right|\I\left\{(1,u)\in E(G(\sZ_n))\right\}\right]\nonumber\\
&\ \leq \E\left[\|g_{Z_{1}} \|_{\cH}\|g_{Z_{N(1)}} - g_{Z_{1}} \|_{\cH}\right] \nonumber\\
&\ \leq\sqrt{\E\left[\|g_{Z_{1}}\|_{\cH}^2\right]}\sqrt{\E\left[\left\|g_{Z_{N(1)}} - g_{Z_{1}}\right\|_{\cH}^2\right]} \label{eq:Hbarlimitdiffbdd}
\end{align}

where the last 2 steps follows by the Cauchy-Schwarz inequality.  From <cit.>,

\begin{align}\label{eq:gZN1Z1}
\E\left[\left\|g_{Z_{N(1)}} - g_{Z_{1}}\right\|_{\cH}^4\right]\lesssim \E\left[\left\|g_{Z_{N(1)}}\right\|_{\cH}^4\right] + \E\left[\left\|g_{Z_{1}}\right\|_{\cH}^4\right]\lesssim \E\left[\left\|g_{Z_{1}}\right\|_{\cH}^4\right].
\end{align}

Now, recalling $g_{Z_{1}} = \E\left[\sfK(\cdot, X_{1}) - \sfK(\cdot, Y_{1})| Z_{1}\right]$, observe that

\begin{align}
\E \left[ \|g_{Z_{1}}\|_{\cH}^4 \right ] & = \E \left [ \left\|\E\left[\sfK(\cdot,X_{1}) - \sfK(\cdot, Y_{1})|Z_{1}\right]\right\|_{\cH}^4  \right] \nonumber \\
& \leq \E \left[ \E\left[\left\|\sfK(\cdot, X_{1}) - \sfK(\cdot, Y_{1})\right\|_{\cH}^4|Z_{1}\right] \right] \tag*{ (by Jensen's inequality [To and Yip, 1975]) } \nonumber \\
& \lesssim \E\left[\left\|\sfK(\cdot, X_{1})\right\|_{\cH}^4\right] + \E\left[\left\|\sfK(\cdot, Y_{1})\right\|_{\cH}^4\right] \nonumber \\
\label{eq:gZ1}  & = \E\left[\sfK(X_{1},X_{1})^2\right] + \E\left[\sfK(Y_{1},Y_{1})^2\right] <\infty ,
\end{align}

where the last equality follows by the reproducible property of $\cH$. Applying  (<ref>) in (<ref>) implies, $\|g_{Z_{N(1)}} - g_{Z_{1}}\|_{\cH}^2$ is uniformly integrable. Recalling that $K = o(n/\log n)$, by <cit.>,
$$\E\left[\left\|g_{Z_{N(1)}} - g_{Z_{1}}\right\|_{\cH}^2\right] = o(1).$$
The proof of Lemma <ref> is completed by recalling (<ref>).

§.§ Proof of Lemma <ref>

The establish the result in (<ref>) we will invoke the Efron-Stein inequality. For this, suppose $(\bW_1', Z_1'), (\bW_2', Z_2'), \ldots, (\bW_n', Z_n')$, where $\bm W_i' = (X_i', Y_i')$, for $1 \leq i \leq n$, are generated i.i.d. from the joint distribution $P_{XYZ}$, independent of the data $(\bW_{1},Z_{1}),\ldots, (\bW_{n},Z_{n})$. Define a new collection

\begin{align}\label{eq:WZs}
\bcW_{n, s} := \left(\bcW_n \setminus\{\bW_{s}\}\right)\bigcup\{ \bW_{s}'\}\text{ and }\cZ_{n,s} := \left(\sZ_n\setminus\{Z_{s}\}\right)\bigcup\{Z_{s}'\} ,
\end{align}

where we replace the $s$-th data point with an independent copy. (Recall that $\bcW_n:= \{\bW_i = (X_i,Y_i):1\leq i\leq n \}$ and $\sZ_n := \{Z_1,Z_2,\ldots,Z_n\}$.) The Efron-Stein inequality [Efron and Stein, 1981] (see also <cit.>) then implies,

\begin{align}\label{eq:EfromStein}
\mathrm{Var}[T_n] & \leq \sum_{s=1}^{n}\E\left[\left(T_n - T_{n, s}\right)^2\right] ,
\end{align}

where, for $1 \leq s \leq n$, $T_{n,s} := T_n(\bcW_{n,s},\cZ_{n,s})$ is the test statistic as defined in (<ref>) with $\bcW_{n}$ and $\cZ_{n}$ replaced by $\bcW_{n,s}$ and $\cZ_{n,s}$, respectively.

We now proceed to bound the individual terms in the sum in (<ref>). To this end, recall that $G(\sZ_n)$ is the directed $K$-NN graph constructed using $\cZ_{n}$. For $1 \leq s \leq n$, denote by $G(\sZ_{n, s})$ the directed $K$-NN graph constructed using $\cZ_{n, s}$. We will now decompose $T_n-T_{n,s}$ by filtering out the common terms. For this, we need the following definitions. Throughout we fix $1 \leq s \leq n$.

* Denote by $\cA_s$ the set of all ordered pairs $(u, v)$, for $1 \leq u \ne v \leq n$, such that either $u=s$ or $v=s$.

* Now, define

\begin{align}\label{eq:defCns}
\cC_{n, s} := \left(E(G(\sZ_n))\bigcap E(G(\sZ_{n, s}))\right)\setminus \cA_{s} ,
\end{align}

which is the collection of edges that are common between $G(\sZ_n)$ and $G(\sZ_{n, s})$ such that none of the end points of the edges are the vertex $s$.

* Also, for $1 \leq v \leq n$, define

\begin{align*}
E_{G(\sZ_n)}(v):= \left\{(v,u): u\in N_{G(\sZ_n)}(v)\right\},
\end{align*}

to be the collection of directed edges in $G(\sZ_n)$ with $Z_v$ as the starting point.

* Finally, let

\begin{align*}
V_{G(\sZ_n)}:= \left\{v\in [n]: E_{G(\sZ_n)}(v)\bigcap\left(E(G(\sZ_n))\bigcap\left(\cC_{n, s} \right)^{c}\right)\neq \emptyset \right\}
\end{align*}

be the indices which have at least one incident edge in $G(\sZ_n)$ that belongs to $\cC_{n,s}^c$. Similarly, for $G(\sZ_{n, s})$, consider

\begin{align*}
V_{G(\sZ_{n, s})}:= \left\{v\in [n]: E_{G(\sZ_{n, s})}(v)\bigcap\left(E(G(\sZ_{n, s}))\bigcap\left(\cC_{n,s}\right)^{c}\right)\neq \emptyset \right\}.
\end{align*}

We now rewrite the difference between $T_n$ and $T_{n,s}$ as follows:

\begin{align*}
T_n - T_{n, s} = \frac{1}{nK}\left[\sum_{u=1}^{n}\sum_{v \in N_{G(\sZ_n)}(u)}\sfH(\bW_{u}, \bW_{v}) - \sum_{u=1}^{n}\sum_{v\in N_{G(\sZ_{n, s})}(u)}\sfH\left(\bW_{u,s},\bW_{v,s}\right)\right] ,
\end{align*}

\begin{align}\label{eq:newnotation}
\bW_{u,s}=
\begin{cases}
\bW_u' & \text{if }u=s\\
\bW_u & \text{otherwise},
\end{cases} \text{ and }
\begin{cases}
Z_u' & \text{if }u=s\\
Z_u & \text{otherwise}.
\end{cases}
\end{align}

Note that the edges in $\cC_{n,s}$ (recall (<ref>)) are not effected when $Z_{s}$ is replaced by $Z_{s}'$. Thus,

\begin{align}\label{eq:diff1stTn}
& T_n - T_{n, s} \nonumber \\
& = \frac{1}{nK} \left[\sum_{u=1}^{n}\sum_{\substack{ v \in N_{G(\sZ_n)}(u)\\ (u, v)\in E(G(\sZ_n))\setminus\cC_{n, s} }} \sfH(\bW_{u},\bW_{v}) - \sum_{u=1}^{n} \sum_{\substack{ v \in N_{G(\sZ_{n, s})}(u) \\ (u, v) \in E(G(\sZ_{n, s}) )\setminus\cC_{n, s} } } \sfH\left(\bW_{u, s}, \bW_{v, s} \right) \right] .
\end{align}

Next, observe that if $u \in V_{G(\sZ_n)}\bigcap V_{G(\sZ_{n, s})}$ is such that $u\neq s$, then the following hold:

$$u\in \left(V_{G(\sZ_n)}\right)^{c} \implies E_{G(\sZ_n)}(u)\subseteq \cC_{n,s} \text{ and } u\in \left(V_{G(\sZ_{n, s})}\right)^{c} \implies E_{G(\sZ_{n, s})}(u)\subseteq \cC_{n,s}.$$
In other words, vertices in $\left(V_{G(\sZ_n)}\right)^c$ and $\left(V_{G(\sZ_{n, s})}\right)^c$ contribute all their edges to $\cC_{n,s}$ and, hence, have no contribution to the difference in (<ref>). Thus,[
For two sequences $\{a_n\}_{ n \geq 1 }$ and $\{b_n\}_{ n \geq 1 }$, denote $a_n \lesssim_{\square} b_n$ to mean $a_n \leq C(\square) b_n$ for $n$ large enough, where $C(\square) > 0$ is a constant that depends on the subscripted parameters. ]

\begin{align}\label{eq:diff2ndTn}
& T_n - T_{n, s} \nonumber \\
& = \frac{1}{nK}\left[\sum_{u\in V_{G(\sZ_n)}}\sum_{\substack{v\in N_{G(\sZ_n)}(u)\\(u,v)\in E(G(\sZ_n))\setminus\cC_{n,s}}}\sfH(\bW_{u},\bW_{v}) - \sum_{u\in V_{G(\sZ_{n, s})}}\sum_{\substack{v\in N_{G(\sZ_{n, s})}(u)\\(u,v)\in E(G(\sZ_{n, s}))\setminus\cC_{n,s}}}\sfH\left(\bW_{u,s},\bW_{v,s}\right)\right] \nonumber \\
& \leq \frac{1}{nK}\left[\left|\sum_{u\in V_{G(\sZ_n)}}\sum_{\substack{v\in N_{G(\sZ_n)}(u)\\(u,v)\in E(G(\sZ_n))\setminus\cC_{n,s}}}\sfH(\bW_{u},\bW_{v})\right| + \left|\sum_{u\in V_{G(\sZ_{n, s})}}\sum_{\substack{v\in N_{G(\sZ_{n, s})}(u)\\(u,v)\in E(G(\sZ_{n, s}))\setminus\cC_{n,s}}}\sfH\left(\bW_{u,s},\bW_{v,s}\right)\right|\right] \nonumber \\
& \lesssim_{d} \frac{1}{n}\left[\max_{1\leq u\neq v\leq n}\left|\sfH(\bW_{u},\bW_v)\right| + \max_{1\leq u\neq v\leq n}\left|\sfH\left(\bW_{u,s},\bW_{v,s}\right)\right|\right] ,
\end{align}

where the last step uses the bounds:

\begin{align}\label{eq:ECn}
\sum_{u\in V_{G(\sZ_n)}}\left|v\in N_{G(\sZ_n)}(u): (u,v)\in E(G(\sZ_n))\setminus\cC_{n,s}\right|
\leq 2\left|E(G(\sZ_n))\setminus\cC_{n,s}\right| \lesssim_d K,
\end{align}

and, similarly,

\begin{align}\label{eq:ECns}
\sum_{u\in V_{G(\sZ_{n, s})}}\left|v\in N_{G(\sZ_{n, s})}(u): (u,v)\in E(G(\sZ_{n, s}))\setminus\cC_{n,s}\right|\leq 2\left|E(G(\sZ_{n, s}))\setminus\cC_{n,s}\right| \lesssim_d K .
\end{align}

The bound in (<ref>) (and similarly (<ref>)) follows by noting that $E(G(\sZ_n))\setminus\cC_{n,s}$ contains the edges in $G(\sZ_n)$ which are either absent in $G(\sZ_{n,s})$ or have $s$ as one of the endpoints; and from the fact that the  maximum total degree (sum of the out-degree and the in-degree) in a $K$-NN graph is bounded by $c(d) K$, for some constant $c(d) > 0$ depending on the dimension $d$ (see, for example, <cit.>).

Combining (<ref>) and the bound from (<ref>) gives,

\begin{align}
\mathrm{Var}[T_n]
& \leq \sum_{s=1}^{n}\E\left[\left(T_n - T_{n, s}\right)^2\right]\lesssim_{d}\frac{1}{n}\E\left[\left(\max_{1\leq u\neq v\leq n}\left|\sfH(\bW_u,\bW_v)\right|\right)^2\right] .
\label{eq:EfronSteinbd}
\end{align}

Now recalling the definition of $\sfH$ from (<ref>), note that,

\begin{align*}
\sfH(\bw,\bw') = \left\langle\sfK(x,\cdot) - \sfK(y,\cdot),\sfK(x',\cdot) - \sfK(y,\cdot)\right\rangle_{\cH}.
\end{align*}

Hence, by the Cauchy-Schwartz inequality,

\begin{align*}
\mathrm{Var}[T_n]  \lesssim_{d} \frac{1}{n}\E\left[\left(\max_{1\leq u\leq n}\left\|\sfK(\cdot, X_{u}) - \sfK(\cdot, Y_{u})\right\|_{\cH}^2\right)^2\right]
\end{align*}

Finally, using the triangle inequality and the reproducing property of $\cH$,

\begin{align}\label{eq:TnETndiffbdd}