From 68fb26626bd39a08b59f2ee45a4146dbfbc05205 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Sun, 26 Oct 2025 21:22:02 +0100 Subject: [PATCH 01/18] Center predictors to simplify the calculation of predictions and the derivates --- src/Loess.jl | 31 +++++++++++++------------------ 1 file changed, 13 insertions(+), 18 deletions(-) diff --git a/src/Loess.jl b/src/Loess.jl index 7e65b74..0e3ec80 100644 --- a/src/Loess.jl +++ b/src/Loess.jl @@ -61,7 +61,6 @@ function _loess( # Initialize the regression arrays us = Array{T}(undef, q, 1 + degree * m) - du1dt = zeros(T, m, 1 + degree * m) vs = Array{T}(undef, q) for vert in kdtree.verts @@ -84,12 +83,22 @@ function _loess( dmax = maximum([ds[perm[i]] for i = 1:q]) dmax = iszero(dmax) ? one(dmax) : dmax + # We center the predictors since this will greatly simplify subsequent + # calculations. With centerting, the intercept is the prediction at the + # the vertex and the derivative is linear coefficient. + # + # We still only support m == 1 so we hard code for m=1 below. If support + # for more predictors is introduced then this code would have to be adjusted. + # + # See Cleveland and Grosse page 54 column 1. The results assume centering + # but I'm not sure, the centering is explicitly stated in the paper. + x₁ = xs[perm[1], 1] for i in 1:q pᵢ = perm[i] w = sqrt(tricubic(ds[pᵢ] / dmax)) us[i, 1] = w - for j in 1:m - x = xs[pᵢ, j] + for j in 1:m # we still only support m == 1 + x = xs[pᵢ, j] - x₁ # center wxl = w for l in 1:degree wxl *= x @@ -99,17 +108,6 @@ function _loess( vs[i] = ys[pᵢ] * w end - # Compute the gradient of the vertex - pᵢ = perm[1] - for j in 1:m - x = xs[pᵢ, j] - xl = one(x) - for l in 1:degree - du1dt[j, 1 + (j - 1)*degree + l] = l * xl - xl *= x - end - end - if VERSION < v"1.7.0-DEV.1188" F = qr(us, Val(true)) else @@ -117,10 +115,7 @@ function _loess( end coefs = F\vs - predictions_and_gradients[vert] = [ - us[1, :]' * coefs; # the prediction - du1dt * coefs # the gradient of the prediction - ] + predictions_and_gradients[vert] = coefs[1:2] end LoessModel(convert(Matrix{T}, xs), convert(Vector{T}, ys), predictions_and_gradients, kdtree) From 1ea668ce36d7f649b23714af50c61ea822b3ab00 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Sun, 26 Oct 2025 21:29:54 +0100 Subject: [PATCH 02/18] Compute the hat matrix --- src/Loess.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/Loess.jl b/src/Loess.jl index 0e3ec80..87fae2e 100644 --- a/src/Loess.jl +++ b/src/Loess.jl @@ -40,6 +40,7 @@ function _loess( end n, m = size(xs) + L = zeros(T, n, n) q = max(1, floor(Int, span * n)) # TODO: We need to keep track of how we are normalizing so we can @@ -108,12 +109,12 @@ function _loess( vs[i] = ys[pᵢ] * w end - if VERSION < v"1.7.0-DEV.1188" - F = qr(us, Val(true)) - else - F = qr(us, ColumnNorm()) - end - coefs = F\vs + F = svd(us) + + @show perm + L₁ = F.V[1:1,:] * (Diagonal(F.S) \ (F.U' * Diagonal(us[:, 1]))) + + coefs = F \ vs predictions_and_gradients[vert] = coefs[1:2] end From 84005cad4df65ad509137df666d20489edc01773 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Thu, 30 Oct 2025 21:16:43 +0100 Subject: [PATCH 03/18] Implement confidence intervals for predictions --- Project.toml | 2 + src/Loess.jl | 165 ++++++++++++++++++++++++++++++++++++++++------- test/runtests.jl | 30 +++++---- 3 files changed, 162 insertions(+), 35 deletions(-) diff --git a/Project.toml b/Project.toml index d536a9e..c190dbe 100644 --- a/Project.toml +++ b/Project.toml @@ -4,12 +4,14 @@ version = "0.6.4" [deps] Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0" [compat] Distances = "0.7, 0.8, 0.9, 0.10" +Distributions = "0.25" Statistics = "1" StatsAPI = "1.1" julia = "1.6" diff --git a/src/Loess.jl b/src/Loess.jl index 87fae2e..eecf2cb 100644 --- a/src/Loess.jl +++ b/src/Loess.jl @@ -4,6 +4,7 @@ import Distances: euclidean import StatsAPI: fitted, modelmatrix, predict, residuals, response using Statistics, LinearAlgebra +using Distributions: Normal, TDist, quantile export loess, fitted, modelmatrix, predict, residuals, response @@ -11,9 +12,18 @@ include("kd.jl") struct LoessModel{T <: AbstractFloat} - xs::Matrix{T} # An n by m predictor matrix containing n observations from m predictors - ys::Vector{T} # A length n response vector - predictions_and_gradients::Dict{Vector{T}, Vector{T}} # kd-tree vertexes mapped to prediction and gradient at each vertex + # An n by m predictor matrix containing n observations from m predictors + xs::Matrix{T} + # A length n response vector + ys::Vector{T} + # kd-tree vertexes mapped to prediction and gradient at each vertex + predictions_and_gradients::Dict{Vector{T}, Vector{T}} + # kd-tree vertexes mapped to generating vectors for prediction (first row) + # and gradient (second row) at each vertex. The values values of the + # predictions_and_gradients field are the values of hatmatrix_generator + # multiplied by ys + hatmatrix_generator::Dict{Vector{T}, Matrix{T}} + # The kd-tree kdtree::KDTree{T} end @@ -40,7 +50,6 @@ function _loess( end n, m = size(xs) - L = zeros(T, n, n) q = max(1, floor(Int, span * n)) # TODO: We need to keep track of how we are normalizing so we can @@ -63,6 +72,9 @@ function _loess( # Initialize the regression arrays us = Array{T}(undef, q, 1 + degree * m) vs = Array{T}(undef, q) + # The hat matrix and the derivative at the vertices of kdtree + F = Dict{Vector{T},Matrix{T}}() + for vert in kdtree.verts # reset perm @@ -109,17 +121,25 @@ function _loess( vs[i] = ys[pᵢ] * w end - F = svd(us) + Fact = svd(us) - @show perm - L₁ = F.V[1:1,:] * (Diagonal(F.S) \ (F.U' * Diagonal(us[:, 1]))) + Ftmp = (Fact \ Diagonal(us[:, 1]))[1:2, :] + FF = zeros(T, 2, n) + FF[:, perm[1:q]] = Ftmp + F[vert] = FF - coefs = F \ vs + coefs = Fact \ vs predictions_and_gradients[vert] = coefs[1:2] end - LoessModel(convert(Matrix{T}, xs), convert(Vector{T}, ys), predictions_and_gradients, kdtree) + LoessModel( + convert(Matrix{T}, xs), + convert(Vector{T}, ys), + predictions_and_gradients, + F, + kdtree, + ) end """ @@ -175,14 +195,14 @@ end # Returns: # A length n' vector of predicted response values. # -function predict(model::LoessModel{T}, z::Number) where T - adjacent_verts = traverse(model.kdtree, (T(z),)) +function predict(model::LoessModel{T}, x::Number) where T + adjacent_verts = traverse(model.kdtree, (T(x),)) @assert(length(adjacent_verts) == 2) v₁, v₂ = adjacent_verts[1][1], adjacent_verts[2][1] - if z == v₁ || z == v₂ - return first(model.predictions_and_gradients[[z]]) + if x == v₁ || x == v₂ + return first(model.predictions_and_gradients[[x]]) end y₁, dy₁ = model.predictions_and_gradients[[v₁]] @@ -190,26 +210,86 @@ function predict(model::LoessModel{T}, z::Number) where T b_int = cubic_interpolation(v₁, y₁, dy₁, v₂, y₂, dy₂) - return evalpoly(z, b_int) + return evalpoly(x, b_int) +end + +struct LoessPrediction{T} + predictions::Vector{T} + lower::Vector{T} + upper::Vector{T} + sₓ::Vector{T} + δ₁::T + δ₂::T + s::T + ρ::T end -function predict(model::LoessModel, zs::AbstractVector) +""" + predict( + model::LoessModel, + x::AbstractVecOrMat = modelmatrix(model); + interval::Union{Symbol, Nothing}=nothing, + level::Real=0.95 + ) + +Compute predictions from the fitted Loess `model` at the values in `x`. +By default, only the predictions are returned. If `interval` is set to +:confidence` then a `LoessPrediction` struct is returned which contains +upper and lower confidence bounds at `level` as well the the quantities +used for computing the confidence bounds. + +When the prodictor takes values at the edges of the KD tree, the predictions +are already computed and the function is simply a look-up. For other values +of the preditor, a cubic spline approximation is used for the prediction. + +When confidence bounds are requested, a matrix of size `n \times n`` is +constructed where ``n`` is the length of the fitted data vector. Hence, +this caluculation is only feasible when ``n`` is not too large. For details +on the calculations, see the cited reference. + +Reference: Cleveland, William S., and Eric Grosse. "Computational methods +for local regression." Statistics and computing 1, no. 1 (1991): 47-62. +""" +function predict( + model::LoessModel, + x::AbstractVecOrMat = modelmatrix(model); + interval::Union{Symbol, Nothing}=nothing, + level::Real=0.95 +) if size(model.xs, 2) > 1 throw(ArgumentError("multivariate blending not yet implemented")) end - return [predict(model, z) for z in zs] -end + if interval !== nothing && interval !== :confidence + if interval === :prediction + throw(ArgumentError("predictions not implemented. If you know how to compute them then please file an issue and explain how.")) + end + throw(ArgumentError("interval must be either :prediction or :confidence but was $interval")) + end -function predict(model::LoessModel, zs::AbstractMatrix) - if size(model.xs, 2) != size(zs, 2) - throw(DimensionMismatch("number of columns in input matrix must match the number of columns in the model matrix")) + if !(0 < level < 1) + throw(ArgumentError("level must be between zero and one but was $level")) end - if size(zs, 2) == 1 - return predict(model, vec(zs)) + predictions = [predict(model, _x) for _x in x] + + if interval === nothing + return predictions else - return [predict(model, row) for row in eachrow(zs)] + # see Cleveland and Grosse 1991 p.50. + L = hatmatrix(model) + L̄ = L - I + L̄L̄ = L̄' * L̄ + δ₁ = tr(L̄L̄) + δ₂ = sum(abs2, L̄L̄) + ε̂ = L̄ * model.ys + s = sqrt(sum(abs2, ε̂) / δ₁) + ρ = δ₁^2 / δ₂ + qt = quantile(TDist(ρ), (1 + level) / 2) + sₓ = [s*sqrt(sum(abs2, _hatmatrix_x(model, _x))) for _x in x] + lower = [_x - qt * _sₓ for (_x, _sₓ) in zip(predictions, sₓ)] + upper = [_x + qt * _sₓ for (_x, _sₓ) in zip(predictions, sₓ)] + return LoessPrediction(predictions, lower, upper, sₓ, δ₁, δ₂, s, ρ) end end @@ -217,6 +297,45 @@ fitted(model::LoessModel) = predict(model, modelmatrix(model)) residuals(model::LoessModel) = response(model) .- fitted(model) +function hatmatrix(model::LoessModel{T}) where T + n = length(model.ys) + L = zeros(T, n, n) + for (i, r) in enumerate(eachrow(model.xs)) + z = only(r) # we still only support one predictor + L[i, :] = _hatmatrix_x(model, z) + end + return L +end + +# Compute elements of the hat matrix of `model` at `x`. +# See Cleveland and Grosse 1991 p. 54. +function _hatmatrix_x(model::LoessModel{T}, x::Number) where T + n = length(model.ys) + + adjacent_verts = traverse(model.kdtree, (T(x),)) + + @assert(length(adjacent_verts) == 2) + v₁, v₂ = adjacent_verts[1], adjacent_verts[2] + + if x == v₁ || x == v₂ + Lx = model.hatmatrix_generator[[x]][1, :] + else + Lx = zeros(T, n) + for j in 1:n + b_int = cubic_interpolation( + v₁, + model.hatmatrix_generator[[v₁]][1, j], + model.hatmatrix_generator[[v₁]][2, j], + v₂, + model.hatmatrix_generator[[v₂]][1, j], + model.hatmatrix_generator[[v₂]][2, j], + ) + Lx[j] = evalpoly(x, b_int) + end + end + return Lx +end + """ tricubic(u) diff --git a/test/runtests.jl b/test/runtests.jl index 9bcfd9c..57077ef 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -58,7 +58,6 @@ end # not great tolerance but loess also struggles to capture the sine peaks @test abs(predict(model, i * π + π / 2 )) ≈ 0.9 atol=0.1 end - end @test_throws DimensionMismatch loess([1.0 2.0; 3.0 4.0], [1.0]) @@ -93,18 +92,25 @@ end @testset "predict" begin # In R this is `predict(cars.lo, data.frame(speed = seq(5, 25, 1)))`. - Rvals = [7.797353, 10.002308, 12.499786, 15.281082, 18.446568, 21.865315, 25.517015, 29.350386, 33.230660, 37.167935, 41.205226, 45.055736, 48.355889, 49.824812, 51.986702, 56.461318, 61.959729, 68.569313, 76.316068, 85.212121, 95.324047] - - for (x, Ry) in zip(5:25, Rvals) - if 8 <= x <= 22 - @test predict(ft, x) ≈ Ry rtol = 1e-7 - else - # The outer vertices are expanded by 0.105 in the original implementation. Not sure if we - # want to do the same thing so meanwhile the results will deviate slightly between the - # outermost vertices - @test predict(ft, x) ≈ Ry rtol = 1e-3 - end + R_fit = [7.797353, 10.002308, 12.499786, 15.281082, 18.446568, 21.865315, 25.517015, 29.350386, 33.230660, 37.167935, 41.205226, 45.055736, 48.355889, 49.824812, 51.986702, 56.461318, 61.959729, 68.569313, 76.316068, 85.212121, 95.324047] + R_se_fit = [7.568120, 5.945831, 4.990827, 4.545284, 4.308639, 4.115049, 3.789542, 3.716231, 3.776947, 4.091747, 4.709568, 4.245427, 4.035929, 3.753410, 4.004705, 4.043190, 4.026105, 4.074664, 4.570818, 5.954217, 8.302014] + preds = predict(ft, 5:25; interval = :confidence) + # The calcalations of the δs in code of Cleveland and Grosse are + # based on a statistical approximation, see section 4 of their 1991 + # paper. We use the exact but expensive definition so the values + # are slightly different + @test preds.s ≈ 15.29496 rtol = 1e-3 + @test preds.ρ ≈ 44.6179 rtol = 5e-3 + + for (i, Ry) in enumerate(R_fit) + # The outer vertices are expanded by 0.105 in the original implementation. Not sure if we + # want to do the same thing so meanwhile the results will deviate slightly between the + # outermost vertices + @test preds.predictions[i] ≈ Ry rtol = (4 <= i <= 18) ? 1e-7 : 1e-3 + @test preds.sₓ[i] ≈ R_se_fit[i] rtol = 1e-3 end + @test_throws ArgumentError predict(ft, 5:25; interval = :x) + @test_throws ArgumentError predict(ft, 5:25; interval = :prediction) end end From e8840fd18b663edd4eb33d4f66efa2af54a2c48e Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Wed, 5 Nov 2025 20:40:09 +0100 Subject: [PATCH 04/18] Reduce benchmark size to avoid running out of memory --- benchmark/benchmarks.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index 7bd21ce..6601b9e 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -4,7 +4,7 @@ const SUITE = BenchmarkGroup() SUITE["random"] = BenchmarkGroup() -for i in 2:6 +for i in 2:4 n = 10^i x = rand(MersenneTwister(42), n) y = sqrt.(x) From b6ca98d528826670d6a37efcc5866c4bf0fa2c4e Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 11 Nov 2025 07:50:15 +0100 Subject: [PATCH 05/18] Update src/Loess.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: David Müller-Widmann --- src/Loess.jl | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/src/Loess.jl b/src/Loess.jl index eecf2cb..7b862a8 100644 --- a/src/Loess.jl +++ b/src/Loess.jl @@ -320,16 +320,10 @@ function _hatmatrix_x(model::LoessModel{T}, x::Number) where T if x == v₁ || x == v₂ Lx = model.hatmatrix_generator[[x]][1, :] else - Lx = zeros(T, n) + Lv₁ = model.hatmatrix_generator[v₁] + Lv₂ = model.hatmatrix_generator[v₂] for j in 1:n - b_int = cubic_interpolation( - v₁, - model.hatmatrix_generator[[v₁]][1, j], - model.hatmatrix_generator[[v₁]][2, j], - v₂, - model.hatmatrix_generator[[v₂]][1, j], - model.hatmatrix_generator[[v₂]][2, j], - ) + b_int = cubic_interpolation(v₁, Lv₁[1, j], Lv₁[2, j], v₂, Lv₂[1, j], Lv₂[2, j]) Lx[j] = evalpoly(x, b_int) end end From ea96be997a70156de6b9d859b141706bc0557b0c Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 11 Nov 2025 07:51:42 +0100 Subject: [PATCH 06/18] Update src/Loess.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: David Müller-Widmann --- src/Loess.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Loess.jl b/src/Loess.jl index 7b862a8..c5eef42 100644 --- a/src/Loess.jl +++ b/src/Loess.jl @@ -97,7 +97,7 @@ function _loess( dmax = iszero(dmax) ? one(dmax) : dmax # We center the predictors since this will greatly simplify subsequent - # calculations. With centerting, the intercept is the prediction at the + # calculations. With centering, the intercept is the prediction at the # the vertex and the derivative is linear coefficient. # # We still only support m == 1 so we hard code for m=1 below. If support From 31654ecf3050e018be489fe10abd0e815b8da3f2 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 11 Nov 2025 07:53:28 +0100 Subject: [PATCH 07/18] Update src/Loess.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: David Müller-Widmann --- src/Loess.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Loess.jl b/src/Loess.jl index c5eef42..ddbbabe 100644 --- a/src/Loess.jl +++ b/src/Loess.jl @@ -4,7 +4,7 @@ import Distances: euclidean import StatsAPI: fitted, modelmatrix, predict, residuals, response using Statistics, LinearAlgebra -using Distributions: Normal, TDist, quantile +using StatsFuns: tdistinvcdf export loess, fitted, modelmatrix, predict, residuals, response From bcaaaac31e4d8037e800862ab6602454393bb245 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 11 Nov 2025 07:53:38 +0100 Subject: [PATCH 08/18] Update src/Loess.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: David Müller-Widmann --- src/Loess.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Loess.jl b/src/Loess.jl index ddbbabe..d13b514 100644 --- a/src/Loess.jl +++ b/src/Loess.jl @@ -242,7 +242,7 @@ When the prodictor takes values at the edges of the KD tree, the predictions are already computed and the function is simply a look-up. For other values of the preditor, a cubic spline approximation is used for the prediction. -When confidence bounds are requested, a matrix of size `n \times n`` is +When confidence bounds are requested, a matrix of size ``n \times n`` is constructed where ``n`` is the length of the fitted data vector. Hence, this caluculation is only feasible when ``n`` is not too large. For details on the calculations, see the cited reference. From 21eecf654158f649a8f982541e346313cffdf6c0 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 11 Nov 2025 07:53:52 +0100 Subject: [PATCH 09/18] Update src/Loess.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: David Müller-Widmann --- src/Loess.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Loess.jl b/src/Loess.jl index d13b514..c6b43fb 100644 --- a/src/Loess.jl +++ b/src/Loess.jl @@ -264,7 +264,7 @@ function predict( if interval === :prediction throw(ArgumentError("predictions not implemented. If you know how to compute them then please file an issue and explain how.")) end - throw(ArgumentError("interval must be either :prediction or :confidence but was $interval")) + throw(ArgumentError(LazyString("interval must be either :prediction or :confidence but was :", interval))) end if !(0 < level < 1) From a95376e7d58f42162dd9a929c8b0e6cf7943a459 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 11 Nov 2025 07:53:59 +0100 Subject: [PATCH 10/18] Update src/Loess.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: David Müller-Widmann --- src/Loess.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Loess.jl b/src/Loess.jl index c6b43fb..283c14b 100644 --- a/src/Loess.jl +++ b/src/Loess.jl @@ -268,7 +268,7 @@ function predict( end if !(0 < level < 1) - throw(ArgumentError("level must be between zero and one but was $level")) + throw(ArgumentError(LazyString("level must be between zero and one but was ", level))) end predictions = [predict(model, _x) for _x in x] From f97d6906dc0c5eede7376fe8f03dee2b0fd46ecd Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 11 Nov 2025 07:55:20 +0100 Subject: [PATCH 11/18] Update src/Loess.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: David Müller-Widmann --- src/Loess.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Loess.jl b/src/Loess.jl index 283c14b..a7c575b 100644 --- a/src/Loess.jl +++ b/src/Loess.jl @@ -285,7 +285,7 @@ function predict( ε̂ = L̄ * model.ys s = sqrt(sum(abs2, ε̂) / δ₁) ρ = δ₁^2 / δ₂ - qt = quantile(TDist(ρ), (1 + level) / 2) + qt = tdistinvcdf(ρ, (1 + level) / 2) sₓ = [s*sqrt(sum(abs2, _hatmatrix_x(model, _x))) for _x in x] lower = [_x - qt * _sₓ for (_x, _sₓ) in zip(predictions, sₓ)] upper = [_x + qt * _sₓ for (_x, _sₓ) in zip(predictions, sₓ)] From 8c2d537c94f145a53caa4812ba2eb319c88f920b Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 11 Nov 2025 11:12:12 +0100 Subject: [PATCH 12/18] Add StatsFuns and remove Distributions in Project.toml And some other cleanup after merging review comments --- Project.toml | 4 ++-- src/Loess.jl | 5 +++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index c190dbe..f35f735 100644 --- a/Project.toml +++ b/Project.toml @@ -4,16 +4,16 @@ version = "0.6.4" [deps] Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" -Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0" +StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" [compat] Distances = "0.7, 0.8, 0.9, 0.10" -Distributions = "0.25" Statistics = "1" StatsAPI = "1.1" +StatsFuns = "1" julia = "1.6" [extras] diff --git a/src/Loess.jl b/src/Loess.jl index a7c575b..6189d6b 100644 --- a/src/Loess.jl +++ b/src/Loess.jl @@ -320,8 +320,9 @@ function _hatmatrix_x(model::LoessModel{T}, x::Number) where T if x == v₁ || x == v₂ Lx = model.hatmatrix_generator[[x]][1, :] else - Lv₁ = model.hatmatrix_generator[v₁] - Lv₂ = model.hatmatrix_generator[v₂] + Lx = zeros(T, n) + Lv₁ = model.hatmatrix_generator[[v₁]] + Lv₂ = model.hatmatrix_generator[[v₂]] for j in 1:n b_int = cubic_interpolation(v₁, Lv₁[1, j], Lv₁[2, j], v₂, Lv₂[1, j], Lv₂[2, j]) Lx[j] = evalpoly(x, b_int) From d5e7b4f51ce37e2d1287310a27005dbcdb77dc23 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 11 Nov 2025 14:38:32 +0100 Subject: [PATCH 13/18] Handle matrix inputs to predict and thereby also the case without new data --- src/Loess.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/Loess.jl b/src/Loess.jl index 6189d6b..272c912 100644 --- a/src/Loess.jl +++ b/src/Loess.jl @@ -259,6 +259,8 @@ function predict( if size(model.xs, 2) > 1 throw(ArgumentError("multivariate blending not yet implemented")) end + # FIXME! Adjust this when multivariate predictors are supported + x_v = vec(x) if interval !== nothing && interval !== :confidence if interval === :prediction @@ -271,7 +273,7 @@ function predict( throw(ArgumentError(LazyString("level must be between zero and one but was ", level))) end - predictions = [predict(model, _x) for _x in x] + predictions = [predict(model, _x) for _x in x_v] if interval === nothing return predictions @@ -286,7 +288,7 @@ function predict( s = sqrt(sum(abs2, ε̂) / δ₁) ρ = δ₁^2 / δ₂ qt = tdistinvcdf(ρ, (1 + level) / 2) - sₓ = [s*sqrt(sum(abs2, _hatmatrix_x(model, _x))) for _x in x] + sₓ = [s*sqrt(sum(abs2, _hatmatrix_x(model, _x))) for _x in x_v] lower = [_x - qt * _sₓ for (_x, _sₓ) in zip(predictions, sₓ)] upper = [_x + qt * _sₓ for (_x, _sₓ) in zip(predictions, sₓ)] return LoessPrediction(predictions, lower, upper, sₓ, δ₁, δ₂, s, ρ) From 730608ddebd73a99bf995d73cc774b31e2683f3a Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 11 Nov 2025 14:49:06 +0100 Subject: [PATCH 14/18] Some memory optimizations --- src/Loess.jl | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/Loess.jl b/src/Loess.jl index 272c912..d1d868a 100644 --- a/src/Loess.jl +++ b/src/Loess.jl @@ -118,17 +118,23 @@ function _loess( us[i, 1 + (j - 1)*degree + l] = wxl # w*x^l end end - vs[i] = ys[pᵢ] * w + vs[i] = ys[pᵢ] end Fact = svd(us) - Ftmp = (Fact \ Diagonal(us[:, 1]))[1:2, :] + # compute the inverse of the singular values with thresholding similar + # to pinv + Sinv = map!(t -> t > eps(one(t)) * q ? inv(t) : zero(t), Fact.S, Fact.S) + # Compute S⁻¹ * Uᵗ * W in place + Fact.U .*= view(us, :, 1) .* Sinv' + # Finalize the pseudo inverse + Ftmp = Fact.V * Fact.U' FF = zeros(T, 2, n) - FF[:, perm[1:q]] = Ftmp + FF[:, view(perm, 1:q)] = view(Ftmp, 1:2, :) F[vert] = FF - coefs = Fact \ vs + coefs = Ftmp * vs predictions_and_gradients[vert] = coefs[1:2] end From 3c3b32ecaf132862515e77407eca6d572489792d Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 11 Nov 2025 14:55:26 +0100 Subject: [PATCH 15/18] Update src/Loess.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: David Müller-Widmann --- src/Loess.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/Loess.jl b/src/Loess.jl index d1d868a..32db670 100644 --- a/src/Loess.jl +++ b/src/Loess.jl @@ -286,7 +286,10 @@ function predict( else # see Cleveland and Grosse 1991 p.50. L = hatmatrix(model) - L̄ = L - I + for i in diagind(L) + L[i] -= 1 + end + L̄ = L L̄L̄ = L̄' * L̄ δ₁ = tr(L̄L̄) δ₂ = sum(abs2, L̄L̄) From 3a9730a8de1c0770fc68eac73038dcb855a4944e Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 11 Nov 2025 15:22:53 +0100 Subject: [PATCH 16/18] Require Julia 1.8 because of LazyStrings. Test on min and lts --- .github/workflows/ci.yml | 5 +++-- Project.toml | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f3f7328..435cb2c 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -16,8 +16,9 @@ jobs: fail-fast: false matrix: version: - - '1.6' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'. - - '1' # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia. + - 'min' + - 'lts' + - '1' - 'pre' os: - ubuntu-latest diff --git a/Project.toml b/Project.toml index f35f735..9be62d2 100644 --- a/Project.toml +++ b/Project.toml @@ -14,7 +14,7 @@ Distances = "0.7, 0.8, 0.9, 0.10" Statistics = "1" StatsAPI = "1.1" StatsFuns = "1" -julia = "1.6" +julia = "1.8" [extras] RDatasets = "ce6b1742-4840-55fa-b093-852dadbb1d8b" From 9ff1e6efbd587a503c9c5a82b3b7460c2873dd68 Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Tue, 11 Nov 2025 16:15:02 +0100 Subject: [PATCH 17/18] Reuse the hat matrix when predicting within sample with uncertainty Also, use norm instead of sqrt(sum(abs2,...)) --- Project.toml | 1 + src/Loess.jl | 11 ++++++----- test/runtests.jl | 6 ++++++ 3 files changed, 13 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index 9be62d2..02ebd23 100644 --- a/Project.toml +++ b/Project.toml @@ -11,6 +11,7 @@ StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" [compat] Distances = "0.7, 0.8, 0.9, 0.10" +LinearAlgebra = "1" Statistics = "1" StatsAPI = "1.1" StatsFuns = "1" diff --git a/src/Loess.jl b/src/Loess.jl index 32db670..9afc5c5 100644 --- a/src/Loess.jl +++ b/src/Loess.jl @@ -286,10 +286,7 @@ function predict( else # see Cleveland and Grosse 1991 p.50. L = hatmatrix(model) - for i in diagind(L) - L[i] -= 1 - end - L̄ = L + L̄ = L - I L̄L̄ = L̄' * L̄ δ₁ = tr(L̄L̄) δ₂ = sum(abs2, L̄L̄) @@ -297,7 +294,11 @@ function predict( s = sqrt(sum(abs2, ε̂) / δ₁) ρ = δ₁^2 / δ₂ qt = tdistinvcdf(ρ, (1 + level) / 2) - sₓ = [s*sqrt(sum(abs2, _hatmatrix_x(model, _x))) for _x in x_v] + sₓ = if model.xs === x + [s * norm(view(L, i, :)) for i in 1:length(x)] + else + [s * norm(_hatmatrix_x(model, _x)) for _x in x_v] + end lower = [_x - qt * _sₓ for (_x, _sₓ) in zip(predictions, sₓ)] upper = [_x + qt * _sₓ for (_x, _sₓ) in zip(predictions, sₓ)] return LoessPrediction(predictions, lower, upper, sₓ, δ₁, δ₂, s, ρ) diff --git a/test/runtests.jl b/test/runtests.jl index 57077ef..4ac773b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -109,6 +109,12 @@ end @test preds.predictions[i] ≈ Ry rtol = (4 <= i <= 18) ? 1e-7 : 1e-3 @test preds.sₓ[i] ≈ R_se_fit[i] rtol = 1e-3 end + + preds_all = predict(ft; interval = :confidence) + @test preds_all.s ≈ 15.29496 rtol = 1e-3 + @test preds_all.ρ ≈ 44.6179 rtol = 5e-3 + @test preds_all.sₓ ≈ [9.885466, 9.885466, 4.990827, 4.990827, 4.545284, 4.308639, 4.115049, 4.115049, 4.115049, 3.789542, 3.789542, 3.716231, 3.716231, 3.716231, 3.716231, 3.776947, 3.776947, 3.776947, 3.776947, 4.091747, 4.091747, 4.091747, 4.091747, 4.709568, 4.709568, 4.709568, 4.245427, 4.245427, 4.035929, 4.035929, 4.035929, 3.753410, 3.753410, 3.753410, 3.753410, 4.004705, 4.004705, 4.004705, 4.043190, 4.043190, 4.043190, 4.043190, 4.043190, 4.074664, 4.570818, 5.954217, 5.954217, 5.954217, 5.954217, 8.302014] rtol = 1e-2 + @test_throws ArgumentError predict(ft, 5:25; interval = :x) @test_throws ArgumentError predict(ft, 5:25; interval = :prediction) end From 640974f961ea58b6f452957c9c1c70fba53aafee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20M=C3=BCller-Widmann?= Date: Sat, 22 Nov 2025 01:18:38 +0100 Subject: [PATCH 18/18] Reduce allocations --- src/Loess.jl | 90 ++++++++++++++++++++++++++++------------------------ src/kd.jl | 27 ++++++++-------- 2 files changed, 62 insertions(+), 55 deletions(-) diff --git a/src/Loess.jl b/src/Loess.jl index 9afc5c5..9ed22c2 100644 --- a/src/Loess.jl +++ b/src/Loess.jl @@ -93,7 +93,7 @@ function _loess( # find the q closest points partialsort!(perm, 1:q, by=i -> ds[i]) - dmax = maximum([ds[perm[i]] for i = 1:q]) + dmax = maximum(ds[perm[i]] for i = 1:q) dmax = iszero(dmax) ? one(dmax) : dmax # We center the predictors since this will greatly simplify subsequent @@ -129,14 +129,12 @@ function _loess( # Compute S⁻¹ * Uᵗ * W in place Fact.U .*= view(us, :, 1) .* Sinv' # Finalize the pseudo inverse - Ftmp = Fact.V * Fact.U' + Ftmp = view(Fact.V, 1:2, :) * Fact.U' FF = zeros(T, 2, n) - FF[:, view(perm, 1:q)] = view(Ftmp, 1:2, :) + FF[:, view(perm, 1:q)] = Ftmp F[vert] = FF - coefs = Ftmp * vs - - predictions_and_gradients[vert] = coefs[1:2] + predictions_and_gradients[vert] = Ftmp * vs end LoessModel( @@ -202,21 +200,22 @@ end # A length n' vector of predicted response values. # function predict(model::LoessModel{T}, x::Number) where T - adjacent_verts = traverse(model.kdtree, (T(x),)) + v₁, v₂ = traverse(model.kdtree, (T(x),)) + + v₁_val = only(v₁) + v₂_val = only(v₂) + if x == v₁_val + return first(model.predictions_and_gradients[v₁]) + elseif x == v₂_val + return first(model.predictions_and_gradients[v₂]) + else + y₁, dy₁ = model.predictions_and_gradients[v₁] + y₂, dy₂ = model.predictions_and_gradients[v₂] - @assert(length(adjacent_verts) == 2) - v₁, v₂ = adjacent_verts[1][1], adjacent_verts[2][1] + b_int = cubic_interpolation(v₁_val, y₁, dy₁, v₂_val, y₂, dy₂) - if x == v₁ || x == v₂ - return first(model.predictions_and_gradients[[x]]) + return evalpoly(x, b_int) end - - y₁, dy₁ = model.predictions_and_gradients[[v₁]] - y₂, dy₂ = model.predictions_and_gradients[[v₂]] - - b_int = cubic_interpolation(v₁, y₁, dy₁, v₂, y₂, dy₂) - - return evalpoly(x, b_int) end struct LoessPrediction{T} @@ -286,7 +285,16 @@ function predict( else # see Cleveland and Grosse 1991 p.50. L = hatmatrix(model) - L̄ = L - I + if model.xs === x + norm_hatmatrix_x = [norm(Li) for Li in eachrow(L)] + else + tmp = Vector{eltype(L)}(undef, length(model.ys)) + norm_hatmatrix_x = [norm(_hatmatrix_x!(tmp, model, _x)) for _x in x_v] + end + L̄ = L + for i in diagind(L̄) + L̄[i] -= 1 + end L̄L̄ = L̄' * L̄ δ₁ = tr(L̄L̄) δ₂ = sum(abs2, L̄L̄) @@ -294,13 +302,9 @@ function predict( s = sqrt(sum(abs2, ε̂) / δ₁) ρ = δ₁^2 / δ₂ qt = tdistinvcdf(ρ, (1 + level) / 2) - sₓ = if model.xs === x - [s * norm(view(L, i, :)) for i in 1:length(x)] - else - [s * norm(_hatmatrix_x(model, _x)) for _x in x_v] - end - lower = [_x - qt * _sₓ for (_x, _sₓ) in zip(predictions, sₓ)] - upper = [_x + qt * _sₓ for (_x, _sₓ) in zip(predictions, sₓ)] + sₓ = lmul!(s, norm_hatmatrix_x) + lower = muladd(-qt, sₓ, predictions) + upper = muladd(qt, sₓ, predictions) return LoessPrediction(predictions, lower, upper, sₓ, δ₁, δ₂, s, ρ) end end @@ -314,29 +318,31 @@ function hatmatrix(model::LoessModel{T}) where T L = zeros(T, n, n) for (i, r) in enumerate(eachrow(model.xs)) z = only(r) # we still only support one predictor - L[i, :] = _hatmatrix_x(model, z) + _hatmatrix_x!(view(L, i, :), model, z) end return L end # Compute elements of the hat matrix of `model` at `x`. # See Cleveland and Grosse 1991 p. 54. -function _hatmatrix_x(model::LoessModel{T}, x::Number) where T +function _hatmatrix_x!(Lx::AbstractVector{T}, model::LoessModel{T}, x::Number) where T + Base.require_one_based_indexing(Lx) n = length(model.ys) + @assert(length(Lx) == n) - adjacent_verts = traverse(model.kdtree, (T(x),)) - - @assert(length(adjacent_verts) == 2) - v₁, v₂ = adjacent_verts[1], adjacent_verts[2] + v₁, v₂ = traverse(model.kdtree, (T(x),)) - if x == v₁ || x == v₂ - Lx = model.hatmatrix_generator[[x]][1, :] + v₁_val = only(v₁) + v₂_val = only(v₂) + if x == v₁_val + copyto!(Lx, view(model.hatmatrix_generator[v₁], 1, :)) + elseif x == v₂_val + copyto!(Lx, view(model.hatmatrix_generator[v₂], 1, :)) else - Lx = zeros(T, n) - Lv₁ = model.hatmatrix_generator[[v₁]] - Lv₂ = model.hatmatrix_generator[[v₂]] + Lv₁ = model.hatmatrix_generator[v₁] + Lv₂ = model.hatmatrix_generator[v₂] for j in 1:n - b_int = cubic_interpolation(v₁, Lv₁[1, j], Lv₁[2, j], v₂, Lv₂[1, j], Lv₂[2, j]) + b_int = cubic_interpolation(v₁_val, Lv₁[1, j], Lv₁[2, j], v₂_val, Lv₂[1, j], Lv₂[2, j]) Lx[j] = evalpoly(x, b_int) end end @@ -403,9 +409,11 @@ Modifies: function tnormalize!(xs::AbstractMatrix{T}, q::T=0.1) where T <: AbstractFloat n, m = size(xs) cut = ceil(Int, (q * n)) - for j in 1:m - tmp = sort!(xs[:,j]) - xs[:,j] ./= mean(tmp[cut+1:n-cut]) + range = (cut - 1):(n - cut) + for xi in eachcol(xs) + copyto!(tmp, xi) + bulk = partialsort!(tmp, range) + rdiv!(xi, mean(bulk)) end xs end diff --git a/src/kd.jl b/src/kd.jl index cd05ea0..828230a 100644 --- a/src/kd.jl +++ b/src/kd.jl @@ -43,10 +43,8 @@ function KDTree( n, m = size(xs) bounds = Array{T}(undef, 2, m) - for j in 1:m - col = xs[:,j] - bounds[1, j] = minimum(col) - bounds[2, j] = maximum(col) + for (j, col) in enumerate(eachcol(xs)) + bounds[1, j], bounds[2, j] = extrema(col) end diam = diameter(bounds) @@ -83,7 +81,7 @@ Returns: """ function diameter(bounds::Matrix) - euclidean(vec(bounds[1,:]), vec(bounds[2,:])) + euclidean(view(bounds, 1, :), view(bounds, 2, :)) end """ @@ -223,7 +221,7 @@ function build_kdtree(xs::AbstractMatrix{T}, rightnode = build_kdtree(xs, view(perm,mid2:length(perm)), rightbounds, leaf_size_cutoff, leaf_diameter_cutoff, verts) - coords = Array{Array}(undef, m) + coords = Vector{Vector{T}}(undef, m) for i in 1:m if i == j coords[i] = [med] @@ -275,19 +273,20 @@ function traverse(kdtree::KDTree{T}, x::NTuple{N,T}) where {N,T} end end - bounds = copy(kdtree.bounds) + lowerbounds = kdtree.bounds[1, :] + upperbounds = kdtree.bounds[2, :] node = kdtree.root - return _traverse!(bounds, node, x) + return _traverse!(lowerbounds, upperbounds, node, x) end -_traverse!(bounds, node::Nothing, x) = bounds -function _traverse!(bounds, node::KDNode, x) +_traverse!(lowerbounds, upperbounds, node::Nothing, x) = lowerbounds, upperbounds +function _traverse!(lowerbounds, upperbounds, node::KDNode, x) if x[node.j] <= node.med - bounds[2, node.j] = node.med - return _traverse!(bounds, node.leftnode, x) + upperbounds[node.j] = node.med + return _traverse!(lowerbounds, upperbounds, node.leftnode, x) else - bounds[1, node.j] = node.med - return _traverse!(bounds, node.rightnode, x) + lowerbounds[node.j] = node.med + return _traverse!(lowerbounds, upperbounds, node.rightnode, x) end end