MMD Two-sample Testing in the Presence of Arbitrarily Missing Data
Authors: Yijin Zeng
Niall M. Adams
Dean A. Bodenham
Venue: Accepted by TMLR ()

Abstract: In many real-world applications, it is common that a proportion of the data may be missing or only partially observed. We develop a novel two-sample testing method based on the Maximum Mean Discrepancy (MMD) which accounts for missing data in both samples, without making assumptions about the missingness mechanism. Our approach is based on deriving the mathematically precise bounds of the MMD test statistic after accounting for all possible missing values. To the best of our knowledge, it is the only two-sample testing method that is guaranteed to control the Type I error for both univariate and multivariate data where data may be arbitrarily missing. Simulation results show that the method has good statistical power, typically for cases where 5% to 10% of the data are missing. We highlight the value of this approach when the data are missing not at random, a context in which either ignoring the missing values or using common imputation methods may not control the Type I error.

Published in Transactions on Machine Learning Research (11/2025)
MMD Two-sample Testing in the Presence of Arbitrarily Missing Data
Yijin Zeng Imperial College London
pfjin..eng20@@mmprriilccc..
Niall Adams Imperial College London
TL. ad.aaams Wwwperaal..aa...m.mm
Dean Bodenham Imperial College London
dean. Sewrdennha@@@ippriiiaaaa
Reviewed OpenReview: https: httpe://openrevieel openreview. net/ forum? id=CfcDellIC5
Abstract
In many real-world applications, it is common that proportion of the data may be missing only partially observed. We develop a novel two-sample testing method based on the Maximum Mean Discrepancy (MMD) which accounts for missing data in both samples, without making assumptions about the missingness mechanism. Our approach is based on deriving the mathematically precise bounds of the MMD test statistic after accounting for all possible missing values. To the best of our knowledge, is the only two-sample testing method that is guaranteed control the Type I error for both univariate and multivariate data where data may be arbitrarily missing. Simulation results show that the method has good statistical power, typically for cases where 5% to 10% of the data are missing, We highlight the value of this approach when the data are missing not at random, a context in which hit ignoring the missing values or using common imputation methods may not control the Type I crror.
1 Introduction
Two-sample hypothesis lesting a fundamental statistical method used to determine two samples of data are different. Numerous two-sample testing methods are available, including Student's 's t-test (Student 1908), the Wilcoxon-Mann-Wien
U test (Mann & Whitney 1947), the Kolmogorov-Smmnne
test (Kolmogorov 1933
the Energy Distance [Rizzo & Székely 2016), and the Maximum Mean Discrepancy (MMD) test (Gretton
et al. 2012a). These methods have proven useful in various fields such as medicine (Tsiatis et al. 2008
O'Brien & Fleming, 1979), finance (Lukić & Milošević 2024 Dave & Korkie 1981), psychology (Pfister &
Janczyk 2013) and machine learning (Gretton et al.) 2012a 2009 Schrab et al. 2023 Nearly all two-sample
testing methods are designed solely for data that are fully observed, with few exceptions discussed
the next section. However, in many real-world datascts, a subset of univariate values may be missing or multivariate values may only be partially observed. For example, in clinical trials it is common that patients drop out during the trial (Bakris et al. 2015), resulting in missing data for those participants.
When data are missing, common practices are either to ignore all missing values or impute these values using some imputation schemes. Following this step, the data are treated as complete, allowing any standard two- sample testing method to be used. However, except in special cases, such as when the missing data are missing completely at random (Rubin| 1976), these practices are often invalid as they risk incrcasing the probability of Type error occurring. Under certain assumptions, such as when the data are missing completely random, missing at random, when the data missing not at random but the missingness mechanisms (the

---

Published in Transactions on Machine Learning Research (11/2025)
samples can be subject to different missingness mechanisms) can be explicitly specified (Schafer & Graham 2002), then certain sophisticated missing data methods exist, such as the expectation-mxxiiiae algorithm and the multiple imputation method (Schafer 1999). However, without these assumptions, relying on these methods is fraught with risk.
Recently a two-sample testing method for data with univariate missing values was proposed which does not rely on a specific imputation procedure (Zeng et al. 2024). This method, based the Wilcoxon-Mann- Soin--hhitte test, rejects the mull hypothesis after accounting for all possible values of the missing data by deriving bounds for the U statistic. This approach avoids ignoring or imputing the missing data and is shown to control the Type I error without making assumptions about the missingness mechanisms, while also having good statistical power when the proportion of missing data is around 10%. However, this approach is restricted to univariate samples and can only detect location shifts since it is based on the Wilcoxon- Mamtoo-----n---- Whitney test.
In this paper, we propose a novel two-sample testing method basod on the MMD test in the presence of missing data that makes no assumption about the missingness mechanisms, following the approach in [Zeng] et al. 2024h. However, our approach can be used for both univariate and multivariate samples, and can detect any distributional change since it is based on the MMD test.
Our approach, named MMD-Miss, is based on deriving the bounds of the unbiased MMD test statistic (Gretton et al. 2012a) in the presence of missing data. To do so, we use the Laplacian kernel, a popular choice of characteristic kernel that enables MMD to detect any distributional shift (Sriperumbudur et al. 2011). To
compute the p-value for the test, we either use Monte Carlo sampling from the permutation distribution of the observed data or, in the case of high-dimensional multivariate data, we use normality approximation recently derived in (Gao & Shao, 2023). We prove that our approach controls the Type crror, regardless of the values of the missing dala. Numerical simulations are presented, showing that our met has good testing power when the proportion of missing data is around 5% to 10%, while controlling the Type I error.
Motivation and contribution. Our work is motivated to ensure valid two-sample testing in the presence of missing data without relying on missingness assumptions. In the context of performing two-sample testing with missing data, a common misunderstanding is that when the proportion of missing data is small, e.g. 5% as suggested by Schafer (1999), and Heymans& & Twisk (2022), the testing result will not be skewed by missing data after either ignoring or imputing these values. While these practices may be justified under certain conditions, such as when the data are missing completely at random (Rubin, 1976), we provide experiments in Figure 1 showing that when the data are missing not at random, these practices can lead to a Type I crror asymplotically equals to even with a small proportion (5%) of missing daía. Our contribution is method based on the MMD two-sample test that is suitable for univariate and multivariate data with missing values. In contrast to other methods, this approach can detect any change in distribution while being guaranteed to control the Type I error, and without making any assumptions about the missing data.
2 Background
This section provides the necessary background for our work. We begin by introducing the two-sample testing problem and then review the different types of missingness mechanisms. Next, we briefly review the MMD two-sample test, and then discuss two standard approaches for implementing the MMD test: the permutation test and the normal approximation. Finally, we highlight related work with a focus on existing approaches to two-sample testing in the presence of missing data as well as extensions of the MMD test.
Two-Sample Testing. A two-sample testing method is used to determine whether two groups of data are statistically significantly different. Consider two samples of obscrvations X = (21,..., and Y = ".....,,,,,,, where each observation is in Rd, ; and suppose that X, and Y are independent random samples from distributions p. and 9, respectively. We define the null hypothesis Ho as p = and the alternative hypothesis H, as p q. The two-sample testing method first computes a statistic basod on the data, and then a p-value based on the statistic. After comparing the p-value to a pre-specified significance threshold o, the null hypothesis Ho is either rejected or fails to be rejected.
2

---

Published in Transactions on Machine Learning Research (11/2025)
A Type I error occurs when the null hypothesis Hg is true but the test incorrectly rejects it. Conversely, Type IT error occurs when Ho is false but the test fails to reject it. A two-sample testing method is usually derived so that, under Ho, the Type I error is not larger than the pre-specified significance threshold a E (0, 1). In this case, we say the test controls the Type crror. For a given a, a two-sample testing method is preferred if it has lower probability of making a Type II error, denoted by 9. The power is defined as - 8, and a two-sample testing method is more desirable if it demonstrates higher power. Two-sample testing mothods are usually assessed on their power, given that the Type I error controlled.
Missingness Mechanisms. Missingness mechanisms are typically classified missing completely at random (MCAR), missing at random (MAR), and missing not at random (MNAR) following Rubin (976)).
Lel Zrom = (.i,...; =²)²² denote the complete data, where each 2 is a R real value. The complete data Zoom can be split into observed parts Zobs and missing parts Zmiss: Let R be a binary indicator matrix matching the dimensions of Zeom, with elements indicating observed (1) or missing (0) values. The core idea of Rubin (1766) to treat R a realization of a random vector following an unknown distribution.
The mechanism is MCAR if P PRR Zeom) = PR)), meaning R is independent of data values, justifying the practice of ignoring missing data. The mechanism is MAR if P (R])eom)= = P Zob) where R depends only on observed data, making certain imputation methods viable. If neither condition holds, the mechanism is MNAR. For known missingness mechanisms, specific imputations might be justified (Rubin 1976 Lachin 1999). However, with unknown mechanisms, ignoring or imputing data can lead to invalid results, as shown in our simulations in Section 5 Notably, our method does not depend on the missing data values, ensuring valid results regardless of the missingness mechanism. This makes our method particularly suitable for the MNAR case.
Maximum Mean Discrepancy. The Maximum Mean Discrepancy (MMD) is a two-sample test used for deciding if two groups of independent observations X = {ZI...., In, and Y = 1y1,,...,nm. are sampled from different distributions p and q. Its test statistic is a kernel-based measure of distance which compares and in a reproducing kernel Hilbert space (Aronszaju, 1950). More formally, let Hi denote a reproducing kernel Hilbert space with kernel function k, the MMD between p and q is defined as (Gretton et al. 2012a) the following statistic (often written as MMD²):
the mean embeddings
( sup
( (Ex= Ex-,(((())) 2
AMDD"H4..P.O/
(
In the context of two-sample testing, an unbiasod estimate of MMD p. is:
R1
n2
ngg
n2
12
n2
MMD = W1 (111111
kz11,1z) + n2(22-1)
k vz)
k( B1).
12
J+
The MMD test. rejects the null hypothesis when MMD Y) exceeds certain threshold, usually chosen to control the probability of Type I error under a pre-specified number (0, 1). It is shown in /Greton et al. (2012a) that when & is characteristic kernel (Sriperumbudur al. 2011 1, MMD? P. q) 0 if and only if p 4, otherwise MMD"(H4, P, 4] > O. Hence, when characteristic kernel is used, the MMD test is able to detect any distribution changes. However, on finite sample sizes, the choices of kernel have significant. impact on the power of MMD test (Biggs et al. 2023). The Laplacian and Gaussian kernels are defined as
ky (x, v) = exp( 3 |- vlli), and ke (a v) = exp(->llx =112)
respectively, where B > 0 and 0 are hyperparameters.
The Laplacian and Gaussian kernels are popular choices for the MMD Gao & Shao, 2023 Sriperumbudur et al. 2011), ; mainly because they are both characteristic kernels (Sriperumbudur et al. 2010 2011). There are other characteristic kernels, such as the Matérn family of kernels (Sriperumbudur et al. 2010). In the
univariate case, the MMD test has been shown perform similarly for the Laplacian and Gaussian kernels [Bodenham & Kawahara 2023
3

---

Published in Transactions on Machine Learning Research (11/2025)
In this work we only use the Laplacian kernel ky because it allows to derive appropriate bounds for our method. Using different kernel with this approach may be possible but would require the derivation of different bounds. While this may seem restrictive, we emphasise that the Laplacian kernel is characteristic and so can be used with the MMD to detect any change in distribution. In the following, we use k(x, y) to denote the Laplacian kernel without causing ambiguity.
Strengths and limitations of the MMD. As discussed, the core principle of MMD is to detect differences between two distributions by measuring the largest difference in expectations over rich class of functions in reproducing kernel Hilbert space. This makes the MMD capable of capturing a wide range of distributional discrepancics, including subtle differences in higher-order moments or tail behaviour, which classical tests such as Sartoon--MMn-WWhhtneee lest may overlook. The MMD is compared to popular univariate testing methods in Bodenham & Kawahara (2023). Furthermore, Since the MMD is based on a kernel, it can be applied to high-dimensional data naturally, while the Wilcoxon- Whitney and Kolmogorov-Smint Smirnoo tests are restricted to univariate data. However, a limitation of the MMD test is that its statistical performance depends on the choice of the kernel and the corresponding parameter value for the kernel. Several mot are proposed for dealing with these issues, which are discussed in the related work (Gretton et al. 2012b Biggs el al. 2023). In this work restrict the MMD kernel to the Laplacian kernel, which is characleristic and so can be used to detect any change in distribution (Sriperumbudur et al. 2010 [001]), and set the kernel paramcter using the median hcuristic, which has been shown to generally perform well (Gretton et al. 2012a
Bodenham & Kawahara 2023).
Permutation Test. The permutation test common numerical procodure for determining the statistical significance of computed two-sample test statistic. It relies on the principle that if the null hypothesis of equal distribution is true, the observed values are exchangeable, and so all permutations are equally likely Good 2005). The permutation test is non-parametric, meaning that it does not assume the distribution of the data. Furthermore, Schrab et al. (2023) has shown that this test produces the p-value non-asymptotically at most the significance level a under the null hypothesis. In contrast, the testing approach based on normal
approximations (Gao & Shao 2023 Shekhar et al. 2022 are only asymptotic true.
To be more specific, let X = {=1..... Fn,) } and Y .. No2): where each sample is Rd real value. The permutation test. begins with uniformly sampling B i.i.d. permutations from (1,..., A1 + n2)- Let us denote (00 (n1 + n2)), i= = L,..., B. Subsequently, define (21,..., 2n,) = (I1,..., )n,) and (2n, (1m,+1...... Amitnol (V1.... Vna), and for any i = 1,..., B, denote
the B samplings as (o(1), ....,, ,.)),,(())),, where each (i))
X =
(000 }, You = { (2 I): (m))))))))))
The test threshold of MMD test, using MMD Y) test statistic with significance level 0 < 1, then computed as the [a(B + 1)]-th largest numbers in the set
(X, U (MMD) (X,
If MMD Y) exceeds this threshold, then MMD test rejects the null hypothesis. Otherwise, the null hypothesis is not rejected.
Normal Approximation. It is well known that for a fixed dimension d, the asymptotic distribution of MMD2 Y) under the null hypothesis (i.e. p = 9) takes the form of an infinite weighted sum of x² random variables with the weights depending on p (Gretton et al., 2012a), which are normally unknown. An estimation of the asymplotic distribution of MMD is established in Gretton et al. (2009). This estimation is proved asymptotically correct regardless of the distribution p under certain conditions. However, empirical experiments indicate this estimation is rather conscrvative.
It was recently shown by Gao & Shao (2023) that the studentized MMD Y) given by equation g convergences to normal distribution when all n, m, + co. This result holds for a wide range of kernels, including Laplacian and Gaussian kernels. This normal approximation is shown empirically to be rather accurate Gaa & Shao 2023) even for values of relatively small n, m, d, e.g. n m 25, d 50, and a higher approximation accuracy is shown for larger values of n, m. d.
4

