MD2sample: R Package for Two-Sample Testing

• MD2sample is an R package for nonparametric multivariate two-sample testing designed to compare distributions in continuous and discrete data.
• It integrates complementary methods—including EDF, distance-based, graph-based, and kernel tests—to robustly detect differences in dependence and marginal shifts.
• The package supports extensive simulation studies, enabling systematic power evaluations and the incorporation of user-defined tests via efficient computational techniques.

Searching arXiv for the supplied MD2sample-related papers and closely related records. MD2sample is an R package for nonparametric multivariate two-sample testing, designed to run and compare a large collection of tests for continuous and discrete data and to support systematic power studies across diverse alternatives (Rolke, 12 May 2026). In the two-sample framework, one observes independent samples $X_1,\dots,X_n \sim F$ and $Y_1,\dots,Y_m \sim G$, with the null hypothesis $H_0: F = G$ against the general alternative $H_1: F \neq G$ (Rolke, 22 Jul 2025). The package is positioned around a central empirical conclusion: for multivariate data, no single method can be relied upon to provide good power uniformly across all alternatives, so MD2sample emphasizes breadth of methodology, comparative evaluation, and short lists of complementary tests (Rolke, 12 May 2026).

1. Statistical setting and motivation

MD2sample addresses the multivariate nonparametric two-sample problem in settings where reliable parametric models for $F$ and $G$ are unavailable or undesirable (Rolke, 22 Jul 2025). The package is intended for continuous data in arbitrary dimension and for discrete data in two dimensions, including binned or histogram representations (Rolke, 12 May 2026). The associated power studies explicitly examine $d=2$ and $d=5$ continuous settings and $2$-dimensional discrete settings (Rolke, 12 May 2026).

The methodological motivation is the heterogeneity of alternatives. Differences between $F$ and $Y_1,\dots,Y_m \sim G$0 may appear in marginals, dependence, covariance structure, local clustering, or mixture composition (Rolke, 22 Jul 2025). The supplied studies distinguish cases with equal marginals, where only dependence changes, from cases with unequal marginals, where mean, variance, or marginal family changes occur (Rolke, 12 May 2026). This distinction matters because univariate marginal tests can completely miss dependence changes, especially in the “D2” and “D5” classes of scenarios described in the studies (Rolke, 22 Jul 2025).

A further motivation is that multivariate methods differ materially in what they detect well. EDF-based procedures may respond strongly to marginal differences; distance- and graph-based methods may be more effective for dependence changes or local structure; chi-square procedures can be strong when binning captures relevant features (Rolke, 12 May 2026). This suggests that MD2sample is less a single test than a comparative framework for selecting and combining tests with complementary sensitivity profiles.

2. Scope of the package and computational design

MD2sample is available on CRAN and is designed to provide a large collection of multivariate two-sample tests in a unified framework (Rolke, 22 Jul 2025). Its stated purposes include two-sample testing for multivariate data, support for power studies and simulation experiments, and the ability to combine multiple tests with multiplicity-corrected $Y_1,\dots,Y_m \sim G$1-values (Rolke, 22 Jul 2025). The package also allows the addition of user-defined tests for comparison in simulation workflows (Rolke, 12 May 2026).

For continuous data, MD2sample implements EDF-based, nearest-neighbor, energy- or distance-based, graph-based, Ball divergence, and binning-based chi-square procedures (Rolke, 22 Jul 2025). For discrete data, it includes discrete analogues of KS, Kuiper, CvM, AD, nearest-neighbor methods, Aslan–Zech, Baringhaus–Franz, and a chi-square test on discrete cell counts (Rolke, 22 Jul 2025). In the broader workflow described in the 2026 paper, MD2sample is also used as the two-sample engine for hybrid goodness-of-fit tests in combination with MDgof, where a Monte Carlo sample from a fitted null model is compared to observed data (Rolke, 12 May 2026).

