Julia: using CurveFit package non linear

Question

why the code below does not work?

xa = [0  0.200000000000000 0.400000000000000  1.00000000000000  1.60000000000000  1.80000000000000  2.00000000000000  2.60000000000000  2.80000000000000  3.00000000000000  3.80000000000000  4.80000000000000  5.00000000000000  5.20000000000000  6.00000000000000  6.20000000000000  7.40000000000000  7.60000000000000  7.80000000000000  8.60000000000000  8.80000000000000  9.00000000000000  9.20000000000000  9.40000000000000  10.0000000000000  10.6000000000000  10.8000000000000  11.2000000000000  11.6000000000000  11.8000000000000  12.2000000000000  12.4000000000000];

ya = [-0.183440428023042  -0.131101157495126  0.0268875670852843 0.300000000120000  0.579048247883555  0.852605831272159  0.935180993484717  1.13328608090532  1.26893326843583  1.10202945535186  1.09201137189664  1.14279083803453  0.811302535321072  0.909735376251797  0.417067545528244  0.460107770989798  -0.516307074859654  -0.333994077331822  -0.504124744955962  -0.945794320817293  -0.915934553082780  -0.975458595671737  -1.09943707404275  -1.11254211607374  -1.29739980589100  -1.23440439602665  -0.953807504156356  -1.12240274852172  -0.609284630192522  -0.592560286759450  -0.402521296049042  -0.510090363150962];

x0 = vec(xa) 
y0 = vec(ya)
fun(x,a) = a[1].*sin(a[2].*x - a[3])
a0 = [1,2,3] 
eps = 0.000001 
maxiter=200 
coefs, converged, iter = CurveFit.nonlinear_fit(x0 , fun , a0 , eps, maxiter ) 
y0b = fit(x0) 
Winston.plot(x0, y0, "ob", x0, y0b, "r-", linewidth=3)

Error: LoadError: MethodError: convert has no method matching convert(::Type{Float64}, ::Array{Float64,1}) This may have arisen from a call to the constructor Float64(...), since type constructors fall back to convert methods. Closest candidates are: call{T}(::Type{T}, ::Any) convert(::Type{Float64}, !Matched::Int8) convert(::Type{Float64}, !Matched::Int16)

while loading In[269], in expression starting on line 8 in nonlinear_fit at /home/jmarcellopereira/.julia/v0.4/CurveFit/src/nonlinfit.jl:75

Answer 1

The fun function has to return a residual value r of type Float64, calculated at each iteration of the data, as follows:

r = y - fun(x, coefs)

so your function y=a1*sin(x*a2-a3) will be defined as:

fun(x,a) = x[2]-a[1]*sin(a[2]*x[1] - a[3])

Where:

  x[2] is a value of 'y' vector  
  x[1] is a value of 'x' vector  
  a[...] is the set of parameters

The fun function has to return a single Float64, so the operators can't be 'dot version' (.*).

By calling the nonlinear_fit function, the first parameter must be an array Nx2, with the first column containing N values of x and the second, containing N values of y, so you must concatenate the two vectors x and y in a two columns array:

xy = [x y]

and finally, call the function:

coefs, converged, iter = CurveFit.nonlinear_fit(xy , fun , a0 , eps, maxiter ) 

Answering to your comment about the returned coefficients are not correct:

The y = 1 * sin (x * a2-a3) is a harmonic function, so the coefficients returning from the function call, depend heavily on the parameter a0 ("initial guess for each fitting parameter") you will send as the third parameter (with maxiter=200_000):

a0=[1.5, 1.5, 1.0]  
coefficients: [0.2616335317043578, 1.1471991302529982,0.7048665905560775]

a0=[100.,100.,100.]  
coefficients: [-0.4077952060368059, 90.52328921205392, 96.75331155303707]

a0=[1.2, 0.5, 0.5]  
coefficients: [1.192007321713507, 0.49426296880933257, 0.19863645732313934]

I think the results you're getting are harmonics, as the graph:

Where:

blue line:  
f1(xx)=0.2616335317043578*sin(xx*1.1471991302529982-0.7048665905560775)

yellow line:
f2(xx)=1.192007321713507*sin(xx*0.49426296880933257-0.19863645732313934)

pink line:
f3(xx)=-0.4077952060368059*sin(xx*90.52328921205392-96.75331155303707)

blue dots are your initial data.

The graph was generated with Gadfly:

plot(layer(x=x,y=y,Geom.point),layer([f1,f2,f3],0.0, 15.0,Geom.line))

_{tested with Julia Version 0.4.3}

Answer 2

from Doc:

we are trying to fit the relationship fun(x, a) = 0

So, if you want to find elements of a in a way that: for each xi,yi in [x0 y0] => a[1].*sin(a[2].*xi - a[3])==yi, then the right way is:

fun(xy,a) = a[1].*sin(a[2].*xy[1] - a[3])-xy[2];
xy=hcat(x0,y0);
coefs,converged,iter = CurveFit.nonlinear_fit(xy,fun,a0,eps,maxiter);

Answer 3

I found the LsqFit package a bit simpler to use, just define first the model and "fit it" with your data:

using DataFrames, Plots, LsqFit

xa         = [0  0.200000000000000 0.400000000000000  1.00000000000000  1.60000000000000  1.80000000000000  2.00000000000000  2.60000000000000  2.80000000000000  3.00000000000000  3.80000000000000  4.80000000000000  5.00000000000000  5.20000000000000  6.00000000000000  6.20000000000000  7.40000000000000  7.60000000000000  7.80000000000000  8.60000000000000  8.80000000000000  9.00000000000000  9.20000000000000  9.40000000000000  10.0000000000000  10.6000000000000  10.8000000000000  11.2000000000000  11.6000000000000  11.8000000000000  12.2000000000000  12.4000000000000];
ya         = [-0.183440428023042  -0.131101157495126  0.0268875670852843 0.300000000120000  0.579048247883555  0.852605831272159  0.935180993484717  1.13328608090532  1.26893326843583  1.10202945535186  1.09201137189664  1.14279083803453  0.811302535321072  0.909735376251797  0.417067545528244  0.460107770989798  -0.516307074859654  -0.333994077331822  -0.504124744955962  -0.945794320817293  -0.915934553082780  -0.975458595671737  -1.09943707404275  -1.11254211607374  -1.29739980589100  -1.23440439602665  -0.953807504156356  -1.12240274852172  -0.609284630192522  -0.592560286759450  -0.402521296049042  -0.510090363150962];
x0         = vec(xa)
y0         = vec(ya)
xbase      = collect(linspace(minimum(x0),maximum(x0),100))
p0         = [1.2, 0.5, 0.5]                                                      # initial value of parameters
fun(x0, p) = p[1] .* sin.(p[2] .* x0 .- p[3])                                     # definition of the model
fit        = curve_fit(fun,x0,y0,p0)                                              # actual fitting job
yFit       = [fit.param[1] * sin(fit.param[2] * x - fit.param[3]) for x in xbase] # building the fitted values
# Plotting..
scatter(x0, y0, label="obs")
plot!(xbase, yFit, label="fitted")

Note that using LsqFit does not solve the problem of dependency from initial conditions highlighted by Gomiero..