---

Lachin (1999) proposes a test method based on Published in Transactions on Machine Learning Research (11/2025)
Related Work.
Various approaches have been proposed to mitigate the issue of missing data for two-sample
hypothesis testing.
When the missing data are missing at random, rank-based two-sample testing met
(Cheung] 2005 Chen et al. 2013) have been proposed which control the Type 1 error and achieve good
Paired two-sample testing with missing data are studied in Martinez-Camblor et al. (2013). 2018), under the assumption that data are missing completely at random and at least. one of the paired samples is observed. When the samples are bounded within an interval, Pan (2000) proposod an algorithm based on multiple imputation (Little & Rubin 2019) by employing approximate Bayesian bootstrap (Rubin & Schenker [986). Shutoh et al. (2010) proposed method for testing the shift in the mean of high-dimensional normal data assuming the two samples share the same missingness mechanism. By considering certain unobserved samples as the "worst case",
testing power.
and Fong et al.
rank called the worst rank test.
All methods above rely on different assumptions about the missing data. One exception of existing methods which makes no missing data assumption, as we mentioned before, is the rank based two-sample testing method proposed in Zeng et. al. (2024). Their method rejects the null hypothesis after accounting for all possible Malcoxon-Manni Whitney test statistics by deriving the bounds of the test statistics. Our method in this work follows the same approach as in Zeng et al. (2024). We make no assumption of the missing data, except assuming they are distinct. The bounds of the MMD test statistic using the Laplacian kernel are derived, and the null hypothesis is rejected when all possible test statistics are significant. The approach in Zeng et al. (2024) is restricted to univariate data for testing a location shift, since it relies on the Wilcoxon- Mann- W Whitney test, while our approach can be used for both univariate and multivariate samples, and can detect any distributional change since it is based on the MMD test.
Our method is deliberately conservative, contrasting to the dominant methods mentioned above, which seek to boost the testing power when certain assumptions of the missing data can be made. Our concern here is whether the testing result is significant when no missing data assumption is considered appropriate. One related concept is known as partial identification (Manskil 2003), which argues that it is not sufficient to know what inferences can be made when certain assumplions are employed. It highlights the importance of making inferences when no assumption is imposed for drawing robust conclusions.
The concern that set of plausible missing data assumptions might lead to different conclusions often motivates the use of sensitivily analysis (Thabane et al. 2013 Mathur 2023 where these assumptions are evaluated to determine whether a consensus can be reached about the testing conclusions. If all evaluated assumptions lead to a consistent conclusion, then the conclusion is considered to be robust and is taken as the final result, as demonstrated by Bakris et al. (2015) Our method can be viewed as a sensitivity analysis that incorporates all possible assumptions. Tf significant testing result is determined by our method, then the result is significant regardless of the values of the missing data, and hence is robust to all possible assumptions on the missing data.
Various authors have proposed extensions of the MMD two-sample test for several different purposes. /retton ct al. (2012b) and Biggs et al. 2023) investigated strategies for optimal kernel selection based on the obscrved data, while Schrab el al. (2023) locused on selecting an appropriate bandwidth for the chosen kernel. To handle complex data modalities such as images, audio, and video in two sample testing, Kirchler et al. (2020) proposed a deep learning-based extension of the MMD test. As we will discuss below, a common implementation of the MMD lest relies on a computationally intensive procedure known as permutation testing to approximate the null distribution and compute p-value. To alleviate the computational burden associated with the permutation procedure, several works (Shekhar et al. 2022 Gao & Shao 2023) have developed permutation-free approaches that significantly reduce computation. Our research in this work explores a different direction by studying MMD-based two-sample testing in the presence of missing data. We propose tests that control the Type I error without requiring assumptions about the underlying missingness mechanism.
3 Bounding MMD with Missing Data
In this section, we provide bounds for the MMD test statistic in the presence of missing data. As mentioned in the Section T. the key idea of the proposed method is to consider all possible test statistics in the presence
5

---

Published in Transactions on Machine Learning Research (11/2025)
of missing data and reject the mull hypothesis only when all test statistics significant. To compute these bounds, we use the Laplacian kernel, which is characteristic kernel (Sriperumbudur et al. 2010 2011 and able to detect any change in distribution.
Univariate Data. If a subset of data in X and Y are unobserved, the MMD2(X; Y) test statistic cannot be directly computed. Howover, the test statistic can be bounded given the missing data. Notice that for any two given real values x and y,, the Laplacian kernel k(x,, y) - exp (-3(x- v/) is bounded by (0, 1] since 3 > and |x- 15-38 y/ > O. Hence, a straightforward way to obtain a lower bound of MMD Y) is to let k(z,, 23) if at least one value of x; and zj is unobserved. Similarly, if at least one value of y and yy is unobserved, let k(y,, vi) = O. On the other hand, to minimize -k yj), we can let -k(2;, My -1 if either 2; yy is unobserved. The upper bound of MMD (X Y) can be obtained following similar manner. This simplistic method, while providing a theoretically correct bound of MMD Y), might yield bounds that are too conservative for effective two-sample testing in practice.
In the following, we construct tighter bounds of MMD Y) order to make useful for the two-sample testing problem. One starts by decomposing MMD (X, Y) into terms which have either none, one or both arguments of the kernel function missing. Without loss of generality, let us assume I1, : 11,..... -yma are not observed, or not fully observed, while Im, +1,..,,. and Um.tl...... Ynz are fully observed. Then, MMD (X,Y) Y) can be decomposed, using the following lemma:
Lemma 1. Suppose X = fx1...., Feyy , Y = fyr,..., &ngt are subsets of Rdd where d 2 1. Suppose that (M, Y) is unbiased MMD test slatistic defined
(X (1
k(rs, 25) m2
k(y: 11
k(ai,1)),
(1)
33
with lo denoting the Laplacian kernel. Then, for any two positive integers mi, m2 such that mi S ni, m2 S ni, MMD Y) can be rewritten as:
(X, = A1 + Ag + As,
2 where 01 (n:--11) C2 m2 1)
2 and
C3
Ai 12222
c Σ
Σ
k(I,,2))
+ 02
k(Ui, 81) -
k( Uj):
1=1 j=
2m1
<=1 j=l
B1-1
R2-1
72
n22
A = 4 Σ
Σ
(1,,115)+ 0
Σ
k(y Hj) C3 Σ
Σ
j=m 1=
i= 1 j=i+1
t= + j=ma+1
A3 =
k(z,, C3
k( Mi);
1=1 1=m1+1
J=
A4 = f=m2+1
k(,,,,))
C3
k(r,,1))
i=m1
The proof of this lemma can be found in Appendix A.1
Lemma I decomposes MMD Y) into four parts: the first part A₁ including only the unobserved data in X and Y; the second part A2 including only the fully observed data in X and Y; the third and the fourth terms Ag and A4 involving data where one sample (either X or Y) contains only unobserved data while the other sample includes only fully observed data.
In order to bound MMD (X,Y), Y), we propose to bound the four terms A,, A2, Ag and A4 separately. Given that Ag includes the fully observed data only, it is thus sufficient to bound A Ag and Aq. For the first term A1, where all data are unobserved, not fully observed, it can be seen that > k(x, v) > 0 for any x, R.
6

---

Published in Transactions on Machine Learning Research (11/2025)
Hence, have mi (mi - 1)/2 Di= Ej=++1 k( 22) > O, m2(m2 1)/2 E Σ Sis+1 k(ys: My)
0,
and m2ma 2 Di=1 1j=1 (r, yj) Then according to the definition of A in Lemma we have
mi(m1 1) m2(m2 1) > A1 11
m2
(2)
(
(2 1)
The challenge of bounding (X, Y) is provide tight bounds for terms Ag and A4, where data can be fully observed in one sample and missing, partially observed in another sample. We focus the term A3 because As can be bounded similarly. The key to bounding Ag is by first rewriting Ag as
A3
C1
Σ
- c3
Σ j=m2+1
k( y/) )
(3)
i=1
j=m +1
2 where C1 II] and
2
Define (z) as
a Σ
112
(1 - c
k(2,0))- -c3 Cs
Σ
kzz,
j=m++1
j=m=+1
Then, we have
A3 Ti
i.e. Ag is the summation of the function values (z) applied to all unobserved samples in X. This equation suggests that in order to bound A3. it is sufficient to bound the values of Ti (z). Therefore, we proceed by bounding (z), using the following result for univariate data: Lemma 2. Let T1...., Be, and y'...., be real values, which are observed. Suppose 0t, bi...., big and B are positive constants. Define
T(2) Σ a, exp(p=Blas 2/) Σ b, exp(-3ly1 2l)
(4)
function of TR. Subsequently, for any given 20 TR, we have
T(z)) 2 min T(r)), (, T(y))),(())))........ T((re))
and
T(zo) S max ((y)),,)).) T(ye,)) Moreover, these bounds for T(zu) are tight, and we have equality when zu E ( Fty ......, Wey}
20
to0.
This lemma is proved in Appendix by computing all the values in the set (T(=) 2 E (#1,..., Fe, Vi,...,yt, and identifying the maximum and minimum values.
A.2 It provides linear time algorithm for computing the bounds of T(z)
After proving Lemma we can provide the bounds of MMD (X,Y) in the presence of missing data, without making missing data. assumptions.
Theorem 1. Suppose X ( Iny} } and V = (Ji,..., Ung} are univoriate real values. Assume
T Tleerey Tmy 01...., Uma are unobserved. Let k denote the Laplacian kernel and define
a Σ j=mi+1
R2
Ti(2) =
k(2,0,) - C3
Σ
k(z, 9;),
j=m2+1
T2(2) = C2 k(2,0,)) - C3
Σ
k(2, 9;),
j=
j=my+1

---

