Multivariate Nonlinear Regression

Specific Domains Optimization (Mathematical)
1

Greetings!

I am trying to perform a multivariate nonlinear regression (2D curve fitting). Let’s assume my data has the following form:

function multimodel(x, p)
    a = x[1]
    t = x[2]
    α = p[1]
    τ = p[2]
    (α .* a) * (τ .* (t.^2)')
end

t = -1.0:0.1:1.0 |> collect
a =  0.0:0.1:1.0 |> collect
p = [0.9,1.1]
x = [a, t]

y = multimodel(x, p) + 0.1 * rand(length(a), length(t))

using PyPlot
mesh(t, a, y)

Which gives the following noisy data:
download-17

I tried fitting using LsqFit:

using LsqFit
fit = curve_fit(multimodel, x, y, [1.0, 1.0])

which results in

MethodError: no method matching isinf(::Array{Float64,1})

I have to admit that that’s not how the LsqFit documentation says how their multivariate regression works. In there it states:

There’s nothing inherently different if there are more than one variable entering the problem. We just need to specify the columns appropriately in our model specification:
@. multimodel(x, p) = p[1]*exp(-x[:, 1]*p[2]+x[:, 2]*p[3])

But I don’t quite understand that because the arrays of the two independent variables must have the same length? And I don’t get a 2D array out of that function.

I also tried modelling my problem with JuMP but I didn’t get very far there.

Does anyone have some experience with this kind of curve fitting?

2

I’m not sure what you’re doing exactly.

Do you want to do non-linear least squares, or do you want to interpolate?

3

it seems that wrong argument is passed to isinf by curve_fit somewhere along the way
try to overload isinf to cover isinf(::Array{Float64,1}) case and see it you get something reasonable out of the fit

4

disregard my suggestion, no use… I’ve tried to overload isinf and isnan but ended up with ERROR: StackOverflowError: with no details

5

The plan is to do nonlinear least squares fitting.

6

In the readme of LsqFit they pass a 2d array as x. I think you passed an array of arrays.

@. multimodel(x, p) = p[1]*exp(-x[:, 1]*p[2]+x[:, 2]*p[3])

7

This is where I don’t get the readme. If I pass a 2D array the first thing is that my independent variable arrays have to have the same length which they don’t have. And secondly that the multimodel function from the readme returns a 1D array instead of a 2D array.

8

maybe try something like this

function multimodel(x, p)
    a = x[:,1]
    t = x[:,2]
    α = p[1]
    τ = p[2]
    @. α * a * τ * t^2
end

a = []
t = []
trange = -1.0:0.1:1.0
arange = 0.0:0.1:1.0
for _t in trange , _a in arange
    push!(a,_a)
    push!(t,_t)
end

p = [0.9,1.1]
x = [a t]

y = multimodel(x, p) + 0.1 * rand(size(a,1))

using Plots
Plots.plot(a,t,y,seriestype=:surface)
2 Likes
9

Yes! This works! That really helps a lot. Thank you! :star_struck:

10

Okay, now I see what you wanted to do. Do you realize why the suggestion worked?

To use LsqFit you need to provide a number (N) of observations. In your case, the observations are chosen by you: you want to sample a the full tensor grid defined by two (in this case) vectors. Then you have to consider each combination/point as one observation. This obviously means that some as are going to be repeated for some ts and vice versa. This is the part I didn’t understand when you said “why do a and t have to be of the same length” - well because you need N pairs!