From d683795bbda0a49fb0d451641baf74bfebff0bc2 Mon Sep 17 00:00:00 2001
From: JamesWrigley <james@puiterwijk.org>
Date: Sun, 28 Jun 2026 23:51:23 +0200
Subject: [PATCH 1/6] Fix `sol.resid` of modified King fits to be in velocity
 space
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The residuals are computed in E² but that's not very intuitive and is not
consistent with the rest of the API.
---
 src/king.jl   | 10 +++++++++-
 src/stats.jl  | 18 ++++++++++++++++++
 test/stats.jl |  2 ++
 3 files changed, 29 insertions(+), 1 deletion(-)

diff --git a/src/king.jl b/src/king.jl
index aaac21f..f6d06b5 100644
--- a/src/king.jl
+++ b/src/king.jl
@@ -105,7 +105,15 @@ function CommonSolve.solve!(cache::ModifiedKingFitCache)
     )
 
     sol = solve(nonlinear_prob, __FallbackNonlinearFitAlgorithm(cache.alg.alg); cache.kwargs...)
-    return CurveFitSolution(cache.alg, sol.u, sol.resid, cache.prob, sol.retcode, sol.original)
+
+    # The inner solve fits E² = A + B·Uⁿ, so its residual `A + B·Uⁿ − E²` is a
+    # difference of E² values. We instead report `residuals` as a difference of
+    # velocities, `U − Û`, so they match `fitted`/`sol(x)` (i.e.
+    # `residuals == prob.y .- fitted`). `vcov` rebuilds the E² residual it needs
+    # from `sol.u`.
+    A, B, n = sol.u
+    resid = @. cache.prob.y - ((cache.prob.x^2 - A) / B)^(1 / n)
+    return CurveFitSolution(cache.alg, sol.u, resid, cache.prob, sol.retcode, sol.original)
 end
 
 function (sol::CurveFitSolution{<:ModifiedKingCurveFitAlgorithm})(x)
