Problème de la Régression Linéaire

L’objectif est de déterminer l’équation d’une courbe de tendance, permettant de modéliser ( et donc de prédire ) le comportement d’un système à partir de données mesurées ( dataset ).

Pour un modèle linéaire à une variable x ( nombre de features n = 1 ),
il faut donc trouver les coefficients a et b pour une droite f(x) = a.x+b ( f(x) est le modèle ) de telle sorte que la distance entre les points réels ( du dataset ) et la droite y soit minimale.

Exemple d’Application :

On veut prédire la distance nécessaire pour qu’une voiture s’arrête complètement en fonction de sa vitesse initiale.

Données collectées :

Vitesse (km/h)	Distance d’arrêt (m)
20	5
40	15
60	25
80	40
100	55

L’objectif est donc de déterminer l’ équation de la droite en rouge faisant le lien entre la vitesse et la distance d’arrêt pour toute vitesse :

Génération du Dataset

Pour aborder le problème du calcul de la régression linéaire, nous allons utiliser un dataset généré à partir de la bibliothèque sklearn.

# Génération du Dataset 
# Nombre de données 
m = 50

import numpy as np
from sklearn.datasets import make_regression
%matplotlib
import matplotlib.pyplot as plt
import random

# la fonction make_regression de sklearn permet de fournir un jeu de valeurs aléatoires 
x, y = make_regression(n_samples = m, n_features = 1, noise = 15)

print('x.shape = ', x.shape)
print('y.shape = ', y.shape)

# Reshape et tri des données
y=y.reshape(m,1)
print('y.shape = ', y.shape)

rd = random.randint(-50,50)
y=y-rd # Ajout d'une ordonnée à l'origine non nulle

plt.scatter(x,y)

for i in range(0,len(x),1):
    print(f'point_{i} :  {x[i][0]} {y[i][0]}')

x.shape =  (50, 1)
y.shape =  (50,)
y.shape =  (50,1)
point_0 :  -0.15189774638632875 -26.05065438545585
point_1 :  -1.4512055894228688 -111.60953597518551
point_2 :  -0.7126960570972625 -61.37126967965305
point_3 :  1.4007134960248657 44.583357456962545
point_4 :  0.28047208495044407 -39.55330043779675
point_5 :  0.42039886367550067 -29.774798254066354
point_6 :  -0.3132497598252263 -68.55046226847585
point_7 :  -0.36037024497739073 -69.80428021450066
point_8 :  0.2625007638567827 -29.447452014898694
....................

png

Fonction coût

Si je définis un modèle au hasard :

# Fonction coût : erreurs entre modèle et points réels

# choix des coefficients a et b au hasard :  
a = 110
b = -20

modele = a*x+b

# Tracé du nuage de points et de la courbe de tendance
plt.figure(figsize=(10, 6))
plt.scatter(x, y, label="Nuage de points", color="blue", alpha=0.6)
plt.plot(x, modele, label="Courbe de tendance y = ax + b", color="red", linewidth=2)

# Ajout des barres entre les points réels et la courbe de tendance
for i in range(len(x)):
    plt.plot([x[i][0], x[i][0]], [y[i][0], modele[i][0]], color="gray", linestyle="--", linewidth=0.8)

plt.xlabel("X")
plt.ylabel("Y")
plt.title("Nuage de points et courbe de tendance avec erreurs")
plt.legend()
plt.grid(True)
plt.show()

Pour estimer l’ensemble des erreurs entre le modèle et les valeurs réelles, nous allons calculer la fonction coût J faisant la somme de ces erreurs ( somme élevée au carré pour considérer des valeurs absolues )

J(a,b) =\frac{1}{2.m}\sum_{i=0}^{m-1}\left(f(x^{(i)})-y^{(i)}\right)^2

# Calcul fonction cout J(a,b)

J = (1/(2*m))*np.sum((modele-y)**2)
print('J = ',J)

J =  306.0313228764071

Afin de minimiser ces erreurs, nous allons chercher un modèle de telle sorte que la fonction coût J soit minimale.
J dépend du choix de a et b

# Représentation de J(a,b) en 3D

# Grille de paramètres a et b
a_values = np.linspace(-2000, 2000, 100)  # Variation des pentes
b_values = np.linspace(-2000, 2000, 100)  # Variation des intercepts
a_grid, b_grid = np.meshgrid(a_values, b_values)

# Calcul de J pour chaque combinaison de a et b
J_grid = np.zeros_like(a_grid)
for i in range(len(a_values)):
    for j in range(len(b_values)):
        a = a_values[i]
        b = b_values[j]
        modele = a * x + b
        J_grid[j, i] = (1 / (2 * m)) * np.sum((modele - y) ** 2)

