numpy.polyfit:如何在估计曲线周围获得1-sigma不确定性? [英] numpy.polyfit: How to get 1-sigma uncertainty around the estimated curve?

问题描述

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

I use numpy.polyfit to fit observations. polyfit gives me the estimated coefficients of the polynomial and can also provides me the error covariance matrix of the estimated coefficients. Fine. Now, I would like to know if there is a way to estimate the +/- 1sigma uncertainty around the estimated curve.

我知道MatLab可以做到( https://stats.stackexchange.com/questions/56596/finding-uncertainty-in-coefficients-from-polyfit-in-matlab )，但我没有找到在python中制作它的方法.

I know MatLab can do it (https://stats.stackexchange.com/questions/56596/finding-uncertainty-in-coefficients-from-polyfit-in-matlab) but I did not found a way to make it in 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可以是单个值或列向量.由于这是一个线性方程，因此协方差矩阵C_p为p，因此值的协方差矩阵为np.dot(tt, np.dot(C_p tt.T)).

If you have enough data points, you can get with the parameter cov=True an estimated covariance matrix from polyfit(). Remember that you can write a polynomial p[0]*t**n + p[1]*t**(n-1) + ... + p[n] as the matrix product np.dot(tt, p) with tt=[t**n, tt*n-1, ..., 1]. t can be either a single value or a column vector. Since this a linear equation, with the covariance matrix C_p of p, the covariance matrix of the values is 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在计算和记忆方面都不是很有效.

Note that calculating the complete matrix C_yi is computationally and memorywise not very efficient.

更新-根据@ oliver-w的要求，对方法进行了说明:

Update - on the request on @oliver-w a couple words on the methodology:

polyfit假定参数x_i是确定性的，而y_i是具有期望值y_i和相同方差sigma的不相关随机变量.因此，这是一个线性估计问题，可以使用普通的最小二乘法.通过确定残基的样本方差，可以近似sigma.可以根据sigma计算pp的协方差矩阵，如 Wikipedia文章中所示在最小二乘法上.这几乎是polyfit()所使用的方法:对于sigma，使用的是更保守的因子S/(n-m-2)而不是S/(n-m).

polyfit assumes that the parameters x_i are deterministic and y_i are uncorrelated random variables with the expected value y_i and identical variance sigma. So it is a linear estimation problem and an ordinary least squares method can be used. By determining the sample variance of the residues, sigma can be approximated. Based onsigma the covariance matrix of pp can be calculated as shown in the Wikipedia article on Least Squares. That is almost the method that polyfit() uses: There for sigma the more conservative factor S/(n-m-2) instead of S/(n-m) is used.

这篇关于numpy.polyfit:如何在估计曲线周围获得1-sigma不确定性?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！