Intermediate sumsquares overflow loses finite variance and Welch results in weightstats

Describe the bug

For the two samples in the code field below, DescrStatsW returns infinite variances and CompareMeans.ttest_ind(usevar="unequal") returns (-0.0, NaN, NaN) after a common power-of-two rescaling. The sample variances and Welch t, p, and degrees of freedom all have finite binary64 approximations. The first range loss occurs when summing finite squared deviations, before division by the sample size.

RuntimeWarnings are emitted. This is an extreme synthetic example, with 4 observations per group; its frequency in practical data is unassessed.

Let s = 2**511. The group means are 0 and s, and both groups have deviations ±s. The input values, means, deviations, and individual squares s**2 = 2**1022 are exactly representable. Each unnormalized sum is 4*s**2 = 2**1024, which exceeds binary64's largest finite value. However, the population variance is s**2, the unbiased sample variance is (4/3)*s**2 = 2**1024/3 (approximately 5.992310449541053e307), and the squared standard error of the difference is (2/3)*s**2 (approximately 2.9961552247705263e307). These are within the finite range. The rational sample variance is not exactly representable, but has a finite rounded binary64 value.

At the tested commit, DescrStatsW.sumsquares evaluates np.dot((self.demeaned ** 2).T, self.weights) at weightstats.py:152. The individual squares are finite; their unit-weighted reduction first becomes infinity. var_ddof, var, and _var then divide this infinity (lines 169, 190, and 199). The calculation is already centered. CompareMeans.std_meandiff_separatevar consumes _var at line 1032, and dof_satt receives infinities and forms inf/inf at lines 1059–1060. _tstat_generic then produces -0.0 and evaluates the tail with NaN df.

I am not expecting sumsquares itself to return an impossible finite value. The exact dot-product sum exceeds binary64 range; the question is whether the variance calculation needs to materialize it before calculating otherwise representable variance and Welch results.

A bounded search on 2026-09-14 of statsmodels issues and PRs, including open and closed records, did not identify an exact duplicate. #9488 concerns broader weightstats modernization; #1131 discusses scaling more generally. SciPy #26169 loses range when squaring already-finite variance contributions in df calculation and then applies a df=1 fallback. Here the first loss is the earlier statsmodels sum of squared deviations, with NaN df and no such fallback. SciPy #26113 and PR #26135 are related range-stability work in one-sample/paired tests.

The convenience weightstats.ttest_ind wrapper shares this path. NumPy var also overflows for these samples. No claim is made that all statsmodels t-tests have this exact cause. The inspected weightstats tests cover ordinary data and ddof behavior; I did not find a regression test for this extreme scale. Would improving the scaling robustness of the variance and Welch path be in scope? This is a numerical-stability question, not a request for a special case or a prescribed implementation. No practical incidence or severity estimate is available.

Code Sample, a copy-pastable example

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])

Expected Output

For either scale, testing a zero difference in means with a two-sided alternative should give finite results:

t  = -sqrt(3/2) = -1.224744871391589049098642037...
df = 6
p  = 0.26656970338006897957779103665614139...

Exact rational moments were calculated from the actual binary64 inputs. Independent arbitrary-precision evaluations of the Student-t tail using the regularized incomplete beta function, a finite polynomial integral for df=6, and direct density quadrature agree. None of these reference paths calls statsmodels.

Actual Output

Development source at `main` commit `77ecdfa3e93ea16ec2179e696753af885759fb01`:


2.5.3 1.16.3 0.15.1.dev41+g77ecdfa3e
0 (np.float64(-1.224744871391589), np.float64(0.266569703380069), np.float64(6.0)) 4.0 4.0
[]
511 (np.float64(-0.0), np.float64(nan), np.float64(nan)) inf inf
['overflow encountered in dot', 'overflow encountered in dot', 'invalid value encountered in scalar divide', 'invalid value encountered in scalar divide']


No exception or traceback is produced. The installed release wheel, statsmodels 0.15.0, produced the same values and warnings in the same environment.

Versions

The release-wheel `statsmodels.api.show_versions()` output follows, with local installation paths omitted and whitespace condensed:


INSTALLED VERSIONS
------------------
Python: 3.12.10.final.0
OS: Windows 11 10.0.26200 AMD64
byteorder: little
LC_ALL: C.UTF-8
LANG: C.UTF-8

statsmodels
===========
Installed: 0.15.0

Required Dependencies
=====================
cython: 3.3.0
numpy: 2.5.3
scipy: 1.16.3
pandas: 3.0.5
dateutil: 2.9.0.post0
patsy: 1.0.3

Optional Dependencies
=====================
matplotlib: Not installed
cvxopt: Not installed
joblib: Not installed

Developer Tools
===============
IPython: Not installed
jinja2: Not installed
sphinx: Not installed
pygments: 2.21.0
pytest: 9.1.1
virtualenv: Not installed


The development-head run used unchanged Python modules directly from the checkout through `PYTHONPATH`, with only the normal generated version file. Sixteen loaded tracked statsmodels source modules were checked against commit `77ecdfa3e93ea16ec2179e696753af885759fb01`, allowing only checkout line-ending conversion. This was not a full compiled installation; the target import path uses no statsmodels compiled extension. A full Windows build was unavailable without MSVC. The release wheel above provides the full installed-package reproduction and `show_versions()` output.

Checklist

• I have searched the issue tracker for a similar issue and did not find one.

• I have confirmed this bug exists on the latest released version of statsmodels, or on the main branch.