# Tracé de la surface 3D
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(a_grid, b_grid, J_grid, cmap='viridis', alpha=0.9)
ax.set_title('Représentation 3D de la fonction coût J(a, b)', fontsize=14)
ax.set_xlabel('Paramètre a (pente)', fontsize=12)
ax.set_ylabel('Paramètre b (intercept)', fontsize=12)
ax.set_zlabel('Coût J', fontsize=12)
plt.show()

png

Une coupe transversale de ce graphe donne :

$\frac{\partial J(a,b) }{\partial a}$ et $\frac{\partial J(a,b) }{\partial b}$ sont les gradients.

Afin de trouver l’erreur minimale, on pourrait résoudre directement les équations $\frac{\partial J(a,b) }{\partial a} = 0$ et $\frac{\partial J(a,b) }{\partial b} = 0$

Si cela est envisageable dans un cas simple comme celui que nous considérons, cela impliquerait des calculs trop exigeants dans des cas plus complexes.

C’est pourquoi nous allons rechercher l’erreur minimale avec l’algorithme de la descente de gradient.

Algorithme de minimisation de la fonction coût

L’algorithme de la descente de gradient consiste à faire évoluer les coefficients a et b, afin de viser le creux du graphique précédent.

a_{i+1}=a_i-\alpha.\frac{\partial J(a_i) }{\partial a}

b_{i+1}=b_i-\alpha.\frac{\partial J(b_i) }{\partial b}

Quel que soit le point de départ ( d’un côté ou de l’autre de la courbe ), le signe du gradient fera que a évoluera toujours vers le minimum de J.

$\alpha$ est le learning rate et détermine la vitesse ou le pas pour atteindre le creux de la courbe.

Calcul du Gradient

J(a,b) =\frac{1}{2.m}\sum_{i=0}^{m-1}\left(f(x^{(i)})-y^{(i)}\right)^2

J(a,b) =\frac{1}{2.m}\sum_{i=0}^{m-1}\left(a.x^{(i)}+b-y^{(i)}\right)^2

On prend :

$f(a,b) = a.x+b$
$g(f) = f^2$

or $(g \circ f)'=f'*g'(f)$

Donc :

\frac{\partial J(a,b) }{\partial a} = \frac{1}{m} \sum_{i=0}^{m-1}\left( x^i.(a.x^i+b-y^i) \right)

\frac{\partial J(a,b) }{\partial b} = \frac{1}{m} \sum_{i=0}^{m-1}\left( a.x^i+b-y^i \right)

a = -100
b = -100
n_iter = 100

# algorithme de Gradient Descent :  

alpha = 0.01 # alpha toujours positif; learning rate
history_a=[]
history_b=[]


for n in range(0,10,1): # plusieurs itération étant donné que a et b sont dépendants
    
    for i in range(0,n_iter,1):
        a = a - alpha*(1/m)*np.sum(x*(a*x+b-y))
        history_a.append(a) 
        
    
    for i in range(0,n_iter,1):
        b = b - alpha*(1/m)*np.sum(a*x+b-y)
        history_b.append(b) 

print('a_final =', a, 'b_final =', b) 

plt.figure(figsize=(12,6))
plt.subplot(2,1,1) # 2 ligne 1 colonne, position 1
plt.title(r'$ a_{i+1}=a_i-\alpha.\frac{\partial J(a_i) }{\partial a} $')
plt.plot(history_a)
plt.grid()
plt.subplot(2,1,2) # 2 ligne 1 colonne, position 2
plt.title(r'$ b_{i+1}=b_i-\alpha.\frac{\partial J(b_i) }{\partial b} $')
plt.plot(history_b)

plt.grid()
plt.subplots_adjust(left=0.1, right=0.9, top=0.9, bottom=0.1, wspace=0.4, hspace=0.4) # Adjust spacing
plt.show()

# Fonction coût : erreurs entre modèle et points réels
modele = a*x+b

plt.figure(figsize=(10, 6))
plt.scatter(x, y, label="Nuage de points", color="blue", alpha=0.6)
plt.plot(x, modele, label="Courbe de tendance y = ax + b", color="red", linewidth=2)
# Ajout des barres entre les points réels et la courbe de tendance
for i in range(len(x)):
    plt.plot([x[i][0], x[i][0]], [y[i][0], modele[i][0]], color="gray", linestyle="--", linewidth=0.8)
plt.xlabel("X")
plt.ylabel("Y")
plt.title("Modèle déterminé par l'algorithme de la descente du gradient")
plt.legend()
plt.grid(True)
plt.show()

# Calcul fonction cout J(a,b)

J = (1/(2*m))*np.sum((modele-y)**2)
print('J=',J)

a_final = 56.53939099648491  
b_final = -38.52047857316444

J= 114.8281245428621

Regression Lineaire Régression Linéaire : Calcul Matriciel