Testing Equality of Several Distributions at High Dimensions: A Maximum-Mean-Discrepancy-Based Approach

Abstract

With the development of modern data collection techniques, researchers often encounter high-dimensional data across various research fields. An important problem is to determine whether several groups of these high-dimensional data originate from the same population. To address this, this paper presents a novel k-sample test for equal distributions for high-dimensional data, utilizing the Maximum Mean Discrepancy (MMD). The test statistic is constructed using a V-statistic-based estimator of the squared MMD derived for several samples. The asymptotic null and alternative distributions of the test statistic are derived. To approximate the null distribution accurately, three simple methods are described. To evaluate the performance of the proposed test, two simulation studies and a real data example are presented, demonstrating the effectiveness and reliability of the test in practical applications.

1. Introduction

Testing whether multiple samples follow the same distribution is a fundamental challenge in data analysis, with wide-ranging applications across diverse fields. Traditional nonparametric tests designed for comparing the distributions of two samples, such as the Wald–Wolfowitz runs test, the Mann–Whitney–Wilcoxon test based on signed ranks, and the Kolmogorov–Smirnov test utilizing the Empirical Distribution Function (EDF), are well established for univariate data [1]. Extending these tests to the multivariate setting in ${\mathcal{R}}^p$ has been the focus of extensive research. This has led to the development of novel tests based on multivariate runs, ranks, EDF, distances, and projections, as pioneered by researchers like [2,3,4,5], among others.

In this paper, our primary focus is on addressing the multi-sample problem for equal distributions in high-dimensional data. In contemporary data analysis, high-dimensional datasets have become increasingly prevalent and easily accessible across various domains. For instance, Section 5 introduces a dataset derived from a keratoconus study involving corneal surfaces. This collaborative effort involves Ms. Nancy Tripoli and Dr. Kenneth L. Cohen from the Department of Ophthalmology at the University of North Carolina, Chapel Hill. The dataset comprises 150 corneal surfaces, each characterized by over 6000 height measurements. These surfaces are categorized into four distinct groups based on their corneal shapes, prompting the question of whether these groups share a common underlying distribution. Consequently, there is a pressing need to develop tests that can assess distributional equality in the context of high-dimensional data.

Mathematically, a multi-sample problem for high-dimensional data can be described as follows. Assume that we have the following k samples of observed random elements in ${\mathcal{R}}^p$:

$${\mathit{y}}_{\alpha 1},\cdots ,{\mathit{y}}_{\alpha n_{\alpha }}\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }{F}_{\alpha },\alpha =1,\dots ,k,$$

(1)

where the data dimension p can be very large and $F_{\alpha },\alpha =1,\dots ,k$, are unknown cumulative distribution functions. Of interest is to test if the k distribution functions are the same:

$$H_0:F_1=\cdots =F_k,\mathrm{vs}.H_1:F_{\alpha }\ne F_{\beta } \mathrm{for} \mathrm{some}\alpha \ne \beta .$$

(2)

When $k=2$, a variety of distance-based tests designed for multivariate data, such as those proposed by [2,3,4,5], can potentially be applied to test (2) in the context of high-dimensional data. However, it has been demonstrated by [6] that these tests may lack the power to effectively detect differences in distribution scales in high dimensions. In response to this challenge, ref. [6] introduced a test based on interpoint distances, while [6] proposed a high-dimensional two-sample run test based on the shortest Hamiltonian path. Nevertheless, ref. [7] showed that these tests, along with a graph-based test by [8], are less effective in detecting differences in location. To address this limitation, ref. [7] introduced two asymptotic normal tests for both location and scale differences, based on interpoint distances. However, it is worth noting that these tests rely on a strong mixing condition and impose a natural ordering on the components of high-dimensional observations, limiting their applicability. Furthermore, these tests involve complex U-statistic estimators, making them computationally intensive. Several other approaches have also been proposed for high-dimensional distribution testing. Ref. [9] presented a high-dimensional permutation test using a symmetric measure of distance between data points, but it comes with high computational costs. Refs. [10,11] proposed tests based on projections, which are generally more effective in detecting low-dimensional distributional differences. Ref. [12] introduced a test based on energy distance and permutation, but it is less powerful in detecting scale differences. Refs. [13,14] proposed kernel two-sample tests based on the Maximum Mean Discrepancy (MMD). Ref. [14] demonstrated the equivalence between the energy-distance-based test and the kernel-based test, showing that the energy-distance-based test can be viewed as a kernel test utilizing a kernel induced by the interpoint distance. The MMD leverages the kernel trick to define a distance between the embeddings of distributions in Reproducing Kernel Hilbert Spaces (RKHSs). It is well suited for checking distribution differences among several high-dimensional samples and is applicable to various data types, such as vectors, strings, or graphs. Recently, there has been further investigation into unbiased and biased MMD-based two-sample tests for high-dimensional data by [15,16], respectively.

On the other hand, when $k>2$, there are limited tests available for testing (2). One such test is the energy test developed by [12]. This test statistic is obtained by directly summing all pairwise energy distances, with its null distribution approximated through permutation. However, this test has some drawbacks, including being time-consuming and yielding p-values that may vary when applied multiple times to the same dataset, as supported by the evidence from Figure 3 in Section 5. Another approach is presented by [17], who extended the idea of MMD from the two-sample problem to the multi-sample problem for equal distributions. They constructed an MMD-based test statistic capable of detecting deviations from distribution homogeneity in several samples. However, the resulting test statistic and its null limit distribution are very complicated in form, which restricts its practical use. To address this limitation, ref. [18] developed a new MMD-based test for testing (2). This test statistic is constructed using U-statistics, making it easy to conduct and yielding accurate results.

In this paper, we maintain our focus on the MMD-based approach for testing (2). However, unlike [18], where a U-statistic technique is employed to construct the test statistic, we take a distinct path by constructing an $L^2$-norm-based test statistic. To achieve this, we first employ a canonical feature map to transform the k original samples (1) into k induced samples (3). Concurrently, we transform the k-sample equal distribution testing problem (2) into a mean vector testing problem (4). This transformation facilitates the straightforward construction of an $L^2$-norm-based test (5) for assessing the mean vector testing problem (4). Leveraging a kernel trick, we derive a formula for computing the $L^2$-norm-based test statistic using the k original samples (1).

Additionally, this paper makes several other significant contributions. Firstly, akin to the work of [17,18], we extend the concept of MMD from two-sample problems to the domain of multi-sample problems with equal distributions. Secondly, we derive the asymptotic null and alternative distributions of the proposed test statistic. Thirdly, we offer three distinct approaches for approximating the null distribution of the MMD-based test statistic, utilizing parametric bootstrap, random permutation, and the Welch–Satterthwaite (W–S) ${\chi }^2$-approximation methods. Lastly, we examine two specific scenarios in our comprehensive simulation studies. In the first scenario, the samples have the same mean vector but different covariance matrices, while in the second scenario, the samples exhibit both distinct mean vectors and covariance matrices. Our simulation results demonstrate that the tests we propose effectively maintain precise control over size in both scenarios. However, in terms of empirical power, they outperform (underperform) the energy test introduced by [12] in the first (second) scenario. In other words, when the primary difference in the distributions of the samples lies in their covariance matrices, the new tests are the preferred choice in terms of statistical power.

The remainder of this paper is organized as follows: Section 2 presents the main results and Section 3 introduces three methods to implement our test. In Section 4, we provide two simulation studies. Section 5 showcases an application to the corneal surface data mentioned earlier. Concluding remarks can be found in Section 6. Technical proofs of the main results are included in Appendix A.

2. Main Results

2.1. MMD for Several Distributions

In this section, we show how the MMD can be defined for several distributions. Let $\mathcal{H}$ be an RKHS associated with a characteristic reproducing kernel $K(·,·)$. For any $\mathit{u}$ and $\mathit{v}$ in $\mathcal{H}$, the inner product and $L^2$-norm of $\mathcal{H}$ are defined as $\langle \mathit{u},\mathit{v}\rangle$ and $\parallel \mathit{u}\parallel =\sqrt{\langle \mathit{u},\mathit{u}\rangle }$, respectively. Let $\varphi (·)$ be the canonical feature mapping associated with $K(·,·)$, i.e., $\varphi \left(\mathit{y}\right)=K(·,\mathit{y})$. Using this feature mapping and setting ${\mathit{x}}_{\alpha i}=\varphi \left({\mathit{y}}_{\alpha i}\right)$ for $i=1,\dots ,n_{\alpha }$, we obtain the following k induced samples in the RKHS $\mathcal{H}$:

$${\mathit{x}}_{\alpha 1},\dots ,{\mathit{x}}_{\alpha n_{\alpha }},\alpha =1,\dots ,k,$$

(3)

which are derived from the k original samples in (1). Define ${\mathit{\mu }}_{\alpha }=E\left({\mathit{x}}_{\alpha 1}\right)$ for $\alpha =1,\dots ,k$, representing the mean embeddings of the k distributions $F_{\alpha }$, where $\alpha =1,\dots ,k$.

Ref. [13] established the MMD for two distributions in a separable metric space (see, e.g., Theorem 5 of [13]). Here, we extend it naturally for multi-sample distributions in ${\mathcal{R}}^p$. According to the MMD of [13], for any $\alpha \ne \beta$, where $\alpha ,\beta =1,\dots ,k$, “testing $H_0:F_{\alpha }=F_{\beta }$ vs. $H_1:F_{\alpha }\ne F_{\beta }$” based on the two samples ${\mathit{y}}_{\alpha i},i=1,\dots ,n_{\alpha }$ and ${\mathit{y}}_{\beta i},i=1,\dots ,n_{\beta }$ is equivalent to “testing $H_0:{\mathit{\mu }}_{\alpha }={\mathit{\mu }}_{\beta }$ vs. $H_1:{\mathit{\mu }}_{\alpha }\ne {\mathit{\mu }}_{\beta }$” based on the two samples ${\mathit{x}}_{\alpha i},i=1,\dots ,n_{\alpha }$ and ${\mathit{x}}_{\beta i},i=1,\dots ,n_{\beta }$. Therefore, testing (2) using the k original samples in (1) is equivalent to testing the following hypothesis using the k induced samples in (3):

$$H_0:{\mathit{\mu }}_1=\cdots ={\mathit{\mu }}_k,\mathrm{vs}.H_1:{\mathit{\mu }}_{\alpha }\ne {\mathit{\mu }}_{\beta } \mathrm{for} \mathrm{some} \alpha \ne \beta .$$

(4)

To test (4), following [19], a natural $L^2$-norm-based test statistic using (3) is given by

$$T_n=\sum _{\alpha =1}^k{n}_{\alpha }{\parallel {\overline{\mathit{x}}}_{\alpha }-\overline{\mathit{x}}\parallel }^2,$$

(5)

where ${\overline{\mathit{x}}}_{\alpha }=n_{\alpha }^{-1}{\sum }_{i=1}^{n_{\alpha }}{\mathit{x}}_{\alpha i}$ and $\overline{\mathit{x}}=n^{-1}{\sum }_{\alpha =1}^k{\sum }_{i=1}^{n_{\alpha }}{\mathit{x}}_{\alpha i}$ denote the group and grand sample means, respectively. Through some simple algebra, as given in Appendix A, we can express $T_n$ as

$$T_n=\sum _{\alpha =1}^k\frac{n_{\alpha }(n-n_{\alpha })}{n}\parallel {\overline{\mathit{x}}}_{\alpha }{\parallel }^2-2\sum _{1\le \alpha <\beta \le k}\frac{n_{\alpha }{n}_{\beta }}{n}\langle {\overline{\mathit{x}}}_{\alpha },{\overline{\mathit{x}}}_{\beta }\rangle =\sum _{1\le \alpha <\beta \le k}\frac{n_{\alpha }{n}_{\beta }}{n}{\parallel {\overline{\mathit{x}}}_{\alpha }-{\overline{\mathit{x}}}_{\beta }\parallel }^2.$$

(6)

Let $\mathit{\pi }={({\pi }_1,\dots ,{\pi }_k)}^{\top }$ denote the weights of the k distributions $F_1,\dots ,F_k$ such that ${\pi }_1,\dots ,{\pi }_k\in (0,1)$ and ${\sum }_{\alpha =1}^k{\pi }_{\alpha }=1$. Then, when we estimate $\mathit{\pi }$ using ${(n_1,\dots ,n_k)}^{\top }/n$, the test statistic $T_n/n$ estimates the following quantity:

$${\mathrm{MMD}}^2(F_1,\dots ,F_k|\mathit{\pi })=\sum _{1\le \alpha <\beta \le k}{\pi }_{\alpha }{\pi }_{\beta }{\parallel {\mathit{\mu }}_{\alpha }-{\mathit{\mu }}_{\beta }\parallel }^2,$$

which can be naturally defined as the MMD of the k distributions $F_1,\dots ,F_k$ with the weight vector $\mathit{\pi }$. It is worth noting that the MMD for multiple distributions presented above is equivalent to the one derived by [18], offering a much simpler alternative compared to the formulation proposed by [17].

It is easy to justify that ${\mathrm{MMD}}^2(F_1,\dots ,F_k|\mathit{\pi })$ is indeed an MMD of the k distributions $F_1,\dots ,F_k$ in the sense that ${\mathrm{MMD}}^2(F_1,\dots ,F_k|\mathit{\pi })=0$ if and only if $F_1=\cdots =F_k$. On the one hand, when ${\mathrm{MMD}}^2(F_1,\cdots ,F_k|\mathit{\pi })=0$, we have ${\mathit{\mu }}_1=\cdots ={\mathit{\mu }}_k$ so that for any $1\le \alpha \ne \beta \le k$, we have ${\mathit{\mu }}_{\alpha }={\mathit{\mu }}_{\beta }$, implying that for any $1\le \alpha \ne \beta \le k$, we have $F_{\alpha }=F_{\beta }$ and hence $F_1=\cdots =F_k$. On the other hand, when $F_1=\cdots =F_k$, for any $1\le \alpha \ne \beta \le k$, we have $F_{\alpha }=F_{\beta }$ so that for any $1\le \alpha \ne \beta \le k$, we have ${\mathit{\mu }}_{\alpha }={\mathit{\mu }}_{\beta }$ and hence ${\mathit{\mu }}_1=\cdots ={\mathit{\mu }}_k$. Therefore, we can use $T_n$ to test the equality of k distributions based on the k induced samples (3).

2.2. Computation of the Test Statistic

Notice that the k induced samples (3) are not directly computable, as the canonical feature mapping $\varphi \left(\mathit{y}\right)$ is explicitly defined through the reproducing kernel. Fortunately, the reproducing kernel $K(·,·)$ and its canonical feature mapping $\varphi (·)$ can be utilized with the following useful kernel trick: $K(\mathit{y},{\mathit{y}}^{\prime })=\langle \varphi \left(\mathit{y}\right),\varphi \left({\mathit{y}}^{\prime }\right)\rangle$. Using this, we can express the inner product as follows:

$$\langle {\mathit{x}}_{\alpha i},{\mathit{x}}_{\beta j}\rangle =\langle \varphi \left({\mathit{y}}_{\alpha i}\right),\varphi \left({\mathit{y}}_{\beta j}\right)\rangle =K({\mathit{y}}_{\alpha i},{\mathit{y}}_{\beta j}).$$

(7)

Let $V_{\alpha \alpha }={\parallel {\overline{\mathit{x}}}_{\alpha }\parallel }^2=\langle {\overline{\mathit{x}}}_{\alpha },{\overline{\mathit{x}}}_{\alpha }\rangle$ and $V_{\alpha \beta }=\langle {\overline{\mathit{x}}}_{\alpha },{\overline{\mathit{x}}}_{\beta }\rangle$. Then, using (7), we have

$$V_{\alpha \alpha }=\frac{1}{n_{\alpha }^2}\sum _{i=1}^{n_{\alpha }}\sum _{j=1}^{n_{\alpha }}K({\mathit{y}}_{\alpha i},{\mathit{y}}_{\alpha j}),\mathrm{and}{V}_{\alpha \beta }=\frac{1}{n_{\alpha }{n}_{\beta }}\sum _{i=1}^{n_{\alpha }}\sum _{j=1}^{n_{\beta }}K({\mathit{y}}_{\alpha i},{\mathit{y}}_{\beta j}).$$

Therefore, using (6), we can compute $T_n$ using any of the following useful expressions:

$$T_n=\sum _{\alpha =1}^k\frac{n_{\alpha }(n-n_{\alpha })}{n}{V}_{\alpha \alpha }-2\sum _{1\le \alpha <\beta \le k}\frac{n_{\alpha }{n}_{\beta }}{n}{V}_{\alpha \beta }=\sum _{1\le \alpha <\beta \le k}\frac{n_{\alpha }{n}_{\beta }}{n}(V_{\alpha \alpha }+V_{\beta \beta }-2V_{\alpha \beta }).$$

(8)

In other words, you can compute the value of $T_n$ using the original k samples (1) and the above expressions.

2.3. Asymptotic Null Distribution

To explore the null distribution of $T_n$, we can rewrite $T_n$ as

$$T_n={\tilde{T}}_n+2S_n+\sum _{\alpha =1}^k{n}_{\alpha }{\parallel {\mathit{\mu }}_{\alpha }-\overline{\mathit{\mu }}\parallel }^2,$$

(9)

where $\overline{\mathit{\mu }}=n^{-1}{\sum }_{\alpha =1}^k{n}_{\alpha }{\mathit{\mu }}_{\alpha }$ represents the weighted average of the mean embeddings of the k distributions. Additionally, we have

$${\tilde{T}}_n=\sum _{\alpha =1}^k{n}_{\alpha }{\parallel ({\overline{\mathit{x}}}_{\alpha }-{\mathit{\mu }}_{\alpha })-(\overline{\mathit{x}}-\overline{\mathit{\mu }})\parallel }^2,\mathrm{and}{S}_n=\sum _{\alpha =1}^k{n}_{\alpha }\langle ({\overline{\mathit{x}}}_{\alpha }-{\mathit{\mu }}_{\alpha })-(\overline{\mathit{x}}-\overline{\mathit{\mu }}),{\mathit{\mu }}_{\alpha }-\overline{\mathit{\mu }}\rangle .$$

