numpy.polyfit vs numpy.polynomial.polynomial.polyfit

Question

Why do numpy.polyfit and numpy.polynomial.polynomial.polyfit produce different plots in the test below?

import numpy as np
from numpy.polynomial.polynomial import polyfit
import matplotlib.pyplot as plt

x = np.linspace(0, 10, 50)
y = 5 * x + 10 + (np.random.random(len(x)) - 0.5) * 5

plt.scatter(x, y,marker='.', label='Data for regression')
plt.plot(x, np.poly1d(np.polyfit(x, y, 1))(x), label='numpy.polyfit')
plt.plot(x, np.poly1d(polyfit(x, y, 1))(x), label='polynomial.polyfit')
plt.legend()
plt.show()

Answer 1

At first glance, the documentation seems to indicate they should give the same result -

numpy.polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False)

Least squares polynomial fit.

Fit a polynomial p(x) = p[0] * x**deg + ... + p[deg] of degree deg to points (x, y). Returns a vector of coefficients p that minimises the squared error in the order deg, deg-1, … 0.

and

numpy.polynomial.polynomial.polyfit(x, y, deg, rcond=None, full=False, w=None)

Least-squares fit of a polynomial to data.

Return the coefficients of a polynomial of degree deg that is the least squares fit to the data values y given at points x. If y is 1-D the returned coefficients will also be 1-D. If y is 2-D multiple fits are done, one for each column of y, and the resulting coefficients are stored in the corresponding columns of a 2-D return. The fitted polynomial(s) are in the form

p(x) = c₀ + c₁ * x + ... + c_n * xⁿ

But the difference is in the order of coefficients returned from the two methods, at least for the use case in question.

• numpy.polyfit returns the coefficients in descending order of degree, according to the generation equation
  p(x) = c_n * xⁿ + c_(n-1) * x^(n-1) + ... + c₁ * x + c₀
• numpy.polynomial.polynomial.polyfit returns the coefficients in ascending order of degree, according to the generation equation
  p(x) = c₀ + c₁ * x + ... + c_(n-1) * x^(n-1) + c_n * xⁿ

though mathematically identical, those two equations are not the same in ndarray representation. This might be obfuscated by the use of different notations in the documentation. For demonstration, consider the following

import numpy as np

x = np.linspace(0, 10, 50)
y = x**2 + 5 * x + 10

print(np.polyfit(x, y, 2))
print(np.polynomial.polynomial.polyfit(x, y, 2))

[ 1.  5. 10.]
[10.  5.  1.]

Both methods get the same result, but in opposite order, the former being what np.poly1d() expects,

print(np.poly1d(np.polyfit(x, y, 2)))
print(np.poly1d(np.polynomial.polynomial.polyfit(x, y, 2)))

   2
1 x + 5 x + 10
    2
10 x + 5 x + 1

and the latter being what the np.polynomial.polynomial.Polynomial() constructor expects.,

print(np.polynomial.polynomial.Polynomial(np.polynomial.polynomial.polyfit(x, y, 2)))
print(np.polynomial.polynomial.Polynomial(np.polyfit(x, y, 2)))

poly([10.  5.  1.])  # 10 + 5 * x + 1 * x**2
poly([ 1.  5. 10.])  # 1 + 5 * x + 10 * x**2

Flipping the result from np.polynomial.polynomial.polyfit before passing it to poly1d() or using a np.polynomial.polynomial.Polynomial will produce the expected result: