mirror of
https://github.com/storytold/gridcraft.git
synced 2026-11-15 14:52:26 +00:00
LINEST: treat only identical raw values as a constant predictor
The 64-ulp spread threshold still classified distinct representable predictors such as [1e15, 1e15 + 1, 1e15 + 2] as constant and dropped them. A column is now constant only when every raw value is the same double; the 1e15 case is added to the offset regression test. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
1 parent
ed15db4101
commit
86ed3d93b4
2 files changed
+13
-7
No files matched your search
@@ -130,10 +130,6 @@ fn norm2(v: &[f64]) -> f64 {
|
||||
m * v.iter().map(|x| (x / m) * (x / m)).sum::<f64>().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<f64>], k: usize, constant: bool) -> R<Fit> {
|
||||
let yc: Vec<f64> = 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<bool> = (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);
|
||||
|
||||
@@ -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);
|
||||
|
||||
Reference in new issue
Block a user