Calculate root mean square of 3D deviation after surface fitting in python

My goal is to determine the 3D deviation (and its RMS) between a set of 3D data points and a fitted paraboloid in Python.

Starting from this: Paraboloid (3D parabola) surface fitting python, I can compute the RMS. If I understand correctly, the error and the RMS are computed along the Z-axis. Is it right?

I tried (without success) to determine the 3D deviation and RMS between the fitted surface and the data points, but I cannot get their. Does anyone have some advices to solve this, please?

import numpy as np
from scipy.optimize import curve_fit

# Initial guess parameters
p0 = [1.5,0.4,1.5,0.4,1]

# INPUT DATA
x = [0.4,0.165,0.165,0.585,0.585]
y = [.45, .22, .63, .22, .63]
z = np.array([1, .99, .98,.97,.96])

# FIT
def paraBolEqn(data,a,b,c,d,e):
    x,y = data
    return -(((x-b)/a)**2+((y-d)/c)**2)+e

data = np.vstack((x,y))
popt, _ = curve_fit(paraBolEqn,data,z,p0)

# Deviation and RMS along Z axis
modelPredictions = paraBolEqn(data, *popt) 
absError = modelPredictions - z
RMSE = np.sqrt(np.mean(np.square(absError))) # Root Mean Squared Error along Z axis
print('RMSE (along Z axis):', RMSE)

# Deviation and RMS in 3D
# ??

Copyright Notice：Content Author:「T_by_PW」，Reproduced under the CC 4.0 BY-SA copyright license with a link to the original source and this disclaimer.
Link to original article：https://stackoverflow.com/questions/59919008/calculate-root-mean-square-of-3d-deviation-after-surface-fitting-in-python