将 3D 数据点拟合到多项式曲面并返回曲面方程答案

【问题标题】：Fitting 3D data points to polynomial surface and getting the surface equation back将 3D 数据点拟合到多项式曲面并返回曲面方程
【发布时间】：2020-08-05 20:01:15
【问题描述】：

我是 Python 3D 拟合和相关优化技术的新手。我试图了解类似的主题并根据最小二乘法找到答案，但我的成功相当有限。

问题：我有一些（大约 400 个）3D 点存储在 np 数组中。

data = np.array([[x1, y1, z1], [x2, y2, z2], [x3, y3, z3], [x4, y4, z4], ... ])

我想将多项式曲面（2 或 3 阶）拟合到这些点，然后获取曲面方程的参数，这样我就可以计算任何给定 (x, y) 值对的 z。例如：

z = (A * x ** 2) + (B * y ** 2) + (C * x * y) + (D * x) + (E * y) + F

我需要得到一个包含 A、B、C、D、E 和 F 参数的输出。有没有一种很好的“Pythonic”方式来做到这一点？

【问题讨论】：

标签： python numpy scipy

【解决方案1】：

curve_fit 接受多维数组作为自变量，但你的函数必须接受相同的：

import numpy as np
from scipy.optimize import curve_fit

data = np.array(
    [[0, 0, 1],
     [1, 1, 2],
     [2, 1, 3],
     [3, 0, 5],
     [4, 0, 2],
     [5, 1, 3],
     [6, 0, 7]]
)

def func(X, A, B, C, D, E, F):
    # unpacking the multi-dim. array column-wise, that's why the transpose
    x, y, z = X.T

    return (A * x ** 2) + (B * y ** 2) + (C * x * y) + (D * x) + (E * y) + F

popt, _ = curve_fit(func, data, data[:,2])

from string import ascii_uppercase
for i, j in zip(popt, ascii_uppercase):
    print(f"{j} = {i:.3f}")

# A = 0.060
# B = 2004.446
# C = -0.700
# D = 0.521
# E = -2003.046
# F = 1.148

请注意，一般情况下，您应该为参数提供初始猜测以获得良好的拟合结果。

【讨论】：

【解决方案2】：

如果您允许 sklearn 依赖，您可以将sklearn.preprocessing.PolynomialFeatures 与sklearn.linear_model.LinearRegression 一起使用：

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

np.random.seed(0)

# Generate fake data 
n,m = 400, 3
data = np.random.randn(n,m)

# Generate polynomial features of desired degree
d = 3
poly = PolynomialFeatures(degree=d, include_bias=False)
X = poly.fit_transform(data[:, :-1]) # X.shape = (400, 9)
y = data[:,-1] # y.shape = (400,)

# Define and fit linear regression 
clf = LinearRegression()
clf.fit(X, y)

# Check results
print(clf.coef_)
[-0.01437971  0.09894586  0.01936384  0.17245758  0.05518938
  0.0239589  -0.09930492 -0.04593238 -0.01588326]

print(clf.intercept_)
-0.12456419670821159

【讨论】：

谢谢！非常有用的一段代码！但是在这种情况下，我需要继续使用 numpy / scipy，因为它们可以使用 Numba 加速（更大项目的一部分）。