Two-sample tests are one of the most frequently used methods for statistical inference. While rooted in classical statistics, the two-sample problem is relevant to numerous cutting-edge applications, including high-energy physics [11], computer vision [14], and genome-wide expression analysis [36].

We begin with two samples X and Y, which may be either univariate or multivariate. In the nonparametric setting, we make minimal assumptions regarding the distributions F and G used to generate X and Y, as we test the null hypothesis F = G. In Section 1.2, we discuss the landscape of existing methods.

While we face no shortage of effective two-sample tests, certain factors may hinder their real-world applicability. For one, we find some nonparametric tests to be more parsimonious than others against the range of potential alternatives. Relatively speaking, a method may excel at detecting location and scale shifts but struggle to catch bimodality when mean and variance are held constant, as one example. We explore this phenomenon in Section 5, showing that the relative performance of well-known univariate tests at detecting a location shift can be reversed by a suitable choice of distribution family. This is unintuitive, as one might expect power against location alternatives to be nearly independent of family.

Furthermore, many tests offer non-transparent rejections of the null hypothesis. While data visualizations and summary statistics offer some degree of explanation for a test, such analyses do not quantify the contribution of various data features to the test’s rejection. For multivariate tests, this problem is compounded, as distributional alternatives may easily exceed human intuition for higher dimensions.

Here, we formulate a new nonparametric two-sample test called AUGUST, an abbreviation of Augmented CDF for Uniform Statistic Transformation. Our method explicitly tests for multiple orthogonal sources of distributional inequality up to a predetermined resolution d, giving it power against a wide range of alternatives. Upon rejection of the null, both resolution control and decomposition into orthogonal signals allow for interpretation of how equality in distribution between X and Y has failed. To promote ease of use, we provide asymptotic theory as well as algorithmic optimizations.

Some well-known rank-based tests are designed for the univariate context, including [13, 27, 33]. Other approaches explicitly refer to a distance between the empirical cumulative distribution functions of X and Y. For instance, [1, 12, 15, 26] are all widely known. The recent test of [16] somewhat combines [1] and [15].

For nonparametric multivariate methods, the range of approaches is quite broad. Tests based on geometric graphs, including [9, 10, 18, 39], have had considerable success [5]. Ball divergence [3, 36] and energy distance [2, 42] are also popular names. Among the family of kernel-based tests are [11, 20, 21, 22, 25, 41, 44], while [4, 24, 35, 40] use generalized ranks. Additional recent work includes [6, 7, 28, 31].

As for interpretable methods, one line of work [25, 44] proposes feature selection in the framework of maximum mean discrepancy [21]. Similar in principle is [34]. The test of [31] inherits the interpretability of the classifier on which the test is based. These methods provide a global type of interpretability, compared to the study of local significant differences [17, 23]. To put AUGUST in context, we propose that our method is more geometrical than feature-selecting tests, but more global than local significant difference methods.

Regarding methodological relatives, the use of subsampling has been considered for the two-sample location and scale problems [32, 37, 38]. In addition, our method builds on the binary expansion framework [45, 46], which has applications to resolution-based nonparametric models of dependency [8]. Using an underlying binary expansion framework, we furnish substantial methodological and algorithmic innovations to yield a practical test in the two-sample context.

We begin by introducing our procedure in the context of univariate data. Given independent samples { X i } i = 1 m and { Y i } i = 1 n , where X i ∼ G and Y i ∼ F, we are interested in testing H 0 : F = G vs. H a : F ≠ G . For our purposes, we will assume that F and G are absolutely continuous functions. We adopt boldface type in { X i } i = 1 m and { Y i } i = 1 n to indicate that these are collections of random quantities. In addition, we use blackboard bold P ( A ) to refer to the probability of an event A.

