How To Include Measurement Errors In Numpy.polyfit

The documentation HERE for numpy.polyfit defines the optional input vector w as the 'weights to apply to the y-coordinates'. This definition is unclear and does not seem the standa

Solution 1:

Short Answer

The relation between w and the 1σ errors in numpy.polyfit is

w = 1/sigma

which is different from what everybody will expect.

Edit: following my comment on Github, the numpy.polyfit 1.16 documentation now explicitly states "use 1/sigma (not 1/sigma**2)" to avoid people thinking there is a typo in the formula.

Explanation

In least-squares fitting one generally defines the weights vector in such a way that the fit minimizes the squared error (see e.g. Wikipedia or NIST)

chi2 = np.sum(weights*(p(x) - y)**2)

In the common situation where the 1σ errors (standard deviations) "sigma" are known one has the familiar textbook relation where the weights are the reciprocal of the variance

weights = 1/sigma**2

However the numpy.polyfit documentation defines the weight as "weights to apply to the y-coordinates". This definition is not quite correct. The weights apply to the fit residuals, not only to the y-coordinates.

More importantly, looking at the math in the Numpy (v1.9.1) code it appears that the Numpy code solves the linear problem below in the last-squares sense, where the w vector does indeed multiply the y-coordinate

(vander*w[:, np.newaxis]).dot(x) == y*w

But solving the above array expression in the least-squares sense is equivalent to minimizing the expression below with w inside the parenthesis

chi2 = np.sum((w*(vander.dot(x) - y))**2)

Or, using the notation of the Numpy documentation

chi2 = np.sum((w*(p(x) - y))**2)

in such a way that the relation between w and the 1σ errors is

w = 1/sigma

which is different from what everybody will expect.