I have made minor changes, please have a look:
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.