Implementation emphasizes computational practicality. The package uses Rcpp and parallel programming to handle computationally demanding distance calculations, graph constructions, permutations, and repeated simulations (Rolke, 12 May 2026). For some procedures, large-sample formulas or asymptotic approximations are used instead of full permutation calibration, which is particularly relevant for very large datasets (Rolke, 22 Jul 2025). Permutation tests are used for $Y_1,\dots,Y_m \sim G$2-value computation for essentially all two-sample methods in the 2026 description unless otherwise stated, while the 2025 description identifies asymptotic or large-sample calibration for chi-square, FR, CF1–CF4, Ball divergence, and NN0 (Rolke, 12 May 2026, Rolke, 22 Jul 2025). This difference in emphasis reflects the package’s dual role as both a simulation platform and an applied testing toolkit.

3. Families of tests implemented in MD2sample

MD2sample collects several distinct classes of nonparametric two-sample procedures, each defined by a different representation of discrepancy between $Y_1,\dots,Y_m \sim G$3 and $Y_1,\dots,Y_m \sim G$4 (Rolke, 22 Jul 2025).

EDF-based tests

For empirical distribution functions $Y_1,\dots,Y_m \sim G$5, $Y_1,\dots,Y_m \sim G$6, and pooled $Y_1,\dots,Y_m \sim G$7, MD2sample implements multivariate versions of Kolmogorov–Smirnov, Kuiper, Cramér–von Mises, and Anderson–Darling (Rolke, 22 Jul 2025). In the pooled-support formulation reported in the studies, the statistics are

$Y_1,\dots,Y_m \sim G$8

$Y_1,\dots,Y_m \sim G$9

$H_0: F = G$0

$H_0: F = G$1

The studies report that these EDF-type tests tend to work well for simple differences in distribution, especially when marginals differ, but can be weaker for pure dependence changes when marginals are equal (Rolke, 12 May 2026).

Nearest-neighbor procedures

Nearest-neighbor procedures in MD2sample include NN1, NN5, and a simple nearest-neighbor test denoted NN0 or NN (Rolke, 22 Jul 2025). For NN1 and NN5, each observation’s $H_0: F = G$2 nearest neighbors in the pooled sample are examined, and the statistic is based on the average number of same-sample neighbors for both samples (Rolke, 22 Jul 2025, Rolke, 12 May 2026). Under $H_0: F = G$3, labels are exchangeable; under $H_0: F = G$4, within-sample clustering increases the statistic (Rolke, 12 May 2026). The simple nearest-neighbor test counts points in the first sample whose nearest neighbor is also in the first sample and uses an approximate binomial distribution under $H_0: F = G$5 for large samples (Rolke, 22 Jul 2025).

These tests are described as particularly sensitive to local clustering differences and dependence changes, and NN5 is repeatedly identified as a stable compromise between locality and variance (Rolke, 12 May 2026).

Energy- and distance-based tests

The package implements the Aslan–Zech, Baringhaus–Franz, and Biswas–Ghosh tests (Rolke, 22 Jul 2025). Using Euclidean distance $H_0: F = G$6, the Aslan–Zech statistic is

$H_0: F = G$7

the Baringhaus–Franz statistic is

$H_0: F = G$8

and the Biswas–Ghosh statistic is based on

$H_0: F = G$9

$H_1: F \neq G$0

with

$H_1: F \neq G$1

These procedures compare cross-sample and within-sample distances in different ways (Rolke, 22 Jul 2025). The studies characterize AZ as very strong in $H_1: F \neq G$2D continuous settings and BG as particularly designed to work well in higher dimensions, including cases where other tests lose power (Rolke, 12 May 2026).

Kernel-based testing via MMD

The 2026 study adds a kernel-based procedure, Maximum Mean Discrepancy, to the MD2sample-centered comparison (Rolke, 12 May 2026). For a positive definite kernel $H_1: F \neq G$3,

$H_1: F \neq G$4

with unbiased empirical estimator

$H_1: F \neq G$5

