Question about LsqFit

Specific Domains Optimization (Mathematical)
,
1

Dear all,

When I using LsqFit.curve_fit, I encountered an error.

the codes are:

using LsqFit

xdata = [[42.7908, 240.168, 32.799],
 [43.262, 346.383, -0.296673],
 [43.842, 332.0, -32.1631]   ,
 [44.184, 324.176, 35.1956]  ,
 [44.8529, 270.56, -21.4661] ,
 [44.8649, 184.976, 11.9061] ,
 [45.1269, 328.77, -10.8362] ,
 [45.4122, 248.579, 37.3079] ,
 [45.4975, 156.928, -76.5363],
 [45.8181, 347.873, -9.41383]];


ydata = [42.37397477812974,
 43.644649744028534,
 43.54054864852959,
 43.93823563647773,
 44.362474498426906,
 45.063235860732235,
 44.30769385086685,
 44.3473803820133,
 45.76719596019439,
 45.74372210459567];

function angular_sepration(l2,b2)
    l1,b1 = 310,-10
    l1 = deg2rad(l1)
    b1 = deg2rad(b1)
    l2 = deg2rad(l2)
    b2 = deg2rad(b2)
    b11 = (pi/2) - b1
    b22 = (pi/2) - b2
    sep = acos( cos(b11) * cos(b22) + sin(b11) * sin(b22) * cos(l1 - l2))
    return rad2deg(sep)
end;

function fit_func(xdatas,p)
    fit_d,fit_l,fit_b = p
    
    return_list = []
    
    for parameter in xdatas
        mu = parameter[1]
        l = parameter[2] # degree
        b = parameter[3] # degree

        angular = angular_sepration(l,b)
        cos_theta = cosd(angular)

        mu_diople = mu * (1 - fit_d * cos_theta)
        append!(return_list,mu_diople)
    end
    return return_list
end

p0 = [0.001,310,-10]
fit = curve_fit(fit_func,xdata,ydata,p0)

The above code runs well. the fit output is

LsqFit.LsqFitResult{Float64,1}(7, [0.00938423, 310.0, -10.0], [0.33999, -0.704888, -0.0545982, -0.0362614, 0.165758, 0.0495903, 0.417554, 0.949986, -0.254515, -0.267582], [-8.18767 0.0 0.0; -34.3382 0.0 0.0; … ; 1.61768 0.0 0.0; -36.4398 0.0 0.0], true, Float64[])

However, when I try to get the standard error,

standard_error(fit)

it raises an error:

LinearAlgebra.LAPACKException(2)

Stacktrace:
 [1] chklapackerror at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/lapack.jl:38 [inlined]
 [2] trtrs!(::Char, ::Char, ::Char, ::Array{Float64,2}, ::Array{Float64,2}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/lapack.jl:3348
 [3] ldiv! at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/triangular.jl:583 [inlined]
 [4] inv at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/triangular.jl:630 [inlined]
 [5] inv(::Array{Float64,2}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/dense.jl:710
 [6] estimate_covar(::LsqFit.LsqFitResult{Float64,1}) at /opt/julia/packages/LsqFit/IsCnQ/src/curve_fit.jl:134
 [7] #standard_error#19(::Float64, ::Int64, ::Function, ::LsqFit.LsqFitResult{Float64,1}) at /opt/julia/packages/LsqFit/IsCnQ/src/curve_fit.jl:150
 [8] standard_error(::LsqFit.LsqFitResult{Float64,1}) at /opt/julia/packages/LsqFit/IsCnQ/src/curve_fit.jl:150
 [9] top-level scope at In[35]:1

I checked the code of standard_error and found this error was caused by

J = fit.jacobian
Q, R = qr(J)
Rinv = inv(R)

where R is singular.

When I using Python (scipy.optimize.curve_fit), scipy.optimize.curve_fit works well for model fit_func.

Sorry for my poor english, could you help me to solve this problem?

Thanks.

1 Like
2

Thanks for providing code and example data. I’ll have a look, this might be fixed on the master branch (in-development version of LsqFit).

Edit: Let me say something up front though. There’s a known issue with the standard errors in LsqFit. They were supposed to be fixed, but were not for various reasons. I’m a bit strapped for time so I can’t promise anything. It’d help if you could provide the full example for scipy as well. Thanks

Edit2: Hmm…

julia> J =fit.jacobian
10×3 Array{Float64,2}:
  -8.18767  0.0  0.0
 -34.3382   0.0  0.0
 -37.9412   0.0  0.0
 -30.0531   0.0  0.0
 -34.5972   0.0  0.0
  26.4194   0.0  0.0
 -42.8008   0.0  0.0
 -12.2369   0.0  0.0
   1.61768  0.0  0.0
 -36.4398   0.0  0.0

Your fit_func doesn’t vary with the last two parameters, that can’t be. Either you need to fix them at their values, or they have to be able to vary.