More stats fixes

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 <pjabardo@gmail.com>, Avik Pal <avikpal@mit.edu> and contributors"]
 
 [deps]

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

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)
 

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

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)
 
@@ -386,7 +473,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

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

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
@@ -133,6 +139,34 @@ using LinearSolve
         @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)