Hi, this is a off-branch topic from this other thread of mine about multiple linear regression. This is not the solution to the original question, but I thought it’s useful, so I’m making a recap here.
Turned out, linear regression is able to fit a paraboloid, but only if the vertex is at (0; 0). In case you wanted to estimate also the vertex as a parameter, like (x_0, y_0), the model becomes non linear. However, @ayushpatnaikgit provided a solution for it. I report it here for convenience:
using DataFrames
using Optim
x = range(-50, 50, length=101);
y = range(-70, 50, length=51);
domain = Iterators.product(x,y) |> collect
β = [1.0, 2.0, 10.0, 20.0]
f((x,y)) = (x, y, (β[1]*(x - β[3])^2 + β[2] * (y - β[4])^2 + 10*randn()))
points = vec(map(f, domain))
df = DataFrame(NamedTuple{(:x, :y, :z)}.(points))
function g(β) # since f is already used.
z = df.z; x = df.x; y = df.y
x₀ = β[3]; y₀ = β[4]
ẑ = β[1]*(x .- x₀).^2 + β[2]*(y .- y₀).^2
return sum((z - ẑ).^2)
end
result = optimize(g, [1.0, 2.0, 10.0, 20.0]) # Actual initial conditions
coefficients = Optim.minimizer(result) # looks like it's the right answer.
result = optimize(g, [1.6, 20.0, 0.0, 0.0]) # Some random initial conditions.
coefficients = Optim.minimizer(result) # Still the correct answer.
Graphing both the dataframe and the fitted surface, you get a result similar to this:
Moreover, @pdeffebach provided useful hints on optimizing the code and making it more robust, suggesting this:
julia> z = df.z; x = df.x; y = df.y;
julia> g = let x = x, y = y, z = z
β -> begin
x₀ = β[3]; y₀ = β[4]
ẑ = β[1]*(x .- x₀).^2 + β[2]*(y .- y₀).^2
return sum((z - ẑ).^2)
end
end;
Kudos to both!