Find uncertainty from polyfit

Question

I use simple polyfit of order 2 to fit a line in sample data:

np.polyfit(x, y, 2)

which returns the coefficients.

Now I want to find uncertainty of the fitted line, and tried to use cov argument, which returns 3x3 covariance matrix:

np.polyfit(x, y, 2, cov=True)

But I'm not sure how to calculate the uncertainty, which according my Google search should be calculated by squaring the diagonal of covariance matrix.

Answer 1

This problem is addressed by "Estimating Errors in Least-Squares Fitting" by P.H. Richter, 1995, TDA Progress Report 42-122.

From the report, this paragraph may already be sufficient to you

The first instance considered above, namely, determining the error of one or more fitting parameters, has a straightforward answer given in terms of the diagonal elements of the covariance matrix of the fit, and is well known.

The diagonal elements you are interested in are for example:

x = linspace(0,1,1000)
# comment and uncomment the last term to see how the fit appears in the figure,
# and how the covariances of the single polynomial coefficients vary in turn.
y = cos(x)*x**2+x+sin(x-1.) #+(x*1.3)**6
p,cov = polyfit(x,y,2,cov=True)
plot(x,y,'b')
plot(x,polyval(p,x),'r')
print sqrt(diag(cov))

More in general, the reference addresses how this error in the polynomial coefficients is also an error of the dependent variable y as a function of the independent variable x. From the report:

It is the purpose of this article to discuss the above errors and, in particular, to present results that will permit one to determine the standard error of the fit as a function of the independent variable, as well as to establish confidence limits for these errors.

Answer 2

For your convenience I made a fully working example for Python 3 based on gg349's answer.

import numpy as np
import matplotlib.pyplot as plt 

x = np.linspace(0,1,1000)
# comment and uncomment the last term to see how the fit appears in the figure,
# and how the covariances of the single polynomial coefficients vary in turn.
y = np.cos(x) * x**2 + x + np.sin(x - 1.) \
#     + (x * 1.3)**6

p, cov = np.polyfit(x, y, 2, cov=True)

plt.plot(x, y)
plt.plot(x, np.polyval(p,x))
plt.show()

print(np.sqrt(np.diag(cov)))