(10)

In the expression for ${\tilde{T}}_n$, you can observe that the mean embeddings of the k distributions have been subtracted. Consequently, ${\tilde{T}}_n$ follows the same distribution as that of $T_n$ under the null hypothesis. Thus, studying the null distribution of $T_n$ is equivalent to studying the distribution of ${\tilde{T}}_n$.

Similar to the proof of (6), we can express ${\tilde{T}}_n$ as follows:

$$\begin{array}{ccc}\hfill {\tilde{T}}_n& =& \sum _{\alpha =1}^k\frac{n_{\alpha }(n-n_{\alpha })}{n}{\parallel {\overline{\mathit{x}}}_{\alpha }-{\mathit{\mu }}_{\alpha }\parallel }^2-2\sum _{1\le \alpha <\beta \le k}\frac{n_{\alpha }{n}_{\beta }}{n}\langle {\overline{\mathit{x}}}_{\alpha }-{\mathit{\mu }}_{\alpha },{\overline{\mathit{x}}}_{\beta }-{\mathit{\mu }}_{\beta }\rangle \hfill \\ & =& \sum _{1\le \alpha <\beta \le k}\frac{n_{\alpha }{n}_{\beta }}{n}{\parallel ({\overline{\mathit{x}}}_{\alpha }-{\mathit{\mu }}_{\alpha })-({\overline{\mathit{x}}}_{\beta }-{\mathit{\mu }}_{\beta })\parallel }^2.\hfill \end{array}$$

(11)

Let $\tilde{K}(\mathit{u},\mathit{v})$ denote the centered version of $K(\mathit{u},\mathit{v})$, defined as

$$\begin{array}{ccc}\hfill \tilde{K}(\mathit{u},\mathit{v})& =& \langle \varphi \left(\mathit{u}\right)-{\mathit{\mu }}_u,\varphi \left(\mathit{v}\right)-{\mathit{\mu }}_v\rangle \hfill \\ \hfill & =& K(\mathit{u},\mathit{v})-E_{{\mathit{v}}^{\prime }}\left[K(\mathit{u},{\mathit{v}}^{\prime })\right]-E_{{\mathit{u}}^{\prime }}\left[K({\mathit{u}}^{\prime },\mathit{v})\right]+E_{{\mathit{u}}^{\prime },{\mathit{v}}^{\prime }}\left[K({\mathit{u}}^{\prime },{\mathit{v}}^{\prime })\right],\hfill \end{array}$$

(12)

where ${\mathit{\mu }}_u=E\left[\varphi \left(\mathit{u}\right)\right]$, ${\mathit{\mu }}_v=E\left[\varphi \left(\mathit{v}\right)\right]$, and ${\mathit{u}}^{\prime }$ and ${\mathit{v}}^{\prime }$ are independent copies of $\mathit{u}$ and $\mathit{v}$, respectively. We can observe two useful properties: when $\mathit{u}=\mathit{v}$, we have

$$E_{\mathit{u}}\left[\tilde{K}(\mathit{u},\mathit{u})\right]=E_{\mathit{u}}{\parallel \varphi \left(\mathit{u}\right)-{\mathit{\mu }}_u\parallel }^2>0.$$

(13)

When $\mathit{u}$ and $\mathit{v}$ are independent, we have

$$E_{\mathit{u}}\left[\tilde{K}(\mathit{u},\mathit{v})\right]=E_{\mathit{v}}\left[\tilde{K}(\mathit{u},\mathit{v})\right]=E_{\mathit{u},\mathit{v}}\left[\tilde{K}(\mathit{u},\mathit{v})\right]=0.$$

(14)

Using (12), we can express

$$\tilde{K}({\mathit{y}}_{\alpha i},{\mathit{y}}_{\beta j})=\langle {\mathit{x}}_{\alpha i}-{\mathit{\mu }}_{\alpha },{\mathit{x}}_{\beta j}-{\mathit{\mu }}_{\beta }\rangle .$$

(15)

Let

$${\tilde{V}}_{\alpha \alpha }={\parallel {\overline{\mathit{x}}}_{\alpha }-{\mathit{\mu }}_{\alpha }\parallel }^2,\mathrm{and}{\tilde{V}}_{\alpha \beta }=\langle {\overline{\mathit{x}}}_{\alpha }-{\mathit{\mu }}_{\alpha },{\overline{\mathit{x}}}_{\beta }-{\mathit{\mu }}_{\beta }\rangle .$$

(16)

It is evident that ${\tilde{V}}_{\alpha \alpha }$ and ${\tilde{V}}_{\alpha \beta }$ can be considered as centered versions of $V_{\alpha \alpha }$ and $V_{\alpha \beta }$, respectively. Consequently, by using (11) and (16), we can express

$${\tilde{T}}_n=\sum _{\alpha =1}^k\frac{n_{\alpha }(n-n_{\alpha })}{n}{\tilde{V}}_{\alpha \alpha }-2\sum _{1\le \alpha <\beta \le k}\frac{n_{\alpha }{n}_{\beta }}{n}{\tilde{V}}_{\alpha \beta }.$$

(17)

Utilizing (15) and performing some straightforward algebraic manipulations, we can express

$$\begin{array}{ccc}\hfill {\tilde{V}}_{\alpha \alpha }& =& \frac{1}{n_{\alpha }^2}\sum _{i=1}^{n_{\alpha }}\sum _{j=1}^{n_{\alpha }}\tilde{K}({\mathit{y}}_{\alpha i},{\mathit{y}}_{\alpha j})=\frac{1}{n_{\alpha }^2}\left[\sum _{i=1}^{n_{\alpha }}\tilde{K}({\mathit{y}}_{\alpha i},{\mathit{y}}_{\alpha i})+2\sum _{1\le i<j\le n_{\alpha }}\tilde{K}({\mathit{y}}_{\alpha i},{\mathit{y}}_{\alpha j})\right],\mathrm{and}\hfill \\ \hfill {\tilde{V}}_{\alpha \beta }& =& \frac{1}{n_{\alpha }{n}_{\beta }}\sum _{i=1}^{n_{\alpha }}\sum _{j=1}^{n_{\beta }}\tilde{K}({\mathit{y}}_{\alpha i},{\mathit{y}}_{\beta j}),\alpha \ne \beta .\hfill \end{array}$$

(18)

Assuming that $\tilde{K}(\mathit{y},{\mathit{y}}^{\prime })$ is square-integrable, i.e., $E\left[{\tilde{K}}^2(\mathit{y},{\mathit{y}}^{\prime })\right]<\infty$, we can express $\tilde{K}(\mathit{y},{\mathit{y}}^{\prime })$ using Mercer’s expansion:

$$\tilde{K}(\mathit{y},{\mathit{y}}^{\prime })=\sum _{r=1}^{\infty }{\lambda }_r{\psi }_r\left(\mathit{y}\right){\psi }_r\left({\mathit{y}}^{\prime }\right),$$

(19)

where ${\lambda }_1,{\lambda }_2,\dots$ represent the eigenvalues of $\tilde{K}(\mathit{y},{\mathit{y}}^{\prime })$, and ${\psi }_1\left(\mathit{y}\right),{\psi }_2\left(\mathit{y}\right),\dots$ are the corresponding orthonormal eigenelements, satisfying

$$\int \tilde{K}(\mathit{y},{\mathit{y}}^{\prime }){\psi }_r\left(\mathit{y}\right)dF\left(\mathit{y}\right)={\lambda }_r{\psi }_r\left({\mathit{y}}^{\prime }\right),\mathrm{and}\int {\psi }_r\left(\mathit{y}\right){\psi }_s\left(\mathit{y}\right)dF\left(\mathit{y}\right)={\delta }_{rs},r,s=1,2,\dots ,$$

(20)

where ${\delta }_{rs}$ equals 1 when $r=s$ and 0 otherwise. Now, let us introduce the following conditions:

C1.

We have ${\mathit{y}}_{\alpha 1},\dots ,{\mathit{y}}_{\alpha ,n_{\alpha }},\alpha =1,\cdots ,k\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }F$.

C2.

As $n\to \infty$, we have $n_{\alpha }/n\to {\tau }_{\alpha }\in (0,1),\alpha =1,\dots ,k$.

C3.

$K(\mathit{y},{\mathit{y}}^{\prime })$ is a reproduced kernel such that $E_{\mathit{y}}\left[{\tilde{K}}^2(\mathit{y},\mathit{y})\right]\le \infty$.

Condition C1 assumes that the null hypothesis is satisfied and the common distribution function is F. Condition C2 is a regularity condition for k-sample problems and it requires that the group sample sizes tend to ∞ proportionally. Condition C3 is required such that $\tilde{K}(\mathit{y},{\mathit{y}}^{\prime })$ is square-integrable and expression (19) is valid.

In fact, under Condition C3, using (19) and the Cauchy–Schwarz inequality, we obtain the following results:

