statsmodels loses finite Welch results when an intermediate sum overflows

Four observations expose a lost finite result

A compact boundary test makes statsmodels return infinite variances and a Welch result of (−0.0, NaN, NaN), even though the variances, t statistic, degrees of freedom, and p-value all have finite float64 approximations. The first overflow occurs in an intermediate sum of squared deviations.

Tasuku Kobayashi reported the example in statsmodels issue #10252 on September 14, 2026. Licklider reproduced it with the statsmodels 0.15.0 Windows wheel and with the Python source at main commit 77ecdfa3e. On September 29, 2026, Kevin Sheppard merged repair PR #10255 into main and closed the issue as completed. The repair is not yet in a statsmodels release.

The merged repair keeps the finite calculation in range

The merged change rescales centered deviations before squaring them, then restores the scale after reduction. It also normalizes the two variance contributions before combining them for pooled and unequal-variance tests. This avoids materializing the overflowing intermediate sums while preserving the statistical formulas.

The upstream regression tests cover the original four-observation example and a nearby two-observation boundary found during review. For the original example, the merged expectations are t = −√(3/2), df = 6, and p ≈ 0.266569703380069 for both pooled and unequal-variance paths. The change was merged as commit 1767776. The implementation was authored by twelfthlabor and merged by Kevin Sheppard.

Scale the inputs without losing their information

Each group contains four values. The second run multiplies every value by the exact power of two 2**511. This changes the magnitude, but it should not change the Welch t statistic, degrees of freedom, or p-value.

import warnings
import numpy as np
import scipy, statsmodels
from statsmodels.stats.weightstats import DescrStatsW, CompareMeans

print(np.__version__, scipy.__version__, statsmodels.__version__)
for k in (0, 511):
    x = np.ldexp(np.array([-1., 1.] * 2), k)
    y = np.ldexp(np.array([0., 2.] * 2), k)
    with warnings.catch_warnings(record=True) as caught:
        warnings.simplefilter('always')
        a, b = DescrStatsW(x), DescrStatsW(y)
        result = CompareMeans(a, b).ttest_ind(usevar='unequal')
    print(k, result, a.sumsquares, b.sumsquares)
    print([str(w.message) for w in caught])

The large-scale call emitted two overflow encountered in dot warnings and two invalid value encountered in scalar divide warnings. It raised no exception. The recorded environment used Python 3.12.10, NumPy 2.5.3, SciPy 1.16.3, and statsmodels 0.15.1.dev41+g77ecdfa3e on Windows AMD64.

The sum overflows before the variance should

Let s = 2**511. The group means are zero and s, and every deviation is either −s or s. The inputs, means, deviations, and each individual squared deviation s² = 2^1022 are exactly representable and finite.

Adding the four squared deviations produces 4s² = 2^1024. That exact unnormalized sum exceeds the largest finite float64 value. The quantities needed after division are still in range: the population variance is 2^1022, and the unbiased sample variance is 2^1024 / 3, approximately 5.992310449541053 × 10^307. The latter is rounded in float64, but finite.

In the inspected source, DescrStatsW.sumsquares calculates np.dot((self.demeaned ** 2).T, self.weights). Its inputs are finite, but the reduction becomes infinity. Later variance calculations divide that infinity, and the Welch degrees-of-freedom path forms inf / inf. This is the cause analysis in the submitted report, separate from a maintainer's confirmation or choice of repair.

Independent references recover the finite answer

Exact rational moments reconstructed from the actual binary64 inputs give t = −√(3/2) and df = 6. Three arbitrary-precision Student-t tail calculations—regularized incomplete beta, a finite polynomial integral for six degrees of freedom, and direct density quadrature—agree on the two-sided result:

t  = -1.224744871391589049098642037...
df = 6
p  = 0.26656970338006897957779103665614139...

These reference calculations do not call NumPy, SciPy, or statsmodels for the moments or probability. The upstream issue contains the complete submitted reproducer, outputs, source trace, and environment.

What this establishes

The example shows that a finite statistical result can be lost when software first materializes a larger intermediate quantity. Exact power-of-two scaling, small inputs, independent references, and warning capture make the first range loss inspectable. Those are the same engineering methods Licklider uses to build scientific verification infrastructure.

This is a deliberately extreme synthetic example. Its frequency in practical data has not been measured, and the merged repair does not establish that every statsmodels variance or t-test path is protected from range loss. The upstream PR explicitly leaves covariance, correlation, and unrelated statistics functions outside this change. The finding and repair do not add a new method to nomue's supported scope.