To begin, imagine a one-sample setting where we test H 0 : X i ∼ F, with F known. It is a well-known result that X i ∼ F if and only if the transformed variables { F ( X i ) : i ∈ [ m ] } follow a Uniform ( 0 , 1 ) distribution. Given this fact, we can test the goodness-of-fit of F by testing of the uniformity of the collection { F ( X i ) : i ∈ [ m ] }. Moreover, examining how { F ( X i ) : i ∈ [ m ] } fails to be uniform indicates why F does not fit the distribution of X.

Returning to the two-sample setting, the same intuition holds true: we might construct transformed variables that are nearly uniform in [ 0 , 1 ] when F = G, and that are not uniform otherwise. When the distributions of the two samples are different, the way that uniformity fails should be informative.

Given the fact that the transformed variables { G ( X i ) : i ∈ [ m ] } follow a uniform distribution, an intuitive choice would be { F ˆ Y ( X i ) : i ∈ [ m ] }, where F ˆ Y is the empirical cumulative distribution function of Y: F ˆ Y ( t ) = 1 n ∑ i = 1 n I ( Y i ≤ t ) . The binary expansion testing framework introduced in [45] provides a way to test { F ˆ Y ( X i ) : i ∈ [ m ] } for uniformity up to a given binary depth d, which is equivalent to testing multinomial uniformity over dyadic fractions { 1 / 2 d , … , 1 }. In particular, we define the random vector P of length 2 d such that, for 1 ≤ i ≤ 2 d , P i = # { k : F ˆ Y ( X k ) ∈ i − 1 2 d , i 2 d } m . That is, P counts the number of transformed observations falling in dyadic intervals of width 1 / 2 d . The associated vector S = H 2 d P is said to contain symmetry statistics, where H 2 d is the Hadamard matrix of size 2 d according to Sylvester’s construction. As the top row of H 2 d contains only ones, the first coordinate of S is always equal to ∑ i = 1 2 d P i = 1, and we may as well restrict our attention to S − 1 , dropping the first component. As shown in [45], S − 1 is a sufficient statistic for uniformity in the one sample setting, and the binary expansion test based on S − 1 achieves the minimax rate in sample size required for power against a wide variety of alternatives.

We can think of S − 1 in a signal-processing context: the Hadamard transform maps the vector of cell probabilities P in the physical domain to the vector of symmetries S − 1 in the frequency domain. This transformation is advantageous since, in the one sample setting, the entries of S − 1 have mean zero and are pairwise uncorrelated under the null. As a result, fluctuations of S − 1 away from 0 2 d − 1 unambiguously support the alternative, and the coordinates of S − 1 are interpretable as orthogonal signals of nonuniformity. Moreover, the vector P always satisfies ∑ i = 1 2 d P i = 1, meaning that the mass of P is constrained to a ( 2 d − 1 )-dimensional hyperplane in R 2 d . In contrast, the vector S − 1 is non-degenerate and summarizes the same information about non-uniformity with greater efficiency. We elaborate on the interpretability of S − 1 in Section 2.4.

One possible choice of test statistic is the quantity S = ‖ S − 1 ‖ 2 2 . A test based on S is essentially a χ 2 test and has decent power at detecting F ≠ G. However, we can substantially improve the power by modifying our construction of  P.

For our testing purposes, recall that we are only interested in the uniformity of { F ˆ Y ( X i ) , i ∈ [ m ] } up to binary depth d. The range of F ˆ Y ( x ) as a function of x comprises n + 1 possible values, namely, 0 , 1 / n , … , 1. However, in our construction of cell counts P, the collection { F ˆ Y ( X i ) , i ∈ [ m ] } is binned across 2 d -many dyadic intervals of depth d. Whenever 2 d < m, some distinct values in the range of F ˆ Y ( · ) correspond to the same dyadic interval by the pigeonhole principle, which indicates that a coarser transformation than F ˆ Y ( · ) should work at least as well, and possibly better. Our approach is to consider transformed variables F ˆ Y ∗ ( X i ) based on a small, random subsample Y ∗ of some size r from Y. The following discussion makes this alternate process explicit. In addition, we comment on the success of this approach in Section 4, and we include empirical power comparisons against the non-subsampled transformation F ˆ Y ( · ) in the supplementary materials.