$$\begin{array}{cc}\hfill E\left[\tilde{K}(\mathit{y},\mathit{y})\right]& =\sum _{r=1}^{\infty }{\lambda }_r<\sqrt{E\left[{\tilde{K}}^2(\mathit{y},\mathit{y})\right]}<\infty ,\hfill \\ \hfill E_{\mathit{y},{\mathit{y}}^{\prime }}{\left[\tilde{K}(\mathit{y},{\mathit{y}}^{\prime })\right]}^2& =\sum _{r=1}^{\infty }{\lambda }_r^2<E_{\mathit{y},{\mathit{y}}^{\prime }}\left[\tilde{K}(\mathit{y},\mathit{y})\tilde{K}({\mathit{y}}^{\prime },{\mathit{y}}^{\prime })\right]=E_{\mathit{y}}^2\left[\tilde{K}(\mathit{y},\mathit{y})\right]<\infty ,\hfill \end{array}$$

(21)

where $\mathit{y},{\mathit{y}}^{\prime }\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }F$. These inequalities hold due to the square-integrability assumption and the properties of Mercer’s expansion. Now, let us state the following theorem that establishes the asymptotic distribution of ${\tilde{T}}_n$.

Theorem 1.

Under Conditions C1–C3, as $n\to \infty$, we have ${\tilde{T}}_n\stackrel{\mathcal{L}}{⟶}\tilde{T}$, where

$$\tilde{T}\stackrel{d}{=}\sum _{r=1}^{\infty }{\lambda }_r{A}_r,A_r\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }{\chi }_{k-1}^2.$$

(22)

It is worth highlighting that the limit null distribution of the proposed test statistic differs from the one derived in [18] (Theorem 1) and it offers a more straightforward alternative compared to the limit null distribution obtained by [17] (Theorem 1). This explains why the limit null distribution presented by [17] is not employed to approximate the null distribution of their test statistic. However, as demonstrated in Section 3, it is indeed feasible to utilize this distribution if desired.

2.4. Mean and Variance of ${\tilde{T}}_n$

Based on (13), (14) and (18), through some simple calculation, under Condition C1, we have

$$\begin{array}{ccc}& & E\left[{\tilde{V}}_{\alpha \alpha }\right]=\frac{1}{n_{\alpha }}E\left[\tilde{K}(\mathit{y},\mathit{y})\right],E\left[{\tilde{V}}_{\alpha \beta }\right]=0,\hfill \\ & & Var\left[{\tilde{V}}_{\alpha \alpha }\right]=\frac{1}{n_{\alpha }^3}Var\left[\tilde{K}(\mathit{y},\mathit{y})\right]+\frac{2(n_{\alpha }-1)}{n_{\alpha }^3}E\left[{\tilde{K}}^2(\mathit{y},{\mathit{y}}^{\prime })\right],\hfill \\ & & Var\left[{\tilde{V}}_{\alpha \beta }\right]=\frac{1}{n_{\alpha }{n}_{\beta }}E\left[{\tilde{K}}^2(\mathit{y},{\mathit{y}}^{\prime })\right],\mathrm{and}Cov\left({\tilde{V}}_{\alpha \alpha },{\tilde{V}}_{\alpha \beta }\right)=0,\alpha \ne \beta ,\hfill \end{array}$$

(23)

where $\mathit{y},{\mathit{y}}^{\prime }\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }F$.

Theorem 2.

Under Condition C1, we have $E\left({\tilde{T}}_n\right)=(k-1)E\left[\tilde{K}(\mathit{y},\mathit{y})\right]$, and

$$Var\left({\tilde{T}}_n\right)=\left[\sum _{\alpha =1}^k\frac{{(n-n_{\alpha })}^2}{n^2{n}_{\alpha }}\right]Var\left[\tilde{K}(\mathit{y},\mathit{y})\right]+2\left[(k-1)-\sum _{\alpha =1}^k\frac{{(n-n_{\alpha })}^2}{n^2{n}_{\alpha }}\right]E\left[{\tilde{K}}^2(\mathit{y},{\mathit{y}}^{\prime })\right],$$

where $\mathit{y},{\mathit{y}}^{\prime }\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }F$.

Note that under Condition C2, we have ${\sum }_{\alpha =1}^k{(n-n_{\alpha })}^2/\left(n^2{n}_{\alpha }\right)\le {\sum }_{\alpha =1}^k{n}_{\alpha }^{-1}\to 0$ as $n\to \infty$. Then as $n\to \infty$, we have $Var\left({\tilde{T}}_n\right)=2(k-1)E\left[{\tilde{K}}^2(\mathit{y},{\mathit{y}}^{\prime })\right][1+o\left(1\right)]$.

2.5. Asymptotic Power

In this subsection, we examine the asymptotic power of the proposed test under the following local alternative hypothesis:

$$H_{1n}:{\mathit{\mu }}_{\alpha }=\mathit{\mu }+n^{-(1/2-\Delta )}{\mathit{h}}_{\alpha },\alpha =1,\dots ,k,$$

(24)

where $0<\Delta \le 1/2$ and ${\mathit{h}}_{\alpha },\alpha =1,\dots ,k$ are constant elements in the RKHS $\mathcal{H}$ such that

$$\begin{array}{c}\hfill 0<\sum _{\alpha =1}^k\frac{n_{\alpha }}{n}{\parallel {\mathit{h}}_{\alpha }-\overline{\mathit{h}}\parallel }^2\le B_h<\infty \mathrm{with} \overline{\mathit{h}}=\frac{1}{n}\sum _{\alpha =1}^k{n}_{\alpha }{\mathit{h}}_{\alpha },\end{array}$$

and

$$\sum _{\alpha =1}^k{\sigma }_{\alpha }^2>0,\mathrm{where}{\sigma }_{\alpha }^2=E{\left[\langle {\mathit{x}}_{\alpha 1}-{\mathit{\mu }}_{\alpha },{\mathit{h}}_{\alpha }-\overline{\mathit{h}}\rangle \right]}^2,\alpha =1,\dots ,k.$$

It is important to note that when $\Delta =1/2$, the local hypothesis (24) simplifies into a fixed alternative hypothesis $H_1:{\mathit{\mu }}_{\alpha }=\mathit{\mu }+{\mathit{h}}_{\alpha }$ for $\alpha =1,\dots ,k$. However, when $0<\Delta <1/2$, Equation (24) represents a strict local alternative hypothesis. In this case, as $n\to \infty$, the strict local hypothesis tends to converge toward the null hypothesis. Detecting a strict local alternative hypothesis becomes exceedingly challenging in such scenarios. A test is typically considered root-n-consistent if it can detect a strict local alternative hypothesis with a probability approaching 1 as $n\to \infty$. A root-n-consistent test is considered effective because it achieves the best possible detection rate for a local alternative hypothesis as n grows.

Through (9) and (24), $T_n$ can be written as

$$T_n={\tilde{T}}_n+2S_n+n^{2\Delta }\sum _{\alpha =1}^k\frac{n_{\alpha }}{n}{\parallel {\mathit{h}}_{\alpha }-\overline{\mathit{h}}\parallel }^2,$$

(25)

where

$$S_n=n^{-(1/2-\Delta )}\sum _{\alpha =1}^k{n}_{\alpha }\langle {\overline{\mathit{x}}}_{\alpha }-{\mathit{\mu }}_{\alpha },{\mathit{h}}_{\alpha }-\overline{\mathit{h}}\rangle .$$

Theorem 3.

Assume that $|K(\mathit{y},{\mathit{y}}^{\prime })|\le B_K$ for all $\mathit{y},{\mathit{y}}^{\prime }\in {\mathcal{R}}^p$ for some $B_K<\infty$. Then, we have $Var\left({\tilde{T}}_n\right)\le 32(k-1)B_K^2$ and $Var\left(S_n\right)=n^{2\Delta }{\sum }_{\alpha =1}^k{n}_{\alpha }{\sigma }_{\alpha }^2/n$, where ${\sum }_{\alpha =1}^k{n}_{\alpha }{\sigma }_{\alpha }^2/n\le 4B_K{B}_h<\infty$.

Theorem 4.

Assume that $|K(\mathit{y},{\mathit{y}}^{\prime })|\le B_K$ for all $\mathit{y},{\mathit{y}}^{\prime }\in {\mathcal{R}}^p$ for some $B_K<\infty$. Then, under Condition C2 and the local alternative hypothesis (24), as $n\to \infty$, we have (a) ${\tilde{T}}_n/\sqrt{Var\left(S_n\right)}\stackrel{P}{⟶}0$; (b) $S_n/\sqrt{Var\left(S_n\right)}\stackrel{\mathcal{L}}{⟶}\mathcal{N}(0,1)$; and (c)

$$\begin{array}{c}\hfill Pr(T_n\ge C_{ϵ})=\Phi \left[\frac{n^{\Delta }{\sum }_{\alpha =1}^k{\tau }_{\alpha }{\parallel {\mathit{h}}_{\alpha }-{\mathit{h}}^{*}\parallel }^2}{2\sqrt{{\sum }_{\alpha =1}^k{\tau }_{\alpha }{\sigma }_{\alpha }^{*2}}}\right][1+o\left(1\right)]⟶1,\end{array}$$