Published in Transactions on Machine Learning Research (11/2025)
where C1 "1(m1-1)1
2
2 2 and C3 min
Further, let
m2(
S {0, Ti (Imi +1), (TTn)) ). Ti (1mm (1mm-1)) ...... (yn2))
S = {0, T2 (mm++))....... T)( T2 (1mm-1))...... (2(Un,)
Then, (X,) defined in equation using Laplacian kernel k is bounded, over all possible values of Fmi y1- .um2 0s follows:
21
m1 1))
m2 1)
+ my max S1 + m2 max S2 + A2 2 MMD (X,Y),
n1 (n1 - 1) n2 - 1)
2
(X. A2 + ma min S1 + m2 min S2
m 72
R1
where A2 is defined in lemma 7
The proof of this theorem is provided in Appendix A.5
Multivariate Data. We now provide the bounds of MMD for multivariate data, where X x; (d));
(a1,..., .nn }. Y (yI...., |nu and X,Y E Rd with d> 1. For each 22 in X, denote 2; = (1) where x1(() E IR is the value the Ith component of o;. Similarly, for each y: in Y, denote y/i (1) 99 (d)). Note that X and Y are assumed to have the same dimension d.
As in the case for univariate data, one starts by decomposing MMD (X, u (X,Y) Y) into A1,A2, As and A4 using Lemma [ Term A2 only includes data that are fully observed. Thus it can be computed directly. Term A can be bounded using simple inequalities, shown in equation 23 in the Appendix. The challenge is to provide tight bounds for term Ag and A4-
To bound term A3 and A4, we first provide notation for which components of a d-dimensional observation are missing and define what it means to impute such a value.
Definition 1. For a value = = that are missing; in other words {z(i) : j U2} are missing.
E TRd, let U2 (1,,, ... denote the set of components of
Definition 2. For a value 2 E Rd, that has missing components, x" is called on imputation of z if 2" E Rd is fully observed and z" (j) = = for all j E .I. U..
We can now exlend Lemma to the d-dimensional case:
Lemma 3. Let .1,..., Ft, ......, E Rdd be values that are fully observed. Suppose 01,...,;4, are positive constants. For = = (z(1), .,,),,,..,,.00) 2(d)) E Idd with mássing components, define
u8, bi,.... beg B
T({(() J :.
Σ a; exp
-B Σ
xx -=:5)=(001
2(1))
Σ
by exp
-B Σ (vs(5) 2(i)
jEU.
==1
JEU,
as a function of the unobserved components of a and let
X= x= : (T((())):)EE j.EU))) =(1) (.. (1), ............. (i),, U.))
Then, for any possible imputation of z, m0n min X3 s <T((>(() j U.}) < maxx), max X}
Lemma 3 is proved in Appendix A.3 It shows that in order to compute the maximum and minimum values of T({= : j E Uz}) for 2 E Rd, it is sufficient to check the imputations of z where its missing components are imputed using the components æi...., Itt.W/.... V/.......... However, computing T({(()) j E U,}) for all possible imputations using Lemma 3 is (4, + (=:1, which is exponential in the number of unobserved components of and impractical to compute. To address this computational challenge, we propose to further bound min X' and max X, using the following result:
8

---

Published in Transactions on Machine Learning Research (11/2025)
Lemma 4. Following the notation and definitions in Lemma
denote
*+(j)= = max (|xx(i) - xn (i)), (ms(1) Te,, xe, (1)l [x(i)) - -n (o,(1) - yes (i))
for any i {1..... ( E U2, and denote
]+(5) = max (lbx (i) - 21 (.)).....,.luu - In, (3)1,11(5))
yc, (i)
for any i E (1):..,(2):] E Uz. Suhsequently, we have
4
E
6
min X 2 Σ exp -8
by max X'
Σ
b,
exp
-B Σ 1500)
JEU,
JEU.
This lemma is proved in Appendix in Lemma g
A.4
it provides the lower and upper bounds for min X' and max X' defined
The bounds given by Lemma 4 can be computed efficiently. For each E (1....,(1) and jE U, |0,()= = 2 according to its dcfinition in Lemma Hence, in order to bound min X, computational times [U + E2) are needed, which is of order (ddn (w) +2))). Similarly, to bound max X, the order of computation complexity is (dng(nz + n2)). This computational cost reduction, compared with Lemma 3 which grows exponentially with [U=], makes bounding MMD (X, possible in practice.
We are now ready to prove the main result. for the multivariate data. In order to state the result, we first make the following two definitions. The first definition defines incomplete Laplacian kernel, which allows us to compute the Laplacian distance between two partially observed data points.
Definition 3. For any 2, y E RRd, let Uza {l,...,d d} be the index of dimensions for which either or y have components that are not observed. Let le denote the Laplacion kernel with parameter B. Then, the incomplele Luplacion kernel is defined us
K(a,y) eepp -3
Σ
|z(i)- (z(1) - s(0)
The second definition defines constants 61, 02,02, C3 and C1, C2, C3 C4 depending on the observed data for bounding the MM, Y) test statistic.
Definition 4. Suppose X = (ZI,..., and Y = (v....., Nrzt are samples of Rd real values. Assume Th...., Imny , "1...., .uz are data that are missing, or not completely observed. Let k denote the Laplacian kernel and K. denote the incomplete Laplacian kernel defined in Definition 3 Further, denote
(1))) max (1): Jmn (mm++1(1)) (1)) (x:(1) Fnn (1)1l,1z; Ymng + +1 (1) -.. (1;(1) Bng (1)
for any i E {m: + ,,...,,1]]... E (l, and denote
(i(()= max (1) 2m (1)) (11(1) - In, (v(()) Um (1)==(())--- (y((() yrez (I)
for any i E (m2 + 1,....,.y m E (1....,d) ,d. Let
mi
N: K(x,,0;)
22 E
G max
0,
C3
K(o,, exp
-3 Σ 81(1)
j=my
j jmmg+1
LEU
mi
n2
C2 min
0, C1
Σ
K( exp
3 Σ
3,(1)
C3
K( U3}
j=m++1
IEU::
j=ma+1
mg
n2
12 Σ
C3=
max
0, C2
K(yi. My) C3
exp
-9 Σ 11
1=
j=mi+1
LEU
m2
C4
min
0, 02
K(yi exp
-B Σ
c3
Σ
K(yi,
i=1
j my
IEUM.
j j=myt +1
9

---

Published in Transactions on Machine Learning Research (11/2025)
where C1 "1 2 (11 12 m2(#2-1) 1)
2 and C3 1ig
Then, using Definition 3 and Definition the main result can be stated as
Theorem 2. Suppose X = {z1...., Ing and Y = fyh...., Yngt are samples of TRd real values. Assume
#1 E1,..., Im "1..... Yozz are data that are missing, or not completely observed. Subsequently, the (X, using Laplacian kernel k bounded, over all possible imputations of :i,...;2mp; ......, 9m2 as follows:
m2g
0,) +
+2 (X.)(X
m
m2
MMD C3
x(o,)))) A2 + C2
where K. is the incomplete Laplacian kernel defined in Definition S Ag is a constant of observed data defined in Lemma 7 and C1,02 C2, C3 and C1, C2 C3, C4 are also constants of the observed data defoned in Defovition 4
The proof of this theorem is included in Appendix A.6
Remark 1. As we discussed previously, our strategy for bounding the MMD test statistic MMD", (X,Y) Y) is to provide the bounds for the terms AI, A2, A3 and A4, separately. Since A includes only the unobserved data and Y, the bounds of can be constructed according to Inequality (2) and cannot be improved. The second term Ag includes only observed data. The third and the fourth term A3, A4 can be decomposed into the sum of functions of T(x), shown in Incquality (3). We show in Lemma 2 that the tight lower and upper bounds of Txx) can be constructed. Hence the bounds of term Ag and A4 are tight. While all the bounds of A,, A2, A3 and A4 are tight, our current bounds for the MMD test statistic are not tight since we provide the bounds of A1, A2, A3, A4 separately, rather than jointly.
4 Two-Sample Testing in the Presence of Arbitrarily Missing Data
In this section, we discuss methods for employing bounds derived in Section 3 to develop valid two-sample lesting moth Besides computing the test statistic, we also need a method for computing y-value from the test statistic.
Bounding y-value using permutations. The permutation lest, as discussed Section 2 is a numerical procedure that can be used to compute p-value for test statistic based on the observed data. The p-value of the MMD test using permutations is computed
1 B+1
(X,o Xxa) MXD
p
1+ I+Σ
(5)
where I(E) is the indicator function such that if the event E occurs, I(E) - 1, otherwise I(E)) The null hypothesis is then rejected if p is smaller or equal to the pre-specified significance level
Using the bounds of the MMD in the presence of missing data, we proceed by providing bounds of the p-value, using the following result:
Theorem 3. Suppose X =
(uI..... Mog) E Rd. Let (o(1), ,.... (o())...,oBB) be B i.i.d. random permutations of {l,...,m, + n2n and denoted as g( = (1), ( + +n))), i= I,...,.. B. Let
21 = E1...., 2ng = Inys Zn1t1 = Y1..... 2antaa = ynn and for any = 1,..., B, denote
X = (2.0)).......................... (0m) boo = (2)00(111))) Zoth (m)+
Define p according to equation Suppose further have (X, and MMD (X,Y) such that
(X,Y < MMD MMD Yow)) S MMD, Yoio),, i= I,.... B.
10

---

Published in Transactions on Machine Learning Research (11/2025)
Define
8
B+1 1+
T(MM) (X,00, (x0,, Yo1) 2 MMD (x, k) (X,
Then, have
The proof of this theorem is included in Appendix A.7 In the presence of missing data, MMD X) cannot be computed directly. However, shown in Section the lower and upper bounds of MMD Y) can be computed when the Laplacian kernel is used. Then, Theorem B can be applied to compute p, the upper bound of p-value. It is proved in Schrab et al. (2023) that the MMD test using the p-value defined in equation 5 controls the Type I error. Therefore, using the upper bound p defined in Theorem 3 also controls the Type error, since the bounds were derived without any assumptions about the missing data mechanisms.
Bounding p-value using normality approximation. We have considered how to bound the p-value when it is computed using permutations. When the values of nu, n2, are large enough, alternative method to compute a p-value is to use the normality approximation (Gao & Shao 2023). It is suggested in Gao &
Shao
2023 that this approximation is effective when 71, n2 2 25, d 2 50.
To more formally describe uormality approximation proposed in Gao Shao 2023 we will follow their notation and definitions the following of this section. Tet us denote the sample sizes for and n and m, respectively, rather than ni, n2.. Denote m N, and let k denote the Laplacian kernel. For any s, < N, define
k(Xo, k(X,, Y(-n) k(X,, Ya-n)
s,
I n+1<s,tN,
k(Y Yi-n)
and Ab; ath,, a ag a. where
a. N 2
a
Σ j=1
d aj*
1
Σ
N
(N 1)(N 2) i==1
Further, define
Vkka n.m
Σ (Att))
(6)
N(N 3)
(N 1)(N 3)
spt
and
2
4
2
Cn,m
(7)
n(n
1)
nm
m(m 1)
Then, it is proved in Gao & Shao (2023) that the studentized MMD (Y,Y), taking the form
(8)
Cr.m
Vik n.m
converges to standard normal distribution under the null hypothesis with some conditions (see Theorem 16 in Gao & Shao (2023) ). In order to use this result, it is necessary to compute an estimation of variance of MMD Y). i.c. Ca.m n.m: which can not be computed directly with missing data. To overcome this problem, propose to also bound p.m- We begin by making the following definitions.
bbs
11

---

Published in Transactions on Machine Learning Research (11/2025)
Definition 5. Suppose X Ci.....,x1):: (Yi,..., ,,,, R and assume X .mi and Yi...., Ymg ore not observed. Let K. denote the incomplete Laplacion kernel defined in Definilion and for S s,t N, denote
K(X,,
1s8,t
Lol { g,t
1ss5m, m1 15 N,
K(X,, Yz-n)
IssSn<ist N, I
i < S n m2
K(X,, yn-n)
otherwise.
K(Y (nn)
s+
Let Az, o,,, e,t,
0 + a
where
N
N
N
k
N
2
3,
N 2 j=l
Caj
(N 1)(N 2)
1=1
i,==
Define
Ki
n,m
N(N - 3)
5,t
(N IN(N 3)
sft
Then, we bound Vik defined in equation 6 and the test statistic defined in equation g using the following results.
Theorem 4. Suppose X
(YI...., Yn.....,nng E Rd, and assume XI...,, Xm, and Yi...., Ying not observed. Let K. denote the incomplete Laplacian hernel defined in Definition 3 Then, (X,Y) (X,)(XXY)) such that
we have Dhe n.m 2 Vkk nim: Further, suppose we have MMD
0 < (X,) S (X,),
have
(M MX,
DrCo n.m
Cr.m Vos
Cr
where Ca,m and Vk+ n,m. are defined in equation and equation
separately, and the am is defined in Definition 5
The proof of this Theorem is included in the Appendix A.8.
(X.Y
In Gao & Shao 2023), the p-value is defined as p = +
where & denotes the cumulative density function of the standard normal distribution. In the presence of missing data, define
/enmm Vie
MMD (X.Y) Cn.m D
MMD 2 0,
F
otherwise.
Then, according Theorem we have p 2 p. Hence, rejecting the null hypothesis when p < a, that is, when P smaller than the significance level, will control the Type I error, provided that the p-value p defined in Gao & Shao (2023) is valid. This follows because p < implics p <
Applying MMD-Miss. Now we give details of applying our proposed method, namely MMD-Miss, using the results we have shown in Section B and Section !
In the presence of missing data, the bounds of "A,D" can be computed using the results in Theorem and Theorem 2 when 1 and d > 1, respectively. Provided the bounds of MMD (X, Y), Theorem gives non-trivial upper bound of p-value of MMD () when the permutation test used, and Theorem gives upper bound of p-value when the normality approximation is used.
12

---

Published in Transactions on Machine Learning Research (11/2025)
Denoting the upper bound of p-value as p, MMD-Miss rejects the mull hypothesis when p5 is larger than equal to the pre-specified significance level a. Further, we refer to the MMD-Miss that computes using the results in Theorem 3 with permutation tests as MMD-Miss-Peres and the MMD-Miss that computes p using the results in Theorem 4 with normality approximation as MMDD-Niss-Normalliy
We note that, since MMDD-iss-Normality (2023), it is valid only if the normality approximation is considered valid, which requires sufficient large sample sizes n, n2 and dimension d. As we discussed before, one rule to use the approximation is when R1, 2 25, d 2 50, suggested by Gao & Shao 2023
relied on the normality approximation proposed in Gao & Shao
5 Experiments
This section investigates the Type I error and statistical power of MMD-Miss proposed in Section 4 and compare it with other missing data methods, including missForest Stekhoven & Bühlmann (2012), the Wilcoxon-Mann- Whnnney test with missing data method Zeng et al. (024) (WMWM), and three common missing data approaches, namely case deletion, mean imputation and hot deck imputation. All these approachos are described in more detail in Appendix [B.3] For MMD-Miss, the parameter |3 in the Laplacian kernel is chosen using the median heuristic, which generally works well (Gretton et al.) 2012a Bodenham &
Kawahara 2023) and is described in Appendix B.2 The number of permutations used for MMD-Perm and
the imputation methods is set to B = 100, as described in Appendix B.2
Asymptotic Type I error for small proportion of missing data (5%). The first experiment investigates the Type I crror of the proposed method MMD-Miss and the common missing data approaches, as the sample size increases, while keeping the proportion of missing data fixed at 5%. In this experiment, the data are missing not at random (MNAR). To assess the Type error, random samples X = {zI,... Zn, } y = fy........ } are independently generated from a d-dimensional normal distribution with mean vector 4 = (0, (0......)" 0)* and covariance matrix equal to the identity matrix. Subsequently, 5% of data will be selected in both X and Y to be labelled missing (for d = 1), or incomplete (for d 1). For multivariate data, when observation is labelled incomplete, 30% of its components will be missing values. We select 5% of data to be missing or incomplete to show that even small proportion of missing data in the data set could lead to invalid testing results, if the missing data are not laken inlo account. Similarly, for multivariate data, the 30% missingness per incomplete observation is chosen to introduce a moderate amount of missing data.
We run this experiment for different dimensions d E 10, 50). In the univariate case, when the dimension 10, d I, observations I; E X with 22 < 0 will be randomly selected to possibly be missing, while only observations y E Y with yi > 0 will be randomly selected to possibly be missing; in this case the data are missing not at random (MNAR), i.e. the data are informatively missing. In the multivariate cases, when d [10, the observations in X that will possibly be partially missing are those with Σ'_ <1(l/Vd < 0.8, and then for cach chosen sample, 30% of its components with values smaller than m2dian(z,(1),... (.). (d)) will be randomly selected to be missing. In other words, only observations with low total values, i.e., with a sum smaller than - have a chance to be incomplete, and within these, components that fall below the median will be missing. This design creates MNAR patterns, where the presence of the values depend on both missing and observed data. On the other hand, only observations in Y with Σi= (1)/v/ > 0.8 will possibly be partially missing, and again for each chosen observation, 30% its components with values larger than modian(ys (1), (d) will be randomly selocted to be missing. We use these missingness mechanisms, since if only the fully observed observations are taken into account, the two samples will appear to be different.
Figure I shows that MMD-Miss- Perm controls the Type I error for dimension d E {1, 10, 50}. MMD-Miss Normality is only suitable for d = 50, and controls the Type I error in this case. WMWM (Zeng et al. 2024) controls the Type I error for d = 1, but is not suitable for multivariate data. All imputation mods fail to control the Type I error. Note that MissForest (Stekhoven & Bühlmann 2012) is only suitable for multivariate data, but similarly to the other imputation methods, it fails to control the Type I error. In particular, Figure shows that even when only 5% of the data are missing, as the sample size increases to 5000, in these experiments the Type I error of all imputation methods approaches 1.
13

---

Published in Transactions on Machine Learning Research (11/2025)
Univariate Samples: d=1
Multivariate Samples: d-10
Multivariate Samples: d-50
1.00
1.00
Type I Enrr
Type Erro
I Error
0.75
0.75
0.75
0.50
0.50
0.50
0.25
0.25
0.25
0.00
0.00
1000
2000
3000
4000
5000
0
1000
2000
3000
4000
5000
1000
2000
5000
4000
5000
aample Size
Sample Size
Sample Size
MMD Miss Peem
Case Deletion
Mean Imputation -MissForest Sig level: 0.05
MMD -Mss -Normalityy
Hot Deck Imputation - WMWM
Figure 1: The Type error of the proposed method and common missing data methods when data are missing not at random. (Left): d = 1; (Middle): d = 10; (Right): d = 50. MMDD-Miss-Normallit is considered when d= 50. For all figures, the significance level is a = 0.05, and 5% of the data are missing or partially observed. The displayed values show the average times of the mull hypothesis rejected over 100 repetitions. The error bars represent one standard error of the mean. Note that WMWM is only applicable for univariate data, while MissForest in only applicable for multivariate data.
Type I error and power for univariate data. The second experiment investigates the Type error and statistical power for a fixed sample size, while the proportion of missing data increases. In this experiment, the data are again missing not at random (MNAR). To assess the Type I error, the data X (.i,...,, Ln,} and Y = {.1...... .n2) } are independently sampled from a standard normal distribution N(0, 1). To assess the statistical power, X are sampled independently from N(0, 1) while Y are sampled independently from either (a) N(1, 1) (b) N(1.5, 1). Then, proportion E {0, 0.01, 0.20} of samples are selected in X and Y to be missing. The missingness mechanisms are the same as for Figure [ when d = 1; only observations E X with H; < and observations y, E Y with > 0 will be randomly selected to possibly be missing.
Figure 2 shows that the Type error is not controlled by the case deletion, hot deck imputation or mean 5%, the Type I error is above 25% for each of these three methods, and when the s 10% the Type error of these three methods is beyond 75%. On the other hand, the proposed method MMD-Miss-Pernn controls the Type I error, regardless of the proportion of missing data. While all the common missing data methods have good statistical power, the power of the proposed method decreases significantly when there are more than 5% missing data for alternative N(0, 1) vs N(1, 1), and when there are more than 10% missing data for alternative N(O, vs N(1.5, 1). This experiment demonstrates that MMD-Miss-PPrmn is useful when the proportion of missing data is in the range 5% to 10%, while the three common missing data approaches fail to control the Type I error.
imputation methods. When the proportion of missing data
We note that WMWM ZZng al.) 2024) controls the Type error in all cases, and has better power than MMD-Miss-Pedd Perm when there is change in the mean. However, it has almost zero power when there is change in variance (bottom row of Figure 2 this is not surprising, since it is based on the Wilcoxon- Mann-Whitna Whitney test, which is designed to test for a change in location. On the other hand, MMD-Miss-Poott cont the Type I error and has good statistical power when 5% to 10% of the data are missing, both when there is a change in mean, and when there is a change in variance.
Remark 2. As reflected in these results, our proposed method is guaranteed to control the Type I error, even in cases where data are missing not at random. However, the cost is that our method is conservative, since it only declares a significant result when the maximum possible p-value, considered over all possible imputations, is significant. When the proportion of missing data increases, the upper and lower bounds of the MMD test statistic increase and decrease, respectively, and so the power decreases.
Type I error and power for multivariate data. The third experiment compares the Type error and power of MMD-Miss and the three common missing data approaches for multivariate observations. For a value R. let. Pa,d (a, a .. a) be mean vector with d components all equal to To assess the
14

---

Published in Transactions on Machine Learning Research (11/2025)
N(0.1) vs ((01))
N(0,1) vs N(1,1)
N(0.1) vs N(1.5,1)
1.00
1.00-
Type EEror
0.75
0.75
0.75
Power
Power
0.50
0:50
0.50
0.25
0.25
0.25
0.00
0.00
0.00
0.05
0.10
0.00
0.05
0.10
0.15
0.20
0.30
0.05
0.10
0.15
Proportion of Missing Data
P Prportion of Missing Data
Proportion of Missing Data
N(0,1) vs N(01))
N(0,1) vs N(0.4)
(O0,1) vs N(0.9)
1.00
1.00
1.00
Type Enor
0:75
0.75
0.75
Power
Power
0.50
0.50
0.50
0.25
0:25
0.25
0.00
0.00
0.00
0.05
0.10
0.15
0.20
0.00
0.05
0,10
6.15
0.20
0.50
0.05
0,30
0.15
0.20
P Proportoon of Missing Data
P Proportion Missing Data
Proportion of Missing Data
MMD M-MM-Misssmssss---- Hot Deck Imputation
WMWM
Case Deletion
M Maan Imputation
Sig level: 0.05
The Type I error and power of the proposed method and common missing data methods for univariate samples when data are missing not at random. (Left): ((0, 1) vs N(O, 1); (Middle): N(O, 1) n2 = 500 are used.
Figure 2:
N(1,)); (Right): 1)0, vs N(1.5, 1). A significance level a = 0.05 and sample sizes The displayed values show the average times of the mull hypothesis is rejected over 100 repetitions. The error bars represent one standard error of the mean. The Type error figure in the left-hand column appears in each row for convenience of comparison.
Type I error, obsorvations X = (x1,..... Ing} } and Y = (yh,...,,un,) are independently sampled from the d-dimensional normal distribution with zero mean vector Pu,d = 0...,....,00 0) and covariance matrix Σ Id the identity matrix. To assess the statistical power; X are independently sampled from N(Po,d; Ta) and Y are samples independently from either (a) N(u1,d, Ia) (b) N(HI.,,d, Ia). Then, proportion E [0. 0.2] of obscrvations from both X and Y will be randomly sclected to be partially missing. The missingness mechanisms follows the same approach used for Figure [ when d > 1. The results displayed in Figure 3 show a similar pattern to the results for univariate case: (i) the three common missing data approaches cannot control the Type error although they all demonstrate good power, while (ii) MMD-Miss controls the Type error, and also have good power when proportion of missing data is around 5% to 10%. Notably, For
in the case where d = 50, MMDD-Miss-Normalll the alternative N((1.5,... 1.5)*, La), MMD-Miss-Normalli approximation has good power when 15% of the data have missing values.
has better statistical power than MMD-Miss-Pern
MNIST dataset. We evaluate the performance of MMD-Miss on real-world data using MNIST images LeCun et al. (1998), with examples shown in Figure 4 in Appendix B.4 Fach image in the MNIST dataset. has dimensions of 28 28 pixels and is labelled from to 9. For our analysis, the pixel values of each image are scaled between 0 and 1. To assess Type I error, datasets X and Y are generated by randomly sampling with replacement from MNIST training set images labelled as 3. For evaluating power, X is generated by randomly sampling with replacement from images labelled as O, while Y continues to be sampled from images labelled as 3. A proportion s E [0, 0.25] of the samples in Y are then randomly selected and labelled as incomplete if there are more than 85 non-zero pixels in the region defined by rows 1 to 14 and columns 8 to 21 (i.c. a sub-block the upper half of each image), which will be marked as missing. In other words, images with more non-zero pixels in the specified region are more likely to be partially observed. Examples of these incompletely observed images are shown in Figure 5 in Appendix B.4 Table shows that the three common
15

---

Published in Transactions on Machine Learning Research (11/2025)
NO lasa) vs N(O", laxel: d=10
NO laxe) vs N(1", laxal; d-10
N(O% laxe) vs (1.5", laxa); d=10
1.00
1.00
1.00
Type Ennor
0.75
Power
0.75
Power
0.75
0.50
0.50
0.50
0.25
0:25
0.25
0.00
0.00
0.00
0.00
0.05
0.10
0.15
3.20
0.05
0.10
0.15
0.20
P roportion of Missing Data
P Proportion of Missing Data
Proportion of Missing Data
N(0 laxa) vs N(O" lasa): d=50
N(0 lasa) N(1 laxa): d-50
N(0 lasa) (1.55% laxal: d=50
Type Emor
Power
Power
0.00
0.00
0.00
0.00
0.05
0.10
0.15
0.20
0.00
0.05
0.10
0.15
0.20
0.50
0.05
0,30
0.15
0.20
roportion of Missing Data
Proportion Missing Data
P Proportion of Missing Data
MMD -Miss Perm
Case Deletion
Mean Imputation
Sig level: 0.05
MMD -Mis -Nrrmalit
Hot Deck Imputation - MissForest
Figure 3: The Type I error and power of MMD-Miss and the three common missing data approaches for multivariate samples when data missing not at random. The dimension of the data is either d = 10 d = 50. When = 50, the SMMD-Mis--Normall ni = n2 = 500 are used. The displayed values show the average times of the null hypothesis is rejected over 100 repetitions. The error bars represent one standard error of the mean.
is used. A significance level = 0.05, and sample sizes
Table 1: Comparison of Type error and power MNIST dataset. NA stands for MMD-Miss MMD-Miss-Normaltty PT for MMD-MMsss-Perm; CD stands for case deletion; MT stands for mean imputation; HD stands for hot deck imputation. A significance level a = 0.05, and sample sizes nh table are the average times the null hypothesis is rejected over 100 repetitions.
500 are used. The values in the
Proportion of Missing Data (s)
Type Error
Power
MMD-Miss
Standard
MMD-Miss
Standard
NA
PT
CD
MI
HD
NA
PT
CD
MI
HD
0.00
0.07
0.07
0.06
0.07
0.07
1.00
1.00
1.00
1.00
1.00
0.05
0.00
0.00
0.08
1.00
1.00
1.00
1.00
1.00
1.00
1.00
0.10
0.00
0.00
0.21
1.00
1.00
1.00
1.00
1.00
1.00
1.00
0.15
0.00
0.00
0.64
1.00
1.00
1.00
0.00
1.00
1.00
1.00
0.20
0.00
0.00
0.91
1.00
1.00
0.96
0.00
1.00
1.00
1.00
0.25
0.00
0.00
1.00
1.00
1.00
0.00
0.00
1.00
1.00
1.00
missing data approaches cannot control the Type error for this task. On the other hand, MMD-Miss controls the Type I error while enjoying good power. For MMDD-Mis-NNrmmalln the power good except for the case when 25% images from Y are missing.
6 Conclusion
The proposed method MMD-Miss performs two-sample testing on both univariate and multivariate data with missing values for any distributional shifts. This method is based deriving precise bounds of the Laplacian
16

---

Published in Transactions on Machine Learning Research (11/2025)
kernel, and rejects the mull hypothesis when all possible missing data lead to significant results. MMD-Miss controls the Type I error, while other common missing data methods such case deletion fail to control the Type I error when the data are missing not at random (MNAR). One limitation of MMD-Miss is that it is restricted to the Laplacian kernel for computational reasons. However, the Laplacian kernel is characteristic, and MMD with this kernel can detect any distributional change. Moreover, this approach will be effective for other kernels if appropriate bounds can be derived and computed efficiently.
Simulation results show that MMD-Miss has good statistical power, typically when 5% to 10% of the data are missing. Our experiments show that in high-dimensional situations, i.e. d 2 50, using normal approximation (Gao & Shao 2023) (MMD-NNrrmalty) improves the statistical power of MMD-Miss based on permutation test HAMD-Permi) Another limitation of MMD-Miss is that it is only suitable when up to 20% of the data is missing. However, our experiments show that using imputation methods when even 5% of the data are missing not at random can lead to the Type error not being controlled. Overall, MMD-Miss is a robust approach to two-sample testing the presence of missing data for detecting changes in univariate or multivariate data
References
Nachman Aronszajn. Theory of reproducing kernels. Tronsactions of the American Mathematical Society, 68 (3):337-400 440,, 1950.
George L Bakris, Rajiv Agarwal, Juliana Chan, Mark Cooper, Ron T Ganscvoort, Hermann Haller, Giuseppe Remuzzi, Peter Rossing, Roland E Schmieder, Christina Nowack, et al. Effect of finerenone on albuminuria in patients with diabetic nephropathy: A randomized clinical trial. JAMA, 314(9))884-894 2015.
Felix Biggs, Antonin Schrab, and Arthur Gretton. Mmd-fuse: Learning and combining kernels for two-sample testing without data splitting. Advances in Neural Information Processing Systems, 36:75151-75188, 2023.
Dean A Bodenham and Yoshinobu Kawahara. eummd: efficiently computing the mnd two-sample test statistic for univariate data. Statistics and Computing, 33(5):110, 2023.
Song Xi Chen, Jing Qin, and Cheng Yong Tang. Mann-whitney test with adjustments to pretreatment variables for missing values and observational study. Journal of the Royal Statistical Society Series B: Statistical Methodology, 75(1):81-102, 2013.
Ying Kuen Cheung. Exact two-sample inference with missing data. Biometrics, (1()):224-531, 2005.
Jobson Dave and Bob M Korkic. Performance hypothesis testing with the sharpe and trcynor measures. Journal of Pinance, 36(4):889 908, 1981
Youyi Fong, Ying Huang, Maria P Lemos, and M Juliana Mcelrath. Rank-based two-sample tests for paired data with missing values. Biostatistics, (881)221-2294 2018
Kenji Fukumizu, Arthur Gretton, Bernhard Schölkopf, and Bharath K Sriperumbudur Characteristic kernels on groups and semigroups. Advances in Neural Information Processing Systems, 21, 2008.
Hanjia Gao and Xiaofeng Shao. Two sample testing in high dimension via maximum mean discrepancy. Journal of Machine Learning Research, 24(304):1-33, 2023.
Phillip Good. Permulation, parametric and bootstrop tests of hypotheses. Springer, 2005.
Arthur Gretton, Kenji Fukumizu, Zaid Harchaoui, and Bharath Sriperumbudur. fast, consistent kernel two-sample test. Advances in Neural Information Processing Systems, 22, 2009.
Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf, and Alexander Smola. kernel (wo-sample test. The Journal of Machine Learning Research, 13(1):723 773, 2012a.
Arthur Gretton, Dino Sejdinovic, Heiko Strathmann, Sivaraman Balakrishnan, Massimiliano Pontil, Kenji Fukumizuu, and Bharath K Sriperunbudur Optimal kernel choice for large-scale two-sample tests. Advances in neural information processing systems, 25, 2012b.
17

---

Published in Transactions on Machine Learning Research (11/2025)
Martijn W Heymans and Jos WR Twisk. Handling missing data in clinical research. Journal of Clinical Epidemiology, 151:185-188 188, 2022.
Matthias Kirchler, Shahryar Khorasani, Marius Kloft, and Christoph Lippert. Two-sample testing using deep learning. In International Conference on Artificial Intelligence and Statistics, pp. 1387-1398. PMLR, 2020.
Andrey Kolmogorov. Sulla determinazione empirica di una legge didistribuzione. Giorn Dell'ina Ital Degli Att, 4:89 91, 1933.
John M Lachin. Worst-rank score analysis with informatively missing observations in clinical trials. Controlled Clinical Trials, 20(5):408-422, 1999.
Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient- learning applied document recognition. Proceedings of the IEEE, 86(11):2278-233 2324, 1998.
Roderick JA Little and Donald B Rubin. Statistical analysis with missing data, volume 793. John Wiley & Sons, 3 edition, 2019.
Žikica Lukić and Bojana Milošević. A novel two-sample test within the space of symmetric positive definite matrix distributions and its application in finance. Annals of the Institute of Statisticul Mathemalics, 76 (5):797-820, 2024.
Henry B Mann and Donald R Whitney, On test of whether one of two random variables is stochastically larger than the other. The Annals of Mathematical Statistics, 18(1):50-60, 1947.
Charles F Manski. Partial identification of probability distributions. Springer Science & Business Media, 2003.
Pablo Martinez-Camblonn Norberto Corral, and Jesus Maria de la Hera. Hypothesis test for paired samples in the presence of missing data. Journal of Applied Statistics, 40(1))76-87, 2013.
Maya B Mathur. The m-value: A simple scnsitivity analysis for bias duc to missing data in treatment effect estimates. American Journal of Epidemiology, 192(4):612 620, 2023.
Peter O'Brien and Thomas Fleming. A multiple testing procodure for clinical trials. Biometrics, 35(3): 549-55 156,556,440 1979.
Wei Pan. A two-sample test with interval censored data via multiple imputation. Statistics in Medicine, 19 (1):1-11, 2000.
Roland Pfister and Markus Janezyk. Confidence intervals for two sample means: Calculation, interpretation, and a few simple rules. Advances in Cognitive Psychology, 9(2):74, 2013.
Maria L Rizzo and Gábor Székely. Energy distance. Wiley Interdisciplinary Reviews- CompuussCompttttionnl Statistics, 8(1):27-38, 2016.
Donald B Rubin. Inference and missing data. Biometrika, 63(3))581-592, 1976.
Donald B Rubin and Nathaniel Schenker. Multiple imputation for interval estimation from simple random samples with ignorable nonresponse. Journal of the American Statistical Association, (1(34))3366331 1986.
Joseph Schafer. Multiple imputation: a primer. Statistical Methods in Medical Research, 8(1):3 15, 1999.
Joseph L Schafer and John W Graham. Missing data: our view of the state of the art. Psychological Methods, 7(2):147, 2002.
Antonin Schrab, Ilmun Kim, Mélisande Albert, Béatrice Laurent, Benjamin Guedj, and Arthur Gretton. Mmd aggregated two-sample test. Journal of Machine Learning Research, 24(194):1-81, 2023.
Shubhanshu Shekhar, Immn Kim, and Aaditya Ramdas. permutation- free kernel two-sample test. Advances in Neural Information Processing Systems, 85:18168-18118 18180, 2022.
18

---

Published in Transactions on Machine Learning Research (11/2025)
Nobumichi Shutoh, Makiko Kusumi, Wataru Morinaga, Shunichi Yamada, and Takashi Sco. Testing equality of mean vectors in two sample problem with missing data. Communications in Slatistics Simulation and Computation, 39(3):487 500, 2010.
Bharath K Sriperumbudur, Arthur Gretton, Kenji Fukumizu, Bernhard Schölkopf, and Gert RG Lanckriet. Hilbert space embeddings and metrics on probability measures. The Journal of Machine Learning Research, 11:1517-155 1561, 2010.
Bharath Sriperumbudur, Kenji Fukumizu, and Gert RG Lanckrict. Universality, characteristic kernels and rkhs embedding of measures. Journal of Machine Learning Research, 12(7), 2011.
Daniel J Stekhoven and Peler Bühlmann. Misslorest Missooreat-mon-paraaaeeee missing value imputation for mixed-type data. Bioinformatics, 28(1): 1182,1111111 2012.
Daniel J Stekhoven and Maintainer Daniel J Stekhoven. Package 'missforest. R packuge version, 1:21, 2013.
Student. The probable error of mean. Biometrika, 6(1):1 25, 1908.
Lehana Thabane, Lawrence Mbuagbaw, Shiyuan Zhang, Zainab Samaan, Maura Marcucci, Chenglin Ye, Marroon Thabane, Lora Giangregorio, Brittany Dennis, Daisy Kosa, et al. A tutorial on sensitivity analyses in clinical trials: the what, why, when and how. BMC Medical Research Methodology, 13(92):1 12, 2013.
Anastasios A Tsiatis, Marie Davidian, Min Zhang, and Xiaomin Lu. Covariate adjustment for two-sample treatment comparisons in randomized clinical trials: a principled yet flexible approach. Statistics in Medicine, 27(23)44658 4677, 2008.
Yijin Zeng, Niall M Adams, and Dean A Bodenham. On two-sample testing for data with arbitrarily missing values. arXiv preprint arXiv:2408. 15327, 2024.
19

---

Published in Transactions on Machine Learning Research (11/2025)
Overview of Appendix
Appendix Mathematical details of results given in Section and Section
Appendix B Details of experiments in Section
Appendix Additional experiments.
A
Mathematical Details
A.1 Decomposition of MMD
Without loss of generality, let assume Fi:...;Tmi V!, M...... Umz data that are not observed, not fully observed. Then, MMD"(X,>) can be decomposed, using the following lemma:
Lemma 5. Suppose X ( Fni fy1.-.., Una) are subsets of RRd where d 2 Suppose that MMD" Y) is unbiased MMD test statistic defined
(X,)) m1
k( 223) + m2(22-1
k(
k(x,,8g),
2m
J
3
with l. denoling the Laplacian kernel. Then, for any two posilive inlegers mi, m2 such that mı S n1, m2 S n2, (X, Y) can be rewritten as:
- - A A1 A4,
where C1
01(n-11) (11 C2
22 1J
C3
and
m22
A = 9 Σ
Σ
13), + 02
U1) C3
k Uj);
1=1 j=it1
1=1
j=1
B1-1
R2
n2
ni
12
A = C1 Σ
211 +
k M))
C3
(z).8;).
i= 1=
i=m2
j=it1
i=mi+ j=m=+11
A3 = C1 Σ
Σ
C3
i=1 j=m+++
J=
A4 = Σ
Σ
ky,,8;)
-3
Σ
Σ
k(
j=m2+1
i=m +=
Proof. To start, for any x, y of Rd real values, the Laplacian kernel is defined as follows
Mx,, y) exp(- -011).
It can be seen that & is symmetric, i.c. k(o, 3) k(y,2) Subsequently, according to Lemma 2 in Bodenham & Kawahara (2023), the MMD? (X,Y) can then be rewritten as
2
2
n2
MMD (X,Y) n2 [n
k( R25) n2(n22 - 1)
Σ
k y,)
j=1
j=1
R1
( yz)-
12g
1=1 1=1
20

---

Published in Transactions on Machine Learning Research (11/2025)
Notice that
21
13
k,);,,,)) Σ
k(
i=1
n1=1
Σ
k() = Σ k(yi,
i=1 1=1
i=1
j=1-1
Thus, we further have
2 -1
2 1
(X
k(I;, x5) +
k(ys, Uj)
m1
n2(n2 - 1)
i=1
j=++1
i=1 1=++1
2
k( Bj)
722
Let m1, m2 be any two positive integers such that my S n2, m2 I n2. Denote C1 m1(11 11 Subsequently,
2
2 R2 -1
n1 1zz;,
m1
R1
R2
R2
c
Σ
Σ
ta
Σ
Σ
k(1,,
M1(n - 1) j=it1
imm 1iji+1
m1
Ci
k(x;; mz)) +
Σ
4z)
1=1+1
1=1 j=mn+1
81
C1
Σ
Σ
2z),
i=my
+=1+111
Also, denote 11 it follows that
2
2
1
k(ys, Uj) = C2 Σ
m2
Σ k(ys, Uj) +
122-1
Σ
k(yi.
12 (n2
1) i=1 1=++1
j=i+1
im +
1 1-1++1
Σ
k(ys, Uj) +
k(y.. Vj)
t=1 1
ng-1
+ 02
Σ
k(y,
1=m-
Denote C3
then
ming
2 11 12
171 ng
ng
k( Mj)
C3
k( U3)
+ C3
Σ
k( Uj)
70.1 772
1=1
(=m
mi
122
mi
12
C3
k( 91)
+
k(0,,)))
jmmg+1
+ C3
Σ
k( Uj) +
Σ
Σ
k( vz).
1 1
+1
i +1 +ij=ma1
21

---

Published in Transactions on Machine Learning Research (11/2025)
Combine the above together, MMD (X, can be rewritten
meg
mi +
MMD) C1
Σ
2
Σ
2z)
j=it1
;=1
j=mi+
+
Σ
Σ
2z) +
k(b..)))
i= 1=++1
1=
fgg
fg
f2
02
k(y: +
Σ
k(ye. Mi)
f=
(=
f=++1
m1
C3
k( U5) - C3 Σ
Σ
k( Uj)
j=1
j j=mz +1
k(
k( Vj)-
+m1+1ja 1=
i j=
By rearranging the above equation, we conclude that
mi
m1
ma
m1g
m
my
XXD = e
k)i, + 02
k(yh: Mb)
C3
k(r; U))
j=1
R2-1
+
Σ
k(I,, 2z) +
k( Us)
j=i+1
i=mn
j=i+1
i= j=ma++
m
81
mi
n2
6
Σ
k(
C3 Σ
U))
j=m1-1
j=m++1
m22
72
k(,,,;)) - 03 Σ
Σ
k(x,,
j=m2-
i=m1 +=
which proves our result.
A.2 Proof of Lemma 2
Lemma 6. Let II...., Ft, and ,1,..., Vt, be univariate real values, that are observed. Suppose a1,... Cty bi....,br bi, and 3 are positive constants. Define
4 by exp(-33lus 21)
T(:)
Σ a4 exeepppBBaa 21)
as a function of 2 E IR. Subsequently, for any given 20 E R,
T(=)) 2 min(0, (())((,,..... T(x)), ((nn), ((,))..., T(ve.))
(9)
On the other hand,
T(zo) < Ttttt
(10)
22

---

Published in Transactions on Machine Learning Research (11/2025)
Proof. We will first prove that inequality 9 holds. Let 21 = min(21..... ...,...,,....................... #21 Mi....... .62). When 21, it follows
T(2)
dg exp(-B((x; - =))
Σ 1=1
by ( 1-((i - =])
-
dg
21 + - -2))=
bg exp(=B((k- 21 + 21 -z))
i=1
4
4
Σ
21 exp(-3((3 - -2)]
Σ
bs exp(px3(1 =1) 2))
8
by exp (-3(y) 21)) }
= exp(- - =))
Σ
exeepp-((e - -21)) -
Σ
1=1
1=1
= exp(-8((21 - -2)) (21).
Notice that, since 8 > 0 and < 21,
< exp(-8(<1 - 2)) < 1.
Hence, if T(21) 20 0,
0 S exp(-8(21 - (
S T(());
0 $ T(2)
if T(z)) <0,
0> - -)) > T(z));
= T()) > T(z)).
Thus, for any given zo E RR if zo E (-00, 21), we have
T(zo) 2 min (z)))
which proves equation 9 when E (-00, 21).
Similarly, let 22- = max(f1,...... Le, .1.....,........................ He When 2 > 22, it follows
T(2)=
Gj exp(-B((2 x:)) (x:))
b, (exp(=B(e - 9:))
4
lg
ag (xxp==((== - + 22 - + -xx))==
Σ
b, eexpppB( - -22 + 22 - -v))
(=1
1=1
4
8
Σ (exp(=8(2)) x:1)) exp(-8(2= - -22)) -
Σ
by exp(-3(22 - v1)) exp(- 0[2=[[]] 22))
1=1
i=1
exp(- B(= =22)
Σ
ai (exp(-3((2 xx))
Σ
by exp((= Ui))
exp(-B((= 222
Notice that, since 8 >>0 0 and 2 > 22,
0 < exp(-3(2 - 22)) 1.
23

---

Published in Transactions on Machine Learning Research (11/2025)
Hence, if T222) 2
0 S xxp(-B(z- - 22))((22) & T(22;;
0 S T())
if T(z)) < <,
0> exp(-8(2= (2)) T(22);
= T(z) T(22).
Thus, for any given 20 E IR, 20 E (22, 20), we have
T(z)) > min(0,
which proves equation g when 2p E (22, 00).
Suppose 21 20 < 22 and 20 $ (.i,,.., 2 : V-2) Then, there must be at least one number in {.i,..., Ft, V1....,......................... ye,} smaller than 2, and at least one number larger than Subsequently, denote the maximum number in (II,..., Ft, .h..... Vt,) smaller than 20; denote as the minimum number in fai...., Fli V1..... :eat larger than 20.
Notice that
f Σ
T(:)=
0: 5xepp(-31rr -2)
by exp( -Blyy - =1)
i=1
a,I[(r; I 23) 2x)(expppppp) expp($flx; - 21) +
a,I((r; > 23) exp(-$la; 21)
t=1
Σ byI((ys $ 2x)(exppp10) =1) Σ
b,I(yy 23) exp(-Bly: =1)-
When 23 < 2 < 24, according to the definition of 2g and 24, for any Fn, V: S 23, it follows F1, Vc- On the other hand, for for any 24, Yi > 23, it follows 24, 31 2 21, which gives < 34 yi-
Hence, when 23 < < 24.
81 T(x)
4
a,I( I 23) expp(BB2z xz)) +
a,I((r, > ex)ex((B(61 - 2))
i=1
l
byI(ys I 23) cx)cxp((BB B(z -v:))
b,I(va > 23) exp( sap(-B(ys - 2))
1=1
1=1
24

---

Published in Transactions on Machine Learning Research (11/2025)
Further,
4 T(²) = Σ
4
azl S 23) expl-xx - 23 2g xx)) +
a;I((r; > 23) exp(-3(2z - 24 + 24 - 2))
i=1
i=l
Σ
bıl S 23) exp(-3(3 23 + Wi))
b,l((y: > ex) + z))
1=1
= -cxp(-3(( 23))
a,I( S (x))cxpp exp(-8(23 xz))
4
+ eex((-3( - =))
Σ
0; a:T > 23) exp( - -3(33 - 24))
1=1
-exp(-8(= - 23)) Σ
by11(11 I 23) exp(- y1)
f2
-xp(-33( - =))
bell(yy > 2)) exp(-3(li 24)).
exp(- 23))
a;ll(x; S ex)e - xxi))
Σ
bill(ni S 22)) 2x)exp((3 - -1:)
1=1
exp(-B(22 2))
a,I(x; > exp((3(( -24))
b,I(y > exp( 2a))
For notation ease, let us denote
A a,I( I 23) exp(-8(23 xx:))
byI(y) I 23) exp(-3(23 -xi))
i=1
a,I( > exp( -3(x; -24))
Σ
b,I(y] > exp( 24))
Thus,
T(:) A Acxp(=((2== 23)) B exp(-8(22 - 2)).
Notice that if B = 0, for any 23 <2<<<< 24, T(=) A expp-AB( -23) is monotonic increasing function
in (23, 24) when A < 0, or monotonic decreasing function in (23, =4) when A >> 0, or a constant function in (2g, 24) when A = 0. Thus, we always have min(T(2s), T(z)) I T(20) for zu E (23,24), which proves inequality directly.
If, however, + 0, let T(:) = 0, we have
A exp(
(-20--(22- - 23)) - B exp(-3(22 - 2))
A
exp( (3(4-3((4 2))
B
exp(- (-((2 2s))
When - < since
exp(-A(zı =2) exp(-D(s
0, then T(z) cannot Lake 0 in (23, 24). If however, - > 0, further
23))
have
A
log
B(2s -2) - f-3(2
23))
B
A
log
-B24 + 32 B2 B2g
1
A
log
+ 24 + 23-
B
B
25

---

Published in Transactions on Machine Learning Research (11/2025)
That is, [(2) = 0 when log (- + 24 + 22 21). Overall, we have shown T(2) takes 0 at most once in (23, 24) when + 0. This result will be used subsequently for proving our final conclusion.
Now, taking derivative of T22) of 2,
(T(=) -AAxxp( 1xAep((-82 za)) + 3B exp(( =z)).
Ō=
0T(=) of 2, 8s
Further, taking derivative
TT())
3 A expepp(3( 23)) + 8 exp(-3(22 2))
02
(2)(
We are going to prove T(20) 2 min
[(za]} when B + using contradiction. Let us assume that T(aa): (i,
T(zo) m0 (0, T(2)), ((=4))
(11)
Then, since B > O, it must have
0T(())
= ==0T(= <.
02
10
For aT(=) 8s
it is either (i)
oT(=) 0s
or
(ii)
[T(=)
0s
2
We are going to discuss the two cases
separately.
Suppose (i) 0T(=)
0. If for any [20, 21),
PT(x)
0. Then, 0T(:) is monotonic decreasing function
82
20 of [20, 24). Hence, for any z E [20; 24),
(T(z)
0((z)
<0
Oz
dz
Thus, T(2) is monotonic decreasing function in 21); giving us
T(z)) > T(za).
which is contradicted to the assumption equation
If, however, for any [20, 24), (*)(2)
does not holds, recall that we have shown T(:) take most
82z
once in 22) when + 0. This then implies there is only one point t [zo: 24) such that T²T(z)
0.
8'
Since 82T(2) is a continuous function and T2)(2) 03
< 0, we have
d
T
0, E [20, t).
82
Thus, 8T(x) is monotonic decreasing function in [Zo, t), giving
0T(=)
0T(=)
<, 2 E (2
(12)
Oz
Oz
However, notice that we also have
8²T(2)
(2T(2)
0>>
82 1
82
120
+ T(t) = 0 > T(2o)
0T(2)
T()) - -T(u)) > 0,
E1/ E (zo, 1) such that
(mean value theorem)
Oz
"
26

---

Published in Transactions on Machine Learning Research (11/2025)
which contradicts equation 12
Suppose (ii)
aT(=) 2 0. The prove method is similar as to when (1) holds.
0s
If for any (z3, =0], any 2 (Z3, 20l:
8'T(2) 0=
< 0. Then, oT(2)
is monotonic decreasing function of in (=3, z0}. Hence, for
0T(2)
0T(z)
2 0.
dz
Oz
Thus, T(2) monotonic increasing function in 2 E (23, 20]: giving
T(2a) < T(z)).
which contradicts to the assumption equation 11
If, however, for any E (23; zol,
@+(())
0 docs not holds, recall that we have shown T(2) taken 0 at most
8
8T(())
8TT(s)
once in (23, 24). This then implics there is only one point (23, 2ol such that
0. Since
is
82z 8
882
a continuous function and TT(())
< <0, we have
82
TT
< 0, 2 E (t, zo]-
82z
8T(2) Thus, 82
is a monotonic decreasing function in 2 E (t, 20], giving us
0T(2)
(T(z)
2 0, (4, so].
(13)
dz
20
However, notice that we also have
0+T(2) =0>
T2)
02 1
82z
Tt)) = -0> T(20)
31' E (t, =0) such that (T(2)
T(t)
T(20) 0. (mean value theorem)
Ōz
20
which contradicts 13 Hence, finish our prove for inequality g
We now prove incquality 10 Notice that
TT() b; exp(-33l: =1)
a; exep(_-8ler) 2().
1=1
Subsequently, using the result of inequality
we have
T(zo) 2 min(0. (o.nn)))T(22)...... -00.......(((p1. -T(e,,) -T(,))
T(zo) <
-T(xe,), - -T(rs))
T(z)) <
which proves inequality 10
A.3
Proof of Lemma
This section proves Lemma
3
extending Lemma
2 into higher dimensions.
27

---

Published in Transactions on Machine Learning Research (11/2025)
Lemma 7. Let 21,...,,2 Fei ,.,..., E Rd be values that are fully observed. Suppose o1....;;ey di,.... ....... beg are posilive constants. For = (z(1),.... ...,, =(d)) E R2d with missing components, define
T(((() 1 J :jE Σ a, exp -B Σ
x:(i) 2(1))
Σ
b, exp
-B E
(vs(s) 2(i)]
jEU,
JEU,
function of the unobserved components of and let
X = : 2(1) ................................ 2e, (2) ve, (1), Uz}.
Then, for any possible imputation 2 of 2,
T((='()) : :EEU,)) 2 min(0,min min X²,
(14)
and
T()*(() j U.)) < ma max X}.
(15)
Proof. We will only prove inequality 14 and inequality 15 can be proved following the same method.
Let prove inequality 14 using mathematical induction. When |U=| = this lemma become Lemma
which has already been proved. Suppose when [U=] = - 1, where is any integer such that 1 < 1, inequality 14 holds. We are going to show when [U&] = 4, inequality 14 still holds.
Without loss generality, let us assume (after relabelling) that U2 = I....,I. I. Subsequently, for any given ("i,...,,7)) E TR', follows
8 b exp
=
Oi exp
-8 Σ (x) - 2z1
3
Lu - 2j1
1=1
J=1
1=1
j=1
a. exp (-B]z;(1) - 2;l)
B
[x((i) - -2]
j=2
Σ by exp (-11yy(1) =11) exp
-
Σ (v+ - 231
/=2
Denole
d = a, exp (-ppmm)))1) - z11), i == 1,...,
6 b, exp (1)) - -2i)), i= 1....,2..
Consider a new function
(2,,,..,,,))
Σ at exp -B (xx(1) 1z1
Σ
b1 exp -B - 23|
j=
2=2
Let
x = (22,......
2)12, (f)..)...... (i). (f)..........)))) V.. ......
Then, using the assumption that when [U&| = - 1, inequality 14 holds, it follows that.
T ( 21) > min min X'),
28

---

Published in Transactions on Machine Learning Research (11/2025)
Notice that
82
Σ
df exp -3 -B
(m((i) - -2yl
U/ exp
-9 Σ
(vs =z]
-=2
j=2
Σ
a; exp (-312x(1) 1) exp
3 Σ
(x,(1)) 2z1
(16)
i=1
Σ
exp (-3lyi(1) = exp
-8 _
1=1
J=2
F2, 22.,,
Hence,
We therefore have
T(E min min X'}.
(17)
Further, let us denote
())
arg min (0) (1) 900,....,.0)),2....... i=)..
(1)))..,))))
(18)
That is,
(1 = min X'
Applying equation 16
it follows
T(2 22,,,,,)) = (') 2) min
(19)
Denote
a. a, exp
3 Σ (x) 2/11
1....,.1
by b, exp
-3 Σ
(vi(i) - =1)
1= ..... 62
j=2
Consider a new function
& Σ
4
(2)
a exp (-8)o,(1) - 21)) -1)
Σ
6y exp(-3llyi11 22)).
1=1
i=1
Let
(T" E (1,, (1) (1)))-
Fz,
Then, using the result that when |U_| = 1, equation 14 holds, it follows
T (zi) > min min X").
(20)
29

---

Published in Transactions on Machine Learning Research (11/2025)
Notice that
4
T (21)
Σ a exp (-81z,(1) 21))
Σ
6g exp (-8lyx =1l
04 exp
-3 Σ
(x,(1) =/11
exp (-3)2,(1) 21l)
j=2
Σ
by exp
-8
Σ
lv) z/11
exp (-p(-A)y<11) - =1)
==2
(21 ,
21)
According to equation 19
(zi) = min X²,
Notice that since equation 20 holds, we have
min X" 2 min min )"),
put which back into equation we further have
T( , 2/) 2 min min X²} 2 min min X").
Notice that
min x" = (z))2 E (1)( It, (1), 21 (1),...., (1)}}
= (2. 2f12 E (1), ... 0e, (1), 31 (1), ..).....3... (1))},
where
o (z), (1) (.. ... ())) 2, 2.....
according to its definition in equation 18
Hence,
min X" 2 min = T( 2/) 2 min(0, min X),
which proves inequality [4
A.4 Proof of Lemma
Lemma 8. Following the notation and definitions in Lemma denote
a+(i) max{(a; x (i), (x:(i) xz, (i)1. xx( 21 .1.... (x:(i ve, (i))
f for any i E (1, (1....,():3 E U2, and denote
(i( max lus (i) x (7)1 I,.))) Fz, (o)). 91 GII. ..l........()))............) 9r, (i)l
for any i (1....,A2)) f2).j E U.. Subsequently,
min X
>
Og exp
-B Σ a,(i)
bi
(21)
1=1
jEU,
1=1
f
8x
max X <
b exp
-3 Σ
9.(1)
(22)
1=1
1=1
jeU.
30

---

Published in Transactions on Machine Learning Research (11/2025)
Proof. We will only prove inequality 21 and inequality 22 can be proved similarly.
Without loss of generality, let us assume (after relabelling) U-= = l....,,, I, where 1 == 1U=1. Subsequently, let
(zi....,,2)) .() (. (0. (ty (1)..))
arg min
That is,
T(E),...,,5)) min
Notice that for any j E {I,...,.
2j E (z). ..,, (f), (2),)).)..., ve,
according to its definition. Hence, for any E {l,..., 11)
1zj - xx(1)) & max((mx(i) - 21 (i))...... (o((()) - Ie, (1) =R,(1),
following which
Σ j=1
Tt(1)= 2 Σ
(n,(i) - 231b i (1....,(1)
/=1
Σ
Tv( -3 Σ
(mt(i) -2]b i (1, .1,....,1)) 81).
==1
Subsequently,
Σ 0g exp
8 Σ (x)(() -zā1
Σ
by exp
- Σ
(weli) - =3]
ag exp
8 Σ ix(1)
Σ
by exp
- 1-8
lv)(t) 23]
a exp
- B Σ a+(
bi.
1=1
1=1
1=1
which proves inequality 21
A.5 Proof of Theorem 1
We now prove Thcorem 1 which bounds the MMD Y) test statistic when data are univariate.
Theorem 5. Suppose
(ZI....,IrI] IngJ } and (yI...., ]ng} are univariate real volues. Assume 21,..., Tmyy 31. .... Vmz are unobserved. Let & denote the Laplacian kernel and define
a
n2
Ti(2) C1
k(2,2,) - c3
k(z,8;),
j=my+1
j=mg+1
ag k(2,0,) - 03
m
T2(2) 02
Σ
k(z,8z),
j=m2+1
f=m++1
where C1
17. = mg(n2-1] 1j and
2
2
2
Further, let
11:11"
S {0, Ti (5mi+1)
Ti (I ), T (Yma+ ....., T (i(2n)) }
S = {0, Tá (1
T2( Th ,)..,
31

---

Published in Transactions on Machine Learning Research (11/2025)
Then, the MX, using Laplacian kernel & is bounded, over all possible values of Fim...;fm1
y1. Mmgz as follows:
my (m² 1) m2 - 1) + my max S1 m2 max S2 + A MMD (12 (12 1)
n1 (11 1)
(X,) > A2 + my min S1 + m2 min S2 71122 m2
where A2 is defined lemma
Proof.
To start, according to Lemma !
MX = A + A2 + A3 + As,
where
m2
mi 122
A1 Σ
Σ
2z) + 02 Σ
k(yl, 93) C3
k Uj).
1=1
11-1
A =
Σ
Σ
2(11,2z) +02
Σ
k(y.. Un)
(z,,8.).
1= +1J=it
1=m
i= 1 j=i+1
t= + j=
A3
k((1))))
k( Mi),
1= j=m++1
J=
12/
12
A4
k(y,, Uj) C3
Σ
Σ
k< Uj),
1=m++1
and Ci
C2 B2( 1) and C3 mar
11] nn Ai
Following Equation equation it provided
mi (m² 1)
m2 -
A >
m, m2.
(n - 1)
n2 - 1)
11
Further, using Lemma 2 it follows
mi max S 2 Ag 2 my min S1,
and
m2 max S2 2 A4 2 mg min S2,
which concludes our proof.
In order to compute the bounds of MMD (X, Y) for univariate samples with missing data using Theorem we need to compute min S1, max S1, min S2 and max Sz. Since |S1| = ni mi and computation complexity is of order O(n + n2).
Sy = n2 -22- The
A.6 Proof of Theorem 2
This subsection proves Theorem which bounds the MMD test statistic for multivariate data. To start, we prove the following result regarding the incomplete Laplacian kernel K. defined in Definition 3
Lemma 5. Suppose a and y are two samples of Iad real values and assume not all dimensions of values of x and y are observed. Denote [d] = d}. Let Uz be set includes all unobserved dimensions in x and U, be set includes all unobserved dimensions in y. Let K. denole the incomplete Laplacian kernel defoned in Definition 3 Subsequently, it follows
K(r,y) 2 k(x,, y) 2 0.
32

