When and how to use Polynomial.fit() as opposed to polyfit()?

Using Python 3.10.0 and NumPy 1.21.4.

I'm trying to understand why Polynomial.fit() calculates wildly different coefficient values from polyfit().

In the following code:

import numpy as np

def main():
    x = np.array([3000, 3200, 3400, 3600, 3800, 4000, 4200, 4400, 4600, 4800, 5000, 5200, 5400, 5600, 5800, 6000, 6200, 6400, 6600, 6800, 7000])
    y = np.array([5183.17702344, 5280.24520952, 5758.94478531, 6070.62698406, 6584.21169885, 8121.20863245, 7000.57326186, 7380.01493624, 7687.97802847, 7899.71417408, 8506.90860692, 8421.73816463, 8705.58403352, 9275.46094996, 9552.44715196, 9850.70796049, 9703.53073907, 9833.39941224, 9900.21604921, 9901.06392084, 9974.51206378])

    c1 = np.polynomial.polynomial.polyfit(x, y, 2)
    c2 = np.polynomial.polynomial.Polynomial.fit(x, y, 2).coef

    print(c1)
    print(c2)

if __name__ == '__main__':

    main()

c1 contains:

[-3.33620814e+03  3.44704650e+00 -2.18221029e-04]

which produces the the line of best fit when plugged a + bx + cx^2 that I predicted while c2 contains:

[8443.4986422  2529.67242075 -872.88411679]

which results in a very different line when plugged into the same formula.

The documentation seems to imply that Polynomial.fit() is the new preferred way of calculating the line but it keeps outputting the wrong coefficients (unless my understanding of polynomial regression is completely wrong).

If I am not using the functions correctly, what is the correct way of using them?

If I am using both functions correctly, why would I use Polynomial.fit() over polyfit(), as the documentation seems to imply I should?

Solution

According to Polynomial.fit() documentation, it returns:

A series that represents the least squares fit to the data and has the domain and window specified in the call. If the coefficients for the unscaled and unshifted basis polynomials are of interest, do new_series.convert().coef.

You can find in https://numpy.org/doc/stable/reference/routines.polynomials.html#transitioning-from-numpy-poly1d-to-numpy-polynomial that

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.

You can check

import numpy as np

def main():
    x = np.array([3000, 3200, 3400, 3600, 3800, 4000, 4200, 4400, 4600, 4800, 5000, 5200, 5400, 5600, 5800, 6000, 6200, 6400, 6600, 6800, 7000])
    y = np.array([5183.17702344, 5280.24520952, 5758.94478531, 6070.62698406, 6584.21169885, 8121.20863245, 7000.57326186, 7380.01493624, 7687.97802847, 7899.71417408, 8506.90860692, 8421.73816463, 8705.58403352, 9275.46094996, 9552.44715196, 9850.70796049, 9703.53073907, 9833.39941224, 9900.21604921, 9901.06392084, 9974.51206378])

    c1 = np.polynomial.polynomial.polyfit(x, y, 2)
    c2 = np.polynomial.polynomial.Polynomial.fit(x, y, 2).convert().coef
    c3 = np.polynomial.polynomial.Polynomial.fit(x, y, 2, window=(x.min(), x.max())).coef

    print(c1)
    print(c2)
    print(c3)

if __name__ == '__main__':

    main()

# [-3.33620814e+03  3.44704650e+00 -2.18221029e-04]
# [-3.33620814e+03  3.44704650e+00 -2.18221029e-04]
# [-3.33620814e+03  3.44704650e+00 -2.18221029e-04]

Another argument for using np.polynomial.Polynomial class is stated here in the docs https://numpy.org/doc/stable/reference/routines.polynomials.package.html