Let Y ∗ be a random subsample from Y of size r = 2 d + 1 − 1. We explain this choice of r momentarily. For any x ∈ R and integer 1 ≤ k ≤ 2 d , let p k ( x ) be the probability, conditional on Y, that either 2 k − 2 or 2 k − 1 elements of Y ∗ are less than or equal to x. The probabilities p k ( x ) are essentially hypergeometric and simple to compute: p k ( x ) = # { i : Y i ≤ x } 2 k − 2 # { i : Y i > x } 2 d + 1 − 1 − ( 2 k − 2 ) n 2 d + 1 − 1 + # { i : Y i ≤ x } 2 k − 1 # { i : Y i > x } 2 d + 1 − 1 − ( 2 k − 1 ) n 2 d + 1 − 1 . Using the scalar function p k ( · ), we define the P : R → R 2 d such that, for each coordinate k, P k ( x ) = p k ( x ) , for 1 ≤ k ≤ 2 d . It holds that F ˆ Y ∗ ( x ) ∈ ( k − 1 ) / 2 d , k / 2 d when exactly 2 k − 2 or 2 k − 1 subsampled elements in Y ∗ are less than or equal to x. Therefore, we could equally say that for 1 ≤ k ≤ 2 d , (2.1) P k ( x ) = P F ˆ Y ∗ ( x ) ∈ k − 1 2 d , k 2 d | Y . It is in precisely this sense that P ( x ) can be considered an augmented cumulative distribution function: instead of mapping x to a single value in the unit interval, x ↦ P ( x ) maps x to a distribution. Moreover, this characterization explains the choice of subsample size r = 2 d + 1 − 1. Any r satisfying r = 2 q − 1, q ≥ d, guarantees that the discrete random variable F ˆ Y ∗ ( x ) has the same number of point masses inside every interval of the form ( k − 1 ) / 2 d , k / 2 d . In Section 4, we give further intuition behind the meaning of q, and in the supplementary materials, we asses our default choice of q = d + 1 empirically.

To collect information about every X i , we define the vector P X to be the average of all P ( X i ): P X = 1 m ∑ i = 1 m P ( X i ) . Given that the formula for p k ( x ) is computed from hypergeometric probabilities, we refer the coordinates of P X as hypergeometric cell probabilities. Just as we expect the distribution of the transformed variables { F ˆ Y ( X i ) : i ∈ [ m ] } to be uniform under the null, we expect the mass of P X to be nearly uniform over its coordinates. The vector of symmetry statistics S X = ( H 2 d P X ) − 1 quantifies non-uniformity in P X .

Importantly, the cell probabilities in P ( x ) are computed with reference to a subsampling procedure, but without actually subsampling. As the discussion above suggests, these probabilities could indeed be approximated by a bootstrap procedure: take many subsamples Y ∗ of size 2 d + 1 − 1 from Y, compute F ˆ Y ∗ ( x ) each time, and bin the results as cell counts at intervals of 1 / 2 d . The exact cell probabilities P ( x ) derived above are the limiting values of this bootstrap procedure as the number of subsamples tends to infinity. The following theorem makes this result explicit.

Theorem 1 shows that the hypergeometric cell probabilities are equivalent to the limiting values of a certain bootstrap procedure. Effectively, one could say that actual subsampling is a valid way to approximate P X . In practice, it is much faster to directly compute the limiting hypergeometric probabilities. While the procedure described in this subsection is more complicated than the original approach from Section 2.1, we achieve superior power using the augmented cumulative distribution function introduced here, which we illustrate empirically in the supplementary materials. In Section 5, we provide comparisons of empirical power against well-known nonparametric tests.