MD2sample is described there as implementing an aggregated MMD using $H_1: F \neq G$6 different bandwidths, usually Gaussian kernels with different $H_1: F \neq G$7, following the “MMD aggregated two-sample test” of Schrab et al. 2023 (Rolke, 12 May 2026). In that study, MMD is identified as the single best general-purpose two-sample test for continuous data (Rolke, 12 May 2026).

Graph-based and ball-based procedures

MD2sample includes the Friedman–Rafsky MST test, Chen–Friedman tests CF1–CF4 via the gTests methodology, and Ball divergence via the Ball package (Rolke, 22 Jul 2025). Friedman–Rafsky constructs the minimum spanning tree on the pooled sample and evaluates how strongly the labels mix across edges; Chen–Friedman procedures use graph counts on MST or $H_1: F \neq G$8-nearest-neighbor-type constructions; Ball divergence compares distributions of distances inside metric balls (Rolke, 22 Jul 2025, Rolke, 12 May 2026). These are omnibus procedures and are reported as competitive in complex shape and dependence changes, with CF2 sometimes strong in certain Dalitz-type or structured alternatives (Rolke, 12 May 2026).

Chi-square and discrete variants

For $H_1: F \neq G$9-dimensional data, MD2sample includes chi-square tests based on equal-size or equal-probability binning, typically $F$0 in the reported studies (Rolke, 22 Jul 2025). For pooled cell counts $F$1, sample totals $F$2, $F$3, and $F$4, the statistic given in the 2026 study is

$F$5

The package merges bins until each contains at least $F$6 observations from the combined sample to satisfy chi-square validity, and the papers note that power is often maximized using a relatively small number of bins, regardless of sample size (Rolke, 22 Jul 2025). Discrete analogues of EDF, nearest-neighbor, and energy tests are also implemented for $F$7D discrete data (Rolke, 12 May 2026).

4. Simulation studies and empirical evaluation

The empirical characterization of MD2sample is grounded in extensive simulation studies (Rolke, 12 May 2026). The 2025 paper reports $F$8 case studies in $F$9D, together with $G$0D discrete and $G$1D continuous scenarios (Rolke, 22 Jul 2025). The 2026 paper broadens the discussion to goodness-of-fit as well as two-sample settings and explicitly states that the studies were carried out using MD2sample and MDgof (Rolke, 12 May 2026).

The simulated alternatives cover multivariate normal and $G$2 families, uniform distributions on rectangles, copulas including Frank, Clayton, Gumbel, Galambos, Joe, Husler–Reiss, and mixtures, as well as Dalitz plot examples motivated by high-energy physics (Rolke, 12 May 2026, Rolke, 22 Jul 2025). The case taxonomy distinguishes equal-marginal dependence changes from unequal-marginal alternatives such as mean shifts, variance changes, or changes in marginal family (Rolke, 22 Jul 2025). Discrete scenarios are produced either intrinsically or by binning continuous scenarios into $G$3 grids (Rolke, 12 May 2026).

The studies target nominal significance level $G$4 (Rolke, 12 May 2026). For the 2026 power comparisons, typical sample sizes are $G$5 or $G$6, with $G$7 used for reported power comparisons, and the alternative parameters are chosen so that the best method reaches approximately $G$8–$G$9 power at $d=2$0 (Rolke, 12 May 2026). Thousands of replications per scenario are used to estimate power reliably, and rejection frequencies are summarized both as mean power across scenarios and as the proportion of scenarios in which a method is within $d=2$1 of the best power for that scenario (Rolke, 12 May 2026, Rolke, 22 Jul 2025).

A central methodological conclusion is repeated across both papers: there is no universally best test (Rolke, 12 May 2026, Rolke, 22 Jul 2025). A method can be excellent for mean shifts and poor for subtle dependence changes, or superb on some copula alternatives and almost powerless on others (Rolke, 22 Jul 2025). MD2sample’s comparative structure is therefore integral to its purpose rather than an accessory feature.

5. Main findings on power and recommended short lists

The reported recommendations depend on data type and dimension.

