NumPy polyfit过0点的拟合问题

Question

假设我有一个带有权重向量wgt的x和y向量。我可以使用np.polyfit拟合一个三次曲线(y = a x^3 + b x^2 + c x + d)，如下所示:

y_fit = np.polyfit(x, y, deg=3, w=wgt)

现在，假设我想进行另一个适配，但这次，我希望适配通过0 (即y = a x^3 + b x^2 + c x，d = 0)，我该如何指定一个特定的系数（即在这种情况下的d）为零？

谢谢

- uday

Answer 1

你可以使用np.linalg.lstsq并手动构建系数矩阵。首先，我将创建示例数据x和y，以及“完美拟合”y0：

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(100)
y0 = 0.07 * x ** 3 + 0.3 * x ** 2 + 1.1 * x
y = y0 + 1000 * np.random.randn(x.shape[0])

现在我将创建一个完整的三次多项式“训练”或“自变量”矩阵，其中包括常数d列。

XX = np.vstack((x ** 3, x ** 2, x, np.ones_like(x))).T

让我们看一下使用这个数据集进行拟合并将其与 polyfit 进行比较时得到的结果：

p_all = np.linalg.lstsq(X_, y)[0]
pp = np.polyfit(x, y, 3)

print np.isclose(pp, p_all).all()
# Returns True

我使用了np.isclose，因为这两个算法确实存在非常小的差异。

你可能会想：'很好，但我还是没有回答问题。从这里开始，强制拟合具有零偏移量与从数组中删除np.ones列是相同的：

p_no_offset = np.linalg.lstsq(XX[:, :-1], y)[0]  # use [0] to just grab the coefs

好的，让我们看一下这个拟合与我们的数据相比如何：

y_fit = np.dot(p_no_offset, XX[:, :-1].T)

plt.plot(x, y0, 'k-', linewidth=3)
plt.plot(x, y_fit, 'y--', linewidth=2)
plt.plot(x, y, 'r.', ms=5)

这将得到以下这个图形： 警告：如果您在数据上使用此方法，而数据实际上没有通过（x，y）=（0,0），您的输出解决方案系数（p）的估计值将会偏差，因为lstsq将试图弥补您的数据中存在偏移量的事实。有点像一个“圆孔方木”的问题。

此外，您还可以通过执行以下操作将数据仅拟合为立方体函数：

p_ = np.linalg.lstsq(X_[:1, :], y)[0]

再次提醒，如果您的数据包含二次、一次或常数项，则三次系数的估计值将存在偏差。在数值算法中，有时候这种情况是有用的，但对于统计学目的而言，包括所有较低阶段是非常重要的。如果测试结果显示低阶项与零没有统计学上的显著差异，那么就没问题了，但为了安全起见，在估计三次项时最好将它们保留。

祝您好运！

- farenorth

---

感谢您提供详细的答案。我还学习了一个新的函数np.isclose，在其他情境下也可能对我有用。 - uday

Answer 2

您可以尝试以下操作：

从scipy导入curve_fit，即：

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

定义曲线拟合函数。在您的情况下，

def fit_func(x, a, b, c):
    # Curve fitting function
    return a * x**3 + b * x**2 + c * x  # d=0 is implied

执行曲线拟合操作，

# Curve fitting
params = curve_fit(fit_func, x, y)
[a, b, c] = params[0]
x_fit = np.linspace(x[0], x[-1], 100)
y_fit = a * x_fit**3 + b * x_fit**2 + c * x_fit

请绘制结果。

plt.plot(x, y, '.r')         # Data
plt.plot(x_fit, y_fit, 'k')  # Fitted curve

它并没有直接回答问题，而是使用了 numpy 的 polyfit 函数通过原点，但它解决了这个问题。

希望有人会发现它有用 :)

- Gerhard

---

不错的解决方案！y_fit 必须从 x_fit 计算，而不是 x。它甚至可以用更 Pythonic 的方式计算，如：y_fit = fit_func(x_fit, *params[0])，其中星号运算符 * 解包数组，但你的解决方案当然非常易于理解。 - Marcelo Bergweiler

谢谢您的更正，@marcelo。确实应该是 x_fit。不好意思犯了这个错误。 - Gerhard