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

我一直找不到在numpy或python中获取系数误差的任何内置方法。我写了一个简单的工具，参考了John Taylor的《误差分析导论》第8.5和8.6节。也许这个工具能满足你的需求（注意，默认返回的是方差，而不是标准差）。在提供的例子中，由于显著的协方差，你可能会得到很大的误差。

def leastSquares(xMat, yMat):
'''
Purpose
-------
Perform least squares using the procedure outlined in 8.5 and 8.6 of Taylor, solving
matrix equation X a = Y

Examples
--------
>>> from scipy import matrix
>>> xMat = matrix([[  1,   5,  25],
                   [  1,   7,  49],
                   [  1,   9,  81],
                   [  1,  11, 121]])
>>> # matrix has rows of format [constant, x, x^2]
>>> yMat = matrix([[142],
                   [168],
                   [211],
                   [251]])
>>> a, varCoef, yRes = leastSquares(xMat, yMat)
>>> # a is a column matrix, holding the three coefficients a, b, c, corresponding to
>>> # the equation a + b*x + c*x^2

Returns
-------
a: matrix
    best fit coefficients
varCoef: matrix
    variance of derived coefficents
yRes: matrix
    y-residuals of fit 
'''
xMatSize = xMat.shape
numMeas = xMatSize[0]
numVars = xMatSize[1]

xxMat = xMat.T * xMat
xyMat = xMat.T * yMat
xxMatI = xxMat.I

aMat = xxMatI * xyMat
yAvgMat = xMat * aMat
yRes = yMat - yAvgMat

var = (yRes.T * yRes) / (numMeas - numVars)
varCoef = xxMatI.diagonal() * var[0, 0]

return aMat, varCoef, yRes

Answer 2

我觉得你可以使用 curve_fit 这个函数，它在 scipy.optimize 里可以找到（文档链接）。这里有个基本的使用示例：

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 会给出：

popt 的估计协方差。对角线提供了参数估计的方差。

这样你就可以计算出系数的误差估计。要得到标准差，你可以对方差开平方。

现在你有了系数的误差，但这只是基于 ydata 和拟合之间的偏差。如果你还想考虑 ydata 本身的误差，curve_fit 函数提供了 sigma 参数：

sigma : None 或 N 长度的序列

如果不是 None，它表示 ydata 的标准差。如果提供了这个向量，它将作为最小二乘问题中的权重。

这是一个完整的示例：

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)