fix(smooth): Standardize x before Loess fit so degree-2 works on Date axes

CHANGELOG.md

@@ -4,6 +4,7 @@
 
 - Analyses (`linear`, `smooth`, `density`, `histogram`, `expectation`, `filled_contours`) now fit in transformed space and back-transform their output when the relevant aesthetic carries a `scale` function set via `scales` (e.g. `scales(Y = (; scale = log10))`), so a fit on log-scaled data is computed in log space. This is distinct from `axis = (; yscale = ...)`, which only transforms the display [#773](https://github.com/MakieOrg/AlgebraOfGraphics.jl/pull/773).
 - `histogram` now accepts a `direction` (`:x` or `:y`) keyword for 1D histograms, which sets the bar orientation and the aesthetic (and thus scale space) the bins are computed in [#773](https://github.com/MakieOrg/AlgebraOfGraphics.jl/pull/773).
+- `smooth` with the default `degree = 2` no longer collapses on `Date` (and other large-magnitude) x axes [#774](https://github.com/MakieOrg/AlgebraOfGraphics.jl/pull/774).
 
 ## v0.12.13 - 2026-06-16
 

src/transformations/smooth.jl

@@ -6,10 +6,19 @@ Base.@kwdef struct SmoothAnalysis
     level::Float64 = 0.95
 end
 
-# Loess loses precision when x lives far from 0 (e.g. `datetime2float(DateTime)` is ~1e11),
-# producing wildly wrong predictions between the endpoints. Recenter x around its mean before
-# fitting; Loess is translation-invariant in x, so this is exact, not approximate.
-_recenter(xn) = (offset = sum(xn) / length(xn); (xn .- offset, offset))
+# Loess loses precision when x lives far from 0 or spans a large magnitude (e.g.
+# `datetime2float(Date)` is ~1e11 and a multi-day range is ~1e9): the degree-2 normal equations
+# then mix terms differing by ~1e18 and the local fit is ill-conditioned, collapsing the prediction.
+# Standardize x (center on the mean, scale by the standard deviation) before fitting. For continuous
+# x this leaves predictions unchanged to machine precision; with few distinct x values Loess's
+# kd-tree becomes mildly scale-sensitive, shifting predictions by a fraction of a percent.
+function _standardize(xn)
+    offset = sum(xn) / length(xn)
+    centered = xn .- offset
+    scale = sqrt(sum(abs2, centered) / length(centered))
+    scale = iszero(scale) ? one(scale) : scale
+    return centered ./ scale, offset, scale
+end
 
 function (l::SmoothAnalysis)(input::ProcessedLayer)
     tx = position_transform(input.axis_transforms, Lines, input.attributes, 2, 1)
@@ -21,10 +30,10 @@ function (l::SmoothAnalysis)(input::ProcessedLayer)
             x, y = p
             xn = to_transformed_numerical(x, tx)
             yn = to_transformed_numerical(y, ty)
-            xn_c, x_offset = _recenter(xn)
-            model = Loess.loess(xn_c, yn; l.span, l.degree)
+            xn_s, x_offset, x_scale = _standardize(xn)
+            model = Loess.loess(xn_s, yn; l.span, l.degree)
             x̂n = collect(range(extrema(xn)..., length = l.npoints))
-            ŷn = Loess.predict(model, x̂n .- x_offset)
+            ŷn = Loess.predict(model, (x̂n .- x_offset) ./ x_scale)
             x̂ = from_transformed_numerical(x̂n, x, tx)
             ŷ = from_transformed_numerical(ŷn, y, ty)
             return (x̂, ŷ), (;)
@@ -37,10 +46,10 @@ function (l::SmoothAnalysis)(input::ProcessedLayer)
             x, y = p
             xn = to_transformed_numerical(x, tx)
             yn = to_transformed_numerical(y, ty)
-            xn_c, x_offset = _recenter(xn)
-            model = Loess.loess(xn_c, yn; l.span, l.degree)
+            xn_s, x_offset, x_scale = _standardize(xn)
+            model = Loess.loess(xn_s, yn; l.span, l.degree)
             x̂n = collect(range(extrema(xn)..., length = l.npoints))
-            pred = Loess.predict(model, x̂n .- x_offset; interval = l.interval, level = l.level)
+            pred = Loess.predict(model, (x̂n .- x_offset) ./ x_scale; interval = l.interval, level = l.level)
             ŷ = from_transformed_numerical(pred.predictions, y, ty)
             lower = from_transformed_numerical(pred.lower, y, ty)
             upper = from_transformed_numerical(pred.upper, y, ty)
@@ -74,8 +83,8 @@ Smaller values result in smaller local context in fitting.
 `degree` is the polynomial degree used in the loess model.
 Use `interval` to specify what type of interval the shaded band should represent,
 for a given coverage `level` (the default `0.95` equates `alpha = 0.05`).
-Valid values of `interval` are `:confidence` (the default), to delimit the uncertainty 
-of the predicted relationship. Use `interval = nothing` to only compute the line fit, 
+Valid values of `interval` are `:confidence` (the default), to delimit the uncertainty
+of the predicted relationship. Use `interval = nothing` to only compute the line fit,
 without any uncertainty estimate.
 `npoints` is the number of points used by Makie to draw the line and shaded band.
 

test/analyses.jl

@@ -609,24 +609,30 @@ end
     insert!(labels, :color, "c")
     @test labels == map(AlgebraOfGraphics.to_label, processedlayer.labels)
 
-    # Also test integer input
+    # Few distinct integer x values make Loess's kd-tree scale-sensitive, so the reference must
+    # standardize x exactly as `SmoothAnalysis` does to match.
     df = (x = rand(1:3, 1000), y = rand(1000), c = rand(["a", "b"], 1000))
     npoints = 150
     layer = data(df) * mapping(:x, :y, color = :c) * smooth(; npoints)
     processedlayers = AlgebraOfGraphics.ProcessedLayers(layer)
     processedlayer = processedlayers.layers[2]
 
+    function loess_standardized(x, y)
+        offset = sum(x) / length(x)
+        centered = x .- offset
+        scale = sqrt(sum(abs2, centered) / length(centered))
+        model = Loess.loess(centered ./ scale, y)
+        x̂ = collect(range(extrema(x)...; length = npoints))
+        return x̂, Loess.predict(model, (x̂ .- offset) ./ scale)
+    end
+
     x1 = df.x[df.c .== "a"]
     y1 = df.y[df.c .== "a"]
-    loess1 = Loess.loess(x1, y1)
-    x̂1 = range(extrema(x1)...; length = npoints)
-    ŷ1 = Loess.predict(loess1, x̂1)
+    x̂1, ŷ1 = loess_standardized(x1, y1)
 
     x2 = df.x[df.c .== "b"]
     y2 = df.y[df.c .== "b"]
-    loess2 = Loess.loess(x2, y2)
-    x̂2 = range(extrema(x2)...; length = npoints)
-    ŷ2 = Loess.predict(loess2, x̂2)
+    x̂2, ŷ2 = loess_standardized(x2, y2)
 
     x̂, ŷ = processedlayer.positional
 
@@ -1044,6 +1050,20 @@ end
         @test all(<=(sqrt(20.0)), ŷ)
         @test issorted(ŷ)
     end
+
+    @testset "smooth (Date) degree 2" begin
+        ts = Date(2024, 1, 1):Day(1):Date(2024, 1, 20)
+        df = (; x = collect(ts), y = sqrt.(1.0:20.0))
+        layer = data(df) * mapping(:x, :y) * smooth(; interval = nothing)
+        pl = AlgebraOfGraphics.ProcessedLayer(layer)
+        ŷ = only(pl.positional[2])
+
+        df_num = (; x = collect(1.0:20.0), y = sqrt.(1.0:20.0))
+        pl_num = AlgebraOfGraphics.ProcessedLayer(data(df_num) * mapping(:x, :y) * smooth(; interval = nothing))
+        ŷ_num = only(pl_num.positional[2])
+
+        @test ŷ ≈ ŷ_num
+    end
 end
 
 @testset "datalimits broadcast vs per-dim" begin