fit numpy polynomials to noisy data

Question

I want to exactly represent my noisy data with a numpy.polynomial polynomial. How can I do that?.

In this example, I chose legendre polynomials. When I use the polynomial legfit function, the coefficients it returns are either very large or very small. So, I think I need some sort of regularization.

Why doesn't my fit get more accurate as I increase the degree of the polynomial? (It can be seen that the 20, 200, and 300 degree polynomials are essentially identical.) Are any regularization options available in the polynomial package?

I tried implementing my own regression function, but it feels like I am re-inventing the wheel. Is making my own fitting function the best path forward?

from scipy.optimize import least_squares as mini
import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(0, 5, 1000)
tofit = np.sin(3 * x) + .6 * np.sin(7*x) - .5 * np.cos(3 * np.cos(10 * x))

# This is here to illustrate what I expected the legfit function to do
# I expected it to do least squares regression with some sort of regularization.
def myfitfun(x, y, deg):
    def fitness(a):
       return ((np.polynomial.legendre.legval(x, a) - y)**2).sum() + np.sum(a ** 2)
    return mini(fitness, np.zeros(deg)).x

degrees = [2, 4, 8, 16, 40, 200]
plt.plot(x, tofit, c='k', lw=4, label='Data')
for deg in degrees:
    #coeffs = myfitfun(x, tofit, deg)
    coeffs = np.polynomial.legendre.legfit(x, tofit, deg)
    plt.plot(x, np.polynomial.legendre.legval(x, coeffs), label="Degree=%i"%deg)
plt.legend()

Answer 1

Legendre polynomials are meant to be used over the interval [-1,1]. Try to replace x with 2*x/x[-1] - 1 in your fit and you'll see that all is good:

nx = 2*x/x[-1] - 1
for deg in degrees:
    #coeffs = myfitfun(x, tofit, deg)
    coeffs = np.polynomial.legendre.legfit(nx, tofit, deg)
    plt.plot(x, np.polynomial.legendre.legval(nx, coeffs), label="Degree=%i"%deg)

Answer 2

The easy way to use the proper interval in the fit is to use the Legendre class

from numpy.polynomial import Legendre as L
p = L.fit(x, y, order)

This will scale and shift the data to the interval [-1, 1] and track the scaling factors.