Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Create Measurements package extension #248

Open
wants to merge 1 commit into
base: master
Choose a base branch
from

Conversation

chriswaudby
Copy link

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
Copy link
Author

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
Copy link
Member

pkofod commented Jan 16, 2024

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
Copy link

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants