numpy polynomial.Polynomial.fit() gives different coefficients than polynomial.polyfit()

Viewed 1124

I do not understand why polynomial.Polynomial.fit() gives coefficients very different from the expected coefficients :

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))
print(np.polynomial.polynomial.Polynomial.fit(x, y, 2))

Gives :

[ 1.  5. 10.]
[10.  5.  1.]
poly([60. 75. 25.])

The two first results are OK, and thanks to this answer I understand why the two arrays are in reversed order.

However, I do not understand the signification of the third result. The coefficients looks wrong, though the polynomial that I got this way seems to give correct predicted values.

2 Answers

The answer is slightly hidden in the docs, of course. Looking at the class numpy.polynomial.polynomial.Polynomial(coef, domain=None, window=None) It is clear that in general the coefficients [a, b, c, ...] are for the polynomial a + b * x + c * x**2 + .... However, there are the keyword parameters domain and window both with default [-1,1]. I am not into that class, so I am not sure about the purpose, but it is clear that a remapping takes place. Now in the case of polynomial.Polynomial.fit() one has a class method that automatically takes the x data as domain, but still makes the mapping to the window. Hence, in the OP [0-10] is mapped onto [-1,1]. This is done by x = x' / 5 - 1 or x' -> 5 * x + 5. Putting the latter in the OP polynomial we get

( 5 x' + 5 )**2 + 5 * ( 5 * x' + 5 ) + 10 = 25 * x'**2 + 75 * x' + 60

Voila.

To get the expected result one has to put

print(np.polynomial.polynomial.Polynomial.fit(x, y, 2, window=[0, 10] ) )

wich gives

poly([10.  5.  1.])

Buried in the docs:

Note that the coefficients are given in the scaled domain defined by the linear mapping between the window and domain. convert can be used to get the coefficients in the unscaled data domain.

So use:

poly.convert()    

This will rescale your coefficients to what you are probably expecting.

Example for data generated from 1 + 2x + 3x^2:

from numpy.polynomial import Polynomial

test_poly = Polynomial.fit([0, 1, 2, 3, 4, 5],
                           [1, 6, 17, 34, 57, 86],
                           2)

print(test_poly)
print(test_poly.convert())

Output:

poly([24.75 42.5  18.75])
poly([1. 2. 3.])
Related