where $C_{ϵ}$ denotes the upper $100ϵ$ percentile of ${\tilde{T}}_n$ with ϵ being the given significance level, ${\tau }_{\alpha }$; $\alpha =1,\dots ,k$ are defined in Condition C2; ${\mathit{h}}^{*}={\sum }_{\alpha =1}^k{\tau }_{\alpha }{\mathit{h}}_{\alpha };{\sigma }_{\alpha }^{*2}=E{\left[\langle {\mathit{x}}_{\alpha 1}-{\mathit{\mu }}_{\alpha },{\mathit{h}}_{\alpha }-{\mathit{h}}^{*}\rangle \right]}^2$; and $\Phi (·)$ denotes the cumulative distribution of $\mathcal{N}(0,1)$.

Theorem 4 shows that the proposed test is indeed a root-n-consistent test.

3. Methods for Implementing the Proposed Test

In this section, we will outline three different approaches for approximating the null distribution of $T_n$ (8) in order to conduct the proposed test. These methods include a parametric bootstrap approach, a random permutation technique, and a ${\chi }^2$-approximation method. We will evaluate and compare their performance in the next section.

3.1. Parametric Bootstrap Method

Theorem 1 reveals that the asymptotic null distribution of $T_n$ takes the form of a chi-squared-type mixture denoted as $\tilde{T}$ (22). The coefficients of this mixture are determined by the unknown eigenvalues ${\lambda }_r,r=1,2,\dots$ of $\tilde{K}(\mathit{y},{\mathit{y}}^{\prime })$, where $\mathit{y},{\mathit{y}}^{\prime }$ are independently and identically distributed according to the common distribution function F representing the k distributions when the null hypothesis is valid. Consequently, in order to estimate the asymptotic null distribution of $T_n$, it is essential to consistently estimate ${\lambda }_r,r=1,2,\dots$. This consistency can be achieved by utilizing the empirical eigenvalues of the centered Gram matrix, as suggested by [20], to construct a reliable estimator for $\tilde{T}$ (22).

Let us recall that $n={\sum }_{\alpha =1}^k{n}_{\alpha }$ represents the total sample size. We pool the k samples (1) and denote it as

$${\mathit{y}}_1,\dots ,{\mathit{y}}_n.$$

(26)

Under the null hypothesis, ${\mathit{y}}_1,\dots ,{\mathit{y}}_n$ are independently and identically distributed from the common distribution F. Let $\mathit{K}$ represent the $n\times n$ Gram matrix, where the $(i,j)$th entry is defined as $K({\mathit{y}}_i,{\mathit{y}}_j)$ for $i,j=1,\dots ,n$. Additionally, let ${\mathbf{1}}_n$ denote a vector of ones with dimensions $n\times 1$, and ${\mathit{I}}_n$ denote the identity matrix of size $n\times n$. Then, the matrix ${\mathit{P}}_n={\mathit{I}}_n-{\mathbf{1}}_n{\mathbf{1}}_n^{\top }/n$ is a projection matrix of size $n\times n$ with rank $n-1$.

Now, define ${\tilde{\mathit{K}}}_n={\mathit{P}}_n\mathit{K}{\mathit{P}}_n$, commonly referred to as the centered Gram matrix. Its $(i,j)$th entry is given by

$${\tilde{K}}_n({\mathit{y}}_i,{\mathit{y}}_j)=K({\mathit{y}}_i,{\mathit{y}}_j)-n^{-1}\sum _{v=1}^{n}K({\mathit{y}}_i,{\mathit{y}}_v)-n^{-1}\sum _{u=1}^{n}K({\mathit{y}}_u,{\mathit{y}}_j)+n^{-2}\sum _{u=1}^n\sum _{v=1}^{n}K({\mathit{y}}_u,{\mathit{y}}_v),$$

(27)

where $i,j=1,\dots ,n$. As n approaches infinity, for any fixed i and j, we can observe that, by the law of large numbers:

$${\tilde{K}}_n({\mathit{y}}_i,{\mathit{y}}_j)\stackrel{\mathcal{L}}{⟶}\tilde{K}(y_i,y_j)=K({\mathit{y}}_i,{\mathit{y}}_j)-E_{{\mathit{y}}^{\prime }}\left[K({\mathit{y}}_i,{\mathit{y}}^{\prime })\right]-E_{\mathit{y}}\left[K(\mathit{y},{\mathit{y}}_j)\right]+E_{\mathit{y},{\mathit{y}}^{\prime }}\left[K(\mathit{y},{\mathit{y}}^{\prime })\right].$$

Let ${\varpi }_1,\dots ,{\varpi }_q$ be all the non-zero eigenvalues of ${\tilde{\mathit{K}}}_n$ that can be obtained via an eigen-decomposition of ${\tilde{\mathit{K}}}_n$. Set ${\widehat{\lambda }}_r={\varpi }_r/n,r=1,\dots ,q$. Then, following [20], we can show that under Condition C3, the distribution of $\tilde{T}$ can be consistently estimated by

$$\widehat{\tilde{T}}=\sum _{r=1}^q{\widehat{\lambda }}_r{A}_r,A_r\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }{\chi }_{k-1}^2.$$

(28)

The parametric bootstrap method can be described as follows. Let us choose a large value for N, for example, $N=1000$. Using expression (28), we can obtain a sample ${\widehat{\tilde{T}}}^{\left(1\right)},\dots ,{\widehat{\tilde{T}}}^{\left(N\right)}$ of $\widehat{\tilde{T}}$ by independently generating $A_r,r=1,\dots ,q$ a total of N times. Now, let $T_{obs}$ represent the observed test statistic calculated using (8) based on the k samples (1). Using the parametric bootstrap method, we can conduct the proposed test by calculating the approximate p-value, which is given by $N^{-1}{\sum }_{i=1}^{N}I\{{\widehat{\tilde{T}}}^{\left(i\right)}\ge T_{obs}\}$, where $I\left\{S\right\}$ is an indicator function that takes 1 when S is a true event and 0 otherwise.

3.2. Random Permutation Method

We can also approximate the null distribution of $T_n$ using a random permutation method. Let ${\ell }_1,\dots ,{\ell }_n$ represent a random permutation of the indices $1,\dots ,n$ from the pooled sample (26). Consequently, the sequence

$${\mathit{y}}_{{\ell }_1},\dots ,{\mathit{y}}_{{\ell }_n},$$

(29)

forms a permutation of the pooled sample (26). To create permuted samples, we utilize the first $n_1$ observations in the permuted pooled sample (29) as the first permuted sample, the next $n_2$ observations as the second permuted sample, and so on, until we obtain k permuted samples. These permuted samples are denoted as

$${\mathit{y}}_{\alpha i}^{*},i=1,\dots ,n_{\alpha };\alpha =1,\dots ,k.$$

(30)

The permutated test statistic, denoted as $T_n^{*}$, is calculated using (8) but with the k samples (1) replaced by the k permuted samples (30).

The random permutation method proceeds as follows. Let N be a sufficiently large number, for instance, $N=1000$. Suppose we repeat the permutation process described above N times, resulting in N permutated test statistics denoted as $T_n^{*\left(i\right)},i=1,\dots ,N$. Then, we can use the empirical distribution of $T_n^{*\left(i\right)},i=1,\dots ,N$ to approximate the null distribution of $T_n$. Recall that $T_{obs}$ represents the test statistic computed using (8) based on the k original samples (1). Following the random permutation method, the proposed test can be conducted by calculating the approximated p-value, given by $N^{-1}{\sum }_{i=1}^{N}I\{T_n^{*\left(i\right)}\ge T_{obs}\}$.

3.3. Welch–Satterthwaite ${\chi }^2$-Approximation Method

The parametric bootstrap method and the random permutation method are effective for controlling size but can be computationally intensive, particularly with large total sample sizes. To address this issue, we can utilize the well-known Welch–Satterthwaite (W–S) ${\chi }^2$-approximation method [21,22]. This method is known to be reliable for approximating the distribution of a chi-squared-type mixture. Theorem 1 demonstrates that the asymptotic null distribution of $T_n$ is a chi-squared-type mixture $\tilde{T}$ (22).

The core concept of the W–S ${\chi }^2$-approximation method is to approximate the null distribution of $T_n$ using that of a random variable of the form: $W\stackrel{d}{=}\beta {\chi }_d^2$, where $\beta$ and d are unknown parameters. These parameters can be determined by matching the means and variances of ${\tilde{T}}_n$ and W, where ${\tilde{T}}_n$, defined in (10), has the same distribution as $T_n$ under the null hypothesis. Specifically, the mean and variance of W are $\beta d$ and $2{\beta }^{2}d$, respectively, while the mean and variance of ${\tilde{T}}_n$ are given in Theorem 2. Equating the means and variances of ${\tilde{T}}_n$ and W, we obtain

$$\beta =\frac{Var\left({\tilde{T}}_n\right)}{2E\left({\tilde{T}}_n\right)}\mathrm{and}d=\frac{2E^2\left({\tilde{T}}_n\right)}{Var\left({\tilde{T}}_n\right)}.$$

