Python - 计算有错误的趋势线

Question

所以我将一些数据存储为两个列表，并使用它们绘制

plot(datasetx, datasety)

然后我设置了一条趋势线

trend = polyfit(datasetx, datasety)
trendx = []
trendy = []

for a in range(datasetx[0], (datasetx[-1]+1)):
    trendx.append(a)
    trendy.append(trend[0]*a**2 + trend[1]*a + trend[2])

plot(trendx, trendy)

但是我有第三个数据列表，这是原始数据集中的错误。我很擅长绘制误差线，但我不知道如何使用它来找出多项式趋势线系数中的误差。

所以说我的趋势线是 5x^2 + 3x + 4 = y，5、3 和 4 值需要有某种错误。

有没有使用 NumPy 的工具可以帮我计算这个？

Answer 1

我认为您可以使用 scipy.optimize ( documentation ) 的函数 curve_fit。用法的基本示例:

import numpy as np
from scipy.optimize import curve_fit

def func(x, a, b, c):
    return a*x**2 + b*x + c

x = np.linspace(0,4,50)
y = func(x, 5, 3, 4)
yn = y + 0.2*np.random.normal(size=len(x))

popt, pcov = curve_fit(func, x, yn)

根据文档，pcov 给出:

The estimated covariance of popt. The diagonals provide the variance of the parameter estimate.

因此，通过这种方式，您可以计算系数的误差估计值。要获得标准偏差，您可以对方差求平方根。

现在您在系数上有一个错误，但它仅基于 y 数据和拟合之间的偏差。如果您还想考虑 ydata 本身的错误，curve_fit 函数提供了 sigma 参数:

sigma : None or N-length sequence

If not None, it represents the standard-deviation of ydata. This vector, if given, will be used as weights in the least-squares problem.

一个完整的例子:

import numpy as np
from scipy.optimize import curve_fit

def func(x, a, b, c):
    return a*x**2 + b*x + c

x = np.linspace(0,4,20)
y = func(x, 5, 3, 4)
# generate noisy ydata
yn = y + 0.2 * y * np.random.normal(size=len(x))
# generate error on ydata
y_sigma = 0.2 * y * np.random.normal(size=len(x))

popt, pcov = curve_fit(func, x, yn, sigma = y_sigma)

# plot
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111)
ax.errorbar(x, yn, yerr = y_sigma, fmt = 'o')
ax.plot(x, np.polyval(popt, x), '-')
ax.text(0.5, 100, r"a = {0:.3f} +/- {1:.3f}".format(popt[0], pcov[0,0]**0.5))
ax.text(0.5, 90, r"b = {0:.3f} +/- {1:.3f}".format(popt[1], pcov[1,1]**0.5))
ax.text(0.5, 80, r"c = {0:.3f} +/- {1:.3f}".format(popt[2], pcov[2,2]**0.5))
ax.grid()
plt.show()

---

然后是其他，关于使用 numpy 数组。使用 numpy 的主要优点之一是可以避免 for 循环，因为对数组的操作是逐元素应用的。因此，您示例中的 for 循环也可以按以下方式完成:

trendx = arange(datasetx[0], (datasetx[-1]+1))
trendy = trend[0]*trendx**2 + trend[1]*trendx + trend[2]

在我使用 arange 而不是 range 的地方，因为它返回一个 numpy 数组而不是列表。在这种情况下，您还可以使用 numpy 函数 polyval:

trendy = polyval(trend, trendx)