import numpy as np

def RegressionPolynomiale(x,y,k):
    # fonction pour la régression polynomiale sur R^2
    # x, y : vecteurs colonne donnant les points (xi,yi)
    # k : ordre de la régression. Le polynôme optimal sera de 
    # degré au plus k
    n = x.shape[0]
    # on crée la matrice M par concaténations successives
    M = np.empty((n,0))
    for j in range(k+1):
        M = np.hstack((M,x**j))
    # on résout le système linéaire M^T M beta = M^T y
    beta = np.linalg.solve( M.T@M , M.T@y )
    # on retourne le vecteur beta donnant les coefficients du polynôme
    return beta

# on définit les points xi et yi sous forme de vecteurs colonne
x = np.random.rand(10,1)
y = -5 + 12*x + np.random.randn(10,1)

# ordre de la régression polynomiale. k=3 signifie qu'on utilise des polynômes
# de degré au plus 3
k = 3

beta = RegressionPolynomiale(x,y,k)

# affichage de la régression
import matplotlib.pyplot as plt
# affichage des points xi, yi
plt.plot(x,y,'*')
# affichage du graphe du polynôme
xx = np.linspace(0,1,1000)
yy = 0
for j in range(k+1):
    yy += beta[j] * xx**j
plt.plot(xx,yy)
plt.show()