numpy.polyfit:如何获得估计曲线周围的 1-sigma 不确定性?

2024-04-15

我使用 numpy.polyfit 来拟合观察结果。 polyfit 为我提供多项式的估计系数,还可以为我提供估计系数的误差协方差矩阵。美好的。现在,我想知道是否有一种方法可以估计估计曲线周围的 +/- 1sigma 不确定性。

我知道 MatLab 可以做到(https://stats.stackexchange.com/questions/56596/finding-uncertainty-in-coefficients-from-polyfit-in-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) with tt=[t**n, tt*n-1, ..., 1]. t可以是单个值或列向量。由于这是一个线性方程,具有协方差矩阵C_p of 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()

gives Polyfit with 1-sigma intervals

请注意,计算完整矩阵C_yi计算和记忆方面的效率不是很高。

Update- 根据@oliver-w 的请求,就方法论说几句话:

polyfit假设参数x_i是确定性的并且y_i是与期望值不相关的随机变量y_i和相同的方差sigma。因此这是一个线性估计问题,可以使用普通的最小二乘法。通过确定残差的样本方差,sigma可以近似。基于sigma的协方差矩阵pp可以计算如下维基百科关于最小二乘的文章 https://en.wikipedia.org/wiki/Least_squares#Least_squares.2C_regression_analysis_and_statistics。差不多就是这个方法polyfit()用途: 用于sigma更保守的因素S/(n-m-2)代替S/(n-m)用来。

本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

numpy.polyfit:如何获得估计曲线周围的 1-sigma 不确定性? 的相关文章

随机推荐