Régression linéaire avec pseudo-inverse

import numpy as np

# 1. Génération des données
np.random.seed(42)
n = 100
x = np.random.uniform(0, 10, n)
vrai_a, vrai_b = 2.0, 1.0
bruit = np.random.randn(n) * 1.5
y = vrai_a * x + vrai_b + bruit

print("Régression linéaire sur données synthétiques")
print(f"Équation théorique : y = {vrai_a:.1f}.x + {vrai_b:.1f} + eps")
print(f"Bruit : sigma = {1.5}")

# 2. Matrice de conception (design matrix)
X = np.column_stack([x, np.ones(n)]) # forme (n, 2)

# 3. Calcul par pseudo-inverse et autres méthodes
# Méthode 1: inverse explicite
beta = np.linalg.inv(X.T @ X) @ X.T @ y

# Méthode 2: pseudo-inverse (plus stable numériquement)
beta_pseudo = np.linalg.pinv(X) @ y

# Méthode 3: résolution du système normal
beta_solve = np.linalg.solve(X.T @ X, X.T @ y)

print(f"\nCoefficients trouvés :")
print(f"  Méthode inverse : a = {beta[0]:.4f}, b = {beta[1]:.4f}")
print(f"  Pseudo-inverse  : a = {beta_pseudo[0]:.4f}, b = {beta_pseudo[1]:.4f}")
print(f"  Résolution system: a = {beta_solve[0]:.4f}, b = {beta_solve[1]:.4f}")

# 4. Comparaison avec np.polyfit
coeffs_polyfit = np.polyfit(x, y, 1)
print(f"  np.polyfit      : a = {coeffs_polyfit[0]:.4f}, b = {coeffs_polyfit[1]:.4f}")

# 5. Prédictions
y_pred = beta[0] * x + beta[1]

# 6. Métriques
residus = y - y_pred
R2 = 1 - np.sum(residus ** 2) / np.sum((y - np.mean(y)) ** 2)
print(f"\nR2 = {R2:.4f}")
print(f"RMSE = {np.sqrt(np.mean(residus**2)):.4f}")

L'exécution de ce script affiche :

Régression linéaire sur données synthétiques
Équation théorique : y = 2.0.x + 1.0 + eps
Bruit : sigma = 1.5

Coefficients trouvés :
  Méthode inverse : a = 2.0453, b = 0.7835
  Pseudo-inverse  : a = 2.0453, b = 0.7835
  Résolution system: a = 2.0453, b = 0.7835
  np.polyfit      : a = 2.0453, b = 0.7835

R2 = 0.9572
RMSE = 1.4785