在没有numpy polyfit的python中拟合二次函数

Posted 2023-03-12

【中文标题】在没有numpy polyfit的python中拟合二次函数

【英文标题】：Fitting a quadratic function in python without numpy polyfit

【发布时间】：2019-10-04 12:12:19

【问题描述】：

我正在尝试将二次函数拟合到某些数据中，并且我正在尝试在不使用 numpy 的 polyfit 函数的情况下做到这一点。

从数学上讲，我尝试关注这个网站https://neutrium.net/mathematics/least-squares-fitting-of-a-polynomial/，但不知何故，我认为我做得不对。如果有人可以帮助我，那就太好了，或者如果您能提出另一种方法，那也很棒。

到目前为止我已经尝试过：

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

ones = np.ones(3)
A = np.array( ((0,1),(1,1),(2,1)))
xfeature = A.T[0]
squaredfeature = A.T[0] ** 2
b = np.array( (1,2,0), ndmin=2 ).T
b = b.reshape(3)

features = np.concatenate((np.vstack(ones), np.vstack(xfeature), np.vstack(squaredfeature)), axis = 1)
featuresc = features.copy()
print(features)
m_det = np.linalg.det(features)
print(m_det)
determinants = []
for i in range(3):
    featuresc.T[i] = b
    print(featuresc)
    det = np.linalg.det(featuresc)
    determinants.append(det)
    print(det)
    featuresc = features.copy()

determinants = determinants / m_det
print(determinants)
plt.scatter(A.T[0],b)
u = np.linspace(0,3,100)
plt.plot(u, u**2*determinants[2] + u*determinants[1] + determinants[0] )
p2 = np.polyfit(A.T[0],b,2)
plt.plot(u, np.polyval(p2,u), 'b--')
plt.show()

如您所见，我的曲线与 nnumpy 的 polyfit 曲线相比不太好。

---

更新： 我检查了我的代码并删除了所有愚蠢的错误，现在它可以工作了，当我尝试将其拟合超过 3 分时，但我不知道如何拟合超过 3 分。

这是新代码：

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

ones = np.ones(3)
A = np.array( ((0,1),(1,1),(2,1)))
xfeature = A.T[0]
squaredfeature = A.T[0] ** 2
b = np.array( (1,2,0), ndmin=2 ).T
b = b.reshape(3)

features = np.concatenate((np.vstack(ones), np.vstack(xfeature), np.vstack(squaredfeature)), axis = 1)
featuresc = features.copy()
print(features)
m_det = np.linalg.det(features)
print(m_det)
determinants = []
for i in range(3):
    featuresc.T[i] = b
    print(featuresc)
    det = np.linalg.det(featuresc)
    determinants.append(det)
    print(det)
    featuresc = features.copy()

determinants = determinants / m_det
print(determinants)
plt.scatter(A.T[0],b)
u = np.linspace(0,3,100)
plt.plot(u, u**2*determinants[2] + u*determinants[1] + determinants[0] )
p2 = np.polyfit(A.T[0],b,2)
plt.plot(u, np.polyval(p2,u), 'r--')
plt.show()

【参考方案1】：

而不是使用克莱默规则，实际上使用最小二乘法求解系统。请记住，只有当您拥有的总点数等于所需的多项式阶数加 1 时，克莱默规则才有效。 如果你没有这个，那么在你试图找到问题的确切解决方案时，克莱默规则将不起作用。如果您有更多点，则该方法不合适，因为我们将创建一个超定方程组。

为了适应更多的点，numpy.linalg.lstsq 会更合适，因为它通过计算向量 x 来解决 Ax = b使用矩阵 A 的欧几里得范数。因此，从特征矩阵的最后一列中删除y的值并求解系数，并使用numpy.linalg.lstsq求解系数：

import numpy as np
import matplotlib.pyplot as plt


ones = np.ones(4)
xfeature = np.asarray([0,1,2,3])
squaredfeature = xfeature ** 2
b = np.asarray([1,2,0,3])

features = np.concatenate((np.vstack(ones),np.vstack(xfeature),np.vstack(squaredfeature)), axis = 1) # Change - remove the y values

determinants = np.linalg.lstsq(features, b)[0] # Change - use least squares
plt.scatter(xfeature,b)
u = np.linspace(0,3,100)
plt.plot(u, u**2*determinants[2] + u*determinants[1] + determinants[0] )
plt.show()

我现在得到了这个图，它与您的图表中的虚线相匹配，也与numpy.polyfit 给您的相匹配：

【讨论】：

哦，哇...我回到我的代码并查看了我的矩阵，心想，如果我正好使用 3 个点，这可能会起作用，瞧，它起作用了。我回到这里阅读您的答案，它准确地说明了我刚刚发现的内容。谢谢哈哈

@AhmadMoussa 啊哈！我的荣幸。我认为如果您尝试添加更多积分，您会遇到麻烦。最小二乘是解决这个问题的方法。很高兴我能帮上忙！

最小二乘拟合的最佳描述，如果你想理解它并使用它，我发现是在 Gene F 的书 Digital Control of Dynamic Systems (3rd Edition) 3rd Edition 中。富兰克林（作者）、J.大卫鲍威尔（作者）、迈克尔 L.沃克曼（作者）

@user50619 很棒的教科书。当我学习数字控制时，这是我的首选书。

@AhmadMoussa 是的！梯度下降找到成本函数的最小值。如果您实际上将问题的矩阵形式扩展为平方项之和，从而为您提供成本函数，您可以计算每个参数的梯度，从而迭代求解系数。成本函数是凸的，因此只有一个唯一的最小值。只要目标是最小化一些成本函数，将多项式拟合到数据的方法并不多。我可以提供给你的一个例子是遗传算法，来自进化学习的方法。