Coding

Julia 101 – 0.1. Least-Square fitting

year = Float64[0]

for i in 1:6
  push!(year, i)
end

year = year .+ 2013
7-element Array{Float64,1}:
 2013.0
 2014.0
 2015.0
 2016.0
 2017.0
 2018.0
 2019.0
year_x = collect(Float64, 2013:1:2020);
h = [2, 2, 2, 3, 4, 6, 8, 12];
using Plots
scatter(year_x, h, 
  xlab = "year", ylab = "h-index", 
  xlim = [2012, 2023], ylim = [0,15], framestyle = :box)

scatter!(year_x[3:end], h[3:end])
using LsqFit

@. model(x, p) = p[1] + x*p[2] + x^2.0*p[3]  
fit = curve_fit(model, year_x, h, [2.0, 0.5, -3.0])
LsqFit.LsqFitResult{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Int64,1}}([1.1832864618761656e6, -1174.946148708822, 0.2916665972735931], [0.20833284547552466, -0.19642864260822535, -0.017856936203315854, -0.2559520348440856, 0.08928606053814292, 0.017857350641861558, 0.5297618352342397, -0.3750004852190614], [1.0000000000082452 2012.9999999932165 4.05216899997484e6; 1.0000000000082452 2014.0000000073403 4.05619599998367e6; … ; 1.0000000000082452 2018.9999999797847 4.076360999975335e6; 1.0000000000082452 2019.9999999939084 4.08039999999393e6], true, Int64[])
T = typeof(fit)
for (name, typ) in zip(fieldnames(T), T.types)
    println("type of the fieldname name istyp")
end
type of the fieldname param is Array{Float64,1}
type of the fieldname resid is Array{Float64,1}
type of the fieldname jacobian is Array{Float64,2}
type of the fieldname converged is Bool
type of the fieldname wt is Array{Int64,1}
?fit
search: fit LsqFit filter filter! curve_fit first firstindex isfinite popfirst!

No documentation found.

fit is of type LsqFit.LsqFitResult{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Int64,1}}.

Summary

struct LsqFit.LsqFitResult{Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Int64,1}} <: Any

Fields

param     :: Array{Float64,1}
resid     :: Array{Float64,1}
jacobian  :: Array{Float64,2}
converged :: Bool
wt        :: Array{Int64,1}
fitted_p = fit.param
3-element Array{Float64,1}:
     1.1832864618761656e6
 -1174.946148708822
     0.2916665972735931
new_x = collect(Float64, 2013:1:2025)
new_y = fitted_p[1] .+ new_x.*fitted_p[2] .+ (new_x.^2).*fitted_p[3] 
plot!(new_x, new_y, ylim = [0,30], seriestype =:line)
fit.resid
8-element Array{Float64,1}:
  0.20833284547552466
 -0.19642864260822535
 -0.017856936203315854
 -0.2559520348440856
  0.08928606053814292
  0.017857350641861558
  0.5297618352342397
 -0.3750004852190614
fit.jacobian
8×3 Array{Float64,2}:
 1.0  2013.0  4.05217e6
 1.0  2014.0  4.0562e6
 1.0  2015.0  4.06022e6
 1.0  2016.0  4.06426e6
 1.0  2017.0  4.06829e6
 1.0  2018.0  4.07232e6
 1.0  2019.0  4.07636e6
 1.0  2020.0  4.0804e6
fit.wt
0-element Array{Int64,1}
fit.converged
true
Astronomy · computer · Note

Least-Squares Fitting of data containing errors

(This is a note personal only! )

There is a powerful tool dealing with data containing both “X” and “Y” errors to make a lease-square fitting in IDL, that is mpfit. This routine allows you to use a self-defined function to do the fitting. Also, you can set up some conditions such as the parameters’ limit or fixed parameters etc. This is fabulous! So beginning with the actual astronomical data have errors both in “X” and “Y”, you need to use a function called “linfitex” (you can see more info from the mpfit website) to both define the function you want to fit and calculate the least-square needed in the fitting. Here is an example of the function:

FUCNTION myfunct, p, x=x, y=y \$
     , sigma_x=sigma_x, sigma_y=sigma_y, _EXTRA=extra
  a = p[0]
  b = p[1]
  f = a*alog10(X)+b
;#########calculate the derivative of the function#########
;##put into resid as (y-f)/sqrt(sigma_y^2+(y'*sigma_x)^2)## 
  resid = (y - f)/sqrt(sigma_y^2 + (a/(x*alog(10))*sigma_x)^2) 
RETURN, resid

In the code above I just defined a function as

    \[y=a \log_{10} (b)+c\]

Then you can read your data. and when you want to fit the data, just using mpfit as:

p      = [1,1,1]
p[1]   = 1         
p[2]   = 2
p[3]   = 3 ;initial guess of the parameters
param  = MPFIT('MYFUNCT', params,\$
      FUNCTARGS={X:x,Y:y,SIGMA_X:x_er,SIGMA_Y:y_er}, PERROR=p_er)

The FUNCTARGS is the place where your data go into the function. And then you can find the results of the parameters from the fitting by print “param” and also the 1σ error from printing p_er. If you only have the error on “Y”, you can use a function like this:

FUNCTION myfunct, p, x=x, y=y, er=er
  a = p[0]     
  b = p[1]
  c = p[2]
  f = a*alog10(b)+c
RETURN, (y-f)/er
END

And then using the same way described above to call “mpfit” and get the results and errors.

Here is also a PYTHON version of this procedure.