I am trying to implement polynomial regression in Python and have it mostly working, but it is not overfitting enough for higher degree polynomials. This is a weird thing to be upset about, I know, but it may point to an underlying problem with my implementation.
In this case I'm using a toy example where I'm sampling uniformly in the region $[0,1]$ and setting output values to $\sin(2\pi x) + \mathcal{N}(0, 0.3)$, i.e. just sin with added Gaussian noise with variance $0.3$.
I am using the normal equation form rather than gradient descent, so solutions should be optimal least-squares solutions. By Lagrange's interpolation theorem, we should have a polynomial of degree n fits n + 1 points exactly. However, for 10 samples and a 9th-degree polynomial (or even higher degrees), my implementation does not fit every point exactly.
You'll notice in my code that I'm not representing all the values of x in the graph exactly (i.e., I may be off after the ten-thousandths place), but the effect should not be as noticeable as it is if that is causing it. Because I'm sampling the x values uniformly in $[0, 1]$, I could have two overlapping x values within my error region which would cause some possible poor fitting in the graph, but this is extremely unlikely for low sample sizes (see my own question and solution to this problem here )
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import math
import random
plt.rcParams['figure.figsize'] = [10, 5]
def sample(n_samples):
''' Samples uniformly from x axis and returns y values given by sin(2pi)
with Gaussian noise'''
x = np.random.uniform(0, 1, size=n_samples)
y = np.sin(2*np.pi*x) + np.random.normal(0, .3, n_samples)
return [x,y]
def poly_regress(data, n_param):
''' Implementation of polynomial regression with n powers. Uses normal
equation rather than gradient descent.
data: 2 by n matrix with data[0] the x values and data[1] the y values.
n_param: degree of returned polynomial
'''
coeff = np.zeros(n_param) # Coefficient vector initialization
# Calculate the vandermonde matrix of input x
vand_x = np.vander(data[0],N = n_param + 1, increasing=True)
# Computes the coefficient vector via the normal equation
coeff = np.linalg.pinv(vand_x.T @ vand_x) @ (vand_x.T @ data[1])
# This is for plotting purposes. If we wanted to predict a particular x
# or range of values, then we'd use an almost identical form.
x = np.arange(0,1, .0001)
y = 0
for i in range(n_param + 1):
y += coeff[i]*np.power(x, i)
# plots the true value of x. This only works for this toy example.
true_x = np.arange(0,1, .01)
true_y = np.sin(2*np.pi*true_x)
plt.plot(true_x, true_y, color='g', label= "True Curve")
plt.plot(x, y, color='b', label="Prediction")
plt.plot(data[0], data[1], 'r+', label="samples")
plt.ylim(-2,2)
plt.legend()
plt.show()
return [x,y]
Any help is greatly appreciated. This is just base case testing before extending to more general forms of regression, so I want it to be functioning ideally at this point.
Cubic case:
High degree case with no noise:
With similar results for lower degree polynomials.
