Nonlinear curve fitting with weights with LsqFit

[achjaj]

Hello,

I need to fit the exponential decay law to some data with weights. I tried to use LsqFit.jl, however the standard_errors function returns too low values. I discovered it by comparing my results with my friend and I also made the fit in R.

[ufechner7]

Well, if you dont’t show us your code nobody will be able to help you. ):

[achjaj]

I am sorry, I did not think it would be useful as I just called the relevant functions, but anyway here it is:

using LsqFit

x = [0, 975, 1767, 2682, 5959]
y = [5490.0, 4546.0, 3679.0, 3134.0, 1451.0]
weights = [109.0, 102.0, 96.0, 89.0, 71.0]

@. model(t, p) = p[1]*(0.5)^(t/p[2])
fit  = curve_fit(model, x, y, weights, [0.1, 80])
err = round(standard_errors(fit)[2], sigdigits=1)
T = fit.param[2]

println("$T +- $err")

Just to clarify, the weights variable contains the stderr of measurement of y.

[ufechner7]

Sorry, please provide a piece of code that is completely defined. I tried your code and get the message:

julia> @. model(t, p) = p[1]*(0.5)^(t/p[2])
model (generic function with 1 method)

julia> fit  = curve_fit(model, x, y, weights, [0.1, 60])
ERROR: UndefVarError: x not defined
Stacktrace:
 [1] top-level scope
   @ REPL[3]:1

What are x, y and the weights?

Please read: make it easier to help you

Welcome to the Julia Discourse! We are enthusiastic about helping Julia programmers, both beginner and experienced. This public service announcement (PSA) outlines best practices when asking for help. Following these points makes it easier for us to help you and more likely you’ll get a prompt, useful answer. Keywords are highlighted to make it easier to refer to specific points. Choose a descriptive title that captures the key part of your question, eg “plots with multiple axes” instead of …

In particular:
"
Do your best to make your example self-contained (“minimal working example”, MWE ), so that it runs (or gets to the error that you want help with) as is. This means including package loading (e.g. using ThatPackage ) and any data that the code operates on. If your data is large or proprietary, generate example data if possible and include that. This StackOverflow answer explains how a DataFrame can be turned into a string representation that can be copy/pasted.
"

[rafael.guerra]

@achjaj

just to clarify, the weights variable contains the stderr of measurement of y.

This is probably the issue. In LsqFit.jl, the weights are assumed to be the inverse of the variances.

[achjaj]

I edited the code in my previous post.

[achjaj]

I tried that also, but the error is still suspiciously low.

[rafael.guerra]

Using the inverse of the variances as weights, we obtain using LsqFit :

p1 = 5541 +/- 86
p2 = 3135 +/- 110

Which from the QC plot below does not seem unreasonable.
What values do the other software packages compute?

[goerch]

Here is my take on it

using LsqFit

x = [0, 975, 1767, 2682, 5959]
y = [5490.0, 4546.0, 3679.0, 3134.0, 1451.0]
weights = [109.0, 102.0, 96.0, 89.0, 71.0]

@. model(t, p) = p[1]*(0.5)^(t/p[2])
fit  = curve_fit(model, x, y, [0.1, 80])
err = round(standard_errors(fit)[2], sigdigits=1)
T = fit.param[2]

println("$T +- $err")

@. model(t, p) = p[1]*(0.5)^(t/p[2])
fit  = curve_fit(model, x, y, weights, [0.1, 80])
err = round(standard_errors(fit)[2], sigdigits=1)
T = fit.param[2]

println("$T +- $err")

@. model(t, p) = p[1]*(0.5)^(t/p[2])
fit  = curve_fit(model, x, y, 1 ./ (weights .* weights), [0.1, 80])
err = round(standard_errors(fit)[2], sigdigits=1)
T = fit.param[2]

println("$T +- $err")

yielding

3158.802386807851 +- 100.0
3170.773189838312 +- 0.1
3134.79231889205 +- 100.0

[rafael.guerra]

It is much more clear to write:

σy = [109., 102., 96., 89., 71.]   # standard deviations of y
w = 1 ./ σy.^2                     # weights for LsqFit

[achjaj]

I see the problem now, I used 1/σ instead of 1/σ^2. Thank you, for your help!