scipy.optimize workd.在这个代码中，估计是线性函数，但它可能是更好的一个。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

X1=[3,3,3.1,3.1,4.2,5.2,6.3,2.3,7.4,8.4,5.4,3.4,3.4,3.4]
X2=[12.1,12.7,18.5,18.3,18.4,18.6,24.2,24.4,24.3,24.5,30.9,30.7,30.3,30.4]
X3=[0.3,9.2,0.3,9.4,0.1,9.8,0.4,9.3,0.7,9.7,18.3,27.4,0.6,9.44]
Y=[-5.890,-5.894,2.888,-3.8706,2.1516,-2.7334,1.4723,-2.1049,0.9167,-1.7281,-2.091,-6.7394,0.8777,-1.7046]

def fitFunc(x, a, b, c, d):
    return a + b*x[0] + c*x[1] + d*x[2]

fitParams, fitCovariances = curve_fit(fitFunc, [X1, X2, X3], Y)
print(' fit coefficients:\n', fitParams)
# fit coefficients:
#  [-6.11934208  0.21643939  0.26186705 -0.33794415]

然后用fitParams[0] + fitParams[1] * x1 + fitParams[2] * x2 + fitParams[3] * x3估计y。

# get single y
def estimate(x1, x2, x3):
    return fitParams[0] + fitParams[1] * x1 + fitParams[2] * x2 + fitParams[3] * x3

将结果与原始y进行比较。

Y_estimated = [estimate(X1[i], X2[i], X3[i]) for i in range(len(X1))]

fig, ax = plt.subplots()
ax.scatter(Y, Y_estimated)

lims = [
    np.min([ax.get_xlim(), ax.get_ylim()]),  # min of both axes
    np.max([ax.get_xlim(), ax.get_ylim()]),  # max of both axes
]

ax.set_xlabel('Y')
ax.set_ylabel('Y_estimated')
ax.plot(lims, lims, 'k-', alpha=0.75, zorder=0)
ax.set_aspect('equal')

参考scipy、stackoverflow-multifit、stackoverflow-plot xy

通常，您可以使用PolynomialFeatures

import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from sklearn import linear_model

x,y,z = np.random.uniform(0,10,(3,100))
values = 10 + 1.5*x -2*y + 3*z
values = values + np.random.normal( 0, 1, values.shape )

poly = PolynomialFeatures(degree=1) 
X_ = poly.fit_transform( np.array( [(x,y,z)] ).T ) #calculate each needed term (x,y,z eventually x², y²...)

clf = linear_model.LinearRegression()
clf.fit(X_, values )
coef_result = clf.coef_  #degree 2: offset, ax, ay,  az .. axax axay axaz.... 
    
estimate_fit = clf.predict( X_ )

#Print
print('PolynomialFeature, deg=1')
print( np.array( [ poly.get_feature_names_out() , np.round(slope_fov_coef,8) ]).T )
print('RMSE is:', np.sqrt( np.nanmean( (slope_estimate_fit - slope_sxsydiff_z)**2 ) ) )

PolynomialFeature，deg=1 [<$'1' 0.0]<$'x 0' 1.44528075]<$'x1' -1.95544098]<$'x2' 3.02024807]] RMSE为：0.956660546736911
请注意，如果需要，可以通过从X_ array中删除某些元素来删除它们。（例如，如果你不想拟合x²的原因.）这将不会使用样条，但它是解决你心目中的那种多项式的最佳线性回归解决方案。