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

我想你可以使用函数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参数:

西格玛:无或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循环,因为数组上的操作应用于elementwise.因此,您的示例中的for循环也可以按如下方式完成:

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

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

trendy = polyval(trend, trendx)

Answer 2

我还没有找到任何方法来获取 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