# Best-fit parameter error bar using Optimization.jl?

**URL:** <https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186>\
**Category:** General Usage\
**Tags:** question, optimization, curve-fitting\
**Created:** [August 25, 2023, 11:02am UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186 "2023-08-25T11:02:40Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![abehersan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abehersan/32/47221_2.png) [@abehersan](https://discourse.julialang.org/u/abehersan)\
**Post date:** [August 25, 2023, 11:02am UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186/1 "2023-08-25T11:02:40Z")

</div>

Specifically in the context of local solvers with and without constraints.

What is the standard way people report error bars on fitting parameters while using `Optimization.jl`?

In principle this can be calculated if the covariance matrix is known but not even that is returned upon solving.

---

<div class="post-metadata">

**Author:** ![maxkapur](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxkapur/32/21208_2.png) [@maxkapur](https://discourse.julialang.org/u/maxkapur)\
**Post date:** [August 25, 2023, 11:32am UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186/2 "2023-08-25T11:32:35Z")

</div>

Need more context. If you are fitting a curve using Optimization.jl, then it should just compute the best parameters according to your objective function; there is no error to plot (assuming the solver converged). Maybe you want to show a confidence interval for the predicted function values?

---

<div class="post-metadata">

**Author:** ![abehersan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abehersan/32/47221_2.png) [@abehersan](https://discourse.julialang.org/u/abehersan)\
**Post date:** [August 25, 2023, 11:39am UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186/3 "2023-08-25T11:39:29Z")

</div>

Precisely, calculating confidence intervals (what I call “error bars”) can be done straightforwardly if a covariance (or non-singular Hessian) is returned by the solver, see e.g. [https://arxiv.org/pdf/1210.3781.pdf](https://arxiv.org/pdf/1210.3781.pdf)

This appears to not be the case in `Optimization.jl`

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 25, 2023, 1:54pm UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186/4 "2023-08-25T13:54:18Z")

</div>

It’s generally not computed since most optimization would be more expensive if that was used. You can just directly calculate it.

---

<div class="post-metadata">

**Author:** ![abehersan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abehersan/32/47221_2.png) [@abehersan](https://discourse.julialang.org/u/abehersan)\
**Post date:** [August 25, 2023, 3:22pm UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186/5 "2023-08-25T15:22:42Z")

</div>

I think it makes sense to add maybe a flag in the `solve` call to include confidence interval calculations in the output.

For a fit to be meaningful (or to be ruled out as meaningless) the best-fit parameters, some goodness-of-fit statistics as well as confidence intervals for the parameters are necessary.

If the computation of the cost function is efficient (it should be) then, the added computational cost in the way of calculating CIs is absolutely worth it.

Is this a feature that makes sense to add? Is it requested by people?  
I imagine it would take `Optimization.jl` closer in utility to `scipy.optimize`’s functions.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 25, 2023, 3:25pm UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186/6 "2023-08-25T15:25:45Z")

</div>

Instead of a flag, maybe a post-processing function?

I don’t think most people would want this confidence interval since it’s not generally correct for nonlinear problems. It’s only correct in the linear approximation so it’s making quite a bit of assumptions in the nonlinear. case. Using a Bayesian approach like Turing.jl would generally give more accurate and interpretable confidence intervals so that would likely be recommended in most cases which are nonlinear.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [August 25, 2023, 5:19pm UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186/7 "2023-08-25T17:19:22Z")

</div>

Yeah, for nonlinear problems, you need to actually have found the global optimum for the estimate to make any sense at all, and if your residuals are non-Gaussian, it still makes no sense.

---

<div class="post-metadata">

**Author:** ![abehersan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abehersan/32/47221_2.png) [@abehersan](https://discourse.julialang.org/u/abehersan)\
**Post date:** [August 25, 2023, 8:03pm UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186/8 "2023-08-25T20:03:37Z")

</div>

> [@ChrisRackauckas](#):
>
> it’s not generally correct for nonlinear problems

True, but in the case where error in user’s data is assumed as Gaussian and uncorrelated then it is a fair approximation to make. Indeed this is the assumption lots of optimization libraries make in order to approximate a confidence interval, this is not the case for `Optimization.jl` as of now.

> [@baggepinnen](#):
>
> you need to actually have found the global optimum for the estimate to make any sense at all

My above comment applies here as well, under assumptions about the nature of the uncertainty in the data. If the errors are a-priori known to be non-Gaussian, other optimization route is of course necessary.

My main stance is that I think that `Optimization.jl` would strongly benefit from offering a comparable utility to general purpose fitting programs such as e.g. `gnuplot`, `lmfit` in python or even classic fitting frameworks such as `Minuit` without needing to resort to `pycall` or C++ bindings.

@ChrisRackauckas you mentioned a post-processing function to estimate confidence-intervals around a found local minimum. Do you see something like that ever making it in `Optimization.jl` for the rough but albeit common gaussian-error approximation?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 26, 2023, 8:04am UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186/9 "2023-08-26T08:04:22Z")

</div>

> [@abehersan](#):
>
> Do you see something like that ever making it in `Optimization.jl` for the rough but albeit common gaussian-error approximation?

Yes, if someone makes a PR I’d accept it. I think it’s a feature worth having and documenting well, though it has enough caveats that I personally won’t be prioritizing implementing it.

---

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [August 26, 2023, 8:18am UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186/10 "2023-08-26T08:18:22Z")

</div>

> [@abehersan](#):
>
> Indeed this is the assumption lots of optimization libraries make in order to approximate a confidence interval, this is not the case for `Optimization.jl` as of now.

What other optimization libraries do this? Optimization.jl isn’t a statistics package.

I don’t think it would be appropriate to have a post-processing function like you recommend in Optimization.jl since it’s easily misused by people (e.g. I’ve seen to many `+/- 2se` type intervals for cases where it makes no sense…) / it doesn’t make sense for the majority of inputs. Wouldn’t it be better off in some other statistics package? Especially since (last I checked?) Optimization.jl doesn’t really have any post-processing functions in it.

---

<div class="post-metadata">

**Author:** ![maxkapur](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxkapur/32/21208_2.png) [@maxkapur](https://discourse.julialang.org/u/maxkapur)\
**Post date:** [August 27, 2023, 5:21pm UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186/11 "2023-08-27T17:21:05Z")

</div>

I agree; Optimization.jl is neither a statistics package nor a plotting package and should do one thing (optimization) and do it well.

That said, I am surprised that OP’s use case isn’t covered by one of the existing nonlinear regression tools elsewhere in the Julia ecosystem. @abehersan , perhaps you can post a minimal version of your model and someone may be able to point you in the right direction.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [August 27, 2023, 5:24pm UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186/12 "2023-08-27T17:24:08Z")

</div>

> [@maxkapur](#):
>
> That said, I am surprised that OP’s use case isn’t covered by one of the existing nonlinear regression tools elsewhere in the Julia ecosystem.

[https://julianlsolvers.github.io/LsqFit.jl/latest/tutorial/#Goodness-of-Fit-1](https://julianlsolvers.github.io/LsqFit.jl/latest/tutorial/#Goodness-of-Fit-1)

---

<div class="post-metadata">

**Author:** ![abehersan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abehersan/32/47221_2.png) [@abehersan](https://discourse.julialang.org/u/abehersan)\
**Post date:** [August 27, 2023, 6:21pm UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186/13 "2023-08-27T18:21:27Z")

</div>

My use case is extremely simple. Peak lineshape fitting to experimental data. Error bars on the data are known. A minimal example of a simple peak fit is shown below when working with `python` and `lmfit`:

```python
"""
lorentzian.py

Generate noisy data that is Lorentzian in nature, fit a Lorentzian to the
data and report parameter uncertainties using LMFIT
"""

import numpy as np
import matplotlib.pyplot as plt
from lmfit import Model
np.random.seed(777)

def lorentzian(x, a, c, s):
    L = (a / np.pi) * (s / ((x - c) **2 + s** 2))
    return L

"""true parameters"""
a = 10.0
c = 0.0
s = 0.02

x_data = np.linspace(-0.5, 0.5, 101)
noise = np.random.normal(0, 7, len(x_data))
y_data = lorentzian(x_data, a, c, s) + noise
dy_data = np.sqrt(y_data - noise) # assumed poissonian noise
data_xye = np.array([x_data, y_data, dy_data]).T
np.savetxt("lorentzian_peak.dat", data_xye, delimiter="\t")

"""fit with lmfit"""
model = Model(lorentzian)
model.nan_policy = "omit"
params = model.make_params(a=17.0, c=0.0, s=0.04)
fit_result = model.fit(y_data, params, x=x_data, weights=1.0/dy_data)
fit_ax = fit_result.plot_fit(numpoints=701, fit_kws={"lw":3})
fit_fig = plt.gcf()
fit_fig.suptitle("Lorentzian fit with LMFIT")
fit_fig.set_size_inches(11.7, 8.3)

"""fit statistics"""
fit_stats = fit_result.fit_report()
fit_ax.text(0.05, 0.5, fit_stats, transform=fit_ax.transAxes)
plt.tight_layout()
plt.savefig("lmfit_fit.pdf", dpi=300)
plt.show()

```

An analogous work flow in Julia is as follows:

```julia
"""
lorentzian.jl

Fit a Lorentzian using Optimization.jl
"""

using Optimization, OptimizationNLopt
using Plots
default(size=(1080, 720))
default(margin=5Plots.mm)
default(linewidth=3.5)

function lorentzian(x::Real, a::Real, c::Real, s::Real)::Float64
    (a / pi) * (s / ((x - c)^2 + s^2))
end

function chisqr_lorentzian(u, p)::Float64
    x_data = p[1]
    y_data = p[2]
    dy_data = p[3]
    y_calc = [lorentzian(x, u...) for x in x_data]
    sum( ((y_data .- y_calc) ./ dy_data) .^ 2 )
end

"""load data generated with python"""
dfile = open("lorentzian_peak.dat", "r")
x_data = Float64[]
y_data = Float64[]
dy_data = Float64[]
for line in eachline(dfile)
    parts = split(line)
    push!(x_data, parse(Float64, parts[1]))
    push!(y_data, parse(Float64, parts[2]))
    push!(dy_data, parse(Float64, parts[3]))
end

u0 = Float64[17, 0, 0.04]
data_xye = [x_data, y_data, dy_data]

optfunc_lor = OptimizationFunction(chisqr_lorentzian)
optprob_lor = OptimizationProblem(optfunc_lor, u0, data_xye)
sol = solve(optprob_lor, NLopt.LN_NELDERMEAD())

s_names = ["amplitude", "center", "sigma"]
sol_str = [n*": "*string(u) for (n, u) in zip(s_names, sol)]
push!(sol_str, "FWHM: $(2*sol[3])")

p = plot(xlabel="x", ylabel="y", title="Lorentzian fit with Optimization.jl")
scatter!(p, x_data, y_data, yerr=dy_data, label="data")
plot!(p, x_data, [lorentzian(x, sol...) for x in x_data], label="fit")
annotate!(p, (-0.3, 100, join(sol_str, "\n")), annotatiohnalign=:left)

savefig(p, "optimization_fit.pdf")

```

The fit results naturally agree. Of course the added benefit of `lmfit` is the calculated CI of the fit parameters from `scipy`´s covariance matrix computation.

For my use case, in principle the Hessian can be calculated since for most lineshapes I work with, analytical expressions for the derivative of the model w.r.t. the fit parameters are straightforward to calculate.

> [@maxkapur](#):
>
> Optimization.jl is neither a statistics package nor a plotting package

I understand `Optimization.jl` is not a statistics package per se, but it strikes me as odd that it optimizes very well but then it does not provide confidence interval for the optimized parameters (even in the Gaussian approximation) nor any goodness-of-fit metric beyond the value of the objective function at the found minimum.

I also understand that for more complex cases pushing a “+/-” on the users’ fit results willy-nilly is dangerous. However I would contend that _most_ use cases are not _all_ use cases.

> [@maxkapur](#):
>
> I am surprised that OP’s use case isn’t covered by one of the existing nonlinear regression tools elsewhere in the Julia ecosystem

This exact surprise is what I felt and what prompted me to start this conversation topic in the first place. And while `LsqFit` is a counter-example, it’s still quite disheartening to have to give up the whole amazing machinery of `Optimization.jl` just because you want an error bar on your parameters post-fit.

---

<div class="post-metadata">

**Author:** ![abehersan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abehersan/32/47221_2.png) [@abehersan](https://discourse.julialang.org/u/abehersan)\
**Post date:** [August 27, 2023, 6:25pm UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186/14 "2023-08-27T18:25:25Z")

</div>

![Screenshot 2023-08-27 at 20.23.10](https://global.discourse-cdn.com/julialang/original/3X/7/6/76cc0b3d613c47349204dada0fb8d61095bea09a.png)

---

<div class="post-metadata">

**Author:** ![maxkapur](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxkapur/32/21208_2.png) [@maxkapur](https://discourse.julialang.org/u/maxkapur)\
**Post date:** [August 27, 2023, 6:41pm UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186/15 "2023-08-27T18:41:46Z")

</div>

The terms “confidence interval” and “goodness of fit” are quite foreign to the domain of optimization. Optimization solvers solve _math problems_ of the form `minimize f(x) subject to g(x) <= 0`. You either find the answer or you don’t.

It would be weird if Julia returned a confidence interval on `sum(x)` just because it is “common” for `x` to be a vector of observations from a statistical distribution, right? 🙂 If you really wanted that, you would do `using StatsBase; mean_and_std(x) .* length(x)` or whatever.

While there would be utility in a separate library that does nonlinear regression using Optimization.jl on the back end and returns the covariance matrices you need—indeed, Julia’s composability means it shouldn’t be terribly hard for someone to build such a library—as someone who solves optimization problems every day, I would much prefer Optimization.jl be kept at a more limited feature and dependency set.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 30, 2023, 7:31pm UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186/16 "2023-08-30T19:31:11Z")

</div>

> [@maxkapur](#):
>
> While there would be utility in a separate library that does nonlinear regression using Optimization.jl on the back end and returns the covariance matrices you need—indeed, Julia’s composability means it shouldn’t be terribly hard for someone to build such a library—as someone who solves optimization problems every day, I would much prefer Optimization.jl be kept at a more limited feature and dependency set.

This wouldn’t require any dependencies. And no, having such functionality would be inline with other parts of the SciML interface. DifferentialEquations has nice solution interface features that people like to use:

> **[Solution Handling · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/)**
>
> Documentation for DifferentialEquations.jl.

> **[Plot Functions · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/basics/plot/)**
>
> Documentation for DifferentialEquations.jl.

Etc. LinearSolve.jl, NonlinearSolve.jl, Integrals.jl, Optimization.jl, etc. should all follow suit. And I think a feature that allows one to do a post-analysis that adds uncertainties to the parameter results under some (well-documented) assumptions is precisely the kind of post-solution analysis functionality that would be good to offer, along with nice pre-defined logging functions, benchmarking tooling, etc.

---

<div class="post-metadata">

**Author:** ![abehersan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abehersan/32/47221_2.png) [@abehersan](https://discourse.julialang.org/u/abehersan)\
**Post date:** [August 30, 2023, 8:40pm UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186/17 "2023-08-30T20:40:26Z")

</div>

Thanks for considering the utility of such post-processing!

I think the challenging part of this extension to Optimization.jl is that different solvers return a solution struct with different fields and is also not type-inferrable, see e.g. [`solve` not type-inferrable · Issue #556 · SciML/Optimization.jl · GitHub](https://github.com/SciML/Optimization.jl/issues/556)

I would be glad to contribute if some pointers are hinted at since I’m not familiar with the whole of SciML suite, maybe with “good first-issue” tags in the repo’s issues(?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 31, 2023, 11:43am UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186/18 "2023-08-31T11:43:33Z")

</div>

Follow up with @Vaibhavdixit02. I recommend joining the Julia slack and the #sciml-bridged channel where a lot of the dev chat takes place. @Vaibhavdixit02 can probably spend a bit of time today tagging some good first issues and helping get you up to speed. There’s a lot of good first issues in that repo.

I think better error throwing on things not supported by a specific algorithm would be a good one. Or another would be to contribute to the benchmark system. We just got Optimization.jl benchmarks up:

[https://docs.sciml.ai/SciMLBenchmarksOutput/stable/Optimizaton/2drosenbrock/](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/Optimizaton/2drosenbrock/)

but it’s really jank right now. Cleaning that up, but also systematizing [Get BlackBoxOptimizationBenchmarking.jl updated to Optimization.jl's interface and run on SciMLBenchmarks · Issue #640 · SciML/SciMLBenchmarks.jl · GitHub](https://github.com/SciML/SciMLBenchmarks.jl/issues/640) and other optimization benchmarks would be good.

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [July 2, 2025, 10:18pm UTC](https://discourse.julialang.org/t/best-fit-parameter-error-bar-using-optimization-jl/103186/19 "2025-07-02T22:18:01Z")

</div>

> [@abehersan](#):
>
> classic fitting frameworks such as `Minuit`

coming from this, I wonder if you found a solution.

> **[Hesse and Minos — iminuit 2.31.1 compiled with ROOT-v6-25-02-9213-g754d22635f...](https://scikit-hep.org/iminuit/notebooks/hesse_and_minos.html)**

pretty much looking for the `HESSE` approximation here, and I found that take Hessian using auto diff package for my likelihood yields poor results ☹

well, I can try to use [GitHub - JuliaHEP/Minuit2.jl: Julia bindings to Minuit2 minimisation package](https://github.com/JuliaHEP/Minuit2.jl) but I want to know if we have something in Julia ecosystem yet.