📊 statistics

Hash-augmented adaptive multilevel splitting Monte Carlo algorithm for accurate estimation of two-sample permutation test p-values

This paper introduces a hash-augmented adaptive multilevel splitting Monte Carlo algorithm, implemented in the Python package `hamstest`, to accurately estimate arbitrarily small p-values for two-sample permutation tests with complex statistics while addressing challenges related to distribution discreteness and ensuring valid confidence intervals.

Original authors: Nikita Golikov, Vladimir Sukhov, Gennady Korotkevich, Alexey Sergushichev

Published 2026-07-15

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Nikita Golikov, Vladimir Sukhov, Gennady Korotkevich, Alexey Sergushichev

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are a detective trying to catch a very rare criminal in a city of millions. You have a list of suspects (your data), and you want to know: "How likely is it that this specific pattern of clues happened just by pure luck?" In the world of statistics, this is called a permutation test. You shuffle the clues around millions of times to see how often a "lucky" pattern appears.

Usually, if the pattern is common, you can just count the lucky shuffles. But what if the pattern is so rare that it only happens once in a trillion tries? That's like looking for a single specific grain of sand on a beach the size of a planet. If you try to find it by randomly picking grains one by one (the old Monte Carlo method), you might spend your whole life picking sand and still never find that one grain. You'd need to pick $10^{10}$ grains just to get a decent guess for a tiny probability like $10^{-10}$ , which is totally impractical.

The Problem: The "Stuck" Elevator

The authors of this paper realized that standard methods hit a wall when dealing with these tiny probabilities, especially because the "sand grains" (the data combinations) aren't all unique. Sometimes, thousands of different shuffles result in the exact same score. It's like an elevator that only stops on floors 1, 10, and 100, but skips 2 through 99. If you are trying to get to floor 99, the elevator just can't stop there because it doesn't exist. This "discreteness" causes the math to get stuck, making it impossible to estimate how rare an event really is.

The Solution: The "Hash" Tag and the Splitting Ladder

The team, led by Nikita Golikov and colleagues, built a new tool called hamstest. Their secret sauce is a clever trick called hash-augmented adaptive multilevel splitting.

Here is how it works, using a fun analogy:

The Ladder (Multilevel Splitting): Instead of trying to jump straight to the top of the mountain (the rare event), they build a ladder. They start at the bottom and ask, "How many people can reach the first rung?" Then, "How many of those people can reach the second rung?" They keep splitting the group into smaller and smaller groups as they climb higher. This turns one impossible jump into a series of easy, manageable steps.
The "Hash" Tag (The Fix for Stuck Elevators): The big problem was that many people were standing on the same rung (the same score), making it impossible to split the group further. To fix this, the authors gave every single person a unique, invisible hash tag (a random number). Even if two people have the exact same score, their hash tags are different. This allows the algorithm to say, "Okay, we can't split by score, but we can split by hash tag." It turns a flat, stuck floor into a smooth, continuous staircase where the algorithm can always find the next step.

