Various parameter/stats fixes

Project.toml

@@ -1,6 +1,6 @@
 name = "CurveFit"
 uuid = "5a033b19-8c74-5913-a970-47c3779ef25c"
-version = "1.9.2"
+version = "1.9.3"
 authors = ["Paulo José Saiz Jabardo <pjabardo@gmail.com>, Avik Pal <avikpal@mit.edu> and contributors"]
 
 [deps]

docs/src/_changelog.md

@@ -7,7 +7,27 @@ CurrentModule = CurveFit
 This documents notable changes in CurveFit.jl. The format is based on [Keep a
 Changelog](https://keepachangelog.com).
 
-## v[1.9.1] - 2026-04-25
+## [v1.9.3]
+
+### Fixed
+- [`LinearCurveFitAlgorithm`](@ref) will now automatically invert the intercept
+  when `yfun` is given to ensure that the returned parameters match the values
+  being fitted. Also affects [`ExpCurveFitAlgorithm`](@ref) and
+  [`PowerCurveFitAlgorithm`](@ref). *This is considered a bugfix rather than a
+  breaking change.*
+- Fixed the parameter handling and Jacobian of [`KingCurveFitAlgorithm`](@ref)
+  ([#112]).
+- Fixed the statistics methods for [`LinearCurveFitAlgorithm`](@ref) when a
+  transform is applied with `yfun` ([#112]). Previously the Jacobian for a
+  linear function would be computed on the residuals stored in the original
+  y-space.
+
+## [v1.9.2] - 2026-06-24
+
+### Changed
+- Various improvements to CI.
+
+## [v1.9.1] - 2026-04-25
 
 ### Fixed
 - Fixed support for `reinit!()`'ing nonlinear fit caches when using the new

src/common_interface.jl

@@ -223,6 +223,10 @@ be converted to a linear fit in a specific function space by choosing appropriat
 `xfun` and `yfun`. The `yfun_inverse` is used to convert the fitted values back to the
 original space (can be specified by defining `InverseFunctions.inverse`).
 
+When `yfun` is not `identity`, `yfun_inverse` is applied to the fitted intercept before
+storing it, so that `sol.u = (a, b)` directly satisfies the original-space formula
+(e.g. `y = b * exp(a*x)` for [`ExpCurveFitAlgorithm`](@ref)).
+
 This algorithm does not support bounds constraints (`lb`/`ub`).
 """
 function LinearCurveFitAlgorithm(;
@@ -291,7 +295,7 @@ ExpCurveFitAlgorithm() = LinearCurveFitAlgorithm(; xfun = identity, yfun = log)
 Represents a king curve fitting problem where `x` and `y` are the data points to
 fit. This algorithm does not support passing weights through `sigma` in
 [`CurveFitProblem`](@ref). This algorithm does not support bounds constraints (`lb`/`ub`).
-We want to solve for `a` and `b` according to original King's law (1910) that represents
+We want to solve for `A` and `B` according to original King's law (1910) that represents
 the relationship between voltage (E) and velocity (U) in a hotwire anemometer:
 
 ```math
@@ -304,7 +308,7 @@ or
 x^2 = A + B y^{1/2}
 ```
 """
-KingCurveFitAlgorithm() = LinearCurveFitAlgorithm(; xfun = abs2, yfun = sqrt)
+struct KingCurveFitAlgorithm <: AbstractCurveFitAlgorithm end
 
 @doc doc"""
     ModifiedKingCurveFitAlgorithm(alg::Union{Nothing, AbstractNonlinearAlgorithm} = nothing)

src/king.jl

@@ -5,9 +5,40 @@ function __king_fun!(resid, p, x)
     return nothing
 end
 
-# Common Solve Interface
+# Common Solve Interface for KingCurveFitAlgorithm
+@concrete struct KingFitCache <: AbstractCurveFitCache
+    prob <: CurveFitProblem
+    alg <: KingCurveFitAlgorithm
+    kwargs
+end
+
+function CommonSolve.init(prob::CurveFitProblem, alg::KingCurveFitAlgorithm; kwargs...)
+    @assert !is_nonlinear_problem(prob) "King's law fitting doesn't work with nlfunc specification."
+    @assert prob.u0 === nothing "King's law fit doesn't support initial guess (u0) specification"
+    sigma_not_supported(prob)
+    bounds_not_supported(prob)
+    return KingFitCache(prob, alg, kwargs)
+end
+
+function CommonSolve.solve!(cache::KingFitCache)
+    # Fit sqrt(U) = b + a*E^2, then recover King's constants A and B:
+    # sqrt(U) = -A/B + (1/B)*E^2  =>  B = 1/a, A = -b/a
+    b, a = __linear_fit_internal(abs2, cache.prob.x, sqrt, cache.prob.y, nothing)
+    B = 1 / a
+    A = -b * B
+    y_pred = ((cache.prob.x .^ 2 .- A) ./ B) .^ 2
+    resid = cache.prob.y .- y_pred
+    return CurveFitSolution(cache.alg, (A, B), resid, cache.prob, ReturnCode.Success)
+end
+
+function (sol::CurveFitSolution{<:KingCurveFitAlgorithm})(x::Number)
+    A, B = sol.u
+    return ((x^2 - A) / B)^2
+end
+
+# Common Solve Interface for ModifiedKingCurveFitAlgorithm
 @concrete struct ModifiedKingFitCache <: AbstractCurveFitCache
-    initial_guess_cache <: Union{Nothing, GenericLinearFitCache}
+    initial_guess_cache <: Union{Nothing, KingFitCache}
     nonlinear_cache
     prob <: CurveFitProblem
     alg <: ModifiedKingCurveFitAlgorithm

src/linfit.jl

@@ -65,14 +65,15 @@ function CommonSolve.solve!(cache::GenericLinearFitCache)
         resid ./= cache.prob.sigma
     end
 
+    b_stored = cache.alg.yfun_inverse(b)
     return CurveFitSolution(
-        cache.alg, (a, b), resid, cache.prob, ReturnCode.Success
+        cache.alg, (a, b_stored), resid, cache.prob, ReturnCode.Success
     )
 end
 
 function (sol::CurveFitSolution{<:LinearCurveFitAlgorithm})(x)
     a, b = sol.u
-    return sol.alg.yfun_inverse.(b .+ a .* sol.alg.xfun.(x))
+    return sol.alg.yfun_inverse.(sol.alg.yfun(b) .+ a .* sol.alg.xfun.(x))
 end
 
 # Polynomial Fit

src/rationalfit.jl

@@ -107,23 +107,19 @@ function CommonSolve.solve!(cache::LinearRationalFitCache)
     )
     cache.linsolve_cache.A = cache.mat
     sol = solve!(cache.linsolve_cache)
-    resid = if sol.resid === nothing
-        # Linear problem: y ≈ p/q => y*q - p ≈ 0 (linearized residual)
-        # But StatsAPI expects y - p/q (nonlinear residual) or linearized?
-        # Standard definition is y - y_pred.
-        # So we should compute y - p(x)/q(x) using the fitted params.
-
-        # We need to construct the RationalPolynomial to eval it.
-        # Helper function from rationalfit.jl isn't easily accessible inside solve! without alloc.
-        # But we can reuse the logic from call:
-        rpoly = RationalPolynomial(
-            view(sol.u, 1:(cache.alg.num_degree + 1)),
-            vcat(one(eltype(sol.u)), view(sol.u, (cache.alg.num_degree + 2):length(sol.u)))
-        )
-        cache.prob.y .- rpoly.(cache.prob.x)
-    else
-        sol.resid
-    end
+    # Linear problem: y ≈ p/q => y*q - p ≈ 0 (linearized residual)
+    # But StatsAPI expects y - p/q (nonlinear residual) or linearized?
+    # Standard definition is y - y_pred.
+    # So we should compute y - p(x)/q(x) using the fitted params.
+
+    # We need to construct the RationalPolynomial to eval it.
+    # Helper function from rationalfit.jl isn't easily accessible inside solve! without alloc.
+    # But we can reuse the logic from call:
+    rpoly = RationalPolynomial(
+        view(sol.u, 1:(cache.alg.num_degree + 1)),
+        vcat(one(eltype(sol.u)), view(sol.u, (cache.alg.num_degree + 2):length(sol.u)))
+    )
+    resid = cache.prob.y .- rpoly.(cache.prob.x)
     return CurveFitSolution(cache.alg, sol.u, resid, cache.prob, sol.retcode)
 end
 

src/stats.jl

@@ -102,11 +102,24 @@ end
 
 function jacobian(sol::CurveFitSolution{<:LinearCurveFitAlgorithm})
     x = sol.prob.x
-    xfun = sol.alg.xfun
-    J = Matrix{eltype(x)}(undef, length(x), 2)
-    J[:, 1] .= xfun.(x) # Slope
-    J[:, 2] .= 1        # Intercept
-    return J
+    alg = sol.alg
+    if alg.yfun === identity
+        # Linear model, compute Jacobian analytically
+        J = Matrix{eltype(x)}(undef, length(x), 2)
+        J[:, 1] .= alg.xfun.(x)
+        J[:, 2] .= 1
+        return J
+    else
+        # When there's a y-transform applied (e.g. power/exp fits) the residuals
+        # are in the original y-space so it would be incorrect to use the linear
+        # analytical Jacobian, instead we use DI.
+        a, b_stored = sol.u
+        u = [a, b_stored]
+        f_pred = u_curr -> begin
+            alg.yfun_inverse.(alg.yfun(u_curr[2]) .+ u_curr[1] .* alg.xfun.(x))
+        end
+        return DifferentiationInterface.jacobian(f_pred, AutoForwardDiff(), u)
+    end
 end
 
 function jacobian(sol::CurveFitSolution{<:PolynomialFitAlgorithm})
@@ -215,6 +228,19 @@ function jacobian(sol::CurveFitSolution{<:ExpSumFitAlgorithm})
     return DifferentiationInterface.jacobian(f_pred, AutoForwardDiff(), u)
 end
 
+function jacobian(sol::CurveFitSolution{<:KingCurveFitAlgorithm})
+    # King's law: U = t² where t = (E²−A)/B, so ∂U/∂A = −2t/B and ∂U/∂B = −2t²/B
+    A, B = sol.u
+    x = sol.prob.x
+    J = Matrix{eltype(x)}(undef, length(x), 2)
+    @inbounds for i in eachindex(x)
+        t = (x[i]^2 - A) / B
+        J[i, 1] = -2t / B    # ∂U/∂A
+        J[i, 2] = -2t^2 / B  # ∂U/∂B
+    end
+    return J
+end
+
 function jacobian(sol::CurveFitSolution{<:ModifiedKingCurveFitAlgorithm})
     # Modified King: E^2 = A + B * U^n
     # x corresponds to E (Voltage) - Model predicts x^2

test/king_law.jl

@@ -13,6 +13,9 @@ using Test
     prob = CurveFitProblem(E, U)
     sol = solve(prob, KingCurveFitAlgorithm())
 
+    @test sol.u[1] ≈ A
+    @test sol.u[2] ≈ B
+
     @testset for val in range(minimum(E), stop = maximum(E), length = 10)
         @test sol(val) ≈ fn(val)
     end

test/linfit_exp.jl

@@ -11,7 +11,7 @@ using Test
     sol = solve(prob, ExpCurveFitAlgorithm())
 
     @test sol.u[1] ≈ 0.8
-    @test sol.u[2] ≈ log(2.0)
+    @test sol.u[2] ≈ 2.0
 
     @testset for val in (0.0, 1.5, 4.5, 10.0)
         @test sol(val) ≈ fn(val)

test/linfit_power.jl

@@ -11,7 +11,7 @@ using Test
     sol = solve(prob, PowerCurveFitAlgorithm())
 
     @test sol.u[1] ≈ 0.8
-    @test sol.u[2] ≈ log(2.0)
+    @test sol.u[2] ≈ 2.0
 
     @testset for val in (0.0, 1.5, 4.5, 10.0)
         @test sol(val) ≈ fn(val)

test/stats.jl

@@ -245,12 +245,7 @@ using LinearSolve
         x_data = collect(1.0:0.5:5.0)
         y_data = 2.0 .* x_data .^ 0.5 .+ 0.1
 
-        # 1. King Fit (Linear alias)
-        # y = A + B*sqrt(x) -> Linear fit with yfun=sqrt? No, King is xfun=abs2, yfun=sqrt?
-        # Wait, King is LinearCurveFitAlgorithm(; xfun = abs2, yfun = sqrt)
-        # Check source: KingCurveFitAlgorithm() = LinearCurveFitAlgorithm(; xfun = abs2, yfun = sqrt)
-        # This implies sqrt(y) = A + B*x^2.
-        # This is strictly linear in transformed space, so vcov should work.
+        # 1. King Fit
         prob_king = CurveFitProblem(x_data, y_data)
         sol_king = solve(prob_king, KingCurveFitAlgorithm())
         @test size(vcov(sol_king)) == (2, 2)
@@ -281,4 +276,30 @@ using LinearSolve
         @test size(vcov(sol_rat_nl)) == (3, 3)
         @test all(stderror(sol_rat_nl) .> 0)
     end
+
+    @testset "Jacobian in original y-space for transformed fits" begin
+        x = collect(1.0:10.0)
+
+        # Power fit: y = b * x^a  =>  ∂y/∂a = y*log(x),  ∂y/∂b = y/b
+        a_true, b_true = 1.5, 3.0
+        y_pow = @. b_true * x^a_true
+        sol_pow = solve(CurveFitProblem(x, y_pow), PowerCurveFitAlgorithm())
+        a, b = sol_pow.u
+        @test [a, b] ≈ [a_true, b_true]
+        y_pred = sol_pow.(x)
+        J_pow = CurveFit.jacobian(sol_pow)
+        @test J_pow[:, 1] ≈ y_pred .* log.(x)
+        @test J_pow[:, 2] ≈ y_pred ./ b
+
+        # Exp fit: y = b * exp(a*x)  =>  ∂y/∂a = y*x,  ∂y/∂b = y/b
+        a_true, b_true = 0.3, 2.0
+        y_exp = @. b_true * exp(a_true * x)
+        sol_exp = solve(CurveFitProblem(x, y_exp), ExpCurveFitAlgorithm())
+        a, b = sol_exp.u
+        @test [a, b] ≈ [a_true, b_true]
+        y_pred = sol_exp.(x)
+        J_exp = CurveFit.jacobian(sol_exp)
+        @test J_exp[:, 1] ≈ y_pred .* x
+        @test J_exp[:, 2] ≈ y_pred ./ b
+    end
 end