numpy - numpy.polyfit: 如何在估计曲线周围获得 1西格玛不确定度?

  显示原文与译文双语对照的内容
0 0

我用 numpy.polyfit 来拟合观察。 polyfit给出了多项式的估计系数,还可以给出估计系数的误差协方差矩阵。 好了,现在我想知道估计估计曲线周围的+/1sigma 不确定度的方法。

我知道MatLab可以做( http://stats.stackexchange.com/questions/56596/finding-uncertainty-in-coefficients-from-polyfit-in-matlab ),但我没有找到一个方法使它在 python 中。

时间: 原作者:

0 0

如果有足够的数据点,可以从 polyfit() 获得参数 cov=True 估计协方差矩阵。 记住,你可以写一个多项式。 p[0]*t**n + p[1]*t**(n-1) +.. . + p[n] 作为矩阵产品 np.dot(tt, p)tt=[t**n, tt*n-1,.. ., 1]t 可以是单个值,也可以是列向量。 由于线性方程的协方差矩阵 C_p,所以值的协方差矩阵是。

所以一个简单的例子


import numpy as np
import matplotlib.pyplot as plt

# sample data:
x = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0])
y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0, -3.0])

n = 3 # degree of polynomial
p, C_p = np.polyfit(x, y, n, cov=True) # C_z is estimated covariance matrix

# Do the interpolation for plotting:
t = np.linspace(-0.5, 6.5, 500)
# Matrix with rows 1, t, t**2,.. .:
TT = np.vstack([t**(n-i) for i in range(n+1)]).T
yi = np.dot(TT, p) # matrix multiplication calculates the polynomial values
C_yi = np.dot(TT, np.dot(C_p, TT.T)) # C_y = TT*C_z*TT.T
sig_yi = np.sqrt(np.diag(C_yi)) # Standard deviations are sqrt of diagonal

# Do the plotting:
fg, ax = plt.subplots(1, 1)
ax.set_title("Fit for Polynomial (degree {}) with $pm1sigma$-interval".format(n))
ax.fill_between(t, yi+sig_yi, yi-sig_yi, alpha=.25)
ax.plot(t, yi,'-')
ax.plot(x, y, 'ro')
ax.axis('tight')

fg.canvas.draw()
plt.show()

Polyfit with 1-sigma intervals

请注意,计算完整矩阵 C_yi 是计算和memorywise不十分有效。

对 @oliver-w 请求的更新方法上的几个单词:

polyfit 假设参数 x_i 是确定性的,而 y_i 是具有期望值 y_i 和相同方差 sigma的非相关随机变量。 这是一个线性估计问题,一般的最小二乘法可以用。 通过确定残留物的样本方差,sigma 可以近似。 基于 sigmapp 协方差矩阵可以计算出最小二乘方差的维基百科文章。 这几乎就是 polyfit() 所使用的方法: 在 sigma 中,使用了更保守的因子 S/(n-m-2) 而不是 S/(n-m)

原作者:
...