Let \(\varvec{X}\in \mathbb R^p, \ p\ge 1\) be a random vector with an absolutely continuous cumulative distribution function (DF) \(F_{\varvec{X}}\), given as \(F_{\varvec{X}}(\varvec{x})=\mathbb P(\varvec{X}\le \varvec{x})\), \(\varvec{x}\in \mathbb {R}^p\). We shall begin our exposition with the simple hypothesis, whereby on the basis of independent copies \({\varvec{\mathcal {X}}}_{n}:=\big \{\varvec{X}_1,...,\varvec{X}_n\big \}\) on \(\varvec{X}\), we consider the problem of goodness-of-fit (GoF) testing with the null hypothesis

against general alternatives, where \(F_0\) is the DF of a completely specified distribution. The discussion is extended to the problem of testing a composite null hypothesis later in Section 4.

GoF tests for arbitrary distributions in the multivariate setting have become particularly important in recent times with the wide availability of data in high dimension. The case of normality is, not surprisingly, special, as there exists a plethora of GoF tests for the multivariate normal distribution; the reader is referred to the surveys of Thode (2002), Henze (2002), Ebner and Henze (2020), Arnastauskaitė et al. (2021), Chen and Genton (2023). On the other hand, outside the Gaussian context, there exist a few tests which however are only tailored for specific multivariate distributions (Meintanis and Hlávka 2010; Meintanis et al. 2015, 2024), while GoF methods for arbitrary multivariate laws are relatively scarce. A few exceptions are the general tests of Jiménez-Gamero et al. (2009), Meintanis et al. (2014), Khmaladze (2016), Ebner et al. (2018), Hallin et al. (2021b). Most of the aforementioned general methods, however, are computationally intensive, either in computing the test statistic itself and/or in getting critical values because of the need for bootstrap approximations. In other words, an exactly, i.e., for each sample size, distribution-free GoF test rendering the bootstrap unnecessary, while at the same time involving a reasonable amount of computational complexity, is still missing in the literature. The subject matter of this paper is to propose a computationally mild general approach for conducting GoF tests for multivariate laws that enjoys the property of exact distributional freeness.

To this end, write \(\varphi _{\varvec{X}}(\varvec{t}) =\mathsf E [\exp ({\texttt {i}}{\varvec{t^\top X}})]\) and \(\varphi _0(\varvec{t})\), \({\varvec{t}} \in \mathbb {R}^p\), for the characteristic function (CF) corresponding to \(F_{\varvec{X}}\) and \(F_0\), respectively. In view of the uniqueness of CFs, a reasonable CF-based criterion for the null hypothesis \({{\mathcal {H}}}_0\) figuring in (1) is the statistic

where

is the empirical CF corresponding to \({\varvec{\mathcal {X}}}_{n}\) and \(w(\cdot )\) is a nonnegative weight function, satisfying

Large values of \(D_{n}\) indicate the violation of the null hypothesis \({{\mathcal {H}}}_0\), so the corresponding test rejects \({{\mathcal {H}}}_0\) if \(D_n\) exceeds a certain threshold computed from the (asymptotic) distribution under the null.

Goodness-of-fit criteria for multivariate observations incorporating test statistics as the one in (2) may be computationally intensive, depending on the null CF \(\varphi _0(\cdot )\), and have a typical but non-standard asymptotic behavior under the null hypothesis \({{\mathcal {H}}}_0\) that involves the DF of \(\varvec{X}_1\). Our proposal herein also involves a CF-based test, but in the interest of computational simplicity, rather than using the formulation figuring in (2) we employ the Monte Carlo approach of casting any GoF test as a two-sample test, suggested by Chen et al. (2022). Moreover, distribution-freeness is achieved by employing, instead of the original observations, the corresponding multivariate ranks based on optimal transport theory recently developed by Hallin et al. (2021a) and Deb and Sen (2023).

The rest of this work is outlined as follows: Section 2 briefly introduces the concept of multivariate ranks. In Sect. 3, we define and compute the test statistic for a simple null hypotheses, while the case of a composite null hypothesis involving unspecified distributional parameters is considered in Sect. 4. The large-sample results under the null hypothesis as well as under alternatives are provided in Sect. 3. The finite-sample behavior of the methods is investigated by a Monte Carlo study in Sect. 6, and a real data application is provided in Sect. 7. The paper concludes with the discussion of results in Sect. 8.

All proofs are deferred to the Appendix, which also includes material related to the simulation study and the real data analysis.

Univariate rank-based methods provide a robust, flexible, and distribution-free approach to statistical analysis. However, extending rank-based inference to a multivariate setting is not a straightforward task due to the absence of canonical ordering in \(\mathbb {R}^p\) for \(p>1\). Various concepts of multivariate ranks have been considered in the literature, for instance, componentwise ranks, spatial ranks and signs, depth-based ranks, or Mahalanobis ranks and signs, see Hallin et al (2021a, Section 1.2) for more details and further references. Recently, ranks and signs based on the optimal measure transport have become popular and proved to be useful in various multivariate statistical problems, see, e.g., Shi et al. (2022); Deb and Sen (2023); Hallin et al. (2023); Huang and Sen (2026); Hlubinka and Hudecová (2024).

Let \((\mu ,\nu )\) be a pair of probability measures on \(\mathbb {R}^p\). The optimal measure transport (OMT) problem, first formulated by Monge (1781), seeks to find a mapping \(\varvec{G}:\mathbb {R}^p\rightarrow \mathbb {R}^p\) that minimizes