To combine information on all forms of asymmetry, we propose the statistic S = − S X T S Y , with S Y defined analogously to S X by reversing the roles of the two samples. First, this choice of statistic has the advantage of treating the X and Y samples symmetrically. This is desirable because it would be counterintuitive for the value of S to change when the roles of X and Y are switched. In addition, this statistic is a continuous function of the concatenated vector ( S X T , S Y T ) T , and in Theorem 3, we state the asymptotic distribution of ( S X T , S Y T ) T in the case of univariate data. For the multivariate AUGUST test, which is described in Section 3.2, we use permutation for p-value calculation.

The negative sign in − S X T S Y comes from the fact that S X and S Y typically have opposite signs in the case of distributional difference, and we wish the critical values of S to be positive. The proposition below gives intuition for this phenomenon in the context of a location shift.

Informally, we could say the following: if X is to the left of Y, then Y is to the right of X, and the symmetry statistic detecting left/right imbalance will be positive in S X and negative in S Y . As shown in Section 5, the negative inner product − S X T S Y gives good power against a wide range of distributional alternatives. In addition, see Section 4 for exploration of the asymptotic properties of ( S X T , S Y T ) T .

One may use entries of S X and S Y to interpret the outcome of the AUGUST test based on S. With respect to S X , the sample Y serves as the reference sample, meaning that information from S X allows us to make statements about how points of X fall relative to the distribution of Y. Each entry in the vector S X describes the non-uniformity of P X with respect to a row of H 2 d . In particular, the largest entries of S X in absolute value tell us the sources of greatest asymmetry in P X .

Before performing the test based on S, we must first choose some resolution d, which determines the scale on which the test will be sensitive. For convenience, let H ˜ 2 d denote the Hadamard matrix of size 2 d according to Sylvester’s construction, without the first row, which is a row of all ones. Now, the depth d = 1 is sensitive primarily left/right imbalance. When d = 1, the 2 1 = 2 entries of P X roughly correspond to the fraction of the X sample falling above or below the median of Y. In this case, the only symmetry statistic is the product of H ˜ 2 = 1 − 1 with P X , namely the difference between the two components of P X .

For d = 2, both 1 1 − 1 − 1 and 1 − 1 − 1 1 are (necessarily orthogonal) rows of H ˜ 4 . When multiplied by P X , the first of these rows produces a statistic for left/right imbalance, similar to the d = 1 case, while the latter row detects differences in scale. Larger values of d detect more granular varieties of imbalance. We use a depth of d = 2 in our real data example, and d = 3 in simulated power comparisons. (See the supplementary materials for an empirical comparison across depths.) In [46], it is shown that a depth of d = 3 is sufficient for a symmetry statistic-based test of independence to outperform both distance correlation and the F-test, which are known to be optimal, in detecting correlation in bivariate normal distributions.

Higher depths d > 3 can be useful for alternatives that are extremely close in the Kolmogorov–Smirnov metric but have densities that are bounded apart in the uniform norm. As one example, we may have X sampled from Uniform ( 0 , 1 ) and Y sampled from a high frequency square wave distribution with the same support. In Section 6, we use symmetry statistics in visualizations of NBA shooting data. As an additional example, a step-by-step interpretation on simulated data is provided in the supplementary materials.

Below, Algorithms 1 and 2 formalize the steps to calculating the AUGUST statistic outlined in earlier sections. In terms of prior notation, Algorithm 1 computes the vector P ( x ), and Algorithm 2 calculates the overall test statistic S = − S X T S Y . Recall that we use H 2 d to refer to the Hadamard matrix of size 2 d according to Sylvester’s construction, and for a matrix M, we use ( M ) − 1 to refer to M without its first row.

The two samples X and Y have sizes m and n, respectively. Treating d as a constant, Algorithm 2 requires O ( m n ) elementary operations. This is due to the calculation of K = # { i : V i ≤ x } in the Algorithm 1, which necessitates iterating over all entries of V each time that Algorithm 1 ( V , d , x ) is called.

However, by first sorting the concatenated X and Y samples, it is possible to reduce the running time to O ( ( m + n ) log ( m + n ) ) operations.

