# Energy and MMD, with blocked Ball Divergence: validation record ## Scope and status This record covers public `equaldist.energy_ksamp` and `equaldist.mmd_2samp`, plus the private Ball Divergence research implementation. Energy and MMD have published statistics, fixed-size label orbits, numerical invariants, and corrected Monte Carlo contracts with executable checks. They test $$ H_0:F_1=\cdots=F_K $$ against general distributional inequality. Unrestricted relabeling requires independent observations and exchangeability of the pooled sample under the null. Equality of means or covariances alone is insufficient. ## Formula ledger ### DISCO energy test For groups $A_j$ of sizes $n_j$, total size $N$, and $0<\alpha<2$, define $$ g_\alpha(A_j,A_k)=\frac{1}{n_jn_k} \sum_{a\in A_j}\sum_{b\in A_k}\lVert a-b\rVert^\alpha. $$ The implementation evaluates the Rizzo--Székely decomposition $$ W_\alpha=\sum_j\frac{n_j}{2}g_\alpha(A_j,A_j), $$ $$ S_\alpha=\sum_{j, Sections 2--3. ### Maximum mean discrepancy For two samples of sizes $n$ and $m$, pySHT reports the unequal-size unbiased estimator $$ \widehat{\mathrm{MMD}}_u^2= \frac{\sum_{i\ne j}k(X_i,X_j)}{n(n-1)}+ \frac{\sum_{i\ne j}k(Y_i,Y_j)}{m(m-1)}- \frac{2\sum_{i,j}k(X_i,Y_j)}{nm}. $$ Supported characteristic kernels are $$ k_{\rm RBF}(x,y)=\exp\{-\lVert x-y\rVert^2/(2\sigma^2)\}, \qquad k_{\rm Lap}(x,y)=\exp\{-\lVert x-y\rVert/\sigma\}. $$ The bandwidth is either a positive scalar in data units or the median of the strictly positive off-diagonal pooled distances. It is selected before labels are permuted and is therefore part of the fixed test statistic. A constant pooled sample is rejected because the median heuristic is undefined. Primary source: Gretton et al., “A kernel two-sample test,” *JMLR* 13 (2012), 723--773, , Section 2 and Lemma 6. ### Ball Divergence (validation-blocked) Let $\delta(u,v,z)=1\{d(u,z)\le d(u,v)\}$. For a ball centered at an $X$ observation with radius determined by another $X$ observation, let $A^X_{ij}$ and $A^Y_{ij}$ be the empirical masses from the two samples; define $C^X_{kl}$ and $C^Y_{kl}$ analogously for balls centered on $Y$. The implemented statistic is $$ \mathrm{BD}_{n,m}=\frac1{n^2}\sum_{i,j=1}^n (A^X_{ij}-A^Y_{ij})^2+ \frac1{m^2}\sum_{k,l=1}^m(C^X_{kl}-C^Y_{kl})^2. $$ The mathematical membership rule is the closed-ball comparison `distance <= radius`. Stable distance orders and tie-end ranks reduce each permuted evaluation from a literal cubic scan to quadratic time, but the prototype is not public. Exact comparison of computed Euclidean distances can split radii that are mathematically equal after a rotation, while a roundoff tolerance can merge radii that are genuinely distinct. On the complete 20-label orbit for `x=[6.66133815e-16,3,-6.66133815e-16]` and `y=[-1,2,2]`, the direct literal statistic is 0.38271604938271603 with 6 exceedances and p-value 0.3. The optimized prototype instead reports 0.4691358024691358 with 4 exceedances and p-value 0.2. No equality rule audited so far passes both this fixture and exact geometric-symmetry fixtures. `ball_divergence_2samp` therefore remains absent from the public API and audit registry until radius equality has a correctness-certified implementation. Primary source: Pan et al., “Ball Divergence: Nonparametric two sample test,” *Annals of Statistics* 46 (2018), 1109--1137, , equations (2.3)--(2.6). ## Randomization and numerical contract For group sizes $(n_1,\ldots,n_K)$, exact calibration enumerates all $N!/\prod_j n_j!$ ordered fixed-size allocations. `calibration="exact"` fails before calculation if that number exceeds `n_resamples`. `"permutation"` enumerates within budget and otherwise samples exactly `n_resamples` allocations. Monte Carlo results use $(b+1)/(B+1)$ and include the conditional MCSE and 95% Clopper--Pearson tail-probability interval; exact results use $b/B$. Coordinates are centered at an overflow-safe per-feature midpoint, divided by a common positive scale, and distances are divided by their maximum. Energy inference is invariant to this scaling. Median-bandwidth MMD is also scale invariant; explicit bandwidths retain their original data units through a log-ratio calculation. A metric-only canonical form incorporates both the pooled distance graph and the unlabeled group partition. Stable distance ranks, colour refinement, and individualization attempt to resolve symmetric geometries within a fixed work budget. Each refinement pass and twin comparison charges $N^2$ times its number of edge layers against a 1,000,000-unit budget. On exhaustion, deterministic coordinate ordering replaces the optional graph search. Exact inference is invariant under row/group reordering, orthogonal maps, feature permutations, and exactly representable common translations. Fixed-seed Monte Carlo replay has the same invariance when canonicalization completes within budget, away from the floating approximate-distance-rank boundary. At that boundary or after the coordinate fallback, isometric coordinates can select a different uniform plan; both plans remain valid draws from the same conditional orbit, but need not give the same finite-$B$ exceedance count. Exact calibration remains invariant. The committed near-regular-octagon boundary fixture applies a rounded rotation to coordinates perturbed at roughly $10^{-15}$. It requires the same statistic, exact exceedance count, and exact p-value, while requiring each forced Monte Carlo result separately to satisfy the corrected-p contract rather than requiring two boundary-dependent seed realizations to coincide. ## Complexity For $N$ pooled observations and $p$ features, distance or Gram construction is $O(N^2p)$ time and $O(N^2)$ storage. Each of $B$ label allocations costs $O(N^2)$ arithmetic, so calibration is $O(N^2p+BN^2)$ time. Dynamically bounded batches keep peak storage $O(N^2)$. Metric canonicalization is fast colour refinement for ordinary data. The bounded refinement/twin work avoids exponential preprocessing on symmetric inputs, including the 64 vertices of a six-dimensional hypercube. Distance-rank construction and edge-token storage remain quadratic, with sorting costs for ranks and neighbor signatures. The budget bounds optional canonicalization work, not distance storage or the requested resampling cost. ## Independent validation and performance `tests/test_equaldist_expanded.py` implements the public statistics independently with SciPy `cdist`, enumerates every labeling of fixed fixtures, and verifies: - two- and three-group DISCO formulas, including constant/separated boundaries; - literal unequal-size $\mathrm{MMD}_u^2$ and fixed bandwidths; - the private Ball prototype's optimized-versus-literal fixtures, including the adversarial radius case that blocks public exposure; - translation, orthogonal, feature, row, group, and common-scaling invariants; - extreme $10^{308}$ and subnormal coordinate scales; - exact budgets, corrected Monte Carlo metadata, deterministic seeds, and isolation from NumPy's global generator. On the release-development machine at $N=50$ and $p=5$, 999 batched permutations took 0.03 seconds for Energy and 0.05 seconds for MMD. Calibration processes at most 512 labelings at once and shrinks the batch as $N$ grows. ## Named-seed release evidence The reproducible runner is `python -m tools.distribution_independence_audits `. Each null row uses 20,000 independently generated null data sets, not repeated draws from one conditional orbit. The null scenario uses two independent $N_2(0,I)$ samples of sizes 4 and 6 and enumerates the complete 210-label orbit for every data set. At levels 0.01, 0.05, and 0.10, the raw counts were: | Method | Seed | at 0.01 | at 0.05 | at 0.10 | |---|---:|---:|---:|---:| | Energy | 2026091101 | 180 / 20,000 | 904 / 20,000 | 1,892 / 20,000 | | MMD | 2026091102 | 175 / 20,000 | 914 / 20,000 | 1,907 / 20,000 | All six calibration points for the two public methods passed the preregistered release tolerance $\max\{0.005,4\sqrt{\alpha(1-\alpha)/20000}\}$. Targeted power uses 300 independent data sets, 199 permutations, and the prespecified alternative $N_3(0,I)$ versus $N_3(1,2.25I)$. These are response checks, not comparative power claims. Counts are recorded by the same runner with `--power`. | Method | Seed | Rejections at 0.01 | at 0.05 | at 0.10 | |---|---:|---:|---:|---:| | Energy | 2026091201 | 202 / 300 | 276 / 300 | 286 / 300 | | MMD | 2026091202 | 193 / 300 | 275 / 300 | 292 / 300 |