Using LsqFit in Julia

user_5476728 0
Viewed 471

This link indicates that the analogue of Matlab's lsqnonlin is LsqFit.jl in Julia. However, matlab uses the notation

 x = lsqnonlin(fun,x0)

But Julia uses notation

 fit = curve_fit(model, xdata, ydata, p0)

So what is the analogue of the x in Matlab notation, if I'm using LsqFit in Julia? In Julia command

 fieldnames(fit)

Gives

 5-element Array{Symbol,1}:
  :dof      
  :param    
  :resid    
  :jacobian
  :converged

where converged is a boolean type, and fit.param corresponds to the parameters vector p.

So where is the independent variable x?Want to solve for x in F(x)=0 where x and F are column vectors and F is a nonlinear function.

julia matlab
2 Answers
fit = curve_fit(model, xdata, ydata, p0)

is more akin to a call like

[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqnonlin(___)

in Matlab. The resulting fit object contains some similar solver metadata:

# fit is a composite type (LsqFitResult), with some interesting values:
#   fit.dof: degrees of freedom
#   fit.param: best fit parameters
#   fit.resid: residuals = vector of residuals
#   fit.jacobian: estimated Jacobian at solution

So, fit.param should be the same as x.

user_13462760

Here is an MWE demonstrating the proper usage for finding y(x) that minimizes squares for function F(y(x)) = 0. It is modified in that the p0 "parameters" are now actually the y-values. The x-values are meant to be the underlying grid (so notation F(x)=0 is not correct for curve_fit), and it should be that we are solving F(y(x))=0 as follows.

using LsqFit
using PyPlot

# example model, could be nlin function or diff eq
# we will solve DE y' - y with init conds y'(0) = 1
# by using function F(y) = y' - y as follows
F(xs, ys) = [1; diff(ys)/xs[2];] - ys; # column vector
# note: function F above assumes xs[1] = 0 so that xs[2] = dx

# xs: independent variables, these values never change
dx = 0.02; # make smaller to get more accurate
xs = [0:dx:2;]; # column vector
len = length(xs);

# ys: initial guess for solution function y(x), will be optimized
ys = exp(xs)+rand(len,1)-0.5; # column vector

# modified usage for sum of squares to minimize F
fit = curve_fit(F,xs,0.0*xs[:],ys[:])
# note: original usage was
# fit = curve_fit(model, xdata, ydata, p0)

figure(1)
scatter(xs,ys,marker = "o") # initial noisy guess in blue
plot(xs,fit.param,color = "red") # optimized solution in red
plot(xs,exp(xs),linestyle = "--",color = "green") # exact sol
axis("tight")

display("sum of squares")
sum(abs(F(xs,fit.param)).^2)
user_54767280
Related