我使用numpy.polyfit来拟合观测数据。polyfit会给出多项式的估计系数,还可以提供估计系数的误差协方差矩阵。现在,我想知道是否有一种方法来估计估计曲线周围的+/- 1sigma不确定性。
我知道MatLab可以做到这一点(https://stats.stackexchange.com/questions/56596/finding-uncertainty-in-coefficients-from-polyfit-in-matlab),但我没有找到在Python中实现的方式。
我使用numpy.polyfit来拟合观测数据。polyfit会给出多项式的估计系数,还可以提供估计系数的误差协方差矩阵。现在,我想知道是否有一种方法来估计估计曲线周围的+/- 1sigma不确定性。
我知道MatLab可以做到这一点(https://stats.stackexchange.com/questions/56596/finding-uncertainty-in-coefficients-from-polyfit-in-matlab),但我没有找到在Python中实现的方式。
cov=True
参数从 polyfit()
获取一个估计的协方差矩阵。请记住,您可以将多项式 p[0]*t**n + p[1]*t**(n-1) + ... + p[n]
写成矩阵乘积的形式 np.dot(tt, p)
,其中 tt=[t**n, tt*n-1, ..., 1]
。 t
可以是单个值或列向量。由于这是一条线性方程,因此对于参数 p
的协方差矩阵 C_p
,值的协方差矩阵为 np.dot(tt, np.dot(C_p tt.T))
。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 $\pm1\sigma$-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()
给出。
请注意,计算完整矩阵C_yi
在计算和内存方面效率不高。
更新 - 根据@oliver-w的要求,对方法进行简述:
polyfit
假定参数x_i
是确定性的,而y_i
是无关的随机变量,具有期望值y_i
和相同的方差sigma
。因此,这是一个线性估计问题,可以使用普通最小二乘法。通过确定残差的样本方差,可以近似估计sigma
。基于sigma
,可以像最小二乘法维基百科文章中所示那样计算pp
的协方差矩阵。这几乎是polyfit()
使用的方法:其中使用更加保守的因子S/(n-m-2)
而不是S/(n-m)
来表示sigma
。