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

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.

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?

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

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

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

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.

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.

1 Like

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.

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.

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?

1 Like

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.

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.

2 Likes

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.

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

1 Like

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`:

``````"""
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:

``````"""
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

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)

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.

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.

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.

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.

3 Likes

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:

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.

1 Like

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

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(?

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/

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 and other optimization benchmarks would be good.

1 Like