Curve fitting with LsqFit

New to Julia
1

Dear All,

My aim is to fit some real data with a pretty complex curve representing a superposition of dumping, oscillating and other simple ‘sub-curves’. I use LsqFit. Here is the code example:

using LsqFit
  
model(t, p) = p[1] * exp.(-p[2] * t)
    model1(t, p) = p[1] * exp.(-p[2] * t) + p[3] * exp.(-p[4] * t) * cos.(p[5] * t + p[6])
    # x data
    tdata = range(0,stop=10,length=20)
    # func data
    ydata = model(tdata, [1.0 2.0]) + 0.01*randn(length(tdata))
    p0 = [0.5, 0.5 ,0.5, 0.5, 0.5, 0.5]
    fit = curve_fit(model1, tdata, ydata, p0)
    param = fit.param

The example with model works just fine, but if now change it to something i want, e.g., model1 - i get something like this:

julia>     fit = curve_fit(model1, tdata, ydata, p0)
ERROR: MethodError: no method matching +(::Array{Float64,1}, ::Float64)
Closest candidates are:
  +(::Any, ::Any, ::Any, ::Any...) at operators.jl:502
  +(::Bool, ::T<:AbstractFloat) where T<:AbstractFloat at bool.jl:112
  +(::Float64, ::Float64) at float.jl:395

Many thanks in advance for any tips/solutions!

2

You just need some more .s:

model1(t, p) = p[1] * exp.(-p[2] * t) .+ p[3] * exp.(-p[4] * t) .* cos.(p[5] * t .+ p[6])

or

model1(t, p) = @.(p[1] * exp(-p[2] * t) + p[3] * exp(-p[4] * t) * cos(p[5] * t + p[6]))
1 Like
3

Thank you, Sebastian!
It works now, but could you, please, however briefly, explain why it works (in both cases preferably).
At first glance my model1 is basically the same as model in examples…
All the best,
K.

4

Take a look at the docs about vectorized functions.

1 Like
5

Thank you! It is much clearer now.
Yet, there is a new problem. Now i am trying to use a complex function with ternary operator:

model1(t, p) = @.(p[1] * exp(-p[2] * t) + p[3] * exp(-p[4] * t) * cos(p[5] * t + p[6])+ p[7] * (t > p[8] ? 0 : (1 - t / p[8])^2))

The error i get is the following:

ERROR: TypeError: non-boolean (BitArray{1}) used in boolean context
Stacktrace:
 [1] model1(::Array{Int64,1}, ::Array{Float64,1}) at .\REPL[184]:1
 [2] top-level scope at none:0

Many thanks in advance for tips.

6

The ternary operator doesn’t particaipate in broadcasting, but you can just use ifelse instead:

model1(t, p) = @.(p[1] * exp(-p[2] * t) + p[3] * exp(-p[4] * t) * cos(p[5] * t + p[6])+ p[7] * ifelse(t > p[8], 0, (1 - t / p[8])^2))

This does have some implications you should be aware of:

help?> ifelse
search: ifelse

  ifelse(condition::Bool, x, y)

  Return x if condition is true, otherwise return y. This differs from ? or if in that it is an ordinary function, so all the arguments are evaluated first. In some cases, using ifelse instead of an if statement
  can eliminate the branch in generated code and provide higher performance in tight loops.
3 Likes