What They Found (and What They Didn't)

The authors tested this new method on two classic statistical tests: the Kolmogorov–Smirnov test and the Mann–Whitney U test.

The Results: In their simulations, the new method was incredibly accurate. When they tried to estimate probabilities as tiny as $10^{-243}$ (that's a 1 followed by 243 zeros!), the method's estimates landed right on the true value. They also calculated confidence intervals (a range where the true answer is likely hiding), and in about 95% of their test runs, the true answer was inside that range.
The "Full Resampling" Rule: They tried a few different ways to run the simulation. They found that a method called "full resampling" (where they shuffle all the samples at every step) was the most reliable and robust. They suggest using a specific setting called $\alpha = 1$ as the default because it worked best in their tests.
What They Ruled Out: They explicitly showed that the old way of doing things (just using the score without the hash tag) fails when the data has "large jumps" or many ties. They proved that without the hash tag, the algorithm can get stuck and give wrong answers. They also noted that while their method works great for one-sided tests (looking for a pattern in one direction), the two-sided version of the Kolmogorov–Smirnov test is tricky because the "elevator" can get disconnected at the very top, requiring special handling.

How Fast Is It?

The team measured how long the algorithm took on a modern computer (an Apple M3 Pro). They found that the time it takes depends mostly on how rare the event is. If you are looking for something extremely rare (like a p-value of $10^{-100}$ ), it takes longer because you have to climb more rungs on the ladder. However, for the Mann–Whitney U test, the time didn't depend much on the size of the data set because the math for that specific test is very efficient to update.

The Bottom Line

The authors haven't "solved" every statistical problem in the universe, but they have built a very powerful, flexible tool that works for any custom test statistic a scientist might invent. They packaged this tool into a free Python library called hamstest.

They suggest that for most people, using the full resampling method with $\alpha = 1$ is the best bet. They also point out that while their method is fast, the exact time it takes depends on the specific math of the test you are running. If you are a researcher dealing with tiny probabilities and messy data, this tool suggests a way to get accurate answers without waiting for the heat death of the universe.

In short: They turned a broken, stuck elevator into a smooth, high-speed escalator that can take you to the very top of the statistical mountain, even when the path is full of potholes.

Technical Summary: Hash-Augmented Adaptive Multilevel Splitting for Permutation Test p-Values

Problem Statement
Nonparametric permutation tests, such as the Kolmogorov–Smirnov and Mann–Whitney U tests, are fundamental for statistical hypothesis testing, particularly when comparing two groups without assuming a specific data distribution. While exact p-value calculation is ideal, it is often algorithmically intractable for complex or custom test statistics due to the exponential number of possible permutations ( $\binom{n+m}{n}$ ). Consequently, researchers rely on Monte Carlo sampling. However, standard Monte Carlo methods suffer from poor relative accuracy when estimating very small p-values (e.g., $10^{-10}$ or smaller), requiring an impractical number of samples to achieve statistical power after multiple hypothesis correction.

Furthermore, applying multilevel splitting—a technique designed to estimate rare event probabilities by decomposing them into a product of conditional probabilities—to permutation tests introduces specific challenges. Permutation tests operate in discrete spaces where the test statistic distribution often exhibits large jumps (discreteness). In such cases, many distinct combinations of data may yield identical statistic values. This discreteness can prevent the adaptive selection of intermediate levels required for multilevel splitting, causing the algorithm to fail or produce biased estimates when the conditional probability of moving to the next level cannot be maintained near a target constant (e.g., 0.5).

Methodology
The authors propose a hash-augmented adaptive multilevel splitting Monte Carlo algorithm to address the discreteness issues inherent in permutation tests. The core methodology involves three key components:

Hash-Augmented Ordering: To resolve the issue of multiple combinations sharing the same statistic value, the authors introduce a strict total order on the combinations. Instead of ordering solely by the test statistic $S(A)$ , they define an order $A \succ_{S,h} B$ based on the statistic value first, and then by a hash value $h(A)$ if the statistics are equal.
- The hash function $h(A)$ maps combinations to integers. A practical implementation uses a sum of independent random integers assigned to each element, allowing for $O(1)$ updates during Markov Chain Monte Carlo (MCMC) steps.
- This augmentation effectively "smoothes" the discrete distribution, ensuring that for any level boundary, a valid next boundary can be found such that the conditional probability remains estimable.
Adaptive Multilevel Splitting with MCMC: The algorithm estimates the p-value $Pr(S \ge \gamma)$ by decomposing it into a product of conditional probabilities over a sequence of nested events defined by the hash-augmented order.
- Adaptive Levels: Starting from a trivial level, the algorithm samples combinations and uses their quantiles (e.g., the median) to define the next level boundary.
- MCMC Sampling: To generate samples for the next level, the algorithm employs a Metropolis-style MCMC procedure. It starts with existing samples that satisfy the new boundary condition and perturbs them (swapping elements) to generate new independent samples.
- Resampling Strategies: The paper evaluates "full resampling" (duplicating existing valid samples and perturbing all of them) versus "partial resampling" (only perturbing a subset). Full resampling is found to be more robust against sampling biases introduced by the discrete nature of the problem.
Log-Space Estimation: The algorithm estimates the logarithm of the p-value to handle extreme values without numerical underflow. It utilizes the digamma ( $\psi$ ) and trigamma ( $\psi_1$ ) functions to compute the mean and variance of the log-estimates, derived from the properties of order statistics of uniform distributions. This allows for the construction of valid confidence intervals.

Key Contributions

Algorithmic Innovation: The introduction of a hash-augmented ordering scheme that enables adaptive multilevel splitting to function correctly in discrete permutation spaces, overcoming the "large jump" problem where standard methods fail.
Robust Implementation: A detailed analysis of resampling strategies (full vs. partial) and MCMC stopping criteria, demonstrating that full resampling with a specific scaling parameter ( $\alpha=1$ ) provides the best balance of convergence speed and bias reduction.
Software Release: The development and open-source release of the Python package hamstest, which implements the proposed algorithm. It allows users to define custom test statistics and provides efficient interfaces (including C++ backends) for high-performance p-value estimation.
Validation: Rigorous comparison against exact algorithms for Kolmogorov–Smirnov and Mann–Whitney U tests, demonstrating that the proposed method accurately estimates p-values as small as $10^{-243}$ with valid confidence intervals.

Results
The authors evaluated the algorithm on both Kolmogorov–Smirnov and Mann–Whitney U tests across various sample sizes ( $n+m$ up to 10,000) and target p-values (down to $10^{-100}$ ).

Accuracy: The estimated p-values closely matched the true values (calculated via exact methods where available), with no visible bias. The median of the estimates aligned with the true p-value across all tested parameters.
Confidence Intervals: The constructed 95% confidence intervals (based on $\log \hat{p} \pm 2\hat{\sigma}$ ) successfully captured the true log p-value in the vast majority of independent runs, validating the theoretical variance estimates.
Performance: Running times were shown to depend primarily on the magnitude of the p-value and the sample size. The algorithm demonstrated linear scalability with respect to the number of levels required, making it feasible to estimate extremely small p-values that are computationally prohibitive for standard Monte Carlo.
Resampling Comparison: The "full-median" variant (full resampling with median boundary selection) consistently outperformed partial resampling variants in terms of convergence speed and robustness to bias, particularly for the Kolmogorov–Smirnov test where distribution discreteness is more pronounced.

Significance
The paper claims that this work provides a major improvement over classical Monte Carlo approaches for p-value estimation in permutation tests. By addressing the specific pitfalls of discrete spaces—namely the inability to define intermediate levels due to tied statistic values—the proposed algorithm enables the accurate and efficient estimation of arbitrarily small p-values. This capability is critical for modern high-throughput applications (e.g., genomics, neuroscience) where thousands of tests are performed simultaneously, and precise control of false discovery rates requires reliable estimation of tail probabilities. The authors note that a previous version of this algorithm (without the hash augmentation) was already successful in Gene Set Enrichment Analysis (GSEA); the current work extends this utility to a broader class of two-sample statistics and resolves the instability issues encountered with specific input distributions. The availability of the hamstest package ensures that these advanced statistical capabilities are accessible for user-defined statistics in practical research workflows.