For continuous $d=2$2-dimensional data, the 2026 study reports mean power $d=2$3 for MMD, $d=2$4 for AZ, $d=2$5 for ES chi-square, $d=2$6 for EP chi-square, $d=2$7 for BF, $d=2$8 for AD, $d=2$9 for NN5, $d=5$0 for CvM, $d=5$1 for Ball, and $d=5$2 for FR, CF1, and CF3 (Rolke, 12 May 2026). The corresponding “percentage close to best” values are $d=5$3 for MMD, $d=5$4 for AZ, $d=5$5 for ES, $d=5$6 for EP, and $d=5$7 for BF (Rolke, 12 May 2026). The same paper therefore characterizes MMD as the single most consistently powerful method in $d=5$8D and AZ as the second-best, while noting that BG remains useful in specific scenarios such as complex Dalitz-type shapes (Rolke, 12 May 2026).

The 2025 study, which predates the inclusion of MMD in the summarized comparison, gives a different ranking for continuous cases: AZ $d=5$9, ES $2$0, EP $2$1, BF $2$2, AD $2$3, NN5 $2$4, CvM $2$5, KS $2$6, and Ball $2$7, with BG at $2$8 mean power but sometimes essential (Rolke, 22 Jul 2025). In that study AZ is within $2$9 of the best power in $F$0 of cases, far above other methods (Rolke, 22 Jul 2025). This suggests that AZ was the most robust performer among the then-compared methods, whereas the 2026 study’s addition of aggregated MMD shifts the overall default recommendation toward MMD for continuous data.

For continuous $F$1-dimensional data, the 2026 study reports mean power $F$2 for MMD, $F$3 for AZ, $F$4 for BF, $F$5 for AD, $F$6 for Ball, $F$7–$F$8 for KS, CvM, K, and NN5, and $F$9 for BG (Rolke, 12 May 2026). Despite BG’s modest mean power, the smallest set of methods such that for every scenario at least one is best is reported as BG + MMD + NN5 (Rolke, 12 May 2026). The interpretation given is that MMD dominates on average, while BG and NN5 are essential for certain pathological or high-dimensional alternatives (Rolke, 12 May 2026).

For equal marginals, where only dependence changes, graph- and distance-based tests such as MMD, AZ, BG, NN5, and FR generally outperform EDF types (Rolke, 12 May 2026). With unequal marginals, EDF tests and binning-based chi-square become more competitive, but MMD and AZ remain very strong (Rolke, 12 May 2026).

For $Y_1,\dots,Y_m \sim G$00-dimensional discrete data, the 2026 study reports mean power $Y_1,\dots,Y_m \sim G$01 for chi-square, $Y_1,\dots,Y_m \sim G$02 for NN, $Y_1,\dots,Y_m \sim G$03 for AD, $Y_1,\dots,Y_m \sim G$04 for Kuiper, $Y_1,\dots,Y_m \sim G$05 for KS and CvM, and $Y_1,\dots,Y_m \sim G$06 and $Y_1,\dots,Y_m \sim G$07 for AZ and BF respectively (Rolke, 12 May 2026). The 2025 study gives closely aligned figures: chi-square $Y_1,\dots,Y_m \sim G$08, NN $Y_1,\dots,Y_m \sim G$09, AD $Y_1,\dots,Y_m \sim G$10, K $Y_1,\dots,Y_m \sim G$11, KS $Y_1,\dots,Y_m \sim G$12, CvM $Y_1,\dots,Y_m \sim G$13, AZ $Y_1,\dots,Y_m \sim G$14, and BF $Y_1,\dots,Y_m \sim G$15 (Rolke, 22 Jul 2025). Both sources therefore identify chi-square as clearly strongest on average for the studied discrete problems (Rolke, 12 May 2026, Rolke, 22 Jul 2025).

