From 3572a9d52d2117d2c70902ab82487b77abc6f5dc Mon Sep 17 00:00:00 2001 From: JamesWrigley Date: Sat, 27 Jun 2026 12:38:31 +0200 Subject: [PATCH 1/3] Use the residual degrees of freedom in margin_error() --- src/stats.jl | 2 +- test/stats.jl | 6 ++++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/stats.jl b/src/stats.jl index 5c7cd4e..808e5a4 100644 --- a/src/stats.jl +++ b/src/stats.jl @@ -386,7 +386,7 @@ See [`vcov`](@ref) for the meaning of `absolute_sigma`. """ function margin_error(sol::CurveFitSolution, alpha = 0.05; absolute_sigma::Bool = false, rtol::Real = NaN, atol::Real = 0) std_errors = stderror(sol; absolute_sigma, rtol, atol) - dist = TDist(dof(sol)) + dist = TDist(dof_residual(sol)) critical_values = quantile(dist, 1 - alpha / 2) return std_errors * critical_values end diff --git a/test/stats.jl b/test/stats.jl index e7df97b..8a4bae0 100644 --- a/test/stats.jl +++ b/test/stats.jl @@ -4,6 +4,7 @@ using StatsAPI using NonlinearSolveFirstOrder using LinearAlgebra using LinearSolve +using Distributions: TDist, quantile @testset "StatsAPI Integration" begin @testset "Linear Fit" begin @@ -92,6 +93,11 @@ using LinearSolve # Just checking structure and non-error @test cis[1][1] < cis[1][2] + # margin_error() should use the residual dof (n - p = 3) rather than the + # parameter dof (p = 2). + t_resid = quantile(TDist(3), 0.975) + @test margin_error(sol) ≈ stderror(sol) .* t_resid + # Test isconverged @test isconverged(sol) end From 22229fe7815ee6dd51ea065e05cdea6e5c475541 Mon Sep 17 00:00:00 2001 From: JamesWrigley Date: Sat, 27 Jun 2026 12:39:50 +0200 Subject: [PATCH 2/3] =?UTF-8?q?Make=20nonlinear=20fits=20compute=20`y=20?= =?UTF-8?q?=E2=88=92=20=C5=B7`=20as=20the=20residual?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/nonlinfit.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/nonlinfit.jl b/src/nonlinfit.jl index 9364c53..3d0fd04 100644 --- a/src/nonlinfit.jl +++ b/src/nonlinfit.jl @@ -20,7 +20,7 @@ end # Out-of-place function (nlf::NonlinearFunctionWrapper{false})(p, X) - resid = nlf.f(p, X) .- nlf.target + resid = nlf.target .- nlf.f(p, X) if !isnothing(nlf.sigma) resid ./= nlf.sigma @@ -32,7 +32,7 @@ end # In-place function (nlf::NonlinearFunctionWrapper{true})(resid, p, X) nlf.f(resid, p, X) - resid .-= nlf.target + resid .= nlf.target .- resid if !isnothing(nlf.sigma) resid ./= nlf.sigma From 942447a3eced7ac8ee0c8ca0271bd60f3f7fe978 Mon Sep 17 00:00:00 2001 From: JamesWrigley Date: Sat, 27 Jun 2026 12:51:22 +0200 Subject: [PATCH 3/3] Use the delta method for the convariance of transformed linear fits Also bumped the version. --- Project.toml | 2 +- docs/src/_changelog.md | 16 +++++++- src/common_interface.jl | 14 +++++++ src/stats.jl | 87 +++++++++++++++++++++++++++++++++++++++++ test/king_law.jl | 21 ++++++++++ test/stats.jl | 28 +++++++++++++ 6 files changed, 166 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 68b80a3..882dd55 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "CurveFit" uuid = "5a033b19-8c74-5913-a970-47c3779ef25c" -version = "1.9.3" +version = "1.9.4" authors = ["Paulo José Saiz Jabardo , Avik Pal and contributors"] [deps] diff --git a/docs/src/_changelog.md b/docs/src/_changelog.md index 13d0218..cdecc1d 100644 --- a/docs/src/_changelog.md +++ b/docs/src/_changelog.md @@ -7,7 +7,21 @@ CurrentModule = CurveFit This documents notable changes in CurveFit.jl. The format is based on [Keep a Changelog](https://keepachangelog.com). -## [v1.9.3] +## [v1.9.4] - 2026-06-27 + +### Changed +- Previously nonlinear fits would compute the residuals as `ŷ − y`, they are now + computed as `y − ŷ` to be consistent with the linear fits ([#114]). + +### Fixed +- Corrected [`margin_error()`](@ref) to use the residual degrees of freedom + rather than the degrees of freedom of the model ([#114]). +- Fixed the covariance calculation in [`vcov()`](@ref) to correctly handle the + uncertainties produced by linear fits of a transformed nonlinear model + (e.g. from [`PowerCurveFitAlgorithm`](@ref), [`ExpCurveFitAlgorithm`](@ref), + and [`KingCurveFitAlgorithm`](@ref)) by using the delta method ([#114]). + +## [v1.9.3] - 2026-06-26 ### Fixed - [`LinearCurveFitAlgorithm`](@ref) will now automatically invert the intercept diff --git a/src/common_interface.jl b/src/common_interface.jl index b9ee82e..095a910 100644 --- a/src/common_interface.jl +++ b/src/common_interface.jl @@ -266,6 +266,13 @@ This is equivalent to a linear fit in log-log space, i.e., ```math \log(y) = a \log(x) + \log(b) ``` + +Because the fit is performed in log space, it minimizes *relative* (rather than +absolute) error: a fixed fractional deviation counts the same at small and large +`y`. This is a different estimator from a nonlinear least-squares fit of the same +model (which minimizes absolute error), so the coefficients and their +uncertainties will generally differ. [`vcov`](@ref) accounts for this when +computing the uncertainty of `(a, b)`. """ PowerCurveFitAlgorithm() = LinearCurveFitAlgorithm(; xfun = log, yfun = log) @@ -286,6 +293,13 @@ This is equivalent to a linear fit in log-linear space, i.e., ```math \log(y) = a x + \log(b) ``` + +Because the fit is performed in log space, it minimizes *relative* (rather than +absolute) error: a fixed fractional deviation counts the same at small and large +`y`. This is a different estimator from a nonlinear least-squares fit of the same +model (which minimizes absolute error), so the coefficients and their +uncertainties will generally differ. [`vcov`](@ref) accounts for this when +computing the uncertainty of `(a, b)`. """ ExpCurveFitAlgorithm() = LinearCurveFitAlgorithm(; xfun = identity, yfun = log) diff --git a/src/stats.jl b/src/stats.jl index 808e5a4..7c5d16c 100644 --- a/src/stats.jl +++ b/src/stats.jl @@ -323,6 +323,26 @@ The covariance is then rescaled by reduced χ² (i.e. `mse(sol)`) so that `sigma acts as a relative weight. Pass `absolute_sigma = true` to skip this rescaling when `sigma` carries absolute physical uncertainties (analogous to scipy's `curve_fit(absolute_sigma=true)`). + +!!! note "Fits estimated in a transformed space" + Some fits estimate their coefficients by ordinary least squares (OLS) in a + *transformed* space rather than in original `y`-space: + + - [`PowerCurveFitAlgorithm`](@ref) and [`ExpCurveFitAlgorithm`](@ref) fit a + straight line in log space, minimizing relative rather than absolute error. + - [`KingCurveFitAlgorithm`](@ref) fits `sqrt(U) = b + a·E²` by OLS and + recovers `(A, B)` from `(b, a)`. [`ModifiedKingCurveFitAlgorithm`](@ref) + is not affected since it's a regular nonlinear fit. + - Any [`LinearCurveFitAlgorithm`](@ref) with `yfun ≠ identity`. + + To keep the reported uncertainty consistent with how the coefficients were + actually estimated, the covariance is computed in that transformed space and + then mapped back to the original parameters via the delta method, rather than + from the original-space Jacobian used above. + + Note that [`residuals`](@ref), [`rss`](@ref), and [`mse`](@ref) remain in + original `y`-space as fit diagnostics, so for these fits `mse(sol)` is *not* + the variance scaling behind `vcov`. """ function StatsAPI.vcov(sol::CurveFitSolution; absolute_sigma::Bool = false) J = jacobian(sol) @@ -351,6 +371,73 @@ function StatsAPI.vcov(sol::CurveFitSolution; absolute_sigma::Bool = false) return covar end +# Covariance for fits whose coefficients are estimated by ordinary least squares +# in a *transformed* space but whose parameters are reported in original +# space. Computing `vcov` from the original-space Jacobian (as the generic +# method does) gives a Gauss–Newton covariance of a different, nonlinear fit, +# not the uncertainty of the estimator that actually produced `sol.u`. Instead +# we form the transformed-space least-squares covariance `σ²·(XᵀX)⁻¹` for the +# straight line `Y = β0 + β1·t` and delta-method it to the stored +# parameterization via `param_map(β0, β1) -> [p1, p2]`. +function _transformed_ols_vcov(t, Y, β0, β1, param_map, absolute_sigma::Bool) + n = length(t) + T = float(promote_type(eltype(t), eltype(Y))) + st = sum(t) + stt = sum(abs2, t) + det = n * stt - st^2 + + # Residual variance in transformed space (sigma is unsupported for these + # fits, so this is the unweighted estimate). absolute_sigma assumes unit + # variance instead, mirroring the generic method. + σ2 = if absolute_sigma + one(T) + else + rss_t = sum(i -> abs2(Y[i] - (β0 + β1 * t[i])), eachindex(t)) + rss_t / (n - 2) + end + + # σ²·(XᵀX)⁻¹ for (β0 intercept, β1 slope) + Σβ = Matrix{T}(undef, 2, 2) + Σβ[1, 1] = σ2 * stt / det # Var(β0) + Σβ[2, 2] = σ2 * n / det # Var(β1) + Σβ[1, 2] = Σβ[2, 1] = -σ2 * st / det # Cov(β0, β1) + + # Delta method: G = ∂(p1, p2)/∂(β0, β1); Cov(p) = G·Σβ·Gᵀ. + g = β -> param_map(β[1], β[2]) + G = DifferentiationInterface.jacobian(g, AutoForwardDiff(), [β0, β1]) + return G * Σβ * G' +end + +# vcov() for linear fits of a transformed model, falls back to the default +# implementation if there's no y-transform. +function StatsAPI.vcov(sol::CurveFitSolution{<:LinearCurveFitAlgorithm}; absolute_sigma::Bool = false) + alg = sol.alg + if alg.yfun === identity + return @invoke StatsAPI.vcov(sol::CurveFitSolution; absolute_sigma) + end + + a, b_stored = sol.u + t = alg.xfun.(sol.prob.x) + Y = alg.yfun.(sol.prob.y) + β0 = alg.yfun(b_stored) # intercept + β1 = a # slope + # Generic in yfun_inverse so custom inverses work. + param_map = (β0, β1) -> [β1, alg.yfun_inverse(β0)] + return _transformed_ols_vcov(t, Y, β0, β1, param_map, absolute_sigma) +end + +# King's law: OLS of `sqrt(U) = b + a·E²` (intercept b, slope a), with stored +# params (A, B) = (-b/a, 1/a). Same transformed-space estimator as power/exp. +function StatsAPI.vcov(sol::CurveFitSolution{<:KingCurveFitAlgorithm}; absolute_sigma::Bool = false) + A, B = sol.u + t = abs2.(sol.prob.x) # E² + Y = sqrt.(sol.prob.y) # sqrt(U) + a = 1 / B # slope + b = -A / B # intercept (A = -b·B ⇒ b = -A/B) + param_map = (β0, β1) -> [-β0 / β1, 1 / β1] # (A, B) + return _transformed_ols_vcov(t, Y, b, a, param_map, absolute_sigma) +end + """ stderror(sol::CurveFitSolution; absolute_sigma = false, rtol = NaN, atol = 0) diff --git a/test/king_law.jl b/test/king_law.jl index e7476aa..92988b1 100644 --- a/test/king_law.jl +++ b/test/king_law.jl @@ -24,3 +24,24 @@ using Test prob_sigma = CurveFitProblem(E, U; sigma = ones(length(U))) @test_throws AssertionError solve(prob_sigma, KingCurveFitAlgorithm()) end + +@testset "King vcov (delta method)" begin + # King fits sqrt(U) = b + a·E² by OLS, then A = -b/a, B = 1/a. vcov must + # describe that transformed-space estimator, not the U-space Jacobian. + E = collect(1.0:0.25:3.0) + A_true, B_true = 1.5, 0.8 + U = ((E .^ 2 .- A_true) ./ B_true) .^ 2 .+ 1.0e-3 .* sin.(1:length(E)) + + sol = solve(CurveFitProblem(E, U), KingCurveFitAlgorithm()) + + # Reference: plain linear OLS on (E², sqrt(U)) gives the (slope a, + # intercept b) covariance; delta-method it through (A, B) = (-b/a, 1/a). + sol_lin = solve(CurveFitProblem(E .^ 2, sqrt.(U)), LinearCurveFitAlgorithm()) + a, b = sol_lin.u # (slope, intercept) + Vlin = vcov(sol_lin) # ordered (a, b) + Σ = [Vlin[2, 2] Vlin[2, 1]; Vlin[1, 2] Vlin[1, 1]] # reordered to (b, a) + G = [-1 / a b / a^2; 0.0 -1 / a^2] # ∂(A, B)/∂(b, a) + Vref = G * Σ * G' + + @test vcov(sol) ≈ Vref +end diff --git a/test/stats.jl b/test/stats.jl index 8a4bae0..cb4a97f 100644 --- a/test/stats.jl +++ b/test/stats.jl @@ -139,6 +139,34 @@ using Distributions: TDist, quantile @test all(stderror(sol_rat) .> 0) end + @testset "Custom yfun_inverse vcov (delta method)" begin + # Non-exp transform: sqrt(y) = β0 + β1*x => y = (β0 + β1*x)^2. + # The delta factor d(b)/d(β0) = 2β0 ≠ b, so this catches an exp-only or + # missing delta-method bug that the power/exp tests would not. + x = collect(1.0:10.0) + β0, β1 = 2.0, 0.5 + noise = [0.05, -0.04, 0.03, -0.02, 0.06, -0.05, 0.01, -0.03, 0.04, -0.01] + Y = β0 .+ β1 .* x .+ noise + y = Y .^ 2 + + alg = LinearCurveFitAlgorithm(; yfun = sqrt, yfun_inverse = z -> z^2) + sol = solve(CurveFitProblem(x, y), alg) + + # Check that the fit succeeded + @test SciMLBase.successful_retcode(sol) + @test sol.u[1] ≈ β1 atol = 1.0e-2 + @test sol.u[2] ≈ β0^2 atol = 1.0e-1 + + # Reference: plain linear OLS on (x, sqrt(y)) gives the transformed-space + # covariance for (slope β1, intercept β0); delta-method it through b = β0^2. + sol_lin = solve(CurveFitProblem(x, sqrt.(y)), LinearCurveFitAlgorithm()) + Vlin = vcov(sol_lin) + d = 2 * sol_lin.u[2] # = 2β0 + Vref = [Vlin[1, 1] d * Vlin[1, 2]; d * Vlin[2, 1] d^2 * Vlin[2, 2]] + + @test vcov(sol) ≈ Vref + end + @testset "Explicit API Coverage (ExpSum)" begin # Test ExpSumFitAlgorithm explicitly to ensure jacobian works # y = k + p*exp(lam*x)