(31)

To implement the W–S ${\chi }^2$-approximation method, we need to consistently estimate $E\left({\tilde{T}}_n\right)$ and $Var\left({\tilde{T}}_n\right)$ based on the pooled sample (26) from the k samples (1). According to Theorem 2, these estimates can be obtained as follows:

$$\begin{array}{ccc}\hfill \widehat{E}\left({\tilde{T}}_n\right)& =& (k-1)\widehat{E}\left[\tilde{K}(\mathit{y},\mathit{y})\right],\mathrm{and}\hfill \\ \hfill \widehat{Var}\left({\tilde{T}}_n\right)& =& \left[\sum _{\alpha =1}^k\frac{{(n-n_{\alpha })}^2}{n^2{n}_{\alpha }}\right]\widehat{Var}\left[\tilde{K}(\mathit{y},\mathit{y})\right]+2\left[(k-1)-\sum _{\alpha =1}^k\frac{{(n-n_{\alpha })}^2}{n^2{n}_{\alpha }}\right]\widehat{E}\left[{\tilde{K}}^2(\mathit{y},{\mathit{y}}^{\prime })\right],\hfill \end{array}$$

where

$$\begin{array}{ccc}\hfill \widehat{E}\left[\tilde{K}(\mathit{y},\mathit{y})\right]& =& n^{-1}\sum _{i=1}^n{\tilde{K}}_n({\mathit{y}}_i,{\mathit{y}}_i),\hfill \\ \hfill \widehat{Var}\left[\tilde{K}(\mathit{y},\mathit{y})\right]& =& {(n-1)}^{-1}\sum _{i=1}^n{\left\{{\tilde{K}}_n({\mathit{y}}_i,{\mathit{y}}_i)-n^{-1}\sum _{j=1}^n{\tilde{K}}_n({\mathit{y}}_j,{\mathit{y}}_j)\right\}}^2,\mathrm{and}\hfill \\ \hfill \widehat{E}\left[{\tilde{K}}^2(\mathit{y},{\mathit{y}}^{\prime })\right]& =& \frac{2}{n(n-1)}\sum _{1\le i<j\le n}{\left[{\tilde{K}}_n({\mathit{y}}_i,{\mathit{y}}_j)\right]}^2,\hfill \end{array}$$

with ${\tilde{K}}_n({\mathit{y}}_i,{\mathit{y}}_j)$ being defined in (27). Substituting these estimators into (31), we obtain

$$\widehat{\beta }=\frac{\widehat{Var}\left({\tilde{T}}_n\right)}{2\widehat{E}\left({\tilde{T}}_n\right)}\mathrm{and}\widehat{d}=\frac{2{\left(\widehat{E}\left({\tilde{T}}_n\right)\right)}^2}{\widehat{Var}\left({\tilde{T}}_n\right)}.$$

Let ${\chi }_{\widehat{d}}^2\left({\alpha }^{*}\right)$ denote the upper $100{\alpha }^{*}$ percentile of ${\chi }_{\widehat{d}}^2$, where ${\alpha }^{*}$ is the given significance level, and let $T_{obs}$ denote the observed test statistic computed using (8) based on the k samples (1). Then, through the W–S ${\chi }^2$-approximation method, the proposed test can be conducted via rejecting the null hypothesis when $T_{obs}>\widehat{\beta }{\chi }_{\widehat{d}}^2\left({\alpha }^{*}\right)$ or when the approximated p-value $Pr\left({\chi }_{\widehat{d}}^2>T_{obs}/\widehat{\beta }\right)$ is less than ${\alpha }^{*}$.

4. Simulation Studies

In this section, we delve into intensive simulation studies aimed at assessing the performance of the test we propose when compared to the energy test introduced by [12], which we denote as $T_{SR}$. Our proposed test employs the parametric bootstrap, the random permutation, and the W–S ${\chi }^2$-approximation methods as described in Section 3. For simplicity, we refer to the resulting tests as $T_{PB}$, $T_{RP}$, and $T_{WS}$, respectively.

For simplicity, we opt for the Gaussian Radial Basis Function (RBF) kernel, denoted as $K(·,·)$, which is defined as follows:

$$K(\mathit{y},{\mathit{y}}^{\prime })=exp\left(-\frac{\parallel \mathit{y}-{\mathit{y}}^{\prime }{\parallel }^2}{2{\sigma }^2}\right).$$

Here, ${\sigma }^2$ is referred to as the kernel width. It is worth noting that the Gaussian RBF kernel described above is bounded by 1, ensuring that Condition C3 is always met. Following the approach outlined in [20], we set $\sigma$ to be equal to the median distance between observed vectors in the pooled sample (26).

We also employ the Average Relative Error (ARE) introduced by [23] to evaluate the overall effectiveness of a test in maintaining its nominal size. The ARE is calculated as follows: $ARE=100M^{-1}{\sum }_{i=1}^M|{\widehat{{\alpha }_i}}^{*}-{\alpha }^{*}|/{\alpha }^{*}$, where ${\widehat{{\alpha }_i}}^{*},i=1,\dots ,M$ represents the empirical sizes observed across M different simulation settings. A smaller ARE value indicates the better performance of a test in terms of size control. In this simulation study, we set the nominal size to ${\alpha }^{*}=5\%$.

4.1. Simulation 1

We set $k=3$ for simplicity. We generate the $k=3$ samples (1) as follows. We set

$$\begin{array}{ccc}\hfill {\mathit{y}}_{1i}& =& \mathit{\mu }+\Gamma ({\mathit{u}}_{1i}+{\delta }_1{\mathit{v}}_{1i}),i=1,\dots ,n_1,\hfill \\ \hfill {\mathit{y}}_{2i}& =& \mathit{\mu }+\Gamma ({\mathit{u}}_{2i}+{\delta }_2{\mathit{v}}_{2i}),i=1,\dots ,n_2,\hfill \\ \hfill {\mathit{y}}_{3i}& =& \mathit{\mu }+\Gamma {\mathit{u}}_{3i},i=1,\dots ,n_3,\hfill \end{array}$$

where ${\mathit{u}}_{\alpha i},i=1,\dots ,n_{\alpha };\alpha =1,2,3\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }{\mathcal{N}}_p(\mathbf{0},{\mathit{I}}_p)$, while ${\mathit{v}}_{\alpha i}={(v_{\alpha ,i1},\dots ,v_{\alpha ,ip})}^{\top }$, $i=1,\dots ,n_{\alpha };\alpha =1,2$ are generated using the following three models:

Model 1.

$v_{\alpha ,ir},i=1,\dots ,n_{\alpha };\alpha =1,2;r=1,\dots ,p\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }\mathcal{N}(0,1)$.

Model 2.

$v_{\alpha ,ir}=u_{\alpha ,ir}/\sqrt{2},i=1,\dots ,n_{\alpha };\alpha =1,2;r=1,\dots ,p$ with $u_{\alpha ,ir},i=1,\dots ,n_{\alpha }$; $\alpha =1,2;r=1,\dots ,p\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }{t}_4$.

Model 3.

$v_{\alpha ,ir}=(u_{\alpha ,ir}-1)/\sqrt{2},i=1,\dots ,n_{\alpha };\alpha =1,2;r=1,\dots ,p$ with $u_{\alpha ,ir},i=1,\dots ,n_{\alpha };$$\alpha =1,2;r=1,\dots ,p\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }{\chi }_1^2$.

It is important to note that ${\delta }_1$ and ${\delta }_2$ play pivotal roles as tuning parameters that govern the similarity of the distributions among the three generated samples. Specifically, when both ${\delta }_1$ and ${\delta }_2$ are set to zero (${\delta }_1={\delta }_2=0$), the three samples generated from the three models exhibit identical distributions. If at least one of ${\delta }_i$, where i can be 1 or 2, is non-zero, the samples still share the same mean vector but differ in their covariance matrices. Furthermore, it is worth mentioning that the test’s power increases as both ${\delta }_1$ and ${\delta }_2$ increase.

We set $\mathit{\mu }$ as $h{(1,\dots ,p)}^{\top }/\sqrt{{\sum }_{i=1}^p{i}^2}$ and $\Gamma$ as $a\left[(1-\rho ){\mathit{I}}_p+\rho {\mathit{J}}_p\right]$, where ${\mathit{J}}_p$ denotes a matrix of ones with dimensions $p\times p$. To assess the performance of the considered tests across a range of dimensionality settings, we examine three cases: $p=10$, $p=100$, and $p=500$. For each of these cases, we consider three sets of sample sizes $(n_1,n_2,n_3)$: $(20,30,40)$, $(80,120,160)$, and $(160,240,320)$. Additionally, we investigate three levels of correlation: $\rho =0.1$, $\rho =0.5$, and $\rho =0.9$. These three values of $\rho$ correspond to samples with varying degrees of correlation, ranging from nearly uncorrelated to moderately correlated and highly correlated. Notably, correlation increases as $\rho$ grows. For simplicity, we set $h=2$ and $a=1.5$ across all three models.

In the case of the parametric bootstrap and random permutation methods, as well as the energy test, we use a total of $N=1000$ replicates for computing the p-values at each simulation run, as described in Section 3.1 and Section 3.2. It is worth noting that the W–S ${\chi }^2$-approximation method, which does not require generating replicates, is the least time-consuming among the methods considered. The empirical sizes and powers are computed based on 1000 simulation runs.

Table 1. Empirical sizes (in %) of Simulation 1.

			$\mathit{\rho }=0.1$				$\mathit{\rho }=0.5$				$\mathit{\rho }=0.9$
Model	$\mathit{p}$	(n1, n2, n3)	$T_{\mathit{PB}}$	$T_{\mathit{RP}}$	$T_{\mathit{WS}}$	$T_{\mathit{SR}}$	$T_{\mathit{PB}}$	$T_{\mathit{RP}}$	$T_{\mathit{WS}}$	$T_{\mathit{SR}}$	$T_{\mathit{PB}}$	$T_{\mathit{RP}}$	$T_{\mathit{WS}}$	$T_{\mathit{SR}}$
		$(20,30,40)$	4.0	4.5	5.2	4.9	4.4	4.0	4.7	4.1	5.7	5.5	6.1	5.9
	10	$(80,120,160)$	4.9	5.0	5.6	5.3	5.8	5.7	6.0	5.2	4.2	4.1	4.3	4.8
		$(160,240,320)$	5.1	5.1	5.2	4.9	4.0	4.2	4.4	3.9	4.6	5.4	4.9	4.3
		$(20,30,40)$	4.8	5.2	6.5	4.8	4.2	4.6	4.3	4.7	5.9	5.6	6.1	5.9
1	100	$(80,120,160)$	4.6	4.4	4.9	3.9	4.8	5.1	5.3	5.3	5.1	5.3	5.3	5.5
		$(160,240,320)$	5.8	5.6	6.9	5.5	4.7	4.5	4.5	4.2	4.6	4.3	5.0	4.5
		$(20,30,40)$	5.0	5.4	5.8	4.9	4.3	4.4	4.7	4.8	3.8	4.4	4.5	4.5
	500	$(80,120,160)$	4.9	5.3	5.8	4.3	5.0	5.2	5.3	5.4	5.0	5.3	5.5	5.3
		$(160,240,320)$	6.2	6.2	7.0	5.7	6.4	6.2	6.4	6.4	4.8	4.4	5.0	5.2
		$(20,30,40)$	4.9	5.5	6.6	5.7	4.4	4.5	5.0	4.9	5.6	5.0	5.7	5.6
	10	$(80,120,160)$	5.6	5.6	6.4	5.8	3.6	3.5	4.0	4.5	5.4	5.3	5.3	4.9
		$(160,240,320)$	5.7	5.8	6.8	6.3	4.9	5.5	5.8	5.2	5.3	5.5	5.7	5.1
		$(20,30,40)$	4.9	5.1	5.9	5.0	5.4	5.5	5.5	5.3	3.9	4.5	4.7	4.3
2	100	$(80,120,160)$	5.5	5.6	6.6	5.5	4.5	4.3	4.9	3.8	5.8	5.6	5.7	5.6
		$(160,240,320)$	5.1	5.3	5.6	4.9	6.4	6.5	6.8	7.0	5.0	4.8	4.8	4.1
		$(20,30,40)$	5.0	5.1	5.4	5.5	5.9	5.8	6.0	6.6	5.0	5.2	5.3	5.5
	500	$(80,120,160)$	4.8	4.3	5.4	5.2	5.7	5.6	5.7	5.1	4.1	3.9	4.0	3.9
		$(160,240,320)$	4.2	4.3	5.0	4.9	4.6	4.4	4.8	4.4	5.4	5.7	5.4	5.9
		$(20,30,40)$	4.2	5.0	5.6	4.8	5.4	6.0	6.5	5.7	4.9	4.9	5.5	4.8
	10	$(80,120,160)$	4.9	5.0	5.5	4.9	4.9	4.4	5.0	5.0	5.2	5.1	5.4	5.4
		$(160,240,320)$	4.5	4.6	5.5	4.9	5.0	5.0	5.4	5.0	6.0	5.9	5.8	5.8
		$(20,30,40)$	4.1	4.4	5.1	4.7	5.0	5.5	5.8	4.5	5.0	5.0	5.4	4.4
3	100	$(80,120,160)$	3.1	3.5	4.0	4.5	4.8	5.2	5.0	4.4	5.0	5.0	5.0	5.0
		$(160,240,320)$	4.7	4.5	5.4	5.0	4.8	4.8	5.1	4.5	4.6	4.6	4.8	3.9
		$(20,30,40)$	4.9	4.7	5.4	6.0	5.0	4.6	5.4	4.5	4.1	4.1	4.5	3.9
	500	$(80,120,160)$	5.6	5.8	6.1	5.5	4.3	4.2	4.5	4.7	4.8	5.0	5.1	5.2
		$(160,240,320)$	5.0	4.8	5.5	6.0	4.4	4.5	4.6	4.5	5.5	5.4	5.6	6.0
ARE			9.04	9.33	16.22	8.67	10.67	12.52	11.56	11.70	9.26	8.74	9.19	11.56

provides an overview of the empirical sizes of $T_{PB}$, $T_{RP}$, $T_{WS}$, and $T_{SR}$, with the last row displaying the associated ARE values for the three different $\rho$ values. Several observations can be made based on this table: Firstly, for nearly uncorrelated samples ($\rho =0.1$), $T_{WS}$ exhibits a slight tendency to be liberal, with an ARE value of $16.22$ that is marginally higher than the ARE values of the other three tests. Secondly, when the generated samples are moderately correlated ($\rho =0.5$) or highly correlated ($\rho =0.9$), all four tests demonstrate fairly similar empirical sizes and ARE values, making them comparable in terms of size control. Finally, it is seen that the influence of sample sizes on size control is relatively minor, even though, in theory, a larger total sample size should result in better size control.

Figure 1 displays the empirical powers of all four tests in scenarios where all three generated samples have the same mean vectors, but they differ from each other in covariance matrices. Several conclusions can be drawn regarding these power values: Firstly, for $\rho =0.1,0.5$, and $0.9$, $T_{PB}$, $T_{RP}$, and $T_{WS}$ exhibit similar empirical powers. This suggests that these three tests perform comparably, regardless of whether the data are nearly uncorrelated, moderately correlated, or highly correlated. Secondly, it is seen that under similar settings, as expected, the empirical powers of the tests generally increase with larger sample sizes. Finally, the empirical powers of $T_{SR}$ consistently rank the lowest among all four tests. This indicates that $T_{SR}$ is less powerful compared to the other three tests in these scenarios.

4.2. Simulation 2

Certainly, it should be acknowledged that the MMD-based tests $T_{PB}$, $T_{RP}$, and $T_{WS}$ may not consistently demonstrate superior performance when compared to the energy-distance-based test $T_{SR}$ as in Simulation 1. To illustrate this point, in the context of this simulation study, we keep the same experimental framework as described in Simulation 1. However, we now introduce a new collection of three models, which are defined as follows:

Model 4.

$v_{\alpha ,ir},i=1,\dots ,n_{\alpha };\alpha =1,2;r=1,\dots ,p\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }\mathcal{N}(0.5,1)$.

Model 5.

$v_{\alpha ,ir}=u_{\alpha ,ir}/\sqrt{2}+0.5,i=1,\dots ,n_{\alpha };\alpha =1,2;r=1,\dots ,p$ with $u_{\alpha ,ir},i=1,\dots ,n_{\alpha };$$\alpha =1,2;r=1,\dots ,p\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }{t}_4$.

Model 6.

$v_{\alpha ,ir}=(u_{\alpha ,ir}-1)/\sqrt{2}+0.5,i=1,\dots ,n_{\alpha };\alpha =1,2;r=1,\dots ,p$ with $u_{\alpha ,ir}$, $i=1,\dots ,n_{\alpha };\alpha =1,2;r=1,\dots ,p\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }{\chi }_1^2$.

Please be aware that $v_{\alpha ,ir}$, where i ranges from 1 to $n_{\alpha }$, $\alpha$ takes values 1 and 2, and r ranges from 1 to p, is adjusted to ensure that $E\left(v_{\alpha ,ir}\right)=0.5$ and $Var\left(v_{\alpha ,ir}\right)=1$ across all three models. When both ${\delta }_1$ and ${\delta }_2$ are set to 0, the three generated samples follow the same distributions as in Simulation 1. Consequently, this implies that the empirical sizes of all four tests in Simulation 2 will be similar to those observed in Simulation 1. However, when at least one of ${\delta }_i$ (where i takes values 1 and 2) is non-zero, the three samples exhibit distinct mean vectors and covariance matrices, differing from those observed in Simulation 1. As a result, our focus should be on calculating the empirical powers of all four tests based on Models 4–6.

Figure 2 presents the empirical powers of all four tests based on Models 4–6, offering several noteworthy insights. First, it is evident that $T_{PB}$, $T_{RP}$, and $T_{WS}$ demonstrate similar empirical powers. This implies that these three tests exhibit comparable performance regardless of whether the generated data follow a normal or non-normal distribution. Second, it is also observed that under similar settings, the empirical powers of the tests generally increase with larger sample sizes. Lastly, when $\rho$ takes on values of 0.5 and 0.9, the empirical powers of $T_{SR}$ surpass those of the other three tests. However, when $\rho =0.1$, the empirical powers of all four tests are generally comparable. This indicates that, in the scenarios under consideration, $T_{SR}$ demonstrates greater effectiveness compared to the other three tests when the correlation coefficient $\rho$ is relatively high.

From these two simulation studies, it can be inferred that the proposed MMD-based tests $T_{PB}$, $T_{RP}$, and $T_{WS}$ may outperform the energy-distance-based test $T_{SR}$ when the differences in distributions are primarily in covariance matrices, while the reverse could be true when the differences in distribution involve both mean vectors and covariance matrices. Notably, the MMD-based test $T_{WS}$ generally requires less computational effort compared to the bootstrap or permutation-based tests $T_{PB}$, $T_{RP}$, and $T_{SR}$.

5. Application to the Corneal Surface Data

The corneal surface data are briefly mentioned in Section 1. They were acquired during a keratoconus study, a collaborative project involving Ms. Nancy Tripoli and Dr. Kenneth L. Cohen from the Department of Ophthalmology at the University of North Carolina, Chapel Hill. This dataset comprises 150 observations with each corneal surface having more than 6000 height measurements. It can be categorized into four distinct groups: a group of 43 healthy corneas (referred to as the normal cornea group), a group of 14 corneas with unilateral suspect characteristics, a group of 21 corneas with suspect map features, and a group of 72 corneas clinically diagnosed with keratoconus. It is important to note that the corneal surfaces within the normal, unilateral suspect, and suspect map groups exhibit similar shapes, but they significantly differ from the corneal surfaces observed in the clinical keratoconus group (refer to Figure 1 in [24] for visualization). In the process of reconstructing a corneal surface, ref. [24] utilized the Zernike regression model to fit the height measurements associated with the corneal surface. The height of the corneal surface at a specific radius r and angle $\theta$ is denoted as $h(r,\theta )$, while $\widehat{h}(r,\theta )$ represents the height estimated through the fitted model within the predefined region of interest. This region of interest spans from $r=0$ to $r=r_0$ and from $\theta =0$ to $\theta =2\pi$, with $r_0$ being a predetermined positive constant. To naturally represent each corneal surface, a feature vector is constructed, consisting of values $\widehat{h}(r_i,{\theta }_j)$, where i ranges from 1 to K and j ranges from 1 to L. These values are obtained by evaluating the fitted corneal surface $\widehat{h}(r,\theta )$ at a grid of points defined as $r_i=r_{0}i/K$ and ${\theta }_j=2\pi j/L$ for $i=1,\dots ,K$ and $j=1,\dots ,L$. For simplicity, we choose to set $K=20$ and $L=100$, resulting in a feature vector with dimensions of 2000 for each corneal surface.

For simplicity, we put the fitted feature vectors for the complete corneal surface dataset collectively into a feature matrix with dimensions of $150\times 2000$. In this matrix, each row corresponds to a feature vector representing a corneal surface. Specifically, the initial 43 rows of the feature matrix correspond to observations from the normal group, sequentially followed by 14 rows from the unilateral suspect group, 21 rows from the suspect map group, and lastly, 72 rows from the clinical keratoconus group.

Our objective is to examine whether there are significant differences in the distributions among various corneal surface groups, referred to as multi-sample problems for the equality of distributions for high-dimensional data, given that the high-dimensional feature vectors represent the observations of the corneal surface data. In this application, we employ $T_{PB}$, $T_{RP}$, $T_{WS}$, and $T_{SR}$ to address these problems. For both the parametric bootstrap and random permutation methods, as well as the energy test, we perform a total of N = 10,000 replicates to compute the associated p-values. For simplicity, we denote the normal, unilateral suspect, suspect map, and clinical keratoconus groups as NOR, UNI, SUS, and CLI, respectively.

Table 2. p-values (in %) for testing the distribution equality of corneal surface groups.

	Method
Case	$T_{\mathit{PB}}$	$T_{\mathit{RP}}$	$T_{\mathit{WS}}$	$T_{\mathit{SR}}$
NOR vs. UNI vs. SUS vs. CLI	0	0	0.00004	0.09
NOR vs. UNI vs. SUS	31.2	32.3	33.3	38.9
NOR vs. UNI vs. CLI	0.02	0	0.0012	0.25
NOR vs. SUS vs. CLI	0	0	0.00001	0.01
UNI vs. SUS vs. CLI	0.01	0	0.0012	0.24

displays the results obtained from the application of four statistical tests, namely $T_{PB}$, $T_{RP}$, $T_{WS}$, and $T_{SR}$, to assess the equality of distributions among different corneal surface groups. A careful examination of these results yields several noteworthy insights. To begin with, when considering the comparison among the corneal surface groups labeled “NOR vs. UNI vs. SUS vs. CLI”, it is evident that all four tests reject the null hypothesis. This rejection signifies that there exists at least one significant difference among the distributions of these four corneal surface groups. Consequently, further investigation is warranted to identify which specific group or groups differ from the others. Secondly, it is important to highlight the outcome for the comparison involving “NOR vs. UNI vs. SUS”. In this case, none of the four tests reject the null hypothesis at the $5\%$ significance level. This outcome indicates that the normal, unilateral suspect, and suspect map groups share a similar distribution pattern. Thirdly, across the remaining three comparisons, all four tests consistently reject the null hypothesis. These results align with the observations depicted in Figure 1 of [24], which illustrates the distinctiveness of corneal surfaces within the clinical keratoconus group when compared to the other three groups. Fourthly, when focusing on the comparison “NOR vs. UNI vs. SUS”, it is worth noting that the p-values obtained from $T_{PB}$, $T_{RP}$, and $T_{WS}$ are quite similar, suggesting their comparable performance. This consistency in p-values is also reflected in the empirical sizes presented in Table 1. Lastly, when analyzing the cases involving CLI, it becomes evident that the p-values generated by $T_{SR}$ consistently exhibit larger values compared to those produced by the other three tests. This discrepancy implies that $T_{SR}$ may have a lower sensitivity in detecting distribution differences when compared to the other tests, indicating the potentially reduced statistical power in this real data example.

Notice that the test $T_{PB}$ is bootstrap-based and the tests $T_{RP}$ and $T_{SR}$ are permutation-based. This means that their p-values are obtained via bootstrapping or permutating numerous random samples to compute the associated p-values as described in Section 3.1 and Section 3.2. Thus, the p-values of these tests are random, i.e., they are different at different instances. However, the p-value of $T_{WS}$ remains fixed. In order to investigate this clearly, we performed 500 iterations of $T_{PB},T_{RP},T_{WS}$, and $T_{SR}$ on the case “NOR vs. UNI vs. SUS”. The boxplots of the corresponding p-values of the four tests are presented in Figure 3. It is evident that the p-values obtained from $T_{WS}$ remain fixed, whereas those derived from $T_{PB},T_{RP}$, and $T_{SR}$ exhibit variability. This contrast underscores the fact that p-values resulting from bootstrap-based or permutation-based tests indeed differ across various instances.

6. Concluding Remarks

Testing whether multiple high-dimensional samples adhere to the same distribution is a common area of research interest. This paper introduces and investigates a novel MMD-based test to address this question. The null distribution of this test is approximated using three methods: parametric bootstrap, random permutation, and the W–S ${\chi }^2$-approximation approach. Results from two simulation studies and a real data application demonstrate that the proposed test exhibits effective size control and superior statistical power compared to the energy test introduced by [12] when the differences among sample distributions are primarily related to covariance matrices rather than mean vectors. Thus, the proposed test is generally well suited for conducting multi-sample equal distribution testing on high-dimensional data. We particularly recommend its use in scenarios where distribution differences are associated with covariance matrices. Conversely, when distribution differences predominantly pertain to means, the energy test is a more powerful choice. However, in practice, determining whether distribution differences are related to means or covariance matrices can be challenging. Therefore, we suggest considering both the new test and the energy test as viable options. Nevertheless, it is important to note that implementing the proposed test comes with certain challenges. Both the parametric bootstrap and random permutation methods can be computationally intensive, leading to variable p-values across different applications. In contrast, the W–S ${\chi }^2$-approximation method offers computational efficiency and produces fixed p-values. However, its accuracy is limited as it solely relies on matching two cumulants of the test statistic under the null hypothesis.

An intriguing question arises naturally: Can we enhance the accuracy of the proposed test by matching three cumulants of the test statistic? Recent work by [18] suggests that this is indeed possible. However, deriving the third cumulant of the test statistic presents a current challenge and requires further investigation. Another aspect to consider is the choice of kernel width. While the paper opts for simplicity by utilizing the median distance between observed vectors in the pooled sample, it is worth exploring the kernel width choice recommended by [18] to potentially enhance the test’s statistical power. These avenues for future research promise exciting developments and warrant further exploration.