The recommendations in the two papers differ in detail. The 2025 study recommends for continuous data the set AZ, chi-square with small number of bins (ES, $Y_1,\dots,Y_m \sim G$16 in $Y_1,\dots,Y_m \sim G$17D), AD, NN5, and BG; for discrete data it recommends chi-square, AZ, CvM, and KS (Rolke, 22 Jul 2025). The 2026 study recommends for two-sample continuous $Y_1,\dots,Y_m \sim G$18 data AZ, BG, and MMD; for continuous $Y_1,\dots,Y_m \sim G$19, BG, MMD, and NN5; and for discrete $Y_1,\dots,Y_m \sim G$20, AZ, AD, Kuiper, and chi-square (Rolke, 12 May 2026). The difference reflects the broader method set and updated comparative evidence in the later paper.

6. Practical use, limitations, and relation to adjacent tools

The studies describe a standard MD2sample workflow. For continuous data, inputs are matrices $Y_1,\dots,Y_m \sim G$21 and $Y_1,\dots,Y_m \sim G$22 of sizes $Y_1,\dots,Y_m \sim G$23 and $Y_1,\dots,Y_m \sim G$24; for discrete $Y_1,\dots,Y_m \sim G$25D data, inputs may be count vectors together with bin centers or cell coordinates (Rolke, 12 May 2026). Tests compute the chosen statistic, generate $Y_1,\dots,Y_m \sim G$26-values by permutation or large-sample calibration as appropriate, and can be embedded in repeated simulation loops for type I error and power estimation (Rolke, 22 Jul 2025, Rolke, 12 May 2026). The package contains built-in case studies and stored power-study results, and it supports test combination with multiplicity correction (Rolke, 22 Jul 2025, Rolke, 12 May 2026).

The reported guidance is pragmatic. If a single default test is needed for continuous data, the 2026 study recommends MMD (Rolke, 12 May 2026). If robustness across alternatives is the goal, the recommendation is to run a small complementary panel such as MMD + BG + AZ in $Y_1,\dots,Y_m \sim G$27D or MMD + BG + NN5 in $Y_1,\dots,Y_m \sim G$28D or higher, optionally combining $Y_1,\dots,Y_m \sim G$29-values (Rolke, 12 May 2026). For discrete $Y_1,\dots,Y_m \sim G$30D data, chi-square is the default, complemented by AD and Kuiper, with AZ as an optional addition (Rolke, 12 May 2026). The 2025 study similarly advises applying a battery of tests rather than relying on one method, especially when local structure, covariance stretch, or dependence changes are plausible (Rolke, 22 Jul 2025).

Several limitations are explicit. The systematic simulations extend only to $Y_1,\dots,Y_m \sim G$31, so performance in much higher dimensions remains unstudied in these papers (Rolke, 12 May 2026). Discrete methods are implemented only for $Y_1,\dots,Y_m \sim G$32D because higher-dimensional binning rapidly becomes impractical (Rolke, 12 May 2026). The discrete conclusions are specific to the chosen $Y_1,\dots,Y_m \sim G$33 equal-size or related binning schemes, and the 2025 paper explicitly cautions that some methods, notably AZ, may have much better power under different binning choices (Rolke, 22 Jul 2025). Sample-size imbalance is not deeply explored, and computational cost remains nontrivial for full distance matrices, MST construction, and repeated permutations on large samples (Rolke, 22 Jul 2025, Rolke, 12 May 2026).

MD2sample is also closely related to MDgof. In the 2026 study, MD2sample serves as the two-sample testing engine for hybrid goodness-of-fit procedures in which an observed sample is compared to a Monte Carlo sample from a fitted null model, with MC sample sizes equal to the data size or five times larger (Rolke, 12 May 2026). This places MD2sample at the intersection of classical two-sample testing and simulation-based goodness-of-fit assessment.

Taken together, the supplied papers define MD2sample as a research-oriented package whose distinguishing feature is not merely the number of implemented tests, but the empirical program accompanying them: a large simulation library, explicit power comparisons, and concrete short-list recommendations that acknowledge the absence of any uniformly dominant multivariate two-sample procedure (Rolke, 12 May 2026, Rolke, 22 Jul 2025).