diff --git a/crates/functions/src/stat/regression.rs b/crates/functions/src/stat/regression.rs index fe9515860..622022ee9 100644 --- a/crates/functions/src/stat/regression.rs +++ b/crates/functions/src/stat/regression.rs @@ -130,10 +130,6 @@ fn norm2(v: &[f64]) -> f64 { m * v.iter().map(|x| (x / m) * (x / m)).sum::().sqrt() } -/// A column whose values differ by no more than this many ulps of its largest magnitude is -/// taken to be constant: the differences are rounding in how it was computed, not data. -const CONSTANT_ULPS: f64 = 64.0; - /// Result of [`qr_solve`]: the least-squares solution on the kept columns, in the caller's /// column scale, and what the standard errors need. struct QrFit { @@ -276,11 +272,12 @@ fn regress(y: &[f64], x: &[Vec], k: usize, constant: bool) -> R { let yc: Vec = y.iter().map(|v| v - ym).collect(); let sstot: f64 = if constant { yc.iter().map(|v| v * v).sum() } else { y.iter().map(|v| v * v).sum() }; // Only centring leaves rounding residue behind; without an intercept a constant column is data. + // Constant means every raw value is the same double: distinct values, however close, vary. let constant_cols: Vec = (0..k) .map(|j| { - let col = x.iter().map(|r| r.get(j).copied().unwrap_or(0.0)); - let (lo, hi, big) = col.fold((f64::INFINITY, f64::NEG_INFINITY, 0.0f64), |(lo, hi, big), v| (lo.min(v), hi.max(v), big.max(v.abs()))); - constant && hi - lo <= CONSTANT_ULPS * f64::EPSILON * big + let mut col = x.iter().map(|r| r.get(j).copied().unwrap_or(0.0)); + let first = col.next(); + constant && first.is_some_and(|f| col.all(|v| v == f)) }) .collect(); let QrFit { coef, dropped, ssresid, diag, quad } = qr_solve(cols, yc, &xm, &constant_cols); diff --git a/crates/functions/src/stat/tests.rs b/crates/functions/src/stat/tests.rs index 1ce6ae1ce..10fd00e09 100644 --- a/crates/functions/src/stat/tests.rs +++ b/crates/functions/src/stat/tests.rs @@ -582,6 +582,15 @@ fn linest_drops_collinear_columns() { let r = linest_rows(&[0.0, 1.0, 2.0], &[vec![1e13, 1e13 + 1.0, 1e13 + 2.0]], true); rel_close(r[0][0], 1.0, 1e-9, "slope on a predictor with a 1e13 offset"); assert_eq!(r[3][1], Some(1.0), "df: the predictor is kept"); + // Same at 1e15, where the spread is only a few ulps of the offset: still distinct doubles. + let x15 = [1e15, 1e15 + 1.0, 1e15 + 2.0]; + let r = linest_rows(&[0.0, 1.0, 2.0], &[x15.to_vec()], true); + rel_close(r[0][0], 1.0, 1e-9, "slope on a predictor with a 1e15 offset"); + rel_close(r[0][1], -1e15, 1e-9, "intercept on a predictor with a 1e15 offset"); + assert_eq!(r[3][1], Some(1.0), "df: the 1e15 predictor is kept"); + assert!(r[4][1].is_some_and(|v| v.abs() < 1e-9), "ssresid: the 1e15 line is exact, got {:?}", r[4][1]); + let fitted = rows_of(&ev("TREND", vec![c(&[0.0, 1.0, 2.0]), c(&x15)])); + close(fitted[2][0].clone(), 2.0); // More predictors than observations: the surplus columns drop, no panic. let r = linest_rows(&[1.0, 2.0], &[vec![1.0, 2.0], vec![2.0, 1.0], vec![0.0, 5.0]], true); assert_eq!(r[0].len(), 4);