[Paper Review] Numerically more stable computation of the p-values for the two-sample Kolmogorov-Smirnov test

Thomas Viehmann|arXiv (Cornell University)|Feb 16, 2021

Statistical Mechanics and Entropy1 references4 citations

TL;DR

This paper proposes a numerically stable method for computing exact p-values in the two-sample Kolmogorov-Smirnov test by reformulating the standard recurrence relation using a complementary probability $ C_{i,j} $, which avoids catastrophic cancellation in floating-point arithmetic. The new recurrence computes $ C_{i,j} $ as a weighted average of neighboring values, ensuring stable error propagation and enabling accurate p-value estimation even for large sample sizes.

ABSTRACT

The two-sample Kolmogorov-Smirnov test is a widely used statistical test for detecting whether two samples are likely to come from the same distribution. Implementations typically recur on an article of Hodges from 1957. The advances in computation speed make it feasible to compute exact p-values for a much larger range of problem sizes, but these run into numerical stability problems from floating point operations. We provide a simple transformation of the defining recurrence for the two-side two-sample KS test that avoids this.

Motivation & Objective

• To address numerical instability in existing implementations of the two-sample Kolmogorov-Smirnov test p-value computation due to floating-point arithmetic.
• To enable accurate computation of exact p-values for larger sample sizes than previously feasible with stable numerical behavior.
• To replace the standard path-counting recurrence that computes $ 1 - P_2 $, which risks underflow or overflow, with a more robust alternative.
• To provide a stable, numerically reliable algorithm for p-value computation that avoids cancellation errors in floating-point operations.

Proposed method

• The paper introduces a new recurrence for $ C_{i,j} $, defined as $ 1 - A_{i,j}/\binom{i+j}{i} $, representing the proportion of paths that exit the $ d $-corridor.
• The recurrence $ C_{i,j} = C_{i-1,j} \frac{i}{i+j} + C_{i,j-1} \frac{j}{i+j} $ is used for interior points within the corridor, ensuring $ C_{i,j} $ remains a convex combination of neighbors.
• Boundary conditions are set as $ C_{i,j} = 1 $ if $ |i/m - j/n| \geq d $, and $ C_{i,j} = 0 $ if $ i=0 $ or $ j=0 $ and within the corridor.
• The algorithm iterates along diagonals, tracking only the current and previous rows within the corridor to reduce memory usage.
• The implementation uses dynamic programming with a sliding window over the diagonal, enabling efficient memory and computation usage.
• The method is implemented in Python, Numba, and C++, with performance comparisons to SciPy’s path-counting method.

Experimental results

Research questions

• RQ1 Can a numerically stable alternative to the standard path-counting recurrence be developed for exact p-value computation in the two-sample KS test?
• RQ2 Does replacing the computation of $ 1 - P_2 $ with a complementary probability $ C_{i,j} $ improve numerical stability in floating-point arithmetic?
• RQ3 To what extent does the new recurrence prevent cancellation errors and maintain accuracy for large sample sizes?
• RQ4 How does the performance of the new method compare to existing implementations like SciPy’s?

Key findings

• The proposed recurrence for $ C_{i,j} $ avoids the numerical instability of computing $ 1 - P_2 $, which can lead to underflow or overflow in floating-point arithmetic.
• The new method ensures that $ C_{i,j} $ is a weighted average of its neighbors, promoting favorable error propagation and preventing catastrophic cancellation.
• The algorithm maintains high accuracy even for large sample sizes, where traditional methods fail due to floating-point precision limits.
• The implementation is approximately four times slower than SciPy’s optimized path-counting method in naive form, but offers superior numerical stability.
• The method supports efficient memory usage by tracking only a sliding window of values along diagonals, reducing storage requirements.
• The approach is extensible and could be adapted to the outside method, though this remains an open direction for future work.