Fits now include xerrors

Project.toml

@@ -9,5 +9,6 @@ BDIO = "375f315e-f2c4-11e9-2ef9-134f02f79e27"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

src/juobs.jl

 module juobs
-using ADerrors, PyPlot, LaTeXStrings, LinearAlgebra, LsqFit
+using ADerrors, PyPlot, LaTeXStrings, LinearAlgebra, LsqFit, Optim
 import Statistics: mean
 
 include("juobs_types.jl")

src/juobs_tools.jl

@@ -249,7 +249,7 @@ The method return an array upar with the best fit parameters with their errors.
 @. model(x,p) = p[1] + p[2] * exp(-(p[3]-p[1])*x)
 fit_routine(model, ydata, param=3)
 """
-function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing )
+function fit_routine(model::Function, xdata::Array{<:Real}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
     uwerr.(ydata)
     yval = value.(ydata)
     yer = err.(ydata)
@@ -268,6 +268,65 @@ function gen_chisq(f::Function, x::Array{<:Real}, err::Vector{Float64})
     chisq(par, dat) = sum((dat .- f(x,par)).^2 ./err.^2)
     return chisq
 end
+
+#TODO: COMBINED FITS, COVARIANCE X-Y
+function fit_routine(model::Function, xdata::Array{uwreal}, ydata::Array{uwreal}, param::Int64=3; wpm::Union{Dict{Int64,Vector{Float64}},Dict{String,Vector{Float64}}, Nothing}=nothing)
+    uwerr.(ydata)
+    yval = value.(ydata)
+    yer = err.(ydata)
+    uwerr.(xdata)
+    xval = value.(xdata)
+    xer = err.(xdata)
+
+    dat = vcat(xval, yval)
+    ddat = vcat(xer, yer)
+    data = vcat(xdata, ydata)
+
+    #guess
+    fit = curve_fit(model, xval, yval, 1.0 ./ yer.^2, fill(0.5, param))
+
+    chisq_full(p, d) = get_chi2(model, d, ddat, p)
+    min_fun(t) = chisq_full(t, dat)
+    sol = optimize(min_fun, vcat(fit.param, xval), method=LBFGS())
+    
+    (upar, chi2_exp) = isnothing(wpm) ? fit_error(chisq_full, Optim.minimizer(sol), data) : fit_error(chisq_full, Optim.minimizer(sol), data, wpm)
+
+    for i = 1:length(upar)
+        uwerr(upar[i])
+        print("\n Fit parameter: ", i, ": ")
+        details(upar[i])
+    end
+    println("Chisq / chiexp: ", sol.minimum, " / ", chi2_exp, " (dof: ", dof(fit),")")
+    return upar
+    
+end
+function get_chi2(f::Function, x, dx, y, dy, par) #full
+    chi2 = 0.0
+    Ndata = length(x)
+    Npar = length(par) - Ndata
+    p = par[1:Npar]
+
+    for k = 1:Ndata
+        xx = par[Npar + k]
+        yy = f(xx, p)
+
+        delta = [x[k] - xx, y[k] - yy]
+        C = [[1/dx[k]^2, 0.0] [0.0, 1/dy[k]^2]]
+        chi2 += delta' * C * delta
+    end
+    return chi2
+end
+
+function get_chi2(f::Function, data, ddata, par) #full
+    Ndata = div(length(data), 2)
+
+    x = data[1:Ndata]
+    y = data[Ndata+1:end]
+    dx = ddata[1:Ndata]
+    dy = ddata[Ndata+1:end]
+    
+    return get_chi2(f, x, dx, y, dy, par)
+end
 #=
 using LsqFit
 @doc raw"""