From 05503e2b82320b13bc97b418798d5a48e8753b66 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc-Andr=C3=A9=20Bombeck?= Date: Tue, 20 Oct 2026 09:43:47 +0200 Subject: [PATCH] refactor(stats): one least-squares fit for every trend line Six readers fitted a line by hand, four in the determinant form that loses digits on a large x axis, and the two charts carried the same trend function twice. linearFit / fitFromSums in lib/stats is now the one implementation (centred sums, so a fit from points and one composed from stored accumulators are the same arithmetic), trendLineByTime serves both charts with the same two-decimal output, and a guard keeps hand-rolled slopes out of the tree. --- ...ar-fit-single-implementation-guard.test.ts | 69 ++++++++++ src/components/charts/health-chart.tsx | 41 +----- src/components/charts/mood-chart.tsx | 41 +----- src/lib/analytics/trends.ts | 27 ++-- src/lib/insights/derived/trajectory.ts | 20 +-- src/lib/insights/graded-series.ts | 21 +-- src/lib/rollups/measurement-read.ts | 17 +-- src/lib/stats/__tests__/linear-fit.test.ts | 44 ++++++ src/lib/stats/linear-fit.ts | 129 ++++++++++++++++++ 9 files changed, 270 insertions(+), 139 deletions(-) create mode 100644 src/__tests__/linear-fit-single-implementation-guard.test.ts create mode 100644 src/lib/stats/__tests__/linear-fit.test.ts create mode 100644 src/lib/stats/linear-fit.ts diff --git a/src/__tests__/linear-fit-single-implementation-guard.test.ts b/src/__tests__/linear-fit-single-implementation-guard.test.ts new file mode 100644 index 000000000..b1ee601f5 --- /dev/null +++ b/src/__tests__/linear-fit-single-implementation-guard.test.ts @@ -0,0 +1,69 @@ +/** + * Structural guard: a least-squares line is fitted in `lib/stats/linear-fit`. + * + * Six readers fitted one by hand, four with the determinant form + * `n·Σxy − Σx·Σy` and two with centred sums, and two charts carried the same + * function twice. The matcher flags either form of the slope outside the + * implementation. Comments are stripped first; the planted cases prove the + * matcher fires. + */ +import { readFileSync } from "node:fs"; +import { join } from "node:path"; +import { describe, expect, it } from "vitest"; + +import { walkSourceFiles } from "./helpers/source-files"; + +const SRC = join(process.cwd(), "src"); + +const ALLOWED: Record = { + "lib/stats/linear-fit.ts": "The implementation.", +}; + +const HAND_ROLLED_SLOPE = + /\bn\s*\*\s*sum_?xy\s*-\s*sum_?x\s*\*\s*sum_?y\b|\bsxy\s*\/\s*sxx\b|\bsumXy\s*-\s*\(\s*sumX\s*\*\s*sumY\s*\)\s*\/\s*n\b/i; + +function stripComments(text: string): string { + return text.replace(/\/\*[\s\S]*?\*\//g, "").replace(/^\s*\/\/.*$/gm, ""); +} + +describe("one least-squares line", () => { + it("the matcher fires on a hand-rolled slope", () => { + expect( + HAND_ROLLED_SLOPE.test( + "const slope = (n * sumXY - sumX * sumY) / denominator;", + ), + ).toBe(true); + expect(HAND_ROLLED_SLOPE.test("const slope = sxy / sxx;")).toBe(true); + expect( + HAND_ROLLED_SLOPE.test("const sxy = sumXy - (sumX * sumY) / n;"), + ).toBe(true); + expect(HAND_ROLLED_SLOPE.test("const fit = linearFit(xs, ys);")).toBe( + false, + ); + }); + + it("no file outside the implementation fits its own", () => { + const files = walkSourceFiles(SRC, { floor: 3000 }).filter( + (file) => + !file.startsWith("generated/") && + !file.includes("__tests__/") && + !/\.test\.tsx?$/.test(file), + ); + const offenders = files.filter( + (file) => + !(file in ALLOWED) && + HAND_ROLLED_SLOPE.test( + stripComments(readFileSync(join(SRC, file), "utf8")), + ), + ); + expect(offenders).toEqual([]); + // The implementation itself matches, so the sweep read real files. + expect( + HAND_ROLLED_SLOPE.test( + stripComments( + readFileSync(join(SRC, "lib/stats/linear-fit.ts"), "utf8"), + ), + ), + ).toBe(true); + }); +}); diff --git a/src/components/charts/health-chart.tsx b/src/components/charts/health-chart.tsx index 640c8db11..96c31223d 100644 --- a/src/components/charts/health-chart.tsx +++ b/src/components/charts/health-chart.tsx @@ -89,6 +89,7 @@ import { DEFAULT_CHART_RANGE, rangeWindowDays, } from "./chart-range-tabs"; +import { trendLineByTime } from "@/lib/stats/linear-fit"; // The range tabs (7 / 30 / 90 / All, a calendar-day window ending now) and // the "All" tab's ten-year fetch span live with the shared @@ -585,44 +586,6 @@ export function resolveAnnotationPositions( return out; } -function buildTrendSeriesByTime(data: Array<{ date: Date; value: number }>): { - points: Array<{ date: Date; value: number }>; - slopePerDay: number; -} | null { - if (data.length < 2) return null; - - const sorted = [...data].sort((a, b) => a.date.getTime() - b.date.getTime()); - const startTime = sorted[0].date.getTime(); - const points = sorted.map((p) => ({ - x: (p.date.getTime() - startTime) / 86400000, - y: p.value, - })); - - const n = points.length; - const sumX = points.reduce((s, p) => s + p.x, 0); - const sumY = points.reduce((s, p) => s + p.y, 0); - const sumXY = points.reduce((s, p) => s + p.x * p.y, 0); - const sumXX = points.reduce((s, p) => s + p.x * p.x, 0); - - const denominator = n * sumXX - sumX * sumX; - if (denominator === 0) return null; - - const slope = (n * sumXY - sumX * sumY) / denominator; - const intercept = (sumY - slope * sumX) / n; - - return { - points: sorted.map((p) => ({ - date: p.date, - value: - Math.round( - (intercept + slope * ((p.date.getTime() - startTime) / 86400000)) * - 100, - ) / 100, - })), - slopePerDay: slope, - }; -} - export function HealthChart({ types, title, @@ -1235,7 +1198,7 @@ export function HealthChart({ })); if (typeData.length >= 2) { - const trend = buildTrendSeriesByTime(typeData); + const trend = trendLineByTime(typeData); if (trend) { for (const point of trend.points) { const pointTimestamp = point.date.getTime(); diff --git a/src/components/charts/mood-chart.tsx b/src/components/charts/mood-chart.tsx index 2b5641cb7..0e3db438a 100644 --- a/src/components/charts/mood-chart.tsx +++ b/src/components/charts/mood-chart.tsx @@ -54,6 +54,7 @@ import { CHART_HEIGHT_PX } from "@/lib/charts/constants"; import { shouldFireDataReady } from "@/lib/charts/data-ready-latch"; import { moodLabelKeyForScore } from "@/lib/mood/labels"; import { apiGet } from "@/lib/api/api-fetch"; +import { trendLineByTime } from "@/lib/stats/linear-fit"; // --- Types --- @@ -197,44 +198,6 @@ function movingAverageByPoints( }); } -function buildTrendLine(data: Array<{ date: Date; value: number }>): { - points: Array<{ date: Date; value: number }>; - slopePerDay: number; -} | null { - if (data.length < 2) return null; - - const sorted = [...data].sort((a, b) => a.date.getTime() - b.date.getTime()); - const startTime = sorted[0].date.getTime(); - const points = sorted.map((p) => ({ - x: (p.date.getTime() - startTime) / 86400000, - y: p.value, - })); - - const n = points.length; - const sumX = points.reduce((s, p) => s + p.x, 0); - const sumY = points.reduce((s, p) => s + p.y, 0); - const sumXY = points.reduce((s, p) => s + p.x * p.y, 0); - const sumXX = points.reduce((s, p) => s + p.x * p.x, 0); - - const denominator = n * sumXX - sumX * sumX; - if (denominator === 0) return null; - - const slope = (n * sumXY - sumX * sumY) / denominator; - const intercept = (sumY - slope * sumX) / n; - - return { - points: sorted.map((p) => ({ - date: p.date, - value: - Math.round( - (intercept + slope * ((p.date.getTime() - startTime) / 86400000)) * - 100, - ) / 100, - })), - slopePerDay: slope, - }; -} - // --- Aggregation helpers (exported for unit tests) --- interface MoodEntryDay { @@ -542,7 +505,7 @@ export function MoodChart({ })); if (scoreData.length >= 2) { - const trend = buildTrendLine(scoreData); + const trend = trendLineByTime(scoreData); if (trend) { for (const point of trend.points) { const pointTimestamp = point.date.getTime(); diff --git a/src/lib/analytics/trends.ts b/src/lib/analytics/trends.ts index c137b62d3..4e3d18327 100644 --- a/src/lib/analytics/trends.ts +++ b/src/lib/analytics/trends.ts @@ -4,6 +4,7 @@ */ import { READING_SPREAD, standardDeviation } from "@/lib/stats/spread"; +import { linearFit } from "@/lib/stats/linear-fit"; export interface DataPoint { date: Date; @@ -48,26 +49,14 @@ export function trendSlope( y: p.value, })); - const n = points.length; - const sumX = points.reduce((s, p) => s + p.x, 0); - const sumY = points.reduce((s, p) => s + p.y, 0); - const sumXY = points.reduce((s, p) => s + p.x * p.y, 0); - const sumXX = points.reduce((s, p) => s + p.x * p.x, 0); - - const denominator = n * sumXX - sumX * sumX; - if (denominator === 0) return null; - - const slope = (n * sumXY - sumX * sumY) / denominator; - const intercept = (sumY - slope * sumX) / n; - - // R² calculation - const meanY = sumY / n; - const ssTotal = points.reduce((s, p) => s + (p.y - meanY) ** 2, 0); - const ssResidual = points.reduce( - (s, p) => s + (p.y - (intercept + slope * p.x)) ** 2, - 0, + const fit = linearFit( + points.map((p) => p.x), + points.map((p) => p.y), ); - const rSquared = ssTotal === 0 ? 0 : 1 - ssResidual / ssTotal; + if (!fit) return null; + const { slope } = fit; + // No variance in y: a flat line explains nothing, as it always read. + const rSquared = fit.r2 ?? 0; const threshold = 0.01; // units per day threshold for "stable" const direction: "up" | "down" | "stable" = diff --git a/src/lib/insights/derived/trajectory.ts b/src/lib/insights/derived/trajectory.ts index 259d3c952..00b60b45e 100644 --- a/src/lib/insights/derived/trajectory.ts +++ b/src/lib/insights/derived/trajectory.ts @@ -53,6 +53,7 @@ import { import { readDayMeanSeries, type BaselineProfile } from "./baseline"; import type { Derived, DerivedProvenanceSource } from "./types"; import { daysBetweenDateKeys, shiftDateKey } from "@/lib/tz/format"; +import { linearFit } from "@/lib/stats/linear-fit"; /** Trailing window the fit reads (days). Matches the baseline engine. */ const DEFAULT_WINDOW_DAYS = 30; @@ -154,22 +155,9 @@ export function fitOls(xs: number[], ys: number[]): OlsFit | null { const n = xs.length; if (n < 3 || ys.length !== n) return null; - const sumX = xs.reduce((s, x) => s + x, 0); - const sumY = ys.reduce((s, y) => s + y, 0); - const meanX = sumX / n; - const meanY = sumY / n; - - let sxx = 0; - let sxy = 0; - for (let i = 0; i < n; i++) { - const dx = xs[i] - meanX; - sxx += dx * dx; - sxy += dx * (ys[i] - meanY); - } - if (sxx === 0) return null; // all x identical → no slope - - const slope = sxy / sxx; - const intercept = meanY - slope * meanX; + const fit = linearFit(xs, ys); + if (!fit) return null; // all x identical → no slope + const { slope, intercept, meanX, meanY, sxx } = fit; let ssRes = 0; let ssTot = 0; diff --git a/src/lib/insights/graded-series.ts b/src/lib/insights/graded-series.ts index 57caf5be3..4f82da485 100644 --- a/src/lib/insights/graded-series.ts +++ b/src/lib/insights/graded-series.ts @@ -60,6 +60,7 @@ import { import { readingsMean } from "@/lib/measurements/day-mean"; import { usesHourlyMeanDay } from "@/lib/measurements/day-statistic"; import { userDayKey } from "@/lib/tz/format"; +import { linearFit } from "@/lib/stats/linear-fit"; const MS_PER_DAY = 24 * 60 * 60 * 1000; @@ -205,21 +206,11 @@ function foldInto(map: Map, key: string, value: number): void { * over the order the readings landed. */ function leastSquaresSlope(values: number[]): number | null { - const n = values.length; - if (n < 2) return null; - let sumX = 0; - let sumY = 0; - let sumXY = 0; - let sumXX = 0; - for (let i = 0; i < n; i++) { - sumX += i; - sumY += values[i]; - sumXY += i * values[i]; - sumXX += i * i; - } - const denom = n * sumXX - sumX * sumX; - if (denom === 0) return null; - return round((n * sumXY - sumX * sumY) / denom, 4); + const fit = linearFit( + values.map((_, i) => i), + values, + ); + return fit ? round(fit.slope, 4) : null; } /** diff --git a/src/lib/rollups/measurement-read.ts b/src/lib/rollups/measurement-read.ts index 82b7c6070..ede087a84 100644 --- a/src/lib/rollups/measurement-read.ts +++ b/src/lib/rollups/measurement-read.ts @@ -38,6 +38,7 @@ import { getSourceLadder, parseSourcePriority, } from "@/lib/validations/source-priority"; +import { fitFromSums } from "@/lib/stats/linear-fit"; /** * v1.11.1 — load a user's source-priority blob for the rollup collapse. @@ -217,18 +218,12 @@ export function composeRegression( if (n < 2) return MISS; - // Mean-centered (corrected-sum) identities. Sxx/Sxy/Syy are the determinant - // terms divided by n — algebraically identical to `n·Σxx − Σx²` etc., but - // they subtract `Σx²/n` from `Σxx` at matched scale, avoiding the - // catastrophic cancellation the determinant form suffers on the ~1e10 - // epoch-day x-axis. See the function header for the full rationale. - const sxx = sumXx - (sumX * sumX) / n; + // The centred fit over the accumulators (`fitFromSums`), the same + // arithmetic every point fit uses. No x variance: no slope, no r². + const fit = fitFromSums({ n, sumX, sumY, sumXx, sumXy, sumYy }); + const slope = fit ? fit.slope : null; + const r2 = fit ? fit.r2 : null; const syy = sumYy - (sumY * sumY) / n; - const sxy = sumXy - (sumX * sumY) / n; - - const slope = sxx === 0 ? null : sxy / sxx; - // REGR_R2 is null when either variance is zero (matches Postgres). - const r2 = sxx === 0 || syy === 0 ? null : (sxy * sxy) / (sxx * syy); // Population variance: Syy/n = Σyy/n − (Σy/n)². Clamp tiny negative values // float rounding can produce when the variance is effectively zero. diff --git a/src/lib/stats/__tests__/linear-fit.test.ts b/src/lib/stats/__tests__/linear-fit.test.ts new file mode 100644 index 000000000..f9702977e --- /dev/null +++ b/src/lib/stats/__tests__/linear-fit.test.ts @@ -0,0 +1,44 @@ +import { describe, expect, it } from "vitest"; + +import { fitFromSums, linearFit, sumsOf, trendLineByTime } from "../linear-fit"; + +describe("linearFit", () => { + it("fits an exact line", () => { + const fit = linearFit([0, 1, 2, 3], [1, 3, 5, 7])!; + expect(fit.slope).toBeCloseTo(2, 12); + expect(fit.intercept).toBeCloseTo(1, 12); + expect(fit.r2).toBeCloseTo(1, 12); + }); + + it("answers null without a slope, and a null r² without y variance", () => { + expect(linearFit([1], [1])).toBeNull(); + expect(linearFit([2, 2, 2], [1, 2, 3])).toBeNull(); + expect(linearFit([0, 1, 2], [5, 5, 5])!.r2).toBeNull(); + }); + + it("keeps its digits on an epoch-day axis", () => { + const xs = [20_000, 20_001, 20_002, 20_003].map((d) => d * 86_400); + const fit = linearFit(xs, [10, 11, 12, 13])!; + expect(fit.slope * 86_400).toBeCloseTo(1, 6); + }); + + it("fits the same from points and from their sums", () => { + const xs = [0, 2, 3, 7, 9]; + const ys = [4, 5, 5, 9, 11]; + expect(fitFromSums(sumsOf(xs, ys))).toEqual(linearFit(xs, ys)); + }); +}); + +describe("trendLineByTime", () => { + it("draws the fitted value at two decimals per point, in time order", () => { + const day = (n: number) => new Date(Date.UTC(2026, 0, 1 + n)); + const line = trendLineByTime([ + { date: day(2), value: 3 }, + { date: day(0), value: 1 }, + { date: day(1), value: 2.5 }, + ])!; + expect(line.points.map((p) => p.date)).toEqual([day(0), day(1), day(2)]); + expect(line.slopePerDay).toBeCloseTo(1, 12); + expect(line.points.map((p) => p.value)).toEqual([1.17, 2.17, 3.17]); + }); +}); diff --git a/src/lib/stats/linear-fit.ts b/src/lib/stats/linear-fit.ts new file mode 100644 index 000000000..2b05f8592 --- /dev/null +++ b/src/lib/stats/linear-fit.ts @@ -0,0 +1,129 @@ +/** + * The least-squares line, in one place. + * + * Six readers fitted a line by hand: two charts, the trend summary, the + * graded series, the trajectory fit and the rollup regression composed from + * stored accumulators. Two of them used the centred sums, four the + * determinant form `n·Σxy − Σx·Σy`, which loses digits on a large x axis + * (epoch days) that the centred form keeps. Every fit now goes through + * {@link fitFromSums}, the centred form over the accumulators, so a fit read + * from points and one composed from stored sums are the same arithmetic. + * + * `src/__tests__/linear-fit-single-implementation-guard.test.ts` keeps the + * hand-rolled forms out of the rest of the tree. + */ + +/** The sums a fit is composed from. */ +export interface FitSums { + n: number; + sumX: number; + sumY: number; + sumXx: number; + sumXy: number; + sumYy: number; +} + +export interface LinearFit { + slope: number; + intercept: number; + /** Null when y has no variance (Postgres REGR_R2 answers null too). */ + r2: number | null; + n: number; + meanX: number; + meanY: number; + /** Centred sums of squares and products. */ + sxx: number; + sxy: number; + syy: number; +} + +/** + * Fit from accumulated sums. Null below two points or when every x is the + * same (no slope). The centred identities `Sxx = Σxx − (Σx)²/n` subtract at + * matched scale, so the cancellation is bounded by the x variance, not the + * x magnitude. + */ +export function fitFromSums(sums: FitSums): LinearFit | null { + const { n, sumX, sumY, sumXx, sumXy, sumYy } = sums; + if (n < 2) return null; + const sxx = sumXx - (sumX * sumX) / n; + const sxy = sumXy - (sumX * sumY) / n; + const syy = sumYy - (sumY * sumY) / n; + if (sxx === 0) return null; + const slope = sxy / sxx; + const meanX = sumX / n; + const meanY = sumY / n; + return { + slope, + intercept: meanY - slope * meanX, + r2: syy === 0 ? null : (sxy * sxy) / (sxx * syy), + n, + meanX, + meanY, + sxx, + sxy, + syy, + }; +} + +/** The sums of paired values. */ +export function sumsOf(xs: readonly number[], ys: readonly number[]): FitSums { + if (xs.length !== ys.length) { + throw new RangeError("linear fit: xs and ys differ in length"); + } + let sumX = 0; + let sumY = 0; + let sumXx = 0; + let sumXy = 0; + let sumYy = 0; + for (let i = 0; i < xs.length; i += 1) { + const x = xs[i]; + const y = ys[i]; + sumX += x; + sumY += y; + sumXx += x * x; + sumXy += x * y; + sumYy += y * y; + } + return { n: xs.length, sumX, sumY, sumXx, sumXy, sumYy }; +} + +/** Fit `ys` against `xs`. */ +export function linearFit( + xs: readonly number[], + ys: readonly number[], +): LinearFit | null { + return fitFromSums(sumsOf(xs, ys)); +} + +/** One day in milliseconds, the x unit of a trend over dated points. */ +const DAY_MS = 86_400_000; + +/** + * The fitted line over dated points, x in days from the first point: the + * chart trend lines draw it, with each fitted value at two decimals. Null + * below two points or when every point shares one instant. + */ +export function trendLineByTime( + data: ReadonlyArray<{ date: Date; value: number }>, +): { + points: Array<{ date: Date; value: number }>; + slopePerDay: number; +} | null { + if (data.length < 2) return null; + const sorted = [...data].sort((a, b) => a.date.getTime() - b.date.getTime()); + const startTime = sorted[0].date.getTime(); + const xOf = (d: Date) => (d.getTime() - startTime) / DAY_MS; + const fit = linearFit( + sorted.map((p) => xOf(p.date)), + sorted.map((p) => p.value), + ); + if (!fit) return null; + return { + points: sorted.map((p) => ({ + date: p.date, + value: Math.round((fit.intercept + fit.slope * xOf(p.date)) * 100) / 100, + })), + slopePerDay: fit.slope, + }; +}