This improved algorithm, named AUGUST+, is stated explicitly in the supplementary materials. The time complexity is asymptotically equivalent to that of an efficient sorting algorithm applied to the concatenated data. The constant factor in this comparison depends on the resolution d, which is assumed constant in Theorem 2. In terms of storage, the AUGUST+ algorithm defines only one array whose length depends on m and n. This array has dimension 2 × ( m + n ), meaning the space requirement is linear in the combined sample size. In Section 6, we use this algorithm to perform our two-sample test on a large data set on the order of 10 6 observations.

With an appropriate transformation, we can extend the univariate test to the problem of multivariate two-sample testing. For the purposes of this subsection, let X = { X i } i = 1 m be an independent sample from multivariate distribution G, and let Y = { Y j } j = 1 n be an independent sample from multivariate distribution F, with F and G defined on R k and k ≥ 2. We adapt an approach that could be appropriately named mutual Mahalanobis distance.

Given a mean μ ∈ R k and invertible k × k covariance matrix Σ, recall that the Mahalanobis distance of x ∈ R k from μ with respect to Σ is M D ( x ; μ , Σ ) = ( x − μ ) T Σ − 1 ( x − μ ) 1 / 2 . Let μ ˆ X and Σ ˆ X be the sample mean and sample covariance matrix of X, where we assume Σ ˆ X is nonsingular. Consider the transformed collections X ˜ ( X ) = M D ( X i ; μ ˆ X , Σ ˆ X ) : 1 ≤ i ≤ m Y ˜ ( X ) = M D ( Y j ; μ ˆ X , Σ ˆ X ) : 1 ≤ j ≤ n , where the superscript ( X ) indicates that means and covariances are estimated using the X sample. If X and Y come from the same multivariate distribution, then the collections X ˜ ( X ) and Y ˜ ( X ) should have similar univariate distributions. As a result, at a given depth d, it is reasonable to test the univariate samples X ˜ ( X ) and Y ˜ ( X ) as an assessment of the distributional equality of the multivariate samples X and Y. Our choice of Mahalanobis distance is also motivated by the wide class of nonparametric tests based on data depth, particularly Mahalanobis depth [5, 30].

From the univariate method, recall that vectors of symmetry statistics quantify regions of imbalance between the univariate samples, as imbalances in distribution appear as non-uniformity in the vector of cell probabilities. Under a Mahalanobis distance transformation, cells in the domain of X ˜ ( X ) and Y ˜ ( X ) correspond to nested elliptical rings centered on μ ˆ X . This principle extends interpretability of symmetry statistics to the multivariate case.

As in the univariate case, it is desirable for the test statistic to be invariant to the transposition of X and Y. To achieve this, we can use the statistic S m u l t i = max ( AUGUST X ˜ ( X ) , Y ˜ ( X ) , d , AUGUST X ˜ ( Y ) , Y ˜ ( Y ) , d ) wherein we use both possible Mahalanobis distance transformations for X and Y, compute two test statistics, and take the maximum. Aside from ensuring transposition invariance, the simultaneous use of the two test statistics is important for detecting some asymmetric alternatives. As a contrived example, take k = 2; suppose X comprises m independent samples of the bivariate standard normal ( Z 1 , Z 2 ) T ∼ N 2 ( 0 , I 2 ), while Y comprises n independent samples of 2 − 1 / 2 ( χ 2 , χ 2 ) T , where χ 2 follows a chi distribution with two degrees of freedom. In this case, both M D ( ( Z 1 , Z 2 ) T ; 0 , I 2 ) and M D ( 2 − 1 / 2 ( χ 2 , χ 2 ) T ; 0 , I 2 ) follow a chi distribution with two degrees of freedom, which indicates that the test based solely on AUGUST X ˜ ( X ) , Y ˜ ( X ) , d is powerless. (In this highly degenerate situation, the transposed statistic AUGUST X ˜ ( Y ) , Y ˜ ( Y ) , d is technically undefined, because Y has no variance in the ( 1 , − 1 ) T direction.)

