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.

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<string, string> = {
+  "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);
+  });
+});

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();

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();

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" =

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;

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<string, Agg>, 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;
 }
 
 /**

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.

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]);
+  });
+});

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,
+  };
+}