diff --git a/src/stats.jl b/src/stats.jl
index 5ec0c39..f6f1ed8 100644
--- a/src/stats.jl
+++ b/src/stats.jl
@@ -443,6 +443,24 @@ function StatsAPI.vcov(sol::CurveFitSolution{<:KingCurveFitAlgorithm}; absolute_
     return _transformed_ols_vcov(t, Y, b, a, param_map, absolute_sigma)
 end
 
+# Modified King fits E² = A + B·Uⁿ by nonlinear least squares, so its covariance
+# is the Gauss–Newton form `σ²·(JᵀJ)⁻¹` built from the E² residual it minimized.
+# `residuals`/`mse` are reported as velocities, so we rebuild that E² residual
+# here instead of reusing `mse(sol)`.
+function StatsAPI.vcov(sol::CurveFitSolution{<:ModifiedKingCurveFitAlgorithm}; absolute_sigma::Bool = false)
+    covar = _vcov_from_jacobian(jacobian(sol))
+
+    if !absolute_sigma
+        A, B, n = sol.u
+        E, U = sol.prob.x, sol.prob.y
+        resid_fitspace = @. A + B * U^n - E^2
+        σ2 = sum(abs2, resid_fitspace) / dof_residual(sol)
+        covar .*= σ2
+    end
+
+    return covar
+end
+
 """
     stderror(sol::CurveFitSolution; absolute_sigma = false, rtol = NaN, atol = 0)
 
diff --git a/test/stats.jl b/test/stats.jl
index 6bd4658..6de7f54 100644
--- a/test/stats.jl
+++ b/test/stats.jl
@@ -208,6 +208,8 @@ using Distributions: TDist, quantile
 
         @test size(vcov(sol_mod_king)) == (3, 3)
         @test all(stderror(sol_mod_king) .> 0)
+        # Residuals are reported in velocity space, consistent with fitted()
+        @test residuals(sol_mod_king) ≈ y_k .- fitted(sol_mod_king)
     end
 
 

From c8fd4e7bb7cc97df5686505afba56245c2d16c6e Mon Sep 17 00:00:00 2001
From: JamesWrigley <james@puiterwijk.org>
Date: Wed, 1 Jul 2026 14:50:00 +0200
Subject: [PATCH 2/6] Remove/cleanup code from ModifiedKing/RationalPolynomial

- These fit types don't properly support `reinit!()` so they always build a cache
  inside of `solve!()`, meaning that their caches `nonlinear_cache` fields are
  never used.
- Delete the custom Jacobian method for King fits. This is unused now that we
  use the delta method for constructing the covariance matrix.
- Fix the modified King fit tests to pass x/y in the right order, previously
  they were swapped.
---
 src/king.jl        | 17 +----------------
 src/rationalfit.jl | 19 +------------------
 src/stats.jl       | 13 -------------
 test/stats.jl      | 17 ++++++-----------
 4 files changed, 8 insertions(+), 58 deletions(-)

diff --git a/src/king.jl b/src/king.jl
index f6d06b5..ea63d0d 100644
--- a/src/king.jl
+++ b/src/king.jl
@@ -39,7 +39,6 @@ end
 # Common Solve Interface for ModifiedKingCurveFitAlgorithm
 @concrete struct ModifiedKingFitCache <: AbstractCurveFitCache
     initial_guess_cache <: Union{Nothing, KingFitCache}
-    nonlinear_cache
     prob <: CurveFitProblem
     alg <: ModifiedKingCurveFitAlgorithm
     kwargs
@@ -66,21 +65,7 @@ function CommonSolve.init(
         init(nobounds_prob, KingCurveFitAlgorithm(); kwargs...)
     end
 
-    nonlinear_cache = init(
-        NonlinearCurveFitProblem(
-            NonlinearFunction{true}(
-                __king_fun!;
-                resid_prototype = similar(prob.x)
-            ),
-            similar(prob.x, 3),
-            stack((prob.x, prob.y); dims = 1),
-            nothing;
-            lb = prob.lb, ub = prob.ub
-        ),
-        __FallbackNonlinearFitAlgorithm(alg.alg);
-        kwargs...
-    )
-    return ModifiedKingFitCache(initial_guess_cache, nonlinear_cache, prob, alg, kwargs)
+    return ModifiedKingFitCache(initial_guess_cache, prob, alg, kwargs)
 end
 
 function CommonSolve.solve!(cache::ModifiedKingFitCache)
diff --git a/src/rationalfit.jl b/src/rationalfit.jl
index 64c3679..adfc6a2 100644
--- a/src/rationalfit.jl
+++ b/src/rationalfit.jl
@@ -48,7 +48,6 @@ end
 
 @concrete struct NonlinearRationalFitCache <: AbstractCurveFitCache
     initial_guess_cache <: Union{Nothing, LinearRationalFitCache}
-    nonlinear_cache
     prob <: CurveFitProblem
     alg <: RationalPolynomialFitAlgorithm
     kwargs
@@ -82,23 +81,7 @@ function CommonSolve.init(
             A, init(LinearProblem(A, prob.y), alg.alg; kwargs...), prob, alg, kwargs
         )
     end
-    nonlinear_cache = init(
-        NonlinearCurveFitProblem(
-            NonlinearFunction{true}(
-                __rational_fit_residual!(alg.num_degree, alg.den_degree);
-                resid_prototype = similar(prob.x)
-            ),
-            similar(prob.x, coeffs_length),
-            prob.x,
-            prob.y;
-            lb = prob.lb, ub = prob.ub
-        ),
-        __FallbackNonlinearFitAlgorithm(alg.alg);
-        kwargs...
-    )
-    return NonlinearRationalFitCache(
-        initial_guess_cache, nonlinear_cache, prob, alg, kwargs
-    )
+    return NonlinearRationalFitCache(initial_guess_cache, prob, alg, kwargs)
 end
 
 function CommonSolve.solve!(cache::LinearRationalFitCache)
diff --git a/src/stats.jl b/src/stats.jl
index f6f1ed8..f3441eb 100644
--- a/src/stats.jl
+++ b/src/stats.jl
@@ -226,19 +226,6 @@ 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
diff --git a/test/stats.jl b/test/stats.jl
index 6de7f54..a76298f 100644
--- a/test/stats.jl
+++ b/test/stats.jl
@@ -193,23 +193,18 @@ using Distributions: TDist, quantile
         @test all(stderror(sol_nc)[2:end] .> 0)
         @test confint(sol_nc)[1] == (0.0, 0.0)
 
-        # Modified King Fit (E^2 = A + B * U^n)
-        # x corresponds to E (Voltage) in Jacobian logic, but input data order is (U, E^2)?
-        # User creates CurveFitProblem(x, y).
-        # King assumes (Voltage, Velocity).
-        # My implementation: x=E, y=U in Jacobian.
-        # But in verify script I passed (U, E^2).
-        # So I will replicate verify script logic here.
+        # Modified King Fit: E² = A + B·Uⁿ, with x = E (voltage), y = U (velocity)
         A_true = 1.0; B_true = 2.0; n_true = 0.5
-        x_k = collect(1.0:0.5:5.0)
-        y_k = @. A_true + B_true * x_k^n_true
-        prob_mod_king = CurveFitProblem(x_k, y_k)
+        U = collect(1.0:0.5:5.0)
+        E = sqrt.(A_true .+ B_true .* U .^ n_true)
+        prob_mod_king = CurveFitProblem(E, U)
         sol_mod_king = solve(prob_mod_king, ModifiedKingCurveFitAlgorithm())
 
+        @test sol_mod_king.u ≈ [A_true, B_true, n_true]
         @test size(vcov(sol_mod_king)) == (3, 3)
         @test all(stderror(sol_mod_king) .> 0)
         # Residuals are reported in velocity space, consistent with fitted()
-        @test residuals(sol_mod_king) ≈ y_k .- fitted(sol_mod_king)
+        @test residuals(sol_mod_king) ≈ U .- fitted(sol_mod_king)
     end
 
 

From e14ecfdedfbf5771b870bd5ebc7575bdb77a9ada Mon Sep 17 00:00:00 2001
From: JamesWrigley <james@puiterwijk.org>
Date: Sun, 28 Jun 2026 23:55:31 +0200
Subject: [PATCH 3/6] Add asserts to reinit!() to validate inputs

---
 Project.toml             |  2 +-
 docs/src/_changelog.md   | 11 +++++++++++
 src/nonlinfit.jl         |  6 ++++++
 test/nonlinfit_reinit.jl | 18 ++++++++++++++++++
 4 files changed, 36 insertions(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 9e4ea76..367e6f4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "CurveFit"
 uuid = "5a033b19-8c74-5913-a970-47c3779ef25c"
-version = "1.10.0"
+version = "1.10.1"
 authors = ["Paulo José Saiz Jabardo <pjabardo@gmail.com>, Avik Pal <avikpal@mit.edu> and contributors"]
 
 [deps]
diff --git a/docs/src/_changelog.md b/docs/src/_changelog.md
index 5ab7f37..55f0749 100644
--- a/docs/src/_changelog.md
+++ b/docs/src/_changelog.md
@@ -7,6 +7,17 @@ CurrentModule = CurveFit
 This documents notable changes in CurveFit.jl. The format is based on [Keep a
 Changelog](https://keepachangelog.com).
 
+## [v1.10.1] - 2026-06-29
+
+### Fixed
+- Fixed the residuals space of [`ModifiedKingCurveFitAlgorithm`](@ref),
+  previously they were in E² space but now they're calculated in velocity space
+  such that they're consistent with the standard definition of `residuals = y -
+  ŷ` ([#116]).
+- The `reinit!()` method for [`NonlinearCurveFitProblem`](@ref)'s will now
+  explicitly throw an error when invalid inputs are given (e.g. arrays of the
+  wrong length) to prevent internal corruption ([#115]).
+
 ## [v1.10.0] - 2026-06-28
 
 ### Added
diff --git a/src/nonlinfit.jl b/src/nonlinfit.jl
index 0b132ea..3fc7e0d 100644
--- a/src/nonlinfit.jl
+++ b/src/nonlinfit.jl
@@ -54,23 +54,29 @@ end
 
 function SciMLBase.reinit!(cache::GenericNonlinearCurveFitCache; u0 = nothing, x = nothing, y = nothing, sigma = nothing, kwargs...)
     if !isnothing(u0)
+        @assert size(u0) == size(cache.u0) "reiniting `u0` must keep the same size"
         kwargs = (; kwargs..., u0)
         copyto!(cache.u0, u0)
     end
 
     # x becomes `p` (parameter) in the NonlinearLeastSquaresProblem
     if !isnothing(x)
+        @assert size(x) == size(_get_cache(cache).p) "reiniting `x` must keep the same size"
         kwargs = (; kwargs..., p = x)
     end
 
     # Update `y` inplace
     wrapper = _unwrap_nonlinear_function(cache.cache.prob.f.f)
     if !isnothing(y)
+        @assert wrapper isa NonlinearFunctionWrapper "cannot reinit `y` for a problem created without a `y`"
+        @assert size(y) == size(wrapper.target) "reiniting `y` must keep the same size"
         copyto!(wrapper.target, y)
     end
 
     # Update `sigma` inplace
     if !isnothing(sigma)
+        @assert wrapper isa NonlinearFunctionWrapper && !isnothing(wrapper.sigma) "cannot reinit `sigma` for a problem created without a `sigma`"
+        @assert size(sigma) == size(wrapper.sigma) "reiniting `sigma` must keep the same size"
         copyto!(wrapper.sigma, sigma)
     end
 
diff --git a/test/nonlinfit_reinit.jl b/test/nonlinfit_reinit.jl
index d4a24ee..e9283d3 100644
--- a/test/nonlinfit_reinit.jl
+++ b/test/nonlinfit_reinit.jl
@@ -49,4 +49,22 @@ using NonlinearSolveBase: NonlinearSolveBase
     @test cache.cache.prob.f.f isa NonlinearSolveBase.AutoSpecializeCallable
     CurveFit.reinit!(cache; u0 = [1.0, 1.0, 1.0], x, y)
     @test solve!(cache).u ≈ a0 atol = 1.0e-7
+
+    # Test reinit!() input validation with an out-of-place function so the
+    # `y=nothing` cache can be constructed.
+    g(a, x) = @. a[1] + a[2] * x^a[3]
+    cache = CurveFit.init(NonlinearCurveFitProblem(g, [0.5, 0.5, 0.5], x, y, sigma))
+
+    # Size changes are rejected
+    @test_throws AssertionError CurveFit.reinit!(cache; u0 = [1.0, 1.0])
+    @test_throws AssertionError CurveFit.reinit!(cache; x = 1.0:5.0)
+    @test_throws AssertionError CurveFit.reinit!(cache; y = y[1:5])
+    @test_throws AssertionError CurveFit.reinit!(cache; sigma = sigma[1:5])
+
+    # Can't reinit a variable the problem was created without
+    cache = CurveFit.init(NonlinearCurveFitProblem(g, [0.5, 0.5, 0.5], x, y))
+    @test_throws AssertionError CurveFit.reinit!(cache; sigma)
+
+    cache = CurveFit.init(NonlinearCurveFitProblem(g, [0.5, 0.5, 0.5], x))
+    @test_throws AssertionError CurveFit.reinit!(cache; y)
 end

From 49bbb6671cc7eb1b4afc6e730536c642fd6869fa Mon Sep 17 00:00:00 2001
From: JamesWrigley <james@puiterwijk.org>
Date: Wed, 1 Jul 2026 15:59:27 +0200
Subject: [PATCH 4/6] Fix polynomial fit residuals

---
 docs/src/_changelog.md    | 6 ++++--
 src/linfit.jl             | 5 +++--
 test/linfit_polynomial.jl | 2 ++
 3 files changed, 9 insertions(+), 4 deletions(-)

diff --git a/docs/src/_changelog.md b/docs/src/_changelog.md
index 55f0749..7b1189b 100644
--- a/docs/src/_changelog.md
+++ b/docs/src/_changelog.md
@@ -7,7 +7,7 @@ CurrentModule = CurveFit
 This documents notable changes in CurveFit.jl. The format is based on [Keep a
 Changelog](https://keepachangelog.com).
 
-## [v1.10.1] - 2026-06-29
+## [v1.10.1] - 2026-07-01
 
 ### Fixed
 - Fixed the residuals space of [`ModifiedKingCurveFitAlgorithm`](@ref),
@@ -16,7 +16,9 @@ Changelog](https://keepachangelog.com).
   ŷ` ([#116]).
 - The `reinit!()` method for [`NonlinearCurveFitProblem`](@ref)'s will now
   explicitly throw an error when invalid inputs are given (e.g. arrays of the
-  wrong length) to prevent internal corruption ([#115]).
+  wrong length) to prevent internal corruption ([#116]).
+- Fixed the residuals and statistics calculations for polynomial fits,
+  previously they were computed from corrupted data ([#116]).
 
 ## [v1.10.0] - 2026-06-28
 
diff --git a/src/linfit.jl b/src/linfit.jl
index 68c7ea8..36701c3 100644
--- a/src/linfit.jl
+++ b/src/linfit.jl
@@ -106,8 +106,9 @@ function CommonSolve.solve!(cache::PolynomialFitCache)
     __vandermondepoly!(cache.vandermondepoly_cache, cache.prob.x, cache.alg.degree)
     cache.linsolve_cache.A = cache.vandermondepoly_cache
     sol = solve!(cache.linsolve_cache)
-    # Always compute residuals manually as LinearSolve may return incorrect residuals
-    resid = cache.prob.y .- cache.vandermondepoly_cache * sol.u
+    # Compute residuals from the fitted polynomial. We can't reuse
+    # vandermondepoly_cache here as the in-place factorization overwrites it.
+    resid = cache.prob.y .- evalpoly.(cache.prob.x, Ref(sol.u))
     return CurveFitSolution(cache.alg, sol.u, resid, cache.prob, sol.retcode)
 end
 
diff --git a/test/linfit_polynomial.jl b/test/linfit_polynomial.jl
index 31cf99b..118080a 100644
--- a/test/linfit_polynomial.jl
+++ b/test/linfit_polynomial.jl
@@ -21,6 +21,8 @@ using LinearSolve
         @test sol(val) ≈ fn(val)
     end
 
+    @test sol.resid ≈ y .- sol.(x)
+
     @testset "ill-conditioned" begin
         true_coeffs = [80.0, -5.0e-18, -7.0e-20, -1.0e-36]
         x1 = 1.0e10 .* (0:0.1:5)

From 14b6ac6b881d5046949f8f2086d5c2b14b97f1be Mon Sep 17 00:00:00 2001
From: JamesWrigley <james@puiterwijk.org>
Date: Thu, 2 Jul 2026 12:55:05 +0200
Subject: [PATCH 5/6] Copy x/y inputs of expsum fits so that they won't be
 modified

This would previously introduce small errors by first scaling and then inverting
the scaling. Also made the residuals return a vector instead of a matrix.
---
 docs/src/_changelog.md |  3 +++
 src/expsumfit.jl       | 24 +++++++++++++-----------
 2 files changed, 16 insertions(+), 11 deletions(-)

diff --git a/docs/src/_changelog.md b/docs/src/_changelog.md
index 7b1189b..1083507 100644
--- a/docs/src/_changelog.md
+++ b/docs/src/_changelog.md
@@ -19,6 +19,9 @@ Changelog](https://keepachangelog.com).
   wrong length) to prevent internal corruption ([#116]).
 - Fixed the residuals and statistics calculations for polynomial fits,
   previously they were computed from corrupted data ([#116]).
+- Solving with [`ExpSumFitAlgorithm`](@ref) will no longer modify the input x/y
+  arrays, and will return a vector of residuals instead of a `(n, 1)` matrix
+  ([#116]).
 
 ## [v1.10.0] - 2026-06-28
 
diff --git a/src/expsumfit.jl b/src/expsumfit.jl
index 6689304..a00dda0 100644
--- a/src/expsumfit.jl
+++ b/src/expsumfit.jl
@@ -137,6 +137,9 @@ end
     prob <: CurveFitProblem
     alg
     kwargs
+    # Scaled working copies of prob.x/prob.y, the inputs are never modified
+    x <: AbstractVector
+    y <: AbstractVector
     Y <: AbstractMatrix
     S <: AbstractMatrix
     A <: AbstractVector
@@ -162,6 +165,8 @@ function CommonSolve.init(
         prob,
         alg,
         kwargs,
+        similar(prob.x),
+        similar(prob.y),
         similar(prob.x, nY, mY),
         similar(prob.x, nY, alg.n - 1),
         similar(prob.x, mY),
@@ -174,28 +179,25 @@ end
 
 # TODO: allocations in this function aren't optimized
 function CommonSolve.solve!(cache::ExpSumFitCache)
-    sc = __expsum_scale!(cache.prob.x, cache.prob.y)
+    copyto!(cache.x, cache.prob.x)
+    copyto!(cache.y, cache.prob.y)
+    sc = __expsum_scale!(cache.x, cache.y)
 
     __expsum_fill_Y!(
-        cache.Y, cache.S, cache.coeff, cache.prob.x, cache.prob.y, cache.alg.n, cache.alg.m
+        cache.Y, cache.S, cache.coeff, cache.x, cache.y, cache.alg.n, cache.alg.m
     )
 
     λ = __expsum_solve_λ(
-        cache.Y, cache.A, cache.Ā, cache.prob.y, cache.alg.n, cache.alg.m
+        cache.Y, cache.A, cache.Ā, cache.y, cache.alg.n, cache.alg.m
     )
 
     X = isreal(λ) ? cache.Xr : cache.Xc
-    __expsum_fill_X!(cache.prob.x, λ, X, cache.alg.n)
+    __expsum_fill_X!(cache.x, λ, X, cache.alg.n)
 
     qrX = qr!(X)
-    p = qrX \ cache.prob.y
+    p = qrX \ cache.y
     isreal(p) && (p = real(p))
 
-    for i in eachindex(cache.prob.x, cache.prob.y)
-        cache.prob.x[i] *= sc.x
-        cache.prob.y[i] *= sc.y
-    end
-
     for i in eachindex(p)
         p[i] *= sc.y
     end
@@ -214,7 +216,7 @@ function CommonSolve.solve!(cache::ExpSumFitCache)
     end
     backing = (; k, p, λ)
 
-    y_pred = k .+ sum(exp.(cache.prob.x * λ') .* p'; dims = 2)
+    y_pred = k .+ vec(sum(exp.(cache.prob.x * λ') .* p'; dims = 2))
     resid = cache.prob.y .- y_pred
 
     return CurveFitSolution(

From 797651378c734aea44527f4d158788d5deab857f Mon Sep 17 00:00:00 2001
From: JamesWrigley <james@puiterwijk.org>
Date: Thu, 2 Jul 2026 13:20:08 +0200
Subject: [PATCH 6/6] Add support for sigma with residual-only fits, and fix
 stats methods

CurveFit supports passing `y=nothing` when creating a NonlinearCurveFitProblem,
which means that we treat the function as returning residuals directly. But
__wrap_nonlinear_function() would silently drop sigma when y is nothing while
leaving it in the original problem, which meant that the statistics methods
would assume sigma was used during the fit and apply it incorrectly. Now we
always use sigma if it's passed.

Also fixed support for `sol(x)` for in-place functions, and made `nobs()` able
to handle solutions where y is nothing.
---
 docs/src/_changelog.md     |  3 +++
 src/common_interface.jl    |  6 ++++--
 src/nonlinfit.jl           | 22 +++++++++++++++++-----
 src/stats.jl               |  2 +-
 test/nonlinfit_weighted.jl | 22 ++++++++++++++++++++++
 test/stats.jl              | 22 ++++++++++++++++++++++
 6 files changed, 69 insertions(+), 8 deletions(-)

diff --git a/docs/src/_changelog.md b/docs/src/_changelog.md
index 1083507..2c5fd53 100644
--- a/docs/src/_changelog.md
+++ b/docs/src/_changelog.md
@@ -22,6 +22,9 @@ Changelog](https://keepachangelog.com).
 - Solving with [`ExpSumFitAlgorithm`](@ref) will no longer modify the input x/y
   arrays, and will return a vector of residuals instead of a `(n, 1)` matrix
   ([#116]).
+- Fixed various statistics methods when using residual-only functions with
+  [`NonlinearCurveFitProblem`](@ref), as well as adding support for `sigma` with
+  such such functions ([#116])
 
 ## [v1.10.0] - 2026-06-28
 
diff --git a/src/common_interface.jl b/src/common_interface.jl
index 095a910..25b0c4c 100644
--- a/src/common_interface.jl
+++ b/src/common_interface.jl
@@ -99,8 +99,10 @@ optimization problem is solved:
 \argmin_u ~ \left\| f(u, x) - y \right\|_2
 ```
 
-If `y` is `nothing`, then it is treated as a zero vector. `f` is a generic Julia function or
-ideally a `NonlinearFunction` from [`SciMLBase.jl`](https://github.com/SciML/SciMLBase.jl).
+If `y` is `nothing`, then it is treated as a zero vector and `f` is expected to
+return residuals directly. `f` is a generic Julia function or ideally a
+`NonlinearFunction` from
+[`SciMLBase.jl`](https://github.com/SciML/SciMLBase.jl).
 
 Lower and upper bounds on the parameters can be passed through `lb` and `ub` keyword
 arguments. These should be arrays with the same length as `u0`.
diff --git a/src/nonlinfit.jl b/src/nonlinfit.jl
index 3fc7e0d..6408092 100644
--- a/src/nonlinfit.jl
+++ b/src/nonlinfit.jl
@@ -11,17 +11,22 @@ _unwrap_nonlinear_function(f::NonlinearSolveBase.AutoSpecializeCallable) = _unwr
 _unwrap_nonlinear_function(f::NonlinearSolveBase.BoundedWrapper) = _unwrap_nonlinear_function(f.f.f)
 _unwrap_nonlinear_function(f) = f
 
-# If `target` is nothing then we can completely ignore sigma
-__wrap_nonlinear_function(f::NonlinearFunction, ::Nothing, _) = f
+__wrap_nonlinear_function(f::NonlinearFunction, ::Nothing, ::Nothing) = f
 function __wrap_nonlinear_function(f::NonlinearFunction, target, sigma)
     internal_f = NonlinearFunctionWrapper{SciMLBase.isinplace(f)}(target, sigma, f.f)
     @set! f.f = internal_f
-    @set! f.resid_prototype = similar(target)
+    @set! f.resid_prototype = similar(isnothing(target) ? sigma : target)
     return f
 end
 
 # Out-of-place
 function (nlf::NonlinearFunctionWrapper{false})(p, X)
+    # `target` is only nothing for residual-style problems (y = nothing), in
+    # which case the wrapper exists purely to apply `sigma`.
+    if isnothing(nlf.target)
+        return nlf.f(p, X) ./ nlf.sigma
+    end
+
     resid = nlf.target .- nlf.f(p, X)
 
     if !isnothing(nlf.sigma)
@@ -34,8 +39,10 @@ end
 # In-place
 function (nlf::NonlinearFunctionWrapper{true})(resid, p, X)
     nlf.f(resid, p, X)
-    resid .= nlf.target .- resid
 
+    if !isnothing(nlf.target)
+        resid .= nlf.target .- resid
+    end
     if !isnothing(nlf.sigma)
         resid ./= nlf.sigma
     end
@@ -68,7 +75,7 @@ function SciMLBase.reinit!(cache::GenericNonlinearCurveFitCache; u0 = nothing, x
     # Update `y` inplace
     wrapper = _unwrap_nonlinear_function(cache.cache.prob.f.f)
     if !isnothing(y)
-        @assert wrapper isa NonlinearFunctionWrapper "cannot reinit `y` for a problem created without a `y`"
+        @assert wrapper isa NonlinearFunctionWrapper && !isnothing(wrapper.target) "cannot reinit `y` for a problem created without a `y`"
         @assert size(y) == size(wrapper.target) "reiniting `y` must keep the same size"
         copyto!(wrapper.target, y)
     end
@@ -139,6 +146,11 @@ function CommonSolve.solve!(cache::GenericNonlinearCurveFitCache)
 end
 
 function (sol::CurveFitSolution{<:__FallbackNonlinearFitAlgorithm})(x)
+    if SciMLBase.isinplace(sol.prob) && x isa AbstractArray
+        out = similar(sol.resid, length(x))
+        sol.prob.nlfunc(out, sol.u, x)
+        return out
+    end
     return sol.prob.nlfunc(sol.u, x)
 end
 
diff --git a/src/stats.jl b/src/stats.jl
index f3441eb..5cbad83 100644
--- a/src/stats.jl
+++ b/src/stats.jl
@@ -53,7 +53,7 @@ end
 Return the number of observations used in the fit.
 """
 function StatsAPI.nobs(sol::CurveFitSolution)
-    return length(sol.prob.y)
+    return length(sol.resid)
 end
 
 # Return the positions in `sol.u` that are held fixed. These coefficients
diff --git a/test/nonlinfit_weighted.jl b/test/nonlinfit_weighted.jl
index c49fdc8..9665959 100644
--- a/test/nonlinfit_weighted.jl
+++ b/test/nonlinfit_weighted.jl
@@ -28,3 +28,25 @@ using SciMLBase
     @test err_weighted < err_no_weight
     @test SciMLBase.successful_retcode(sol_weighted)
 end
+
+@testset "Weighted residual-style problem (y = nothing)" begin
+    x = collect(1.0:10.0)
+    y = @. 2.0 + 3.0 * x + 0.1 * sin(x)
+    sigma = collect(range(0.5, 2.0; length = length(x)))
+    u0 = [1.0, 1.0]
+
+    pred(u, x) = @. u[1] + u[2] * x
+    sol_pred = solve(NonlinearCurveFitProblem(pred, u0, x, y, sigma))
+
+    # Residual-style problems (y = nothing) must apply sigma the same way as
+    # the equivalent prediction-style problem.
+    resid_oop(u, x) = y .- pred(u, x)
+    sol_oop = solve(NonlinearCurveFitProblem(resid_oop, u0, x, nothing, sigma))
+    @test sol_oop.u ≈ sol_pred.u
+    @test sol_oop.resid ≈ sol_pred.resid
+
+    resid_iip(resid, u, x) = resid .= y .- pred(u, x)
+    sol_iip = solve(NonlinearCurveFitProblem(resid_iip, u0, x, nothing, sigma))
+    @test sol_iip.u ≈ sol_pred.u
+    @test sol_iip.resid ≈ sol_pred.resid
+end
diff --git a/test/stats.jl b/test/stats.jl
index a76298f..031da6a 100644
--- a/test/stats.jl
+++ b/test/stats.jl
@@ -285,6 +285,28 @@ using Distributions: TDist, quantile
         @test vcov(sol) ≈ vcov(sol; absolute_sigma = true) .* mse(sol)
     end
 
+    @testset "Residual-style problems (y = nothing)" begin
+        x = collect(1.0:10.0)
+        y = @. 2.0 + 3.0 * x + 0.1 * sin(x)
+
+        f_oop(u, x) = @. u[1] + u[2] * x - $y
+        f_iip(resid, u, x) = @. resid = u[1] + u[2] * x - $y
+        fn_iip = NonlinearFunction{true}(f_iip; resid_prototype = similar(x))
+
+        # Smoke test that the stats functions handle problems without a `y`
+        for f in (f_oop, fn_iip)
+            sol = solve(NonlinearCurveFitProblem(f, [1.0, 1.0], x))
+
+            @test nobs(sol) == length(x)
+            @test dof_residual(sol) == length(x) - 2
+            @test residuals(sol) == sol.resid
+            @test predict(sol) ≈ residuals(sol)
+            @test mse(sol) ≈ rss(sol) / dof_residual(sol)
+            @test length(stderror(sol)) == 2
+            @test all(lo < c < hi for (c, (lo, hi)) in zip(coef(sol), confint(sol)))
+        end
+    end
+
     @testset "Extended Algorithm Test Coverage" begin
         x_data = collect(1.0:0.5:5.0)
         y_data = 2.0 .* x_data .^ 0.5 .+ 0.1