In practice, we calculate p-values for the multivariate statistic using permutation. While this current necessity sacrifices the computational advantage of the asymptotic result Theorem 3, the multivariate method may nonetheless take advantage of Theorem 2, as it is built upon the univariate procedure. As we show in Section 5, a depth of d = 2 is sufficiently large to detect common multivariate alternatives with empirical power comparable to existing tests.

Here, we compare AUGUST to a sampling of other nonparametric two-sample tests: Kolmogorov–Smirnov distance [26], Wasserstein distance [15], energy distance [42], and the recent DTS [16]. For these simulations, we use a sample size of n = m = 128, and for our resolution-based test, we set a depth of d = 3. For all tests, we use a p-value cutoff of α = 0.05. Simulation results are graphed in Fig. 1.

The first two plots of Fig. 1 correspond to normal and Laplace location alternatives, situations where differences in the first distributional moment are most diagnostic. Third, we have a symmetric beta versus asymmetric beta alternative, and fourth, we include a Laplace scale family. The last two plots of Fig. 1 focus on families with identical first and second moments: normal versus mean-centered gamma in the fifth position, and standard normal versus symmetric, variance-scaled normal mixture in the sixth. For this final alternative distribution, samples are generated by first drawing from a symmetric mixture of normals with unit variance and then dividing by the theoretical standard deviation of the mixture distribution.

For the location alternatives, the power of each method depends on the shape of the distribution. DTS, Wasserstein, and energy distance tests perform slightly better than ours for normal and beta distributions, and ours in turn outperforms Kolmogorov–Smirnov. In contrast, for a Laplace location shift, Kolmogorov–Smirnov outperforms every test, with our test in second place and DTS last. For the Laplace scale family, Kolmogorov–Smirnov performs poorly, with DTS and our test leading. DTS has the edge on the gamma skewness family, while we outperform all other tests at detecting normal versus symmetric normal mixture.

As expected, no single test performs best in all situations. Even for simple alternatives such as location families, the precise shape of the distribution is highly influential as to the tests’ relative performance. In fact, the performance rankings of DTS, Wasserstein, energy distance, and Kolmogorov–Smirnov in the Laplace location trials are exactly reversed compared to the normal location trials. We theorize that because the symmetry statistics S X and S Y are weighted equally in every coordinate, AUGUST is very parsimonious towards the range of potential alternatives. In contrast, other univariate methods show relatively greater sensitivity to location and scale shifts, but may be less robust against more obscure alternatives.

In Fig. 2, we compare our multivariate resolution-based test at depth d = 2 to some other well-known nonparametric multivariate two-sample tests. We perform these simulations in a low-dimensional context with k = 2, using sample size n = m = 128 and cutoff α = 0.05. In particular, we again consider the energy distance test of [42], as well as the classifier test of [31], the generalized edge-count method of [9], and the ball divergence test of [36], where the choice of these comparisons is inspired by [41]. For [9], we use a 5-minimum spanning tree based on Euclidean interpoint distance.

We consider a variety of alternatives. In order:

In Fig. 2, we see that the energy and ball divergence tests dominate the other methods when mean shift is a factor (i.e. in the normal location and log-normal families). On a scale alternative, AUGUST has the best power, with ball divergence at a close second. In contrast, for correlation, rotation, and multimodal alternatives, the edge-count test has superior power, with ball divergence and energy distance coming at or near last place.

Overall, our test is robust against a wide range of possible alternatives, and it has particularly high performance against a scale alternative, where it outperforms all other methods considered. We theorize that, in part, this is because some of the other methods rely heavily on interpoint distances. The scale alternative does not result in good separation between X and Y, meaning that interpoint distances are not as diagnostic as they would be in, say, a location shift.

In the supplementary materials, we include additional comparisons with k = 5 and k = 20, keeping the sample size fixed with m = n = 128. Performance follows the same general pattern as when k = 2: we lag against location alternatives but are very strong at scale, with no universal winner across all scenarios.