---

Published in Transactions on Machine Learning Research (11/2025)
Proof. According to the definition of the Laplacian kernel, we have
k(x,y) = exp -8 8 > [x(1)) - v(11)
2 20.
Also, notice that
E
(x(1) 1 - Σ [z(1) v()).
Then, since 8 > O,
d
B
Σ
[x(1) - (1))-(011 2 -3 Σ x+)) - -()
(U,
According to the definition of incomplete Laplacian kernel, therefore have
K(x, 3) = exp
-3
Σ
[x(1)) 3(01
(E(d)(()))e)))
d
2 exp
-3 Σ
[x(() - v(()
1=1
b(r,y);
which completes our proof.
Recall that in Lemma it is shown that MMD Y) can be decomposed into four parts A2, Ag and A4 with part A₁ including incomplete samples only, A2 including complete samples only, and A3, A4 mixed with both complete and incomplete samples.
Lemma 5 allows us bound term Ai by bounding any two incompletely observed samples separately. Let Uz. be set including all unobserved dimensions in 22, where 1 < i S my. Similarly, let Uy, be set including
all unobserved dimensions in yi, where < S m2. Subsequently, according to Lemina 5 it can be seen that
Em 2j=++1 K(z;,2z) 22) 2 >
ii= -==-1 k(x,, > 0, Em Ej K(y;; VU) 2
mi
my
-= 2j=it1 k(y: U₁) > 0,
mg
and E E K(I;, Ui) 2) D) M( 03) 0. Then, according to the definition of A in Lemma 1 it
follows
K(r,, 13) + 02
Σ
K(y,, U3) > A >
K(ai, Us),
(23)
1=1
where C1 mi n2 and C3 niNg We now bound the terms A3 and A4. Recall that the term Ag can be rewritten as equation 3 in both univariate and multivariate cases. For any de where S mi, let Ug, be set including all unobserved dimensions in x; and let introduce the function
Ti T((:())lEEa)) Ua,) e Σ k(r, ca Σ
k(a,, i (I,...,,.)
(24)
j=mi+1
j=m2+1
where Cj (i( 11 and C3 nina
Then, Ag can be computed as
mi
A
Σ n(a,()))e U2.).
Subsequently, by combining Temma 7 and Lemma have our main result:
33

---

Published in Transactions on Machine Learning Research (11/2025)
Theorem 6. Suppose X {21....,, Zni and Y (....... Uuaa are samples of B real values. Assume #1... .1..... Fmay 1 "h,,.... Umz are data that are missing, not completely observed. Subsequently, the MMD (X.Y) using Laplacian kernel k is bounded, over all possible imputations of #1,..., 2mm M1..... Mm as follows:
K(),, +
K(yi, + A2 + C1 +3 (XY).
Myg
m2g
MMD
C3
K(r,,,,))+++ A2 C2
where K. is the incomplete Loplacian kernel defined in Definition S A2 is constant of observed data defined in Lemma 7 and C1, C2, Cg and C₁, C2, C3, C4 are also constants of the observed data defined in Definition 4
Proof. In Lemmna it is proved that
(X,)) A A2 A3 + A4,
where
12
mgg
12
mi m2
A1
15 + 02
k(y;,
C3
k Vj);
B1-1 A = Σ
k(I +
Σ
k(y:. Un)
k(x,, (1),
i= +1j=
1= 1 j=
i= + j=
A3
k(11,,z))
k( Uj):
1=l j=m++1
1=1 J=
172
112
R1
m2
A4 02
k(y,, Mj) C3
Σ
k(
m2
and Ci
and C3 nine
A1(m1 -1)
m2(n2-1)
According to incquality
23
m2
K(x,, x3) + Σ
Σ
K(y,8,)>Ar C3
K( Uz).
(25)
1=
Notice that
Ag = Σ
C1 Σ
k(I;,2z)
Σ
k( U3)
iml
jmmi+ +1
jmma+1
Let us introduce the function
R2 Σ
n2
(a; E Un) =c
,10,,1;) Cg
Σ
k(z,, i fl....,m,
j=m1+1
j=mz+1
Subsequently, Ag be computed
m4
Σ T(6:00))E
(26)
Denote
Xx, (Ti(Z,(1;; U2)) E (1mm+1(0),.......... Ini (1, (1)+1 ,n,, (0), Uz.)
34

---

Published in Transactions on Machine Learning Research (11/2025)
for any = 1,.... my. Let [d] {l,...,,e d). Notice that
=3,011
T (1;;1)); EU2)) =
Σ
dij exp
exp
-1 Σ
[.,(1)
j=m++1
LEUx
N2
Σ
B exp
-3 Σ
[x:(1) - y5(00)
jmmg+1
LEUN;
where
d.j C1 exp
-8 Σ
(x,(1) 2z(01
b/g cy exp
E LE|d/Us,
(2,(1) -111)-9(011
Then, according to Lemma 3 follows
max max X 2 T((,(0))EEE) 2 min(0, min Nx 1, 1,..., mi
(27)
According to Lemma 4
min Xx 2
de
exp
-3 Σ
Zill)
Σ Ulas
j=m++1 +1
LEUn;
j=
max Xz.
d
Σ
B exp
-3 Σ
35(1)
j=m1+1
i=m2+1
LEUx,
Notice that according to the definition of incomplete kernel in Definition
3
d, = c1 exp
-B Σ
[x:(1) - x0((1
= =.((();))))
WE(d)
8/1 03 exp
B Σ
[x((1) - y,(1
= c,,((1;
LE(dWW,
Hence, for any i 1, 1,..., m1,
min Xx, 2
aK( exp
-3 Σ
2111)
Σ
CaK((zi
j=mp-1
LEU
j=ma+1
max Xx
Σ
CaK((z, ex) exp
-8 Σ 81(1)
j=mit
j=mz+1
LEU
According equation 27 we then have
T (x,())); Un,) 2 min min X
m
n2
min
Σ
aK((,, exp
-B Σ 3.(l)
Σ
caK(Iz)15)
j=mi+1
LEUs,
j=mg+1
35

---

Published in Transactions on Machine Learning Research (11/2025)
and
T (((;(1);EEz)) S ma max 22.}
S max 0,
Σ ()K(z1)(z)) +=
Σ
(() exp
-
I
IEU
Following equation 26 then
C 2 A3 2 Cz.
Using the similar method, it can be prove that
C3 2 1
Combining with equation 27 conclude our result.
For computing bounds of MMD (X,Y) using Theorem we need to compute 5;()) for any i E {mi (...,,71., 1a (I,...,dd, d}, and compute (I) for any E {ma + I, m2). / (I....,d}. Notice that for any T,(l), we have
(1=;(1) - Zm, (1,) - 2n, (x,))) - yo (o)+++(111-000. (D)...., - yn yr. (DIH ny - m + 1/2 112.
Hence, computing all (1) is of computational order at most O (m n2)). Similarly, computing all (l) is of computational complexity order at most O(d22(11 + n2)). After computing 2,(() for any (1....,d} d], C2,C2, C's and Ca can be computed with computational complexity order at most 0((1 + n2)). Overall, the computational complexity order of computing bounds of MMD (X, (X,Y) using Theorem 6 is O(m + +2))).
E {mi + 1, 1,,,,, E (1...., d), and ū:()) for any i E {m2 + 1....,
A.7 Proof of Theorem 3
This section provides bounds of p-value of MMD test statistic based on the bounds of MMD derived in Scction 8 and the permutation test.
Theorem 2. Suppose X = (II,...,In,I and Y = fyh...., ynat are samples of Rdd real value. Let
(1 ,..,, (Bi) be B i.i.d. random samplings, each being random permulaation of (1...., ni + n2) and denoted as (() (1)..., o + R2)), L....,, Subsequently, let 21 ; Zni Fny 2
V1.... Znitng Vo2 and for any I...., B, denote
X = = (1). (o1) }. Yoi = { (2000m+11)--- (20
and define p according to equation
5 Suppose further
MMD (X,Y) < (X,", MMDi Yom) S MMD (o.o), i= 1=,..., B.
Define
1+ M (XMM Ya(o) 2 MMDD((X)
p B+1
Then, we must have 2 p.
Proof. Notice that for any i E (1....,B), B), we have
T(MMD (x,,,,
Yo, Ya6) 2 MMD =1
MMD (ow) > MMD? (X,Y). (X,).
36

---

Published in Transactions on Machine Learning Research (11/2025)
Since
(X,) MX, MMD (X Yoo)) 2 (X You), i= 1...... B,
it then follows
MMD (X,,, Ya0) 2 MMD", (oo) 2 MM,"
L(MMD (X,,0, x2(0 2 (X, =
for any E [l,...,B, If, for any E BI,...,B),
T(MMD (X,o, Yaa) 2 MMD ()) = 0,
then we have
I(MMD (Xoro, Yo1o (X))) or 1.
Hence, for any i E {I.................. B}. we have
IXMMD ,0,, YO)) 2 MMD 2 (MMM))( (x())))()) MX)
According the definitions of p and p, we then must have p 2 p.
A.8 Proof of Theorem
This section provides bounds value of MMD test statistic, based on the bounds MMD derived Section and the normality approximation proposed Gao & Shao (2023).
Theorem 3. Suppose X = (X1....,XX,, N, = (Yi,...,YY,, E TRd and assume X,...., X and YI...., Ymg are not observed. Let denote the incomplete Laplacion kernel defined in Definition Then, have am 2 Vik Further, suppose have MMD2 (X,Y) such that
0 < MMD S (X,))
have
MMD (X.Y)
(M
V Cn.m n,m
Cn Die
where Ca,m and Vkr are defined in equation and equation separately, and Whe is defined in Definition 5
Proof. To start, recall that it is proved in Lemma that
K(o.v) k(x, v) 0.
where x, y are potentially incompletely observed Rdd samples. Hence, for any 1 S 8, S N, we have
a 2 a 2 at...
which then gives
The 2, > Ak.
37

---

Published in Transactions on Machine Learning Research (11/2025)
where Alg is defined in Proposition 10 in Gao & Shao (2023)
k 8, = d, a. a a, + a.
Subsequently,
D Vike n,m 2
n,m
where Vie is defined in Proposition 10 in Gao & Shao (2023) as
Vaka n.m
(Ah:))
N(N 3)
-N- 1)N 3)
Further, suppose
s (X,) & (X,)
Then, we must have
MMD
(X,))
V Cnam DiC
Cr.m Vkk
n.m
which completes our proof.
38

---

Published in Transactions on Machine Learning Research (11/2025)
B Simulation Details
B.1 Computing Environment and Resources
The experiments were run an high performance computing cluster with 325 compute nodes, cach equippod with 2x AMD EPYC 7742 processors (128 cores, 1TB RAM per node). Each job utilized 32 cores and 64 GB of memory, with compute times varying from a couple of hours to 2-3 days, based on the samples sizes and dimensions of data of the job.
B.2 Simulation Settings
Median Hcuristic. All experiments Scction use median heuristic (Gretton et al., 2012a for the parameter B the Laplacian kernel. This heuristic first calculates the median of all pairwise differences, i.e.
m = med[]] - 2 2,2 { Iny V1. ........ yn. * 2').
and then sets 8 = = 1/m. For the Laplacian kernel, the norm is the 1-norm. This is a popular method of setting 3, and it generally works well (Gretton et al. 2012a). For case deletion and the proposed methods, only completely observed data will be used for computing For mean and hot deck imputation, the 3 will be computed using imputed data.
Permutations. Permutation tests Good (0055 Schrab et al. (2023) used for deletion, mean imputation, and hot deck imputation, with the number of permutations set to B = 100. MMD-Perm also uses B = 100 permutations.
B.3 Other missing data methods
MissForest. MissFooret Stekhoven & Bühlmann 2012 is nonparametric imputation method based on random forests. It iteratively imputes missing values by training random forest on observed data for each variable with missing entries, using the other variables as predictors. The process repeats until the imputed values converge or a stopping criterion is met. In our experiments, we evaluate this mothod using the default parameters recommended in Stekhoven & Stekhoven (2013).
Wilcoxon-Mann- Whitney test with missing data method (WMWM). WMWM Zeng (2024) iss two sample vesting with missing data method based the Bectxoon-Mann-Wiie test statistic. This method rejects the null hypothesis when all possible imputed test statistics are significant. Since WMWM is based on the Wilcoxon-Mann- -in--hhhinn test statistic, it is applicable only for univariate data and can is restricted to detect mean value shift.
Case deletion. For univariate data, case delction uses only the obscrved data for testing. For multivariate data, only data that are completely observed are used.
Mean imputation. For univariate data, mean imputation replaces missing values with the mean of the observed samples in the same group, X Y. For multivariate data, the missing components of incomplete vector are imputed with the mean of the observed components of that same vector.
Hot deck imputation. For univariate data, hot deck imputation imputes missing values with randomly selected observed values (with replacement) from the same group, X or Y. For multivariate data, the missing components of an incomplete vector are imputed with randomly selected observed components (with replacement) from that same vector.
39

---

Published in Transactions on Machine Learning Research (11/2025)
B.4 Examples of MNIST dataset.
In this section, present examples of MNIST images used in Section
2
3
3
3
Figure 4: Examples of MNIST images LeCun ct al. (1998). Each image has dimensions of d x 28. The labels for the images, ordered from left to right and top bottom, 5, 0, 1, 9, 1, 3, 1, 3, 3, 1, and 7.
40

---

Published in Transactions on Machine Learning Research (11/2025)
3
2
3
3
5
Figure 5: Examples of incompletely observed MNIST images, where white pixels represent unobserved values. Each image has dimensions of d = 28 x 28 and is labelled as 3. The first column shows fully observed images; the second column shows images with upper half pixels missing (rows to 14 and columns 8 to 21); the third column shows images after mean imputation; and the fourth column shows images after hot deck imputation. Each row presents different variations of the same original image.
41

---

Published in Transactions on Machine Learning Research (11/2025)
C Additional Experiments
In this section, additional experiments are provided for investigating the power of proposed method. Different sample sizes and alternatives are considered compared to those in Section
C.1
Univariate Data
The missingness mechanism for producing the following two figures is the same the mechanisms used to produce Figure 2
Sample Size: 50
Sample Size: 100
Sample Size: 200
1.00
1.00
1.00
0.75
0.75
0.75
Power
Power
Power
0.50
0.00
0.00
0.00
0.10
0.00
0.05
0:10
c.15
3.20
0.50
0.05
0.10
0.15
0.20
P roportion of Missing Data
P Proportion of Missing Data
Proportion of Missing Data
MMD Miss-
Hot Deck Imputation Mean Imputation
WMWM
Case Deletion
Sig level: 0.05
Figure 6: The power of MMD-Miss and the three common missing data approaches for univariate samples when data are missing not random (MNAR). (Left): Sample size n1 n2 = 50; (Middle): Sample size n1 12 100; (Right): Sample size n1 n2 200. For all figures, significance level 0.05, and alternative hypothesis (0, 1) vs 1), are used. The displayed values are the average times of rejecting the null hypothesis over 100 repetitions.
Sample Size: 50
Sample Size: 100
Sample Size: 200
1.00
1.00
1.00
Power
0.75
Power
0.75
Power
0.75
0.50
0.50
0.50
0.00
0.00
0.00
0.00
0.05
0.10
0.15
0.20
0.00
0.05
0.10
6.15
0.20
0.50
0.05
0.10
0.15
0.20
roportion of Missing Data
Proprtion Missing Data
Proportion of Missing Data
MMD MMDD-Mis--aae -erm
Hot Deck Imputation M Mean Imputation
WMWM
Case Deletion
Sig level: 0.05
Figure 7: The power of MMD-Miss and the three common missing data approaches for univariate samples when data are missing not at random (MNAR). (Left): Sample size n, n2 50; (Middle): Sample size n2 100; (Right): Sample size n1 = 200. For all figures, significance level a = 0.05, and alternative hypothesis (0, 1) vs NO,, 3) are used. The displayed values are the average times of rejecting the null hypothesis over 100 repetitions.
42

---

Published in Transactions on Machine Learning Research (11/2025)
C.2
Multivariate Data
The missingness mechanism for producing the following two figures is the same the mechanisms used to produce Figure g
Sample Size: 50
Sample Size: 100
Sample Size: 200
1.00
1.00
1.00
Power
0.75
Power
0.75
Power
0.75
0.50
0.50
0.00
up qpe Gp
0.00
0.00
0.00
0.05
0.10
0.15
0.20
0.00
0.05
0.10
6.15
0.20
0.50
0.05
0.10
0.15
0.20
P Proportoon of Missing Data
P Proportion of Missing Data
Proportion of Missing Data
MMD MMD-Miss-Perm
Case Deletion
Mean Imputation MissForest
MMD-Miss-Normalle
Hot Deck Imputation
Figure 8: The power of MMD-Miss and the three common missing data approaches for multivariate samples (d = 50) when data missing not at random (MNAR). (Left): Sample size n1 n2 50; (Middle): Sample size n1 n2 100; (Right): Sample size n1 n2 200. For all figures, significance level 0.05, and alternative hypothesis N((O.... 0)² Tso) vs average times of rejecting the null hypothesis over 100 repetitions.
1)" Tso) are used. The displayed values are the
Sample Size: 50
Sample Size: 100
Sample Size: 200
Power
0.75
Power
0.75
Power
0.75
0.50
0:50
0.50
0.00
0.00
0.00
0.00
0.05
0.10
0.15
0.20
0.00
0.05
0.10
0.15
0.20
0.30
0.55
0.10
0.15
0.20
roportion of Missing Data
P Poportion Missing Data
Proportion of Missing Data
IMD-Miss-Pern MMD MMD-Miss-Normahlee
Case Deletion
Mean Imputation MissForest Forest
Hot Deck Imputation
Figure 9: The power of MMD-Miss and the three common missing data approaches for multivariate samples (d = 50) when data are missing not at random (MNAR). (Left): Sample size nı n2 50; (Middle): Sample size nı = 100; (Right): Sample size n2 = 200. For all figures, significance level 0.05, and alternative hypothesis N ((O,...,00) 0)² I5) vs N((O...,, 0)2 9/5) are used. The displayed values are the average times of rejecting the null hypothesis over 100 repetitions.
43


Reviews:
[{"review_id": "2ExHwbOdlr", "reviewer": "Reviewer_EBYb", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "S: Overall, I think that the paper is well-written and the contents are scientifically valid, sufficiently.\n\nS: The proposed method is independent from the mechanism of the missing value generation and I think this is a strong advantage in practice.\n\nS: The derivation of the bound is based on simple reasonable approaches, by which a rationale behind why the proposed method works without assuming missing mechanism is clear.\n\nW: The tightness of the bound is not clear. Although it would be tighter than a naive bound, I guess the bound seems loose. It is unclear how the bound looseness affects the decrease of the power of hypothesis test.\n\nW: Only the Laplacian kernel can be used in the kernel function in MMD, which may be a bit restrictive (e.g., well-known Gaussian kernel is not discussed).", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper considers two-sample test by MMD (Maximum Mean Discrepancy), in which the dataset contains missing values. The basic idea is to derive the bound of MMD values for any miss values, by which the bound of the p-value can also be derived. The proposed method can provides hypothesis test that can appropriately control the Type I error without depending on the mechanism of the missing value generation (such as well-known MCAR, MAR, and MNAR). Two approaches to the null distribution calculation are provided based on the permutation and asymptotic normality, respectively. The empirical evaluation on the synthetic and MNIST datasets show that the proposed method can inhibit Type I error less than significance level.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "About the asymptotic normality, since (8) is seemingly non-negative, it is intuitively a bit weird that this quantity converges to standard normal distribution N(0,1). For example, the value of p = 1 - \\Psi(MMD^2 / \\sqrt{c V}) becomes p = 0.5 when MMD^2 = 0, but this seemingly indicates that p cannot becomes larger than 0.5 because MMD^2 >= 0. However, usually, p-value should be the uniform distribution in [0,1]. How is this inconsistency explained?\n\nAbout the tightness of the bound, for example, the tightest upper bound of MMD^2_u(X,Y) should be \\max_{all missing values} MMD^2_u(X,Y). Although I guess it is computationally difficult, could you provide some discussion about it or a relation between the tightest bound and the current bound? \n\nIn page 23, the authors derived dT(z)/dz. I guess the zero point of this derivative can be calculated analytically, i.e, dT(z)/dz = 0 can be solved wrt z. Then, if the extreme point (solution of dT(z)/dz = 0) exists in (z_3,z_4), it would be the minimum or maximum. Otherwise, The minimum or maximum should exist at z_3 or z_4. Is this logic correct? If so, by calculating it for all the intervals (the same logic would be applicable to any intervals if my guess is correct), can the exact maximum (or minimum) be derived?\n\nIn the case of the permutation approach, the inside of the indicator can be written as MMD^2_u(X_\\sigma(i), Y_\\sigma(i)) - MMD^2_u(X, Y) >= 0. Then, can directly considering the lower bound of the left hand side MMD^2_u(X_\\sigma(i), Y_\\sigma(i)) - MMD^2_u(X, Y) derive a tighter bound?\n\nThe authors first mentioned that the Laplacian and Gaussian kernel are standard, but the technical discussion is only for the Laplacian kernel throughout the paper. The reason why focusing only on Laplacian is not revealed. I think explicitly providing the reason would be beneficial for readers. \n\nIn Figure 1, Type I error of the proposed method is much less than 0.05. Does this indicate the conservative property of the proposed method? It would be better to explicitly discuss why it is not close to 0.05.\n\nMinor: \n\n- I cannot find how \\beta in the Laplacian kernel is set.\n\n- I cannot find the setting of the number of permutation B in MMD-Miss-Perm. \n\n- In the first definition of MMD^2 in page 3, the left had side should not be squared? (or right hand side should be squared). \n\n- In the plots of the power, the line at 0.05 is shown. What is it for? If there's no specific reason, it should be removed to avoid confusion with the Type I error.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "fkW8F1FHkt", "reviewer": "Reviewer_u2qJ", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n- Good paper structure\n- Deals with an interesting and practical subject\n- clear objective and positioning of the paper in relation to the state of the art\n- easy-to-follow mathematical demonstration\n- Satisfactory performance\n\n**Weaknesses**\n- Some methodological choices are not properly justified (choice of kernel, etc.).\n- Principles, advantages and disadvantages need to be highlighted, particularly in relation to literature.\n- Contributions not highlighted\n- experiments are difficult to follow\n- The paper could be improved in terms of clarity and transitions between sections, in order to better guide the reader through the narrative.\n- A few clarifications and modifications to be made (see requested changes section)", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes MMD-Miss, a novel two-sample testing method based on the Maximum Mean Discrepancy (MMD) that is designed to handle missing data in both univariate and multivariate settings, without making assumptions on the missingness mechanism (e.g., MCAR, MAR, MNAR). Unlike traditional approaches that rely on data imputation or deletion, MMD-Miss computes mathematically rigorous bounds on the MMD test statistic by accounting for all possible values of missing data. It extends recent work that derived bounds for the Wilcoxon-Mann-Whitney test in the univariate setting, enabling broader application to any distributional shift and multivariate data.\nThe proposed method:\n* Uses the Laplacian kernel to ensure sensitivity to general distributional differences;\n* Provides a permutation-based test or asymptotic approximation (via Gao & Shao, 2023) to compute p-values;\n* Is provably valid in terms of Type I error control, even when the data are missing not at random (MNAR);\n* Demonstrates strong statistical power (particularly for 5–10% missing data) through numerical simulations;\n* Highlights the failure of standard imputation or deletion techniques to control Type I error under MNAR scenarios, even with low missingness.\nTo the best of the authors' knowledge, MMD-Miss is the first MMD-based two-sample test with formal Type I error guarantees under arbitrary missingness, including in high-dimensional multivariate data.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- Introduction: \"*Nearly all two-sample testing methods are designed solely for data that are fully observed.*\": Which tests are applicable to data with missing values? Have such methods been included in the benchmark? This statement should be substantiated with more detail or appropriate citations.\n- Introduction: paragraph \"two-sample testing ...\": This paragraph in question appears somewhat misplaced. Typically, introductions conclude with a summary of the main contributions and, ideally, an outline of the paper’s structure. Consider moving this paragraph to a more appropriate section (e.g., related work or discussion), or rephrasing it to fit more naturally at the end of the introduction if it is meant to serve as motivation.\n- Related work: While the paper is well-situated within the existing literature on two-sample testing with missing data, a few relevant and recent works seem to be missing from the discussion. In particular, \"Two-sample Testing Using Deep Learning\" [1] and \"A Permutation-Free Kernel Two-Sample Test\" [2] may offer complementary perspectives or contrasting approaches that would strengthen the related work section. Including these references would help position the contribution more clearly with respect to the state of the art.\n\n[1] Kirchler, M., Khorasani, S., Kloft, M., & Lippert, C. (2020, June). Two-sample testing using deep learning. In International Conference on Artificial Intelligence and Statistics (pp. 1387-1398). PMLR.\n\n[2] Shekhar, S., Kim, I., & Ramdas, A. (2022). A permutation-free kernel two-sample test. Advances in Neural Information Processing Systems, 35, 18168-18180.\n- Background: The current structure of the paper presents subsections somewhat in isolation, without sufficient transitions or narrative flow to guide the reader. As this is a scientific article rather than a technical note, it is important to articulate clearer motivations and provide smoother transitions between the different components of the method and experiments. For example, it would be beneficial to briefly explain the purpose of each section and how it builds upon the previous one, making the contributions and rationale behind the methodological choices more accessible and coherent. Reorganizing or enriching the structure with guiding commentary would significantly improve readability and clarity.\n- Maximum Mean Discrepancy:The Maximum Mean Discrepancy (MMD) test is introduced in a technical manner, but the paper would benefit from a more accessible presentation of its core principle, strengths, and limitations, particularly in comparison to classical alternatives such as the Wilcoxon-Mann-Whitney or Kolmogorov-Smirnov tests.\n( MMD² formula : \"n1\" for the second sum\n- \"al., 2024). Two popular choices of characteristic kernels are Laplacian and Gaussian kernels\" : could you justify this point (with references if possible) ?\n- \"Our method is based on the MMD with the Laplacian kernel\" : why? this choice must be more justified (beyond technical reasons for demonstrations)\n- Permutation Test: Please explain the principles of the test and its advantages (non-parametric distribution, etc.).\n- Bounding MMD with Missing Data: The section on bounding the MMD statistic in the presence of missing data is technically sound, but would benefit from a clearer explanation of the rationale behind this approach. Specifically, it is important to guide the reader by explicitly stating why bounding is necessary\n- Lemma1 : \"n1\" for the second sum\n- Multivariate data : Is it possible to have different dimensions for the variables? If so, this is a limitation of the test, which must be indicated.\n- Bounding p-value using normality approximation: In which cases is it preferable to use this approximation?\n- Experiments : \n    - the XP protocols are not very clear and difficult to follow.\n    - Would it be possible to add sensitivities with a missing data rate varying from 0% to 5%? For % rates, it would be interesting to compare the proposed approach with conventional tests.\n    - Why is the method based on the Wilcoxon-Mann-Whitney test (work by Zeng et al.) not added as a competitor in the univariate case?\n    - From my understanding, only independent variables are simulated, in order to test type I error. Would it be possible to play scenarios with dependent variables (in a complex way) to analyze the type II error? d'autant plus qu'il est dit en background \"a two-sample testing method is preferred if it has a lower probability of making a Type II error\"\n    - The imputation methods are not very efficient for evaluating the method. Would it be possible to add more efficient methods (missForest, for example)?\n    - MNIST dataset: \"To assess Type I error, datasets X and Y are generated by randomly sampling with replacement from MNIST training set images labeled as 3.\" Why 3 and not another? Is it possible to carry out the experiment with other labels?\n- Conclusion: The conclusion is far too short. It should include limitations (e.g. choice of Laplace kernel, assumptions in the demonstrations, etc.) as well as perspectives (possibility of extending this approach to the Gaussian kernel? data in large dimensions, etc.).", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "ask0tLWGmt", "reviewer": "Reviewer_uzLL", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n- Unlike prior work which is only univariate and detects only location shifts, this method supports multivariate data and identifies any distributional differences\n- The work has theoretical rigour, with precise bounds on the empirical MMD statistic derived and claims justified theoretically\n- The assumption-free approach is valuable for real-world applications when missingness mechanisms are unknown, and the conservatism could ensure reliability in sensitive domains \n\n\nWeaknesses\n- The authors have relied solely on the Laplacian kernel. Would it be possible to justify the selection and discuss this in the context of other kernels? \n- The simulations could be strengthened with real-world datasets with natural missingness like clinical or survey data\n- Minor typos like section 2  Shutoh et al. (2010) proposed  _an_ method\n\nQuestion to the authors: Is it possible to derive or cite how the bounds ensure Neyman Pearson optimality in the context of Neyman Pearson classification, and/or to evaluate on Neyman Pearson classification benchmarks e.g. loan default or disease diagnosis datasets with induced or natural missingness to further demonstrate practical utility?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors have introduced MMD-Miss, an extension of the Maximum Mean Discrepancy (MMD) two-sample test to handle datasets with arbitrarily missing values, without imposing assumptions on the missingness mechanisms and by leveraging the Laplacian kernel. They derive lower and upper bounds of the MMD statistic with missing values present for univariate and multivariate data. They use these bounds to control for Type I error under arbitrary missingness, employing upper bounds of p-values computed using permutations and with normality approximations. Their simulations on synthetic Gaussian data and MNIST images demonstrate that MMD-Miss controls Type I error under missing data compared to baselines and it maintains good power for 5-10% of missingness.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "See above", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)