From 32e2c1efabcc083237cc6420b09d7508c898ea4f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Petr=20Balv=C3=ADn?= Date: Sun, 27 Sep 2026 19:26:04 +0200 Subject: [PATCH] fix(stats): keep regression inference alive when squared deviations underflow --- CHANGELOG.md | 4 + stats/regression.go | 189 +++++++++++++++++++++++++++++- stats/regression_edge_pin_test.go | 109 +++++++++++++++++ 3 files changed, 296 insertions(+), 6 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ecb1c0b..f00fc46 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -33,6 +33,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 p-values and `NoncentralTCDF` at extreme t answer the tail the format still holds instead of a silent zero or an error naming a NaN. +- `LinearRegression` and `WeightedLinearRegression` keep their + inference alive when the squared deviations underflow to zero: the + standard errors, t and F statistics read factored sums of squares + instead of reporting infinity beside zero evidence. ## [1.0.0] - 2026-09-03 diff --git a/stats/regression.go b/stats/regression.go index cb7c9ae..f229e2b 100644 --- a/stats/regression.go +++ b/stats/regression.go @@ -177,6 +177,7 @@ func LinearRegression(x, y *core.Array) (*LinearRegressionResult, error) { rss := 0.0 tss := 0.0 uncentred := 0.0 + maxRes, maxDev, maxY := 0.0, 0.0, 0.0 mean := 0.0 if fy != nil { for _, v := range fy[:n] { @@ -212,26 +213,85 @@ func LinearRegression(x, y *core.Array) (*LinearRegressionResult, error) { rss += res * res tss += (yv - mean) * (yv - mean) uncentred += yv * yv + if a := math.Abs(res); a > maxRes { + maxRes = a + } + if a := math.Abs(yv - mean); a > maxDev { + maxDev = a + } + if a := math.Abs(yv); a > maxY { + maxY = a + } } if !hasConstant { // The null model is y = 0, so the uncentred total is what the // model has to beat, and it carries n degrees of freedom. tss = uncentred } + // The factored sums of squares: a response on a scale whose squared + // deviations fall below the subnormal floor reads as a zero sum while + // its deviations are live, and the statistics below would report the + // evidence backwards (an exact fit the t statistics cannot support, + // an F of zero beside them). Each pair keeps the largest deviation as + // the scale and the scaled sum as the unit, so scale²·unit is the + // true sum wherever the plain product underflows; the unit stays 1 + // whenever the plain sum already holds. + rssScale, rssUnit := 1.0, rss + if rss == 0 && maxRes > 0 { + rssScale, rssUnit = maxRes, 0.0 + for _, res := range out.Residuals { + d := res / maxRes + rssUnit += d * d + } + } + tssScale, tssUnit := 1.0, tss + if tss == 0 { + devScale := maxDev + if !hasConstant { + devScale = maxY + } + if devScale > 0 { + tssScale = devScale + tssUnit = 0.0 + for r := range n { + var yv, dev float64 + if fy != nil { + yv = fy[r] + } else { + yv = y.FloatAt(r) + } + if !hasConstant { + dev = yv + } else { + dev = yv - mean + } + d := dev / devScale + tssUnit += d * d + } + } + } dof := n - p out.ResidualVariance = rss / float64(dof) - if tss == 0 { + tssDOF := n - 1 + if !hasConstant { + tssDOF = n + } + if tss == 0 && tssScale == 1 { // A constant response reproduced exactly: R² is 1 by the // perfect-fit convention, not the 1 − 0/0 NaN every consumer // would propagate. The same guard the F statistic below has. out.RSquared = 1 out.AdjustedRSquared = 1 + } else if tss == 0 { + // The total underflowed while the response varies: the ratio of + // the factored forms, the scale factors divided out one at a + // time. Both R² measures round back to 1 here, but the F + // statistic below reads the same factored pieces and does not. + ratio := rssUnit / tssUnit * (rssScale / tssScale) * (rssScale / tssScale) + out.RSquared = 1 - ratio + out.AdjustedRSquared = 1 - ratio*float64(tssDOF)/float64(dof) } else { out.RSquared = 1 - rss/tss - tssDOF := n - 1 - if !hasConstant { - tssDOF = n - } out.AdjustedRSquared = 1 - (rss/float64(dof))/(tss/float64(tssDOF)) } out.DModel = p - 1 @@ -271,6 +331,24 @@ func LinearRegression(x, y *core.Array) (*LinearRegressionResult, error) { } out.PValues[j] = pv case v == 0: + // The residual sum of squares may have underflowed while the + // residuals live: the factored standard error is representable + // where the squared one is not, and the t test then reports + // the evidence it actually holds instead of an unearned + // infinity. + if rssScale != 1 && inv[j][j] > 0 { + se := rssScale * math.Sqrt(rssUnit*inv[j][j]/float64(dof)) + if se > 0 { + out.StandardErrors[j] = se + out.TStatistics[j] = beta[j] / se + pv, err := twoSidedT(out.TStatistics[j], dof) + if err != nil { + return nil, base.Errf("%s: %w", name, err) + } + out.PValues[j] = pv + continue + } + } // An exact fit: the coefficient is infinitely many standard // errors from zero, and the evidence is total. Reporting // t = 0 next to p = 0 would contradict itself. A zero @@ -300,6 +378,19 @@ func LinearRegression(x, y *core.Array) (*LinearRegressionResult, error) { explained = 0 // rounding only, and a negative F is meaningless } out.FStatistic = explained / float64(out.DModel) / out.ResidualVariance + if (math.IsInf(out.FStatistic, 0) || math.IsNaN(out.FStatistic)) && (rssScale != 1 || tssScale != 1) { + // An underflowed sum of squares drove the quotient to Inf or + // 0/0 while the factored pieces live: F from the factored + // forms, every scale factor applied one division at a time so + // no intermediate leaves the representable range before the + // answer does. rssUnit 0 is the exact fit, whose F is + // genuinely infinite. + out.FStatistic = (tssUnit*tssScale/rssScale/rssScale*tssScale - rssUnit) * + float64(dof) / (float64(out.DModel) * rssUnit) + if out.FStatistic < 0 { + out.FStatistic = 0 + } + } if math.IsNaN(out.FStatistic) { // 0/0: a response with no variation at all, reproduced // exactly by the fit. There is no evidence of a model, so @@ -484,6 +575,7 @@ func WeightedLinearRegression(x, y, w *core.Array) (*LinearRegressionResult, err // regression and reports the uncentred conventions whenever the // weights vary, so the four model-level fields are overwritten here. sumW, sumWY, rssW := 0.0, 0.0, 0.0 + maxResW := 0.0 for r := range n { var wr float64 if fw != nil { @@ -500,8 +592,12 @@ func WeightedLinearRegression(x, y, w *core.Array) (*LinearRegressionResult, err sumW += wr sumWY += wr * yv rssW += wr * out.Residuals[r] * out.Residuals[r] + if a := math.Abs(out.Residuals[r]); a > maxResW { + maxResW = a + } } tssW := 0.0 + maxDevW := 0.0 if hasConstant { meanW := sumWY / sumW for r := range n { @@ -518,6 +614,9 @@ func WeightedLinearRegression(x, y, w *core.Array) (*LinearRegressionResult, err } d := yv - meanW tssW += wr * d * d + if a := math.Abs(d); a > maxDevW { + maxDevW = a + } } } else { // Without an intercept the null model is zero, so Σw·y² is the @@ -535,16 +634,85 @@ func WeightedLinearRegression(x, y, w *core.Array) (*LinearRegressionResult, err yv = y.FloatAt(r) } tssW += wr * yv * yv + if a := math.Abs(yv); a > maxDevW { + maxDevW = a + } + } + } + // The weighted sums of squares carry the same factored form the + // unweighted fit keeps: a response scale whose weighted squared + // deviations fall below the subnormal floor reads as a zero sum + // while the deviations live, and the F below would report zero + // evidence beside the t statistics' infinity. + rssScaleW, rssUnitW := 1.0, rssW + if rssW == 0 && maxResW > 0 { + rssScaleW = maxResW + rssUnitW = 0.0 + for r := range n { + var wr float64 + if fw != nil { + wr = fw[r] + } else { + wr = w.FloatAt(r) + } + d := out.Residuals[r] / maxResW + rssUnitW += wr * d * d + } + } + tssScaleW, tssUnitW := 1.0, tssW + if tssW == 0 && maxDevW > 0 { + tssScaleW = maxDevW + tssUnitW = 0.0 + if hasConstant { + meanW := sumWY / sumW + for r := range n { + var wr, yv float64 + if fw != nil { + wr = fw[r] + } else { + wr = w.FloatAt(r) + } + if fy != nil { + yv = fy[r] + } else { + yv = y.FloatAt(r) + } + d := (yv - meanW) / maxDevW + tssUnitW += wr * d * d + } + } else { + for r := range n { + var wr, yv float64 + if fw != nil { + wr = fw[r] + } else { + wr = w.FloatAt(r) + } + if fy != nil { + yv = fy[r] + } else { + yv = y.FloatAt(r) + } + d := yv / maxDevW + tssUnitW += wr * d * d + } } } tssDOF := n - 1 if !hasConstant { tssDOF = n } - if tssW == 0 { + if tssW == 0 && tssScaleW == 1 { // Constant weighted response, exact fit: 1, as above. out.RSquared = 1 out.AdjustedRSquared = 1 + } else if tssW == 0 { + // The weighted total underflowed while the weighted response + // varies: the factored ratio, both R² measures rounding back + // to 1 while the F below reads the same pieces and does not. + ratio := rssUnitW / tssUnitW * (rssScaleW / tssScaleW) * (rssScaleW / tssScaleW) + out.RSquared = 1 - ratio + out.AdjustedRSquared = 1 - ratio*float64(tssDOF)/float64(out.DResidual) } else { out.RSquared = 1 - rssW/tssW out.AdjustedRSquared = 1 - (rssW/float64(out.DResidual))/(tssW/float64(tssDOF)) @@ -559,6 +727,15 @@ func WeightedLinearRegression(x, y, w *core.Array) (*LinearRegressionResult, err explained = 0 // rounding only, and a negative F is meaningless } out.FStatistic = explained / float64(out.DModel) / out.ResidualVariance + if (math.IsInf(out.FStatistic, 0) || math.IsNaN(out.FStatistic)) && (rssScaleW != 1 || tssScaleW != 1) { + // The factored F, as in the unweighted path: every scale + // factor divided out one step at a time. + out.FStatistic = (tssUnitW*tssScaleW/rssScaleW/rssScaleW*tssScaleW - rssUnitW) * + float64(out.DResidual) / (float64(out.DModel) * rssUnitW) + if out.FStatistic < 0 { + out.FStatistic = 0 + } + } if math.IsNaN(out.FStatistic) { // 0/0, as in the unweighted path: nothing to test, p = 1. out.FStatistic = 0 diff --git a/stats/regression_edge_pin_test.go b/stats/regression_edge_pin_test.go index 187d42e..d01f1f4 100644 --- a/stats/regression_edge_pin_test.go +++ b/stats/regression_edge_pin_test.go @@ -239,3 +239,112 @@ func mustMatrix(t *testing.T, vals []float64, r, c int) *core.Array { } return a } + +// TestLinearRegressionTinyScaleInference pins the inference of a response +// on a scale whose squared residuals fall below the subnormal floor: the +// plain residual and total sums of squares read zero there, and the fit +// used to report the evidence backwards, R² of 1 with an infinite t and +// p = 0 beside an F of zero with p = 1. The response y = [0, 0, e] over +// x = 1, 2, 3 keeps every least-squares quantity exactly representable +// while both sums of squares underflow: slope e/2, residual sum e²/6, +// total 2e²/3, (XᵀX)⁻¹₁₁ = 1/2, so SE(slope) = e/(2√3), t = √3, +// p = 1/3, F = 3 and R² = 3/4, all closed fractions. +func TestLinearRegressionTinyScaleInference(t *testing.T) { + const e = 1e-200 + x := mustMatrix(t, []float64{1, 1, 1, 2, 1, 3}, 3, 2) + y := mustFloats(t, []float64{0, 0, e}, 3) + res, err := LinearRegression(x, y) + if err != nil { + t.Fatalf("LinearRegression: %v", err) + } + if math.Abs(res.Coefficients[1]-e/2) > 1e-12*e/2 { + t.Fatalf("slope = %.17g, want %.17g", res.Coefficients[1], e/2) + } + wantSE := e / (2 * math.Sqrt(3)) + se := res.StandardErrors[1] + if !(se > 0) || math.IsInf(se, 0) { + t.Fatalf("slope standard error = %g beside nonzero residuals, want %.17g", se, wantSE) + } + if math.Abs(se-wantSE) > 1e-12*wantSE { + t.Fatalf("slope standard error = %.17g, want %.17g", se, wantSE) + } + if math.Abs(res.TStatistics[1]-math.Sqrt(3)) > 1e-12 { + t.Fatalf("t = %.17g, want √3", res.TStatistics[1]) + } + if math.Abs(res.PValues[1]-1.0/3) > 1e-12 { + t.Fatalf("p = %.17g, want 1/3", res.PValues[1]) + } + if math.Abs(res.RSquared-0.75) > 1e-12 { + t.Fatalf("R² = %.17g, want 3/4", res.RSquared) + } + if math.Abs(res.AdjustedRSquared-0.5) > 1e-12 { + t.Fatalf("adjusted R² = %.17g, want 1/2", res.AdjustedRSquared) + } + if math.Abs(res.FStatistic-3) > 1e-11 { + t.Fatalf("F = %.17g, want 3", res.FStatistic) + } + if math.Abs(res.FPValue-1.0/3) > 1e-11 { + t.Fatalf("F p-value = %.17g, want 1/3", res.FPValue) + } + + // Unit weights are the same fit, weighted statistics included. + w := mustFloats(t, []float64{1, 1, 1}, 3) + wres, err := WeightedLinearRegression(x, y, w) + if err != nil { + t.Fatalf("WeightedLinearRegression: %v", err) + } + if math.Abs(wres.TStatistics[1]-math.Sqrt(3)) > 1e-12 { + t.Fatalf("weighted t = %.17g, want √3", wres.TStatistics[1]) + } + if math.Abs(wres.RSquared-0.75) > 1e-12 { + t.Fatalf("weighted R² = %.17g, want 3/4", wres.RSquared) + } + if math.Abs(wres.FStatistic-3) > 1e-11 { + t.Fatalf("weighted F = %.17g, want 3", wres.FStatistic) + } + if math.Abs(wres.FPValue-1.0/3) > 1e-11 { + t.Fatalf("weighted F p-value = %.17g, want 1/3", wres.FPValue) + } + + // A response with a live scale beside the tiny spread keeps the same + // behaviour: the slope's own rounding leaves residuals near its last + // ulp, and the standard error must stay representable and the t + // finite rather than answer an exact fit the residuals contradict. + const a = 1e-160 + step := math.Nextafter(3*a, math.Inf(1)) - 3*a + y2 := mustFloats(t, []float64{a, 2 * a, 3*a + step}, 3) + res2, err := LinearRegression(x, y2) + if err != nil { + t.Fatalf("LinearRegression: %v", err) + } + if se2 := res2.StandardErrors[1]; !(se2 > 0) || math.IsInf(se2, 0) { + t.Fatalf("slope standard error = %g beside nonzero residuals, want a representable value", se2) + } + if math.IsInf(res2.TStatistics[1], 0) { + t.Fatalf("t = %g beside nonzero residuals, want a finite statistic", res2.TStatistics[1]) + } + if res2.FPValue > 1e-10 || res2.PValues[1] > 1e-10 { + t.Fatalf("p = %g, F p = %g, want the far tail both", res2.PValues[1], res2.FPValue) + } +} + +// TestLinearRegressionTinyScaleExactLine pins the fully underflowed +// corner: an exact line at a scale where both the residual and the total +// sums of squares fall below the subnormal floor. The t statistics +// already answer the exact fit with infinite evidence; the F test must +// agree with them instead of reporting zero evidence. +func TestLinearRegressionTinyScaleExactLine(t *testing.T) { + const a = 1e-300 + x := mustMatrix(t, []float64{1, 1, 1, 2, 1, 3}, 3, 2) + y := mustFloats(t, []float64{a, 2 * a, 3 * a}, 3) + res, err := LinearRegression(x, y) + if err != nil { + t.Fatalf("LinearRegression: %v", err) + } + if res.FPValue > 1e-10 { + t.Fatalf("F p-value = %g on an exact line at a tiny scale, want the far tail beside the infinite t", res.FPValue) + } + if res.PValues[1] > 1e-10 { + t.Fatalf("slope p-value = %g, want the exact-fit report", res.PValues[1]) + } +}