JuliaNLSolvers / LsqFit.jl

Simple curve fitting in Julia
Other
317 stars 78 forks source link

Create Measurements package extension #248

Open chriswaudby opened 1 year ago

chriswaudby commented 1 year ago

This PR creates a simple extension for the Measurements package, allowing data with uncertainties to be passed directly for fitting, e.g.:

using LsqFit, Measurements
x = [1.0, 2.0, 3.0]
y = [1 ± 0.1, 2.2 ± 0.2, 2.9 ± 0.3]
model(x, p) = p[1] * x
p0 = [1.]
curve_fit(model, x, y, p0)

The extension creates a simple wrapper for LsqFit.curve_fit that dispatches on ydata::AbstractArray{Measurement{T}} where T. Values and uncertainties are extracted from the input, then weights are calculated as the inverse square of the uncertainties and passed to the main curve_fit function.

A unit test checks that the extension gives the same results as when weights are passed manually.

An __init__ function is defined in LsqFit.jl to load the extension for Julia versions <1.9 (see discussion at https://discourse.julialang.org/t/package-extensions-for-julia-1-9/93397).

chriswaudby commented 11 months ago

Hi @pkofod, haven't made a pull request before so don't know the etiquette, but it's been a few weeks since I submitted this - are you able to take a look? Let me know if there's anything else I should be doing!

pkofod commented 9 months ago

Thanks, not it's not bad etiquette to bump :) Sorry for not seeing this.

I think sometimes there can be corner cases with automatic differentiation where things can go wrong, so if we are to merge this PR, I think the tests have to be significantly extended to cover the various ways which the user-facing functions can be called: different combinations of inputs, different choices for AD, etc.

TheFibonacciEffect commented 1 month ago

This works very well, thanks a lot :D

using LsqFit, Plots, Measurements
@. gauss(x,p) = p[1] + p[2] * exp(-(x-p[3])^2/p[4]^2)
x = range(-1,1,100)
yerr = gauss(x,[0,1,0,0.5]) + 0.1*randn(length(x)) .± 0.1
fit = curve_fit(gauss, x,yerr,Float64[0,1,0,1])
begin
    scatter(x,yerr)
    plot!(x,gauss(x,fit.param);ribbon = gauss(x,fit.param.+stderror(fit)) - gauss(x,fit.param.-stderror(fit)))
end
savefig("fit.png")

fit