Question: 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?

Answer: 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