[PSA]: curve_fit from Optim.jl is moving to LsqFit.jl

[Blake Johnson]

The curve fitting functionality in Optim.jl is being moved into its own package: LsqFit.jl:

https://github.com/JuliaOpt/LsqFit.jl

Installable via Pkg.add("LsqFit").

[Robert Feldt]

Can LsqFit.jl be used to fit multivariate model? I tried this but must have missed something (maybe xv and Y need to have same length?):

using LsqFit

# We want to try to use curve_fit to fit the parameters of this multivariate model:

mvmodel(x, p) = begin

  p[1] .* ((x[1,:] .^ p[2]) ./ (x[2,:] .^ p[3]))

end

N = 10

M = 2

X = randn(M, N)

Params = [1.0, 2.0, 2.0] # Actual params that we are looking for

# Dependent var with small error term

Y = mvmodel(X, Params)[:] + 0.01*randn(N)

# We need a vector to submit to curve_fit:

xv = X[:]

# Reshape to get back the original matrix.

orig_matrix(xv, M) = reshape(xv, (M, int(length(xv)/M)))

# The model function we supply to curve_fit needs to take a vector...

model(xv, p) = mvmodel(orig_matrix(xv, M), p)

fit = curve_fit(model, xv, Y, [0.5, 0.5, 0.5])

[Robert Feldt]

Actually this might be a bug. I really don't understand this behavior (latest 0.3.0-rc3):

feldt:~/tmp$ julia03

               _

   _       _ _(_)_     |  A fresh approach to technical computing

  (_)     | (_) (_)    |  Documentation: http://docs.julialang.org

   _ _   _| |_  __ _   |  Type "help()" for help.

  | | | | | | |/ _` |  |

  | | |_| | | | (_| |  |  Version 0.3.0-rc3 (2014-08-10 02:54 UTC)

 _/ |\__'_|_|_|\__'_|  |  Official http://julialang.org/ release

|__/                   |  x86_64-apple-darwin13.3.0

julia> -0.11292333772901704 ^ 1.0000060554544523

-0.11292184633488864

julia> x = -0.11292333772901704

-0.11292333772901704

julia> p = 1.0000060554544523

1.0000060554544523

julia> x ^ p

ERROR: DomainError

 in ^ at math.jl:252

[Hans W Borchers]

Robert,

your first expression is really

    julia> -(0.11292333772901704 ^ 1.0000060554544523)
    -(0.11292333772901704 ^ 1.0000060554544523)

while the second one changes it to

    julia> (-0.11292333772901704) ^ 1.0000060554544523

    ERROR: DomainError
     in ^ at math.jl:252

which is indeed mathematically undefined.

[Robert Feldt]

Yes, sorry, seems I'm still on vacation... ;) Ruby parses this another way so I confused the two languages.

Actually, by changing my model in the LsqFit example I now also get the multivariate curve fitting to work; code below. I guess in general convergence might be a problem unless the initial guess is relatively close.

If there are faster/simpler ways to accomplish this (with LsqFit) please share.

Regards, Robert

using LsqFit

# Try to use curve_fit for this multivariate model:

mvmodel(x, p) = begin

  n = abs(x[1,:]) .^ p[2]

  d = abs(x[2,:]) .^ p[3]

  p[1] .* (n ./ d)

end

N = 10

M = 2

X = randn(M, N)

Params = [1.0, 2.0, 2.0] # Actual params that we are looking for

errors = 0.01*randn(N)

# Dependent var with small error term

Y = mvmodel(X, Params)[:] .+ errors

# We need a vector to submit to curve_fit:

xv = X[:] # reshape(X, (N*M, 1))[:]

# Reshape to get back the original matrix.

orig_matrix(xv, M) = reshape(xv, (M, int(length(xv)/M)))

# The model function we supply to curve_fit needs to take a vector...

model(xv, p) = mvmodel(orig_matrix(xv, M), p)[:]

fit = curve_fit(model, xv, Y, [0.5, 1.0, 1.0])

julia> fit.param

3-element Array{Float64,1}:

 0.999694

 1.99672

 1.99926

--

Best regards,

/Robert Feldt
--
Tech. Dr. (PhD), Professor of Software Engineering
Blekinge Institute of Technology, Software Engineering Research Lab, and

Chalmers, Software Engineering Dept

Explanea.com - Igniting your Software innovation
robert.feldt (a) bth.se    or    robert.feldt (a) chalmers.se    or    robert.feldt (a) gmail.com
Mobile phone: +46 (0) 733 580 580
http://www.robertfeldt.net

[Robert Feldt]

So as not to confuse others I just wanted to clarify that curve_fit does not care what you supply as x as long as what is returned from calling your model is a vector of same length as what you supply as y. Example:

using LsqFit

mvmodel(x, p) = begin

  n = abs(x[1,:]) .^ p[2]

  d = abs(x[2,:]) .^ p[3]

  (p[1] .* (n ./ d))[:]

end

N = 25

M = 100

X = randn(M, N)

Params = [1.0, 2.0, 2.0]

errors = 0.01*randn(N)

Y = mvmodel(X, Params) .+ errors

# curve_fit does not care about what is in X as long as what is returned

# when calling model(X, p) is vector of float of same size as Y.

fit = curve_fit(mvmodel, X, Y, rand(3))

julia> Params .- fit.param

3-element Array{Float64,1}:

  1.91519e-5

 -3.94951e-7

 -5.9904e-6