for a given cost function \(\mathcal {C}:\mathbb {R}^p\times \mathbb {R}^p \rightarrow [0,\infty )\). We write \(\varvec{G}\#\mu =\nu \) when \(\varvec{G}(\varvec{Z})\sim \nu \) for \(\varvec{Z}\sim \mu \), and say that \(\varvec{G}\) pushes \(\mu \) to \(\nu \). The most common choice is \(\mathcal {C}(\varvec{x},\varvec{y})=\Vert \varvec{x}-\varvec{y}\Vert ^2\); therefore, we consider only this cost function in everything what follows. The concept of OMT has been recently used to define multivariate ranks and quantiles (Hallin et al. 2021a; Deb and Sen 2023) and data depth (Chernozhukov et al. 2017).

Let \(\varvec{Z}\) be a random vector with an absolutely continuous distribution \(\mathcal {P}_{\varvec{Z}}\) on \(\mathbb {R}^p\), and let \(\nu \) be some specified reference measure on \(\mathbb {R}^p\). If \(\varvec{Z}\) has finite second moments, then there exists a unique solution \(\varvec{G}^*\) to (5) with \(\mathcal {C}(\varvec{x},\varvec{y})=\Vert \varvec{x}-\varvec{y}\Vert ^2\). More generally (without the assumption \(\textsf{E}\Vert \varvec{Z}\Vert ^2<\infty \)), the mapping \(\varvec{G}^*\) can be defined as described in Remark 1. In any case, \(\varvec{G}^*\) is in some sense unique mapping that pushes \(\mathcal {P}_{\varvec{Z}}\) to \(\nu \), see Remark 1 for details. Deb and Sen (2023) take \(\nu \) as the uniform measure on \([0,1]^p\) and call \(\varvec{G}^*\) the population rank function. Hallin et al. (2021a) specify \(\nu \) as the spherical distribution of a random vector \(U \varvec{S}\), where U is uniformly distributed on [0, 1] and independent of \(\varvec{S}\) with a uniform distribution on the unit sphere \( \mathcal {S}_{p}=\{\varvec{x} \in \mathbb {R}^p; \Vert \varvec{x}\Vert = 1\}\). The optimal mapping \(\varvec{G}^*\) from (5) is then called the multivariate center-outward distribution function.

Let \(\varvec{\mathcal {Z}}_N: = \big \{\varvec{Z}_1,\dots ,\varvec{Z}_N\big \}\) be a random sample from \(\mathcal {P}_{\varvec{Z}}\). Moreover, let \(\mathcal {G}_N =\{\varvec{g}_1,\dots ,\varvec{g}_N\}\) be a given grid of points from the support of \(\nu \). The empirical version \(\widehat{\varvec{G}}_N\) of \(\varvec{G}^*\) is the solution to (5) with \(\mu \) and \(\nu \) being the empirical distributions on \(\varvec{\mathcal {Z}}_N\) and \(\mathcal {G}_N\), respectively. In other words,

where the \(\arg \min \) is computed among all bijections \(\varvec{G}:\varvec{\mathcal {Z}}_N \rightarrow \mathcal {G}_N\). Consequently, \(\widehat{\varvec{G}}_N\) is used to define the multivariate ranks of \(\varvec{Z}_1,\dots ,\varvec{Z}_N\), but the terminology here slightly varies: Deb and Sen (2023) call multivariate rank of \(\varvec{Z}_i\) directly \(\widehat{\varvec{G}}_N(\varvec{Z}_i)\), while Hallin et al. (2021a) define univariate ranks and multivariate signs as simple functions of \(\widehat{\varvec{G}}_N(\varvec{Z}_i)\). For simplicity, we call \(\varvec{R}_i= \widehat{\varvec{G}}_N(\varvec{Z}_i)\) the multivariate rank of \(\varvec{Z}_i\) in what follows. Remark that computation of the discrete OMT \(\widehat{\varvec{G}}_N\) is a standard optimization task that can be formulated as a linear program for which efficient algorithms are available, see Peyré et al (2019, Chapter 3).

Notice that the multivariate ranks depend on the choice of the grid \(\mathcal {G}_N\). Selection of the grid points should be driven by the choice of the theoretical reference measure \(\nu \) in a sense that the uniform measure on \(\mathcal {G}_N\) converges weakly to \(\nu \) as \(N\rightarrow \infty \). In that case, it is possible to show that

see Deb and Sen (2023, Theorem 2.1). Under slightly stronger conditions, see Hallin et al (2021a, Proposition 2.4), the convergence is uniform

Due to their choice of \(\nu \), Deb and Sen (2023) take the grid points \(\mathcal {G}_N\) as a subset of \([0,1]^p\), while Hallin et al. (2021a) consider grid points in the unit ball in \(\mathbb {R}^p\). In general, if the support of \(\nu \) is \(S\subset \mathbb {R}^p\), then the grid points should be distributed in S so that the uniform measure on \(\mathcal {G}_N\) well approximates \(\nu \). There is not a unique approach to achieve this, and a construction of \(\mathcal {G}_N\) often makes use of low-discrepancy sequences, e.g., Halton sequences (Halton 1960) or GLP sets (Fang and Wang 1994, Section 1.3).

If the reference measure \(\nu \) is uniform on \([0,1]^p\) as in Deb and Sen (2023), then the grid set can be taken simply as a Halton sequence \(\{\varvec{x}_i\}_{i=1}^N\) in \([0,1]^p\) (Halton 1960; Dutang and Savicky 2024). Such grid will be referred to as a rectangular grid, abbreviated as \(\mathcal {G}_N^R\).

If \(\nu \) is the distribution of the spherically uniform distribution of the random vector \(U\varvec{S}\), then Hallin et al. (2021a) recommend to take \(\mathcal {G}_N\) as a collection of \(n_0\) replications of \(\varvec{0}\) and a set of points \(\{\varvec{g}_{ij}\}_{i=1,j=1}^{n_R,n_S}\), where \(\varvec{g}_{ij} =\frac{i}{n_R+1} \varvec{s}_j\), where \(n_0\), \(n_R\), \(n_S\) are such that \(n_R\cdot n_S + n_0 = N\) and \(\varvec{s}_1,\dots ,\varvec{s}_{n_S}\) are directional vectors chosen uniformly as possible from the unit sphere \(\mathcal {S}_p\). However, this approach requires the choice of \(n_R\) and \(n_S\). Hlávka et al. (2025) suggest to take \(\{\varvec{g}_i\}_{i=1}^N\) such that

where \(\{\varvec{x}_i\}_{i=1}^N\) is a Halton sequence in \([0,1]^p\), \(\varvec{x}_i=(x_{i,1},\dots ,x_{i,p})^\top \), so \(x_{i,1}\) is the first element of \(\varvec{x}_i\), and \(\varvec{s}_i\) is a directional vector from the unit sphere in \(\mathbb {R}^p\) computed as \(\varvec{s}_i = \tau (x_{i,2},\dots ,x_{i,p})\), where \(\tau \) is a mapping from \([0,1]^{p-1}\) to the unit sphere \(\mathcal {S}_p\) such that \(\tau (\varvec{U})\) has a uniform distribution on \(\mathcal {S}_p\) whenever \(\varvec{U}\) has a uniform distribution on \([0,1]^{p-1}\), see Fang and Wang (1994, Section 1.5.3). Such grid will be referred to as a spherical grid, abbreviated as \(\mathcal {G}_N^S\).

A graphical illustration of \(\mathcal {G}_N^R\) and \(\mathcal {G}_N^S\) is provided in Fig. 2 in Appendix I.

Hallin et al. (2021a) define so-called center-outward distribution function \(\varvec{F}_{\pm }\) more generally, without the assumption of finite second-order moments of \(\varvec{Z}\), as the unique gradient of a convex mapping that pushes \(\mathcal {P}_{\varvec{Z}}\) to \(\nu \), where \(\nu \) is the distribution of the random vector \(U \varvec{S}\) specified above. The uniform convergence in (7) then holds for \(\varvec{G}^*\) replaced by \(\varvec{F}_{\pm }\). If \(\textsf{E}\Vert \varvec{Z}\Vert ^2 < \infty \), then it follows from Brenier’s theorem (Villani 2003, Theorem 2.32) that \(\varvec{F}_{\pm }=\varvec{G}^*\), where \(\varvec{G}^*\) is the solution to the Monge’s problem in (5), see (Villani 2003, Chapter 3).

Hence, if one is not willing to make an assumption about finiteness of the second-order moments of \(\varvec{Z}\), then it is possible to define \(\varvec{G}^*\) as \(\varvec{F}_{\pm }\). Such \(\varvec{G}^*\) still pushes \(\mathcal {P}_{\varvec{Z}}\) to \(\nu \) and enjoys a uniqueness property, since it is the unique gradient of a convex mapping with this property. Moreover, it is also the optimal mapping with this property if \(\textsf{E}\Vert \varvec{Z}\Vert ^2 < \infty \).

Consider a random sample \({\varvec{\mathcal {X}}}_{n}\) from the absolutely continuous distribution \(\mathcal {P}_{\varvec{X}}\) with DF \(F_{\varvec{X}}\) and recall that the aim is to test \(\mathcal {H}_0\) in (1) for a given absolutely continuous DF \(F_0\). Let m be some integer, the choice of which will be discussed later, and let \({\varvec{\mathcal {X}}}^{(0)}_{m}:=\big \{\varvec{X}^{(0)}_{1},...,\varvec{X}^{(0)}_{m}\big \}\) be a random sample drawn from the distribution with DF \(F_0\), independent of \({\varvec{\mathcal {X}}}_{n}\). Furthermore, let \(\mathcal {G}_N\) be a specified grid of \(N=n+m\) points.

Denote as \(\varvec{\mathcal {Z}}_N = {\varvec{\mathcal {X}}}_{n} \cup {\varvec{\mathcal {X}}}^{(0)}_{m}\) the union of the two samples. For \(\varvec{\mathcal {Z}}_N\) and \(\mathcal {G}_N\), we can compute the OMT \(\widehat{\varvec{G}}_N\) from (6). Let \({\varvec{\mathcal {R}}}_{n}:=\big \{\varvec{R}_{i}\big \}_{i=1}^{n}\), \(\varvec{R}_i= \widehat{\varvec{G}}_N(\varvec{X}_i)\), be the collection of rank vectors associated with \({\varvec{\mathcal {X}}}_{n}\), and similarly \({\varvec{\mathcal {R}}}^{(0)}_{m}:=\big \{\varvec{R}^{(0)}_{j}\big \}_{j=1}^{m}\), \(\varvec{R}_j^{(0)} = \widehat{\varvec{G}}_N(\varvec{X}_j^{(0)})\) be the collection of rank vectors associated with \({\varvec{\mathcal {X}}}^{(0)}_{m}\). Then, a rank-based analogue of the test statistic \(D_{n}({\varvec{\mathcal {X}}}_{n})\) from (2) is given by

where \(\widehat{\phi }_n\) and \(\widehat{\phi }^{(0)}_{m}\) are computed as in (3) for \({\varvec{\mathcal {R}}}_{n}\) and \({\varvec{\mathcal {R}}}^{(0)}_{m}\), respectively, that is

The test statistic \(D_{n,m}\) in (8) is a two-sample test statistic that measures the distance between the empirical CF of the ranks \({\varvec{\mathcal {R}}}_{n}\) of the observed data \({\varvec{\mathcal {X}}}_{n}\) and the corresponding empirical CF of the ranks \({\varvec{\mathcal {R}}}^{(0)}_{m}\) obtained from an artificial random sample \({\varvec{\mathcal {X}}}^{(0)}_{m}\) drawn from the reference distribution with DF \(F_0\). The idea of casting any GoF test as a two-sample test goes back to Friedman (2003), and has been used very effectively in Chen et al. (2022) and Karling et al. (2023) for CF-based GoF testing, with promising results. However, while these works employ the original data here we propose to use the corresponding ranks. As it will be discussed below, the use of appropriate multivariate ranks leads to an exactly distribution-free test while, additionally, basing the test on empirical CFs retains the computational simplicity and, of course, the consistency of earlier CF-based tests.

Under the null hypothesis \(\mathcal {H}_0\), both samples \(\varvec{\mathcal {X}}_n\) and \(\varvec{\mathcal {X}}^{(0)}_m\) come from the same distribution with DF \(F_0\) and with a population OMT function \(\varvec{G}^*\) such that both \(\varvec{G}^*(\varvec{X}_i)\) and \(\varvec{G}^*(\varvec{X}_j^{(0)})\), \(i=1,\dots ,n\), \(j=1,\dots ,m\), have the reference distribution \(\nu \). It follows from (7) that for m, n large, both \(\widehat{\phi }_n\) and \(\widehat{\phi }_m^{(0)}\) will be close to the CF of \(\nu \) and the test statistic \(D_{n,m}\) is expected to be small. On the other hand, large values of \(D_{n,m}\) indicate that the distributions of \({\varvec{\mathcal {R}}}_{n}\) and \({\varvec{\mathcal {R}}}^{(0)}_{m}\) differ, indicating that the two samples \({\varvec{\mathcal {X}}}_{n}\) and \({\varvec{\mathcal {X}}}^{(0)}_{m}\) do not come from the same distribution. Therefore, the null hypothesis \(\mathcal {H}_0\) is rejected if

where \(c_{n,m,\alpha }\) is a critical value such that the test keeps the prescribed level \(\alpha \). Remark that \(c_{n,m,\alpha }\) depends not only on the sample sizes n, m and level \(\alpha \), but also on the choice of the grid points \(\mathcal {G}_N\), the weight function w, and the dimension p. Note that the choice of the set \(\mathcal {G}_N\) implicitly involves the choice of the reference measure \(\nu \), i.e., the type of the multivariate ranks. On the other hand, the critical value \(c_{n,m,\alpha }\) does not depend on the distribution \(\mathcal {P}_{\varvec{X}}\), as justified in the following lemma.

Under the null hypothesis \(\mathcal {H}_0\), the distribution of \(D_{n,m}\) is free of \(\mathcal {P}_{\varvec{X}}\).

The significance of the CF-based test statistic (2), computed directly from the original data, is often evaluated via resampling or permutation techniques because its finite as well as asymptotic distribution typically depends on the unknown distribution \(\mathcal {P}_{\varvec{X}}\) in a non-trivial way. In contrast to this, the finite-sample distribution-freeness of our test statistic \(D_{n,m}\) allows to avoid these resampling techniques as the critical value \(c_{n,m,\alpha }\), for specified m, n, \(\mathcal {G}_N\) and w, needs to be computed only once. Table 1 provides an example of critical values \(c_{n,m,\alpha }\) for \(p=2\). Alternatively, the critical values for \(D_{n,m}\) can be calculated from its asymptotic distribution provided in Sect. 3.2, but it is typically easier to compute \(c_{n,m,\alpha }\) using Monte Carlo simulations.

Using the well-known routine calculations for the test statistic in (8) with a weight function satisfying (4) yields the equivalent expression

where

If w is selected as a density of a spherical distribution on \(\mathbb {R}^p\), then \(C_w({\varvec{x}})\) is its characteristic function, and the formula (9) provides a closed-form expression for \(D_{n,m}\) that avoids computation of multivariate integrals. Specifically, if w is a density of a spherical stable distribution, then

If w is a density of a generalized spherical Laplace distribution, then \(C_w(\varvec{x})=(1+\Vert \varvec{x}\Vert ^2)^{-\gamma }, \gamma >0\). Remark that in order to gain extra flexibility, one can compute \(D_{n,m}\) from (9) with \(C_w\) enhanced by an additional scale parameter \(a >0\). For instance, for the function from (11) this leads to

which is for \(\gamma =2\) the characteristic function of the normal distribution \(\textsf{N}(\varvec{0}, 2a^2 \varvec{I})\).

Remark that our test statistic is related to the celebrated energy criterion of Székely and Rizzo (2013), when applied to the ranks \(\varvec{\mathcal {R}}_n\) and \(\varvec{\mathcal {R}}_m^{(0)}\). If w is chosen such that \(C_w\) is given as in (11) and one uses the approximation \(\mathrm{{e}}^x\approx 1+x\), then \(D_{n,m}\approx E_{n,m}\), where

It should be pointed out that the energy statistic based on the original observations is applicable only for \(\gamma \in (0,2)\) and if \(\mathsf E\Vert X_t\Vert ^{\gamma }<\infty \), while the statistic in (2) or (8) works for \(\gamma \in (0,2]\) (stable density) or for all \(\gamma >0\) (Laplace density), and is free of moment assumptions. Because we use ranks though, the energy statistic is expected to work even with heavy-tailed data.

Deb and Sen (2023) proposed a rank energy statistic based on the OMT ranks (with a reference measure \(\nu \) being uniform on \([0,1]^p\)) for a two-sample problem for testing the equality of two continuous distributions. Their test statistic, for testing the equality of distributions of \(\varvec{\mathcal {X}}_n\) and \(\varvec{\mathcal {X}}^{(0)}_m\), takes the form (13) with \(\gamma =1\). Hence, our procedure can be seen also as a generalization of the test recommended by Deb and Sen (2023) for a two-sample problem.

The idea of using an artificial sample to test some multivariate GoF hypothesis has been already used in Chen et al. (2022). Their simulations indicate that efficiency can be improved by repeating the procedure M times, for some chosen \(M>0\), and using the average over these M repetitions as the final test statistic. A similar idea could potentially be applied to the test based on \(D_{n,m}\). However, this would require computing the discrete optimal transport in (6) M-times. As a result, the test would become computationally demanding for n, m, N large. Note also that the critical values would need to be recalculated accordingly.

In the following, we assume that \(\nu \) is a specified absolutely continuous reference measure on \(\mathbb {R}^p\) with a compact support \(\mathcal {S} \subset \mathbb {R}^p\), \(N=n+m\), and the weight function w satisfies (4).

Let \(\{\mathcal {G}_N\}\) be a sequence of grids in \(\mathcal {S}\) such that the uniform measure on \(\mathcal {G}_N\) converges weakly to \(\nu \) as \(n\rightarrow \infty \), and let \(n/m\rightarrow \lambda \in (0,\infty )\) as \(n\rightarrow \infty \). Then, under \(\mathcal {H}_0\),

as \(n\rightarrow \infty \), where \(\{Z(\varvec{t}), \varvec{t}\in \mathbb {R}^p\}\) is a centered Gaussian process with a covariance function

where \(\varvec{Y}\) is a random vector with distribution \(\nu \).

Theorem 1 shows that the test based on \(D_{n,m}\) is also asymptotically distribution-free, and the limiting distribution of \(D_{n,m}\) is the same as the distribution of \( D=\sum _{k=1}^{\infty } \lambda _k U_k^2\), where \(U_k\) are iid standard normal variables and \(\lambda _1\ge \lambda _2\ge \dots \) are constants that depend on the weight function w and the reference measure \(\nu \). For a specified w, the asymptotic distribution from Theorem 1 can be simulated and, subsequently, one can obtain asymptotic critical values \(c_{\alpha }\), as described in Algorithm 2 in the Appendix. However, this procedure is numerically demanding for larger dimensions, so we rather recommend computing the critical values \(c_{n,m,\alpha }\) by Monte Carlo simulations.

Table 1 compares the asymptotic critical values \(c_{\alpha }\) with the finite-sample critical values \(c_{n,m,\alpha }\) for dimension \(p=2\), significance level \(\alpha =0.05\), and various values of sample sizes n, m. The weighting corresponds to \(C_{a,\gamma }\) from (12) with \(\gamma =2\). It is visible that \(c_{n,m,\alpha }\) are close to \(c_{\alpha }\) even for rather small values of n.

The condition \(n/m \rightarrow \lambda \in (0, \infty )\) excludes non-standard cases in which the two sample sizes grow to infinity at different rates, while imposing no restriction on practical applications with a finite-sample size n. The effect of the choice of m on the power of the test is explored in the Monte Carlo simulation study in Sect. 6.

Recall that \(\mathcal {P}_{\varvec{X}}\) stands for the distribution of \(\varvec{X}_i\), and let \(\mathcal {P}_0\) be the distribution corresponding to \(F_0\). Consider now the alternative \(\mathcal {H}_1: F_{\varvec{X}} \ne F_0\). In this case, \( \mathcal {P}_{\varvec{X}} \ne \mathcal {P}_0\) and the pooled data \(\varvec{\mathcal {Z}}_N =\varvec{\mathcal {X}}_n \cup \varvec{\mathcal {X}}_m^{(0)}\) consist of two independent samples from two different distributions \(\mathcal {P}_{\varvec{X}}\) and \(\mathcal {P}_0\). The existence and uniqueness of the mapping \(\widehat{\varvec{G}}_N\) are still guaranteed (by standard arguments). Asymptotically, \(\widehat{\varvec{G}}_N\) behaves as the solution to (5) for \(\mu \) being a mixture of \(\mathcal {P}_{\varvec{X}}\) and \(\mathcal {P}_{0}\). The next theorem describes the asymptotic behavior of the test statistic \(D_{n,m}\) in this situation.

We still assume that \(\nu \) is absolutely continuous with a compact support \(\mathcal {S} \subset \mathbb {R}^p\), \(N=n+m\), and both \(\mathcal {P}_{\varvec{X}}\) and \(\mathcal {P}_0\) are absolutely continuous.

Let both \(\mathcal {P}_{\varvec{X}}\) and \(\mathcal {P}_0\) have finite second moments. Let \(\{\mathcal {G}_N\}\) be a sequence of grids on a compact set S such that the uniform measure on \(\{\mathcal {G}_N\}\) converges weakly to \(\nu \) as \(n\rightarrow \infty \) and \(n/m\rightarrow \lambda \in (0,\infty )\). Let the weight function w satisfy (4) and \(\int _{\mathbb {R}^p} \Vert \varvec{t}\Vert w(\varvec{t}) \textrm{d} \varvec{t} <\infty \). Let \(\varvec{G}\) be the solution to (5) for \(\mu = \frac{1}{1+\lambda } \mathcal {P}_{\varvec{X}} + \frac{\lambda }{1+\lambda } \mathcal {P}_0\). Then as \(n\rightarrow \infty \),

In particular, if \(\mathcal {P}_{\varvec{X}} \ne \mathcal {P}_0\) and w is positive on an open neighborhood of \(\varvec{0}\), then the right-hand side of (16) is positive, and the test based on \(D_{n,m}\) is consistent for \(n \rightarrow \infty \).

The quantity on the right-hand side of (16) is a positive multiple of the weighted \(L_2\) comparison of the characteristic functions of \(\varvec{G}(\varvec{X}_1)\) and \(\varvec{G}(\varvec{X}_1^{(0)})\). The assumptions of Theorem 2 ensure that if \(\mathcal {P}_{\varvec{X}} \ne \mathcal {P}_0\), then the same holds for \(\varvec{G}(\varvec{X}_1)\) and \(\varvec{G}(\varvec{X}_1^{(0)})\). This phenomenon is illustrated by the empirical transform \(\widehat{\varvec{G}}_N\) in Fig. 3 in Appendix I, which graphically compares \(\widehat{\varvec{G}}_N(\varvec{X}_i)\) and \(\widehat{\varvec{G}}_N(\varvec{X}_j^{(0)})\) under various alternatives \(\mathcal {P}_{\varvec{X}} \ne \mathcal {P}_0\).

Note that it may be possible to weaken the assumption of finite second-order moments similarly to Hallin et al. (2021a), but we do not explore this issue further.

This section addresses the problem of testing GoF to parametric families of distributions. Let \( \mathcal {F} = \{ F_{\varvec{\vartheta }}, \varvec{\vartheta }\in \Theta \} \) be a family of distributions indexed by a parameter \(\varvec{\vartheta }\) taking values in \(\Theta \subseteq \mathbb {R}^q, \ q\ge 1\). We wish to test the composite null hypothesis

against a general alternative \(\mathcal {H}_1^C: F_{\varvec{X}} \not \in \mathcal {F}\).

The idea of the goodness-of-fit criterion constructed as a two-sample comparison of the original data \(\varvec{\mathcal {X}}_n\) and an artificial sample can be extended to this composite null hypothesis \(\mathcal {H}_0^C\). However, as the true value of the parameter \(\varvec{\vartheta }\) is unknown, the artificial sample is simulated from a distribution with parameter that is estimated from the data \(\varvec{\mathcal {X}}_n\). Such approach is typical in various parametric bootstrap procedures.

Let \(\widehat{\varvec{\vartheta }}=\widehat{\varvec{\vartheta }}_n\) be an estimator of \(\varvec{\vartheta }\) constructed from the sample \({\varvec{\mathcal {X}}}_n\). Let \({\varvec{\mathcal {X}}}_{m}^{(\widehat{\varvec{\vartheta }})}:=\big \{\varvec{X}_1^{(\widehat{\varvec{\vartheta }})},\dots ,\varvec{X}_m^{(\widehat{\varvec{\vartheta }})}\big \}\) be a Monte Carlo sample simulated from \(F_{\widehat{\varvec{\vartheta }}}\). We propose to evaluate the composite hypothesis figuring in (17) via a test statistic from (8) computed for \({\varvec{\mathcal {X}}}_n\) and \({\varvec{\mathcal {X}}}_{m}^{(\widehat{\varvec{\vartheta }})}\). Consider the pooled sample \( {\varvec{\mathcal {X}}}_n \cup {\varvec{\mathcal {X}}}_{m}^{(\widehat{\varvec{\vartheta }})}\), and let \({\varvec{\mathcal {R}}}_{n}\) and \({\varvec{\mathcal {R}}}_{m}^{(\widehat{\varvec{\vartheta }})}\) be the collection of multivariate ranks associated with observations in \({\varvec{\mathcal {X}}}_{n}\) and \({\varvec{\mathcal {X}}}_{m}^{(\widehat{\varvec{\vartheta }})}\), respectively, in this pooled sample. Define

where \(\widehat{\phi }^{(\widehat{\varvec{\vartheta }})}_{m}\) is the empirical CF computed from \({\varvec{\mathcal {R}}}_{m}^{(\widehat{\varvec{\vartheta }})}\). In other words, \(\widehat{\phi }^{(\widehat{\varvec{\vartheta }})}_{m}(\cdot )\) is the analogue of \(\widehat{\phi }^{(0)}_m(\cdot )\) figuring in (8), with extra randomness induced by estimation of the parameter \(\varvec{\vartheta }\). A closed-form expression, analogous to that in (9), is again available if w is chosen appropriately.

As before, large values of \(\widetilde{D}_{n,m}\) indicate the violation of the null hypothesis. An important difference is that the test based on \(\widetilde{D}_{n,m}\) is no longer distribution-free, because the distribution of the artificial sample \(\varvec{\mathcal {X}}_m^{(\widehat{\varvec{\vartheta }})}\) depends on the distribution of \(\widehat{\varvec{\vartheta }}\). This issue is discussed in more detail in Remark 3. Nevertheless, the null hypothesis \(\mathcal {H}_0^C\) is rejected if

where \(\widetilde{c}_{n,m,\alpha }\) is such that the test keeps the prescribed level \(\alpha \). Namely, the critical value \(\widetilde{c}_{n,m,\alpha }\) can be obtained using the bootstrap procedure described in Algorithm 1.

The procedure is designed such that the artificial observations \(\varvec{\mathcal {X}}_{m}^{(\widehat{\varvec{\vartheta }})} =\{\varvec{X}_1^{(\widehat{\varvec{\vartheta }})},\dots ,\varvec{X}_m^{(\widehat{\varvec{\vartheta }})}\}\) are mutually conditionally independent given \(\widehat{\varvec{\vartheta }}\). Unconditionally, they are dependent on the original data \(\varvec{X}_1,\dots ,\varvec{X}_n\), so generally unconditionally dependent. Furthermore, \(\varvec{X}_j^{(\widehat{\varvec{\vartheta }})}\) is distributed according to \(F_{\widehat{\varvec{\vartheta }}}\), given \(\widehat{\varvec{\vartheta }}\), and so the unconditional distribution of \(\varvec{X}_j^{(\widehat{\varvec{\vartheta }})}\) is not \(F_{\varvec{X}}\). Hence, the pooled sample

is a collection of variables that are not independent nor identically distributed. Moreover, the latter causes that the joint distribution of the pooled sample is not exchangeable. Consequently, the pooled multivariate ranks \(\varvec{\mathcal {R}}_n \cup \varvec{\mathcal {{R}}}_m^{(\widehat{\varvec{\vartheta }})}\) are not exchangeable under the composite null hypothesis \(\mathcal {H}_0^C\). This causes that the test based on \(\widetilde{D}_{n,m}\) is not completely distribution-free for finite n, m. The proof of Theorem 1 relies on the permutation central limit theorem, which is now not applicable due to the fact that the null distribution of \(\varvec{\mathcal {R}}_n \cup \varvec{\mathcal {{R}}}_m^{(\widehat{\varvec{\vartheta }})}\) is not exchangeable (i.e., permutation invariant).

If the estimator \(\widehat{\varvec{\vartheta }}\) is consistent, then under the null hypothesis, the pooled empirical measure on \(\varvec{\mathcal { X}}_n \cup \varvec{\mathcal {X}}_m^{(\widehat{\varvec{\vartheta }})}\) converges in probability to \(F_{\varvec{X}}\). This means that the pooled sample behaves asymptotically as a random sample from \(F_{ \varvec{X}}\), and the empirical optimal transport \(\widehat{\varvec{G}}_N\) should be close to the true transport \(\varvec{G}^*\) of \(F_{\varvec{X}}\) for large n, m. However, it is expected that the distribution of \(\widetilde{D}_{n,m}\) is generally affected by the distribution of \(\widehat{\varvec{\vartheta }}\), as this is often the case in analogous goodness-of-fit procedures. A formal treatment of this would require limit theorems for the empirical transport of a collection of dependent and non-identically distributed data, that are, however, currently unavailable. Therefore, we leave the derivation of the asymptotic distribution of \(\widetilde{D}_{n,m}\) as an open problem for future work.

A special case of (17) is obtained for \(\mathcal {F}\) being a family of elliptical distributions on \(\mathbb {R}^p\) with a specified radial distribution, see Babić et al (2021, Section 1.2) for an extensive list of related tests to this problem. We say that a random vector \(\varvec{X}\) has an elliptical distribution \(\mathcal {E}_p(\varphi ,\varvec{\mu }, \varvec{\Sigma })\) if it can be represented as

where \(\varvec{\mu }\in \mathbb {R}^p\), \(\varvec{A}\) is a \(p\times p\) matrix such that \(\varvec{A} \varvec{A}^\top = \varvec{\Sigma }\), \(\varvec{S}\) is a random vector uniformly distributed on the unit sphere in \(\mathbb {R}^p\) and \(R\ge 0\) is a radial random variable with density \(\varphi \), independent of \(\varvec{S}\).

Let \(\varphi \) be specified. Then, the aim is to test that \(F_{\varvec{X}}\) is a DF of \(\mathcal {E}_p(\varphi ,\varvec{\mu }, \varvec{\Sigma })\) for some (unspecified) \(\varvec{\mu }\) and \(\varvec{\Sigma }\). A very important special case is obtained for \(\varphi \) being the density of \(\chi \) distribution with p degrees of freedom, which leads to the class of p-variate normal distributions, and testing (17) is testing normality.

Following the procedure from Sect. 4, the artificial sample \({\varvec{\mathcal {X}}}_{m}^{(\widehat{\varvec{\vartheta }})}\) could be taken as a sample simulated from \(\mathcal {E}_p(\varphi ,\widehat{\varvec{\mu }},\widehat{\varvec{\Sigma }})\) for \(\widehat{\varvec{\mu }}\) and \(\widehat{\varvec{\Sigma }}\) being some suitable estimates of \(\varvec{\mu }\) and \(\varvec{\Sigma }\), respectively. The test based on \(\widetilde{D}_{n,m}\) then proceeds as described above. For the special case of the normal distribution, one typically takes \(\widehat{\varvec{\mu }}\) and \(\widehat{\varvec{\Sigma }}\) as the sample mean and sample covariance matrix of \(\varvec{\mathcal {X}}_n\), respectively, but some more robust estimators of \(\varvec{\mu }\) and \(\varvec{\Sigma }\) can also be considered.

However, for elliptical families, we propose also an alternative method to testing (17) that removes the additional randomness in \(\widetilde{D}_{n,m}\) arising from sampling the reference sample \({\varvec{\mathcal {X}}}_{m}^{(\widehat{\varvec{\vartheta }})}\). Namely, instead of simulating independent random observations from \(\mathcal {E}_p(\varphi ,\widehat{\varvec{\mu }},\widehat{\varvec{\Sigma }})\), we propose to take the reference data as a grid that is fixed and non-random for given \(\widehat{\varvec{\mu }}\) and \(\widehat{\varvec{\Sigma }}\) and that corresponds to the structure of \(\mathcal {E}_p(\varphi ,\widehat{\varvec{\mu }},\widehat{\varvec{\Sigma }})\). The construction of the grid is based on the representation (19). Let \(H_{p}^{-1}\) be the quantile function of R (which is known because the density \(\varphi \) of R is assumed to be known). Let U be uniformly distributed on [0, 1] and independent of \(\varvec{S}\) from (19). Then, the random vector \(\varvec{Y} = U \varvec{S}\) has spherically uniform distribution in the unit ball, and \(H_p^{-1}(U)\) has the same distribution as R, so

has \(\mathcal {E}_p(\varphi ,\varvec{\mu }, \varvec{\Sigma })\) distribution.

Let \(\{\varvec{y}_i\}_{i=1}^m\) be some fixed points from the unit ball such that the uniform distribution on \(\{\varvec{y}_i\}_{i=1}^m\) is close to spherically uniform. For instance, one can take \(\{\varvec{y}_i\}_{i=1}^m\) as the points from a spherical grid \(\mathcal {G}_m^S\) described in Sect. 2. Define

It follows from the above discussion that \({\widetilde{\varvec{\mathcal {X}}}}_{m}^{(\widehat{\varvec{\vartheta }})}=\big \{\widetilde{\varvec{X}}_1^{(\widehat{\varvec{\vartheta }})},...,\widetilde{\varvec{X}}_m^{(\widehat{\varvec{\vartheta }})}\big \}\) should under the null hypothesis mimic a data set from \(\mathcal {E}_p(\varphi ,\varvec{\mu }, \varvec{\Sigma })\). Subsequently, the test statistic in (18) can be computed for \(\varvec{\mathcal {X}}_n\) and \({\widetilde{\varvec{\mathcal {X}}}}_{m}^{(\widehat{\varvec{\vartheta }})}\) and their OMT ranks \(\varvec{\mathcal {R}}_n\) and \(\widetilde{\varvec{\mathcal {R}}}_m^{(\widehat{\varvec{\vartheta }})}\). The resulting statistic is denoted as \(\widetilde{D}_{n,m}^\star = D_{n,m}(\varvec{\mathcal {R}}_n,\widetilde{\varvec{\mathcal {R}}}_m^{(\widehat{\varvec{\vartheta }})})\). Large values of \(\widetilde{D}_{n,m}^\star \) indicate the violation of the null hypothesis \(\mathcal {H}_0^C\).

Even though the set \({\widetilde{\varvec{\mathcal {X}}}}_{m}^{(\widehat{\varvec{\vartheta }})}\) is non-random for given \(\widehat{\varvec{\mu }}\) and \(\widehat{\varvec{\Sigma }}\), the randomness in the estimation of the parameters \(\varvec{\mu }\) and \(\varvec{\Sigma }\) and the consequent issues from Remark 3 remain present. Hence, the test statistic \(\widetilde{D}_{n,m}^{\star }\) is not distribution-free and its significance needs to be evaluated using a bootstrap procedure analogous to that in Algorithm 1. Nevertheless, compared to \(\widetilde{D}_{n,m}\), the test statistic \(\widetilde{D}_{n,m}^\star \) depends purely on the original data set \(\varvec{{\mathcal {X}}}_n\). Consequently, \(\widetilde{D}_{n,m}^\star \) is fully reproducible and may be more appealing in practical applications. Moreover, the simulation study presented in Sect. 6.2 indicates that the test based on \(\widetilde{D}_{n,m}^\star \) even seems to be slightly more powerful compared to \(\widetilde{D}_{n,m}\) for some settings.

A graphical comparison of a random sample drawn from \(\textsf{N}_2(\widehat{\varvec{\mu }},\widehat{\varvec{\Sigma }})\) and a set of points obtained by (20) is provided in Fig. 4 in Appendix I.

The performance of the proposed GoF test is explored in a simulation study. Both types of problems, simple and composite null hypothesis, are considered. We focus on tests of multivariate normality in dimension \(p\in \{2,3,4\}\), as this allows for a comparison between our approach and existing multivariate normality tests. However, the optimal transport of GoF test based on \(D_{n,m}\) or \(\widetilde{D}_{n,m}\) can be applied for testing any specified multivariate distribution in any dimension, as illustrated later in Sect. 7.

For given sample sizes n, m, the computation of the test statistic requires a choice of the grid set \(\mathcal {G}_N\) and the weight function w, or directly the function \(C_w\) from (9). We present results for \(C_w=C_{a,\gamma }\) from (12) for \(\gamma =2\) and various choices of \(a>0\). The grid is taken either as rectangular \(\mathcal {G}_N^R\) or spherical \(\mathcal {G}_N^S\), where both types are introduced in Sect. 2.

For a simple null hypothesis, the distribution function \(F_0\) needs to be completely specified, and \(D_{n,m}\) is distribution-free. The critical values \(c_{n,m,\alpha }\) were computed for fixed w and grid \(\mathcal {G}_N\) only once from \(10\,000\) independent replications. Remark that these Monte Carlo critical values are very close to the approximate asymptotic critical values, as already demonstrated in Table 1.

We consider \(F_0\) as the DF of the standard normal distribution \(\textsf{N}_p(\varvec{0},\varvec{I})\). Under the alternative, the data are simulated either from the uniform distribution on a specified set \(A \subset \mathbb {R}^p\), denoted as \(\textsf{U}_A\), or from a t distribution with 3 degrees of freedom. In the following, we use the notation \(\textsf{T}_p (\varvec{\mu },\varvec{\Sigma },d)\) for a p-variate t distribution with d degrees of freedom, location vector \(\varvec{\mu }\) and scale matrix \(\varvec{\Sigma }\).

The obtained results for \(p=3\) are summarized in Table 2, while results for \(p=2\) and \(p=4\) are provided in Tables 7 and 8 in Appendix K. The empirical level and power are computed from 1 000 simulations for the nominal significance level \(\alpha =0.05\). The results show that the empirical level of the test is close to the nominal value. Concerning the empirical power, the choice \(a \in \{3,4\}\) can be recommended for the rectangular ranks (grid \(\mathcal {G}_N^R\)) and \(a\in \{2,3\}\) seems to be reasonable for the spherical ranks (grid \(\mathcal {G}_N^S\)), but the optimal value of a depends on the alternative. Clearly, the spherical ranks achieve larger power than the rectangular ranks for all alternatives considered in Tables 2, 7, and 8. Therefore, we further focus solely on tests with spherical ranks \(\mathcal {G}_N^S\) in the next section.

In order to investigate the test for a composite null hypothesis \(\mathcal {H}_0^C\) in (17), we specify \(\mathcal {F}\) as a family of p-variate normal distributions. It follows from Sect. 4 that for this case, the reference sample can be generated either as

• (R)    \({\varvec{\mathcal {X}}}_{m}^{(\widehat{\varvec{\vartheta }})}\) a sample from \(\textsf{N}_p(\widehat{\varvec{\mu }}, \widehat{\varvec{\Sigma }})\), or

• (G)    \({\varvec{\widetilde{\mathcal {X}}}}_{m}^{(\widehat{\varvec{\vartheta }})}\) set of points defined in (20) as appropriately shifted and rescaled spherically uniform grid points,

leading to test statistic \(\widetilde{D}_{n,m}\) or \(\widetilde{D}_{n,m}^\star \), respectively.

The empirical size of the test is computed for samples generated from normal distributions with various mean vectors and variance matrices, while the power is investigated for data simulated from the uniform distribution on \([-1,1]^p\) and from the centered t-distribution with 3 degrees of freedom and identity scale matrix. The dimension is chosen as \(p\in \{2,3,4\}\). The obtained results are summarized in Tables 3 and 4, and in Table 9 in Appendix K.

First, we aim to investigate differences between the two approaches (R) and (G) for generating the reference sample and differences between weight functions, represented by the parameter a. For the composite hypothesis, it is necessary to use a resampling procedure to calculate the critical values, but, unfortunately, running Monte Carlo simulations together with the bootstrap approximation described in Algorithm 1 is prohibitively computationally expensive. At the same time, preliminary simulation results (not presented here) suggest that the dependency of the critical values \(\widetilde{c}_{n,m,\alpha }\) on the vector \(\varvec{\vartheta }\) is not very strong, and \(\widetilde{c}_{n,m,\alpha }\) are close to the critical values obtained for the standard normal distribution provided that \(\varvec{\Sigma }\) is not ‘extremely far’ from the identity matrix. Therefore, to obtain the power results in Table 3 for dimension \(p=2\), the bootstrap critical values were replaced by simpler critical values corresponding to the standard bivariate Normal distribution \(\mathsf {\textit{N}}_2(\varvec{0},\varvec{I})\), listed in Table 6 in Appendix K. We believe that for this initial comparison of the tuning settings, this simplified procedure is satisfactory. Indeed, Table 3 suggests that this approach works quite well for randomly generated reference sample in (R), whereas the empirical level obtained by using the grid points in (G) seems to be somewhat inflated for the non-standard Normal distribution \(\mathsf {\textit{N}}_2 (\varvec{1},\varvec{\Sigma }_{21})\), where \(\Sigma _{21}=\big ({\begin{smallmatrix} 2 & 1\\ 1 & 1 \end{smallmatrix}}\big )\), so for (G), the critical values should be computed rather from the proper bootstrap, which is used in Tables 4 and 9 for \(p=2\) and \(p\in \{3,4\}\), respectively.

Comparing the empirical power in Table 3 observed for \(\mathsf {\textit{U}}_{(-1,1)^2}\) and \(\mathsf {\textit{T}}_2(\varvec{0},\varvec{I},3)\) distributions, we conclude that a larger reference sample size (\(m=1\,000\)) is to be recommended. For this choice, the powers for the two approaches (R) and (G) are similar. The choice of the tuning parameter a seems to be crucial for the \(\mathsf {\textit{U}}_{(-1,1)^2}\) alternative where, for instance, for the random grid with \(m=200\) and \(n=80\), the power ranges from 8.0 % (\(a=0.5\)) to 71.7 % (\(a=3.5\)). For the \(\mathsf {\textit{T}}_2(\varvec{0},\varvec{I},3)\) alternative, the differences between powers for different a are smaller. In general, a from the interval [2, 3] seems to be a reasonable choice for practical usage.

Finally, Tables 4 and 9 in Appendix K compare the test based on \(\widetilde{D}_{n,m}^\star \) (with a grid type reference sample) to a set of standard multivariate normality goodness-of-fit tests for dimensions \(p=2\) and \(p\in \{3,4\}\), respectively. The results for \(\widetilde{D}_{n,m}^\star \) are computed from the proper bootstrap-based estimates of the critical values \(\widetilde{c}_{n,m,\alpha }\). Since the full bootstrap is computationally too demanding, the empirical level and power were computed using the so-called warp-speed method, see Giacomini et al. (2013). The reference normality tests are computed using R library MVN (Korkmaz et al. 2014). Namely, we present results of Mardia’s multivariate skewness and kurtosis coefficients, denoted as M3 and M4, respectively, see Mardia (1974), Henze–Zirkler’s (abbreviated as H-Z) multivariate normality test from Henze and Zirkler (1990), Royston’s multivariate test of Royston (1992), Doornik–Hansen’s (abbreviated as D-H) multivariate normality test from Doornik and Hansen (2008), and Energy multivariate normality test of Székely and Rizzo (2013).

It is discussed in Sect. 3.1 that the test statistic \({D}_{n,m}\) is, for some specific weight functions, approximately related to the two-sample energy statistic for the OMT ranks, considered for a comparison of two independent samples and rectangular grid in Deb and Sen (2023). Hence, Table 4 provides also power results for a test based on \(\widetilde{E}_{n,m}^\star = E_{n,m}(\varvec{\mathcal {R}}_n, \widetilde{\varvec{\mathcal {R}}}_m^{(\widehat{\varvec{\vartheta }})})\), where \(E_{n,m}\) is defined by (13) for \(\gamma =1\), denoted in the table as D-S. The significance was computed analogously as for the test based on \(\widetilde{D}_{n,m}^\star \), using \(m=1000\). In contrast to Deb and Sen (2023), we use spherical ranks that provided better power in Sect. 6.1. The results indicate that the power of the D-S test is comparable to \(\widetilde{D}_{n,m}^\star \) with the same sample sizes n, m.

Table 4 for \(p=2\) shows that the empirical level of most of the considered tests is close to the nominal value \(\alpha =0.05\), with lower values observed for the M4 test and slightly higher values for the Royston test. Under the considered alternatives, the results for the candidate tests are rather comparable for the \(\mathsf {\textit{T}}_2(\varvec{0},\varvec{I},3)\) alternative, while more visible differences are observed for the uniform alternative \(\mathsf {\textit{U}}_{(-1,1)^2}\). The results indicate that the Royston test yields the largest empirical power in most of the cases. Our approach with \(a=2\) has a similar power as the Doornik–Hansen and energy tests for the \(\mathsf {\textit{T}}_2(\varvec{0},\varvec{I},3)\) distribution. For the uniform \(\mathsf {\textit{U}}_{(-1,1)^2}\) distribution, the M3 test completely fails to detect the alternative. The remaining tests (including our approach) have similar power for \(n=80\) observations. The Royston test clearly outperforms all other tests for \(n=50\), and it achieves similar power as our approach (with \(a\in \{2.5,3\}\)) for \(n=20\). Similar conclusions can be made for dimensions \(p=3\) and \(p=4\) based on results from Table 9 in Appendix K.

Altogether, our approach seems to be similarly powerful as the standard multivariate normality tests and, in terms of empirical power, it is significantly outperformed only by the (slightly oversized) Royston test.

Understanding the joint distribution of socioeconomic variables is a central problem in econometrics. In particular, the relationship between income and expenditures plays a fundamental role in analyzing economic behavior and inequality. Examining the joint distribution of these two variables allows us to capture important features such as nonlinear dependence, asymmetry, and tail behavior, which are not revealed by simple correlation analysis or linear regression models.

As an example, we analyze a two-dimensional data set consisting of expenditures per student and district average income (both in USD 1000) in Californian schools (Stock and Watson 2007) in 1998–1999, plotted in Fig. 1. The data are available in the R library AER (Kleiber and Zeileis 2008). For illustrative purposes, we consider both the complete data set and a subset consisting of the first 100 observations. As can be seen from Fig. 1, for a sample size of 100, the heavy tails of the distribution are not yet fully apparent, whereas they become clearly visible in the plot of the complete data set. The first 100 observations tend to come from lower income districts, and therefore, they may not capture all features present in the entire data. Consequently, one may expect the goodness-of-fit test to yield different results in the two cases.

Figure 5 in Appendix K shows histograms of the two variables. They suggest that the marginal distribution of income exhibits positive skewness relative to the normal distribution, while the skewness of expenditures is more modest. This feature is already apparent even for the small sample size \(n=100\).

To get flexible models for the data, we define the tested class of distributions \(\mathcal {F}\) by specifying parametric models for the two marginals and a parametric copula. The marginals are considered as normal, lognormal, and Gamma, while the copula is chosen from the normal, Frank, or Gumbel families. Since it is reasonable to assume the same marginal model for both variables, this results in nine candidate models for the bivariate distribution.

For each specified \(\mathcal {F}\), the goodness-of-fit test is conducted as described in Algorithm 1. In Steps 1 and 6, the parameter vector \(\varvec{\vartheta }\) is estimated by the two-step Inference Functions for Margins (IFM) method (Joe 2005). More precisely, the marginal distribution parameters are estimated first by maximum likelihood using the function fitdistr() from the R package MASS (Venables and Ripley 2002). In the second step, the copula parameter is estimated from the resulting pseudo-observations using the function fitCopula() from the R package copula (Hofert et al. 2025; Yan 2007; Kojadinovic and Yan 2010). The artificial data are generated from the fitted distribution, in steps 2, 5, and 7, by standard procedures. Figure 1 (right panel) compares the data with contours of the distribution with lognormal marginals and a normal copula, fitted by this procedure.

The tuning parameters of the GoF test are chosen based on their good performance in the simulation study. In particular, \(m=200\), the grid is spherical and the weight parameter is \(a=2.5\), with \(B=1000\) bootstrap replications. As mentioned above, for illustrative purposes, we present results not only for the full data set of \(n=420\) observations, but also for a subset of the first \(n=100\) records. The aim is to compare the results of the test in these two cases, since features such as heavy tails and extreme observations are less pronounced in the smaller sample, whereas they become more evident when the entire data set is considered. Comparing the two cases also illustrates the behavior of the goodness-of-fit test with respect to the amount of available data. The resulting p-values are summarized in Table 5.

For the first \(n=100\) observations, all the three bivariate distributions with normal marginals provide rather a poor fit. This is consistent with the conclusions drawn from histograms in Fig. 5. In particular, the null hypothesis is rejected at the test level \(\alpha =0.05\) for the normal copula (p-value 0.004), where in combination with normal marginals, \(\mathcal {F}\) is the class of bivariate normal distributions. A borderline p-value is obtained for the remaining two copula families combined with normal marginals. On the other hand, lognormal and Gamma marginals lead to an acceptable fit for the first \(n=100\) observations when combined with all three copulas, with p-values between 0.078 (lognormal marginals with Frank copula) to 0.936 (lognormal marginals with Frank copula).

When testing the entire dataset of \(n=420\) observations, the model with lognormal marginals and normal copula seems to provide the best fit, see also the right panel of Fig. 1. A bivariate distribution with lognormal marginals and Frank copula would still be acceptable, while the other models are rejected at level \(\alpha =0.05\).

Based on the results, the distribution with lognormal marginals and a normal copula seems to be an appropriate model for the considered data set. The results showed that with only 100 observations, the goodness-of-fit test was not powerful enough to distinguish between lognormal and Gamma marginal distributions, or between normal, Frank, and Gumbel copulas. This illustrates that the power of the goodness-of-fit test may deteriorate in smaller samples, where important distributional features, such as heavy tails and extreme observations, are not yet fully manifested.

This paper deals with goodness-of-fit testing for a specified multivariate distribution and proposes a test statistic that makes the use of multivariate ranks derived from the optimal measure transport theory. We show that the test of a simple null hypothesis is distribution-free, while a composite null hypothesis requires a bootstrap approximation of the critical values.

The empirical results presented in Sect. 6 show that the proposed test achieves power comparable to that of the most powerful existing normality tests. However, the optimal transport-based goodness-of-fit test is significantly more general, as it can be applied to testing arbitrary multivariate distributions with only minor modifications. This flexibility is demonstrated in Sect. 7 for a real data example utilizing separate models for marginal distributions and a copula. The test requires only a procedure for parameter estimation and a method for generating samples from the distribution under the null hypothesis.

Authors

1. Zdeněk Hlávka
2. Šárka Hudecová
3. Simos G. Meintanis

$$\begin{aligned} g_{i,1} = x_{i,1} \cos (2\pi x_{i,2}), \quad g_{i,2} = x_{i,1} \sin (2\pi x_{i,2}), \quad i=1,\dots ,N, \end{aligned}$$

(A.1)

$$\begin{aligned} \widetilde{\varvec{X}}_i = \varvec{\mu }+ \varvec{\Sigma }^{-1/2} H_2^{-1}(\Vert \varvec{y}_i\Vert )\frac{\varvec{y_i}}{\Vert \varvec{y}_i\Vert }, \quad i=1,\dots ,m, \end{aligned}$$

(A.2)

$$ \frac{1}{n}\sum _{i=1}^{n} f({\pi _i}) - \frac{1}{N-n}\sum _{i=n+1}^N f({\pi _i}) = \frac{N}{n (N-n)} \left[ \sum _{i=1}^n f(\pi _i) - \frac{n}{N}\sum _{i=1}^N f(i)\right] . $$

$$\begin{aligned} \frac{1}{N-n}\sum _{i=n+1}^{N} f({\pi _i})&= \frac{1}{N-n}\left[ \sum _{i=1}^{N} f({\pi _i}) - \sum _{i=1}^n f({\pi _i})\right] \\&= \frac{1}{N-n}\left[ \sum _{i=1}^{N} f({i}) - \sum _{i=1}^n f({\pi _i})\right] . \end{aligned}$$

$$ \phi ^{\varvec{g}}(\varvec{t}) = \frac{1}{N}\sum _{j=1}^N \exp (\texttt {i}\varvec{t}^\top \varvec{g}_j). $$

$$ \widehat{\phi }_n(\varvec{t}) - \widehat{\phi }_m^{(0)}(\varvec{t}) = \frac{N}{m}\big [\widehat{\phi }_n(\varvec{t})-\phi ^{\varvec{g}}(\varvec{t})\big ]. $$

$$ \textsf{P}\left( \widehat{\phi }_n(\varvec{t}) = \frac{1}{n}\sum _{j=1}^n \exp ({\texttt {i}\varvec{t}^\top \varvec{g}_{\pi _j}})\right) = \frac{1}{{N \atopwithdelims ()n}} $$

$$ S_N(\varvec{t}) = \sum _{i=1}^N d_{\pi _i} \bigl [\cos (\varvec{t}^\top \varvec{g}_{i})+\sin (\varvec{t}^\top \varvec{g}_{i})\bigr ], $$

$$ d_j = {\left\{ \begin{array}{ll} \sqrt{\frac{N-n}{n N}},& i\le n,\\ -\sqrt{\frac{n}{(N-n)N}},& i>n. \end{array}\right. } $$

$$\begin{aligned} \int _F S_{N}^2(\varvec{t}) w(\varvec{t}) \textrm{d} \varvec{t} {\mathop {\rightarrow }\limits ^{D}} \int _F Z(\varvec{t})^2 w(\varvec{t}) \textrm{d} \varvec{t}, \end{aligned}$$

(J.1)

$$\begin{aligned} \overline{d}_N&:= \frac{1}{N}\sum _{i=1}^N d_i = \frac{1}{N} \sqrt{\frac{N}{nm}} \left[ n \frac{m}{N} - m \frac{n}{N}\right] =0,\\ \sigma ^2_d&:= \frac{1}{N}\sum _{i=1}^N (d_i - \overline{d}_N)^2 = \frac{1}{N}\frac{N}{nm}\left[ n \left( \frac{m}{N}\right) ^2 + m \left( \frac{n}{N}\right) ^2\right] = \frac{1}{N},\\ \max _{1\le i\le N} (d_i - \overline{d}_N)^2&= \frac{1}{N}\max \left\{ \frac{n}{m},\frac{m}{n}\right\} . \end{aligned}$$

$$ \frac{1}{N} \sum _{i=1}^N h(\varvec{g}_i) \rightarrow \int _{\mathbb {R}^p} h(\varvec{x}) \textrm{d}\nu (\varvec{x}) = \textsf{E}h(\varvec{Y}) $$

$$\begin{aligned} \overline{a}_{\varvec{t},N}&:= \frac{1}{N}\sum _{i=1}^N a_{\varvec{t}}(i) \rightarrow \textsf{E}W(\varvec{t}), \\ \sigma ^2_{a,\varvec{t}}&:=\frac{1}{N}\sum _{i=1}^N[a_{\varvec{t}}(i) - \overline{a}_{\varvec{t},N}]^2 \rightarrow \textsf{var}[ W(\varvec{t})], \end{aligned}$$

$$ \max _{1\le i\le N}[a_{\varvec{t}}(i) - \overline{a}_{\varvec{t},N}]^2 \le \max _{1\le i\le N} 2[a_{\varvec{t}}(i)^2 + \overline{a}_{\varvec{t},N}^2]\le 2(4+4)=16. $$

$$\begin{aligned} \textsf{var}S_N(\varvec{t}) = \frac{N^2}{N-1} \sigma ^2_{a,\varvec{t}} \sigma ^2_d \rightarrow \textsf{var}\, W(\varvec{t}). \end{aligned}$$

(J.2)

$$\begin{aligned} N&\frac{\max _{1\le i\le N}[a_{\varvec{t}}(i) - \overline{a}_{\varvec{t},N}]^2}{\sum _{i=1}^N [a_{\varvec{t}}(i) - \overline{a}_{\varvec{t}, N}]^2} \frac{\max _{1\le i\le N}[d_i - \overline{d}_N]^2}{\sum _{i=1}^N [d_i - \overline{d}_N]^2}\\&\quad \quad \le \frac{16}{ \frac{1}{N}\sum _{i=1}^N[a_{\varvec{t}}(i) - \overline{a}_{\varvec{t},N}]^2} \frac{1}{N} \max \left\{ \frac{n}{m},\frac{m}{n}\right\} \rightarrow 0, \end{aligned}$$

$$ \sum _{j=1}^K \lambda _j S_N(\varvec{t}_j) = \sum _{i=1}^N d_{\pi _i} \underbrace{\sum _{j=1}^K \lambda _j a_{\varvec{t}_j}(i)} _{b(i)} = \sum _{i=1}^N d_{\pi _i}b(i). $$

$$ \overline{b}_N := \frac{1}{N} \sum _{i=1}^N b (i ) = \sum _{j=1}^K \lambda _j \overline{a}_{\varvec{t}_j, N} $$

$$\begin{aligned} \sigma ^2_b&:= \frac{1}{N}\sum _{i=1}^N [b(i)-\overline{b}_N]^2 = \frac{1}{N}\sum _{i=1}^N \left[ \sum _{j=1}^K \lambda _j [a_{\varvec{t}_j}(i) - \overline{a}_{\varvec{t}_j,N}]\right] ^2 \\&=\sum _{j=1}^K \sum _{l=1}^K \lambda _j \lambda _l \frac{1}{N} \sum _{i=1}^N [a_{\varvec{t}_j}(i) - \overline{a}_{\varvec{t}_j,N}][a_{\varvec{t}_l}(i) - \overline{a}_{\varvec{t}_l,N}]\\&\rightarrow \sum _{j=1}^K \sum _{l=1}^K\lambda _j \lambda _l \textsf{cov}\big (W(\varvec{t}_j),W(\varvec{t}_l)\big ) = \sum _{j=1}^K \sum _{l=1}^K\lambda _j \lambda _l R(\varvec{t}_j,\varvec{t}_l). \end{aligned}$$

$$\begin{aligned} [b(i) -\overline{b}_N]^2 \le 16 K \sum _{j=1}^K \lambda _j^2, \end{aligned}$$

$$ N \frac{\max _{1\le i\le N}[b(i) - \overline{b}_{N}]^2}{\sum _{i=1}^N [b(i) - \overline{b}_{N}]^2} \frac{\max _{1\le i\le N}[d_i - \overline{d}_N]^2}{\sum _{i=1}^N [d_i - \overline{d}_N]^2}\rightarrow 0, $$

$$ \sum _{j=1}^K \lambda _j S_N(\varvec{t}_j) {\mathop {\rightarrow }\limits ^{D}} \textsf{N}\left( 0, \sum _{j=1}^K \sum _{l=1}^K\lambda _j \lambda _l R(\varvec{t}_j,\varvec{t}_l)\right) . $$

$$\begin{aligned} \sup _N \textsf{E}\int _F S_N(\varvec{t})^2 w(\varvec{t}) \textrm{d} \varvec{t} <\infty \end{aligned}$$

(J.3)

$$\begin{aligned} \sup _N \textsf{E}|S_N^2(\varvec{t}_1) - S_N^2(\varvec{t}_2)|\le C \Vert \varvec{t}_1-\varvec{t}_2\Vert ^\kappa , \end{aligned}$$

(J.4)

$$ \textsf{E}S_N(\varvec{t})^2 = \textsf{var}S_N(\varvec{t}) = \frac{N}{N-1} \sigma ^2_{a,\varvec{t}} \le M_1 $$

$$\begin{aligned} \textsf{E}|S_N^2(\varvec{t}_1) - S_N^2(\varvec{t}_2)&|\le \sqrt{\textsf{E}[S_N(\varvec{t}_1) - S_N(\varvec{t}_2)]^2} \sqrt{\textsf{E}[S_N(\varvec{t}_1) + S_N(\varvec{t}_2)]^2} \\&= \sqrt{\textsf{var}Q_N(\varvec{t}_1,\varvec{t}_2)} \sqrt{\textsf{var}\widetilde{Q}_N(\varvec{t}_1,\varvec{t}_2)}, \end{aligned}$$

$$\begin{aligned} Q_N(\varvec{t}_1,\varvec{t}_2)&= \sum _{i=1}^N d_{\pi _i} [a_{\varvec{t}_1}(i) -a_{\varvec{t}_2}(i)],\\ \widetilde{Q}_N(\varvec{t}_1,\varvec{t}_2)&= \sum _{i=1}^N d_{\pi _i}[a_{\varvec{t}_1}(i) +a_{\varvec{t}_2}(i)]. \end{aligned}$$

$$\begin{aligned} |a_{\varvec{t}_1}(i) -a_{\varvec{t}_2}(i)|&\le |\cos (\varvec{t}_1^\top \varvec{g}_i) - \cos (\varvec{t}_2^\top \varvec{g}_i) | + |\sin (\varvec{t}_1^\top \varvec{g}_i) - \sin (\varvec{t}_2^\top \varvec{g}_i) | \\&\le 2 \Vert \varvec{t}_1 - \varvec{t}_2\Vert \Vert \varvec{g}_i\Vert . \end{aligned}$$

$$\begin{aligned} \textsf{var}Q_N(\varvec{t}_1,\varvec{t}_2)&= \frac{1}{N-1}\sum _{i=1}^N[a_{\varvec{t}_1}(i) -a_{\varvec{t}_2}(i)]^2 - \frac{N}{N-1}[\overline{a}_{\varvec{t}_1,N} - \overline{a}_{\varvec{t}_2,N}]^2 \\&\le \frac{N}{N-1} 4 \Vert \varvec{t}_1 - \varvec{t}_2\Vert ^2 \frac{1}{N} \sum _{i=1}^N \Vert \varvec{g}_i\Vert ^2 < M_2 \Vert \varvec{t}_1 - \varvec{t}_2\Vert ^2, \end{aligned}$$

$$ \int _{\mathbb {R}^p} \cos (\varvec{t}^\top \varvec{x})\sin (\varvec{t}^\top \varvec{y}) w(\varvec{t})\textrm{d} \varvec{t} = 0. $$

$$ D_{n,m} = \int _{\mathbb {R}^p} Z_{n,m}^2(\varvec{t}) w(\varvec{t})\textrm{d}\varvec{t} $$

$$\begin{aligned} Z_{n,m} (\varvec{t})&= \sqrt{\frac{nm}{n+m}} \left[ \textrm{Re}(\widehat{\phi }_n)(\varvec{t}) +\textrm{Im}(\widehat{\phi }_n)(\varvec{t}) - \textrm{Re}(\widehat{\phi }_m^{(0)})(\varvec{t}) -\textrm{Im}(\widehat{\phi }_m^{(0)})(\varvec{t}) \right] \\&= \sqrt{\frac{nm}{n+m}}\Bigg \{ \frac{1}{n}\sum _{i=1}^n \Big [\cos \big (\varvec{t}^\top \widehat{\varvec{G}}(\varvec{X}_i)\big ) +\sin \big (\varvec{t}^\top \widehat{\varvec{G}}(\varvec{X}_i)\big ) \Big ]\\&\quad - \frac{1}{m}\sum _{i=1}^m \Big [\cos \big (\varvec{t}^\top \widehat{\varvec{G}}(\varvec{X}_i^{(0)})\big ) +\sin \big (\varvec{t}^\top \widehat{\varvec{G}}(\varvec{X}_i^{(0)})\big )\Big ] \Bigg \}. \end{aligned}$$

$$ \sqrt{\frac{nm}{n+m}}\Bigg \{ \frac{1}{n}\sum _{i=1}^n \bigl [\cos (\varvec{t}^\top \varvec{g}_{\pi _i}) +\sin (\varvec{t}^\top \varvec{g}_{\pi _i}) \bigr ]- \frac{1}{m}\sum _{i=n+1}^{N} \bigl [\cos (\varvec{t}^\top \varvec{g}_{\pi _i}) +\sin (\varvec{t}^\top \varvec{g}_{\pi _i})\bigr ] \Bigg \}. $$

$$ \sqrt{\frac{n+m}{nm}} \sum _{i=1}^N c_i \bigl [\cos (\varvec{t}^\top \varvec{g}_{\pi _i})+\sin (\varvec{t}^\top \varvec{g}_{\pi _i})\bigr ] $$

$$\begin{aligned} & Z_{n,m}^{(0)}(\varvec{t}) = \sqrt{\frac{n+m}{nm}} \sum _{i=1}^N c_{\pi _i} \bigl [\cos (\varvec{t}^\top \varvec{g}_{i})+\sin (\varvec{t}^\top \varvec{g}_{i})\bigr ]\nonumber \\ & = \sum _{i=1}^N d_{\pi _i} \bigl [\cos (\varvec{t}^\top \varvec{g}_{i})+\sin (\varvec{t}^\top \varvec{g}_{i})\bigr ]. \end{aligned}$$

(J.5)

$$ \int _F [Z_{n,m}^{(0)}(\varvec{t})]^2 w(\varvec{t}) \textrm{d} \varvec{t} {\mathop {\rightarrow }\limits ^{D}} \int _F Z^2(\varvec{t}) w(\varvec{t}) \textrm{d} \varvec{t} $$

$$ |C_w(\varvec{x}_1 -\varvec{y}_1) - C_w(\varvec{x}_2 -\varvec{y}_2)| \le K \left( \Vert \varvec{x}_1 - \varvec{x}_2\Vert + \Vert \varvec{y}_1 -\varvec{y}_2\Vert \right) . $$

$$ \left| \cos (a) - \cos (b) \right| = \left| 2 \sin \left( \frac{a+b}{2}\right) \sin \left( \frac{a-b}{2}\right) \right| \le 2 \left| \sin \left( \frac{a-b}{2}\right) \right| \le |a-b|. $$

$$\begin{aligned} |C_w(\varvec{x}) - C_w(\varvec{y})|&\le \int _{\mathbb {R}^p} |\cos (\varvec{t}^\top \varvec{x}) - \cos (\varvec{t}^\top \varvec{y})| w(\varvec{t}) \textrm{d} \varvec{t} \\&\le \int _{\mathbb {R}^p} |\varvec{t}^\top (\varvec{x} -\varvec{y})| w(\varvec{t}) \textrm{d} \varvec{t} \\&\le \Vert \varvec{x} -\varvec{y} \Vert \int _{\mathbb {R}^p}\Vert \varvec{t}\Vert w(\varvec{t}) \textrm{d} \varvec{t} = K \Vert \varvec{x} -\varvec{y} \Vert . \end{aligned}$$

$$ \frac{1}{n} \sum _{i=1}^n \big \Vert \widehat{\varvec{G}}_N(\varvec{X}_i) - \varvec{G}(\varvec{X}_i) \big \Vert \rightarrow 0, \quad \frac{1}{m}\sum _{j=1}^m \big \Vert \widehat{\varvec{G}}_N(\varvec{X}_i^{(0)}) - \varvec{G}(\varvec{X}_i^{(0)}) \big \Vert \rightarrow 0 $$

$$ \frac{1}{N}\left[ \sum _{i=1}^n \Vert \widehat{\varvec{G}}_N(\varvec{X}_i) - {\varvec{G}}(\varvec{X}_i) \Vert ^2 + \sum _{j=1}^m \Vert \widehat{\varvec{G}}_N(\varvec{X}_i^{(0)}) - {\varvec{G}}(\varvec{X}_i^{(0)}) \Vert ^2 \right] \rightarrow 0 $$

$$ \frac{D_{n,m}}{n+m} = \underbrace{\frac{m}{n (n+m)^2} S_1}_{J_1} + \underbrace{\frac{n}{m (n+m)^2} S_2}_{J_2}- 2\cdot \underbrace{\frac{1}{(n+m)^2} \cdot S_3}_{J_3}, $$

$$\begin{aligned} S_1&= \sum _{j=1}^n\sum _{k=1}^n C_w(\varvec{R}_{j}-\varvec{R}_{k}),\\ S_2&= \sum _{j=1}^{m}\sum _{k=1}^m C_w(\varvec{R}^{(0)}_{j}-\varvec{R}^{(0)}_{k}), \\ S_3&= \sum _{j=1}^n \sum _{k=1}^{m} C_w(\varvec{R}_{j}-\varvec{R}^{(0)}_{k}). \end{aligned}$$

$$\begin{aligned}&|C_w(\varvec{R}_j - \varvec{R}_k ) - C_w(\varvec{G}(\varvec{X}_i) - \varvec{G}(\varvec{X}_k))| \le K\left( \Vert \widehat{\varvec{G}}_N(\varvec{X}_j) - \varvec{G}(\varvec{X}_j) \Vert \right. \\&\left. + \Vert \widehat{\varvec{G}}_N(\varvec{X}_k) - \varvec{G}(\varvec{X}_k) \Vert \right) , \end{aligned}$$

$$ S_1 \ge \sum _{j=1}^n\sum _{k=1}^n C_w(\varvec{G}(\varvec{X}_j) - \varvec{G}(\varvec{X}_k)) - 2K\cdot n\sum _{j=1}^n \Vert \widehat{\varvec{G}}_N(\varvec{X}_j) - \varvec{G}(\varvec{X}_j) \Vert . $$

$$\begin{aligned} \liminf _{n\rightarrow \infty } J_1&\ge \liminf _{n \rightarrow \infty } \frac{m}{n (n+m)^2} \sum _{j=1}^n\sum _{k=1}^n C_w(\varvec{G}(\varvec{X}_j) - \varvec{G}(\varvec{X}_k)) \\&\quad -2K \limsup _{n\rightarrow \infty } \frac{m n}{n (n+m)^2} \sum _{j=1}^n \Vert \widehat{\varvec{G}}_N(\varvec{X}_j) - \varvec{G}(\varvec{X}_j) \Vert \\&= \frac{\lambda }{(1+\lambda )^2} \textsf{E}C_w(\varvec{G}(\varvec{X}_j) - \varvec{G}(\varvec{X}_k)), \end{aligned}$$

$$ S_1 \le \sum _{j=1}^n\sum _{k=1}^n C_w(\varvec{G}(\varvec{X}_j) - \varvec{G}(\varvec{X}_k)) + 2K\cdot n\sum _{j=1}^n \Vert \widehat{\varvec{G}}_N(\varvec{X}_j) - \varvec{G}(\varvec{X}_j) \Vert $$

$$ J_1 {\mathop {\rightarrow }\limits ^{\textsf{P}}} \frac{\lambda }{(1+\lambda )^2} \textsf{E}C_w(\varvec{G}(\varvec{X}_i) - \varvec{G}(\varvec{X}_k)) $$

$$ J_2 {\mathop {\rightarrow }\limits ^{\textsf{P}}} \frac{\lambda }{(1+\lambda )^2} \textsf{E}C_w(\varvec{G}(\varvec{X}_j^{0)}) - \varvec{G}(\varvec{X}_k^{(0)})) $$

$$ J_3 {\mathop {\rightarrow }\limits ^{\textsf{P}}} \frac{\lambda }{(1+\lambda )^2} \textsf{E}C_w(\varvec{G}(\varvec{X}_j^{0)}) - \varvec{G}(\varvec{X}_k)). $$

$$\begin{aligned} \frac{D_{n,m}}{n+m}&{\mathop {\rightarrow }\limits ^{\textsf{P}}} \frac{\lambda }{(1+\lambda )^2} \left[ \textsf{E}C_w(\varvec{U}_1- \varvec{U}_2)+ \textsf{E}C_w(\varvec{V}_1 - \varvec{V}_2) - 2 \textsf{E}C_w(\varvec{U}_1 - \varvec{V}_2)\right] \nonumber \\&= \frac{\lambda }{(1+\lambda )^2} \int _{\mathbb {R}^p}\left| \exp (\texttt {i}\varvec{t}^\top \varvec{U}_1) - \exp (\texttt {i}\varvec{t}^\top \varvec{V}_1)\right| ^2 w(\varvec{t}) \textrm{d} \varvec{t},\nonumber \\&= \frac{\lambda }{(1+\lambda )^2} \int _{\mathbb {R}^p}\left| \varphi _{\varvec{U}}(\varvec{t}) - \varphi _{\varvec{V}}(\varvec{t})\right| ^2 w(\varvec{t}) \textrm{d} \varvec{t}, \end{aligned}$$

(J.6)

Full size table

Full size table

Full size table

Full size table

Hlávka, Z., Hudecová, Š. & Meintanis, S.G. Optimal transport-based multivariate goodness-of-fit tests. TEST (2026). https://doi.org/10.1007/s11749-026-01026-7

Download citation

• Received: 06 June 2025

• Accepted: 25 May 2026

• Published: 23 July 2026

• Version of record: 23 July 2026

• DOI: https://doi.org/10.1007/s11749-026-01026-7