# Quadratic (paraboloid) model fitting

**URL:** <https://discourse.julialang.org/t/quadratic-paraboloid-model-fitting/88007>\
**Category:** Statistics\
**Tags:** regression, fit, curve-fitting\
**Created:** [September 29, 2022, 7:37pm UTC](https://discourse.julialang.org/t/quadratic-paraboloid-model-fitting/88007 "2022-09-29T19:37:47Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![bertulli](https://avatars.discourse-cdn.com/v4/letter/b/ecc23a/32.png) [@bertulli](https://discourse.julialang.org/u/bertulli)\
**Post date:** [September 29, 2022, 7:37pm UTC](https://discourse.julialang.org/t/quadratic-paraboloid-model-fitting/88007/1 "2022-09-29T19:37:47Z")

</div>

Hi, this is a off-branch topic from [this other thread of mine about multiple linear regression](https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692). 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](https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692/12). I report it here for convenience:

```julia
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:  
 ![quadratic_non_linear_plot](https://global.discourse-cdn.com/julialang/original/3X/8/9/89aa4bcb79419501f20dd7e764e35e8d9d7a7d68.png)

Moreover, @pdeffebach provided [useful hints on optimizing the code and making it more robust](https://discourse.julialang.org/t/how-can-i-tell-julia-what-are-the-parameters-in-a-linear-model/87692/15), suggesting this:

```julia
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!