Régression Bayésienne : Guide complet — Principes, Exemples et Implémentation Python

Résumé — La régression bayésienne est une approche probabiliste de la régression linéaire qui applique le théorème de Bayes pour estimer une distribution complète sur les poids du modèle, plutôt qu’un simple point fixe. Contrairement à la régression classique qui retourne une unique valeur pour chaque coefficient, la régression bayésienne fournit une distribution a posteriori permettant de quantifier l’incertitude de chaque prédiction. Ce guide couvre la théorie fondamentale, l’intuition géométrique, une implémentation complète en Python avec scikit-learn, et des cas d’usage concrets.

Principe mathématique de la régression bayésienne

Le théorème de Bayes comme fondement

Au cœur de la régression bayésienne se trouve le théorème de Bayes, qui permet d’inverser une probabilité conditionnelle. Dans le contexte de la régression, on cherche à connaître la distribution des poids w sachant les données observées D = (X, y) :

P(w | D) = P(D | w) × P(w) / P(D)

où :

• P(w | D) est la distribution a posteriori — ce que l’on croit sur les poids après avoir observé les données.
• P(D | w) est la vraisemblance — la probabilité d’observer les données étant donné des poids particuliers.
• P(w) est la distribution a priori (ou prior) — ce que l’on croit sur les poids avant de voir les données.
• P(D) est la probabilité marginale (ou évidence), servant de constante de normalisation.

Modèle de vraisemblance gaussienne

On suppose que les observations suivent un modèle linéaire avec du bruit gaussien :

y = Xw + ε, où ε ~ N(0, α⁻¹ I)

La vraisemblance s’écrit alors comme le produit de N gaussiennes indépendantes, chacune centrée sur la prédiction linéaire wᵀxᵢ avec une variance α⁻¹.

Prior gaussien sur les poids (régularisation Ridge)

Le choix le plus courant pour la distribution a priori est une gaussienne centrée en zéro :

P(w | λ) = N(w | 0, λ⁻¹ I)

Ce prior gaussien a une conséquence majeure : le Maximum A Posteriori (MAP) sous ce prior est exactement équivalent à la régression Ridge. Cela signifie que la régression Ridge classique peut être interprétée comme une estimation ponctuelle MAP dans un cadre bayésien avec prior gaussien. Le paramètre de régularisation λ de Ridge correspond à l’inverse de la variance du prior sur les poids.

Posterior analytique

Grâce à la conjugaison entre la vraisemblance gaussienne et le prior gaussien, la distribution a posteriori admet une forme analytique fermée :

p(w | X, y) = N(w | m_N, S_N)

avec :

S_N = (α XᵀX + λ I)⁻¹

m_N = α S_N Xᵀy

Ce résultat est fondamental : on obtient non seulement une estimation centrale des poids (m_N), mais aussi une matrice de covariance (S_N) qui mesure l’incertitude sur chaque coefficient et leurs corrélations.

Distribution prédictive avec variance épistémique

Pour une nouvelle entrée x_*, la distribution prédictive s’obtient en intégrant sur la posterior des poids :

p(y_ | x_, X, y) = ∫ p(y_ | x_, w) · p(w | X, y) dw

Ce qui donne :

p(y_ | x_, X, y) = N(y_ | m_Nᵀ x_, σ*²)

où la variance prédictive se décompose en deux termes :

σ² = α⁻¹ + x_ᵀ S_N x_*

• Le premier terme (α⁻¹) représente le bruit inhérent aux données (variance aléatorique).
• Le second terme (x_ᵀ S_N x_) représente l’incertitude épistémique — celle qui diminue lorsqu’on ajoute des données.

C’est cette capacité à séparer les deux sources d’incertitude qui rend la régression bayésienne si puissante.

Intuition : pourquoi estimer une distribution plutôt qu’un point ?

En régression classique (moindres carrés ordinaires ou Ridge), on obtient un unique vecteur de poids ŵ. C’est comme si on disait : « Les coefficients sont exactement ces valeurs-là, point final. » Mais cette approche a un défaut majeur : elle ne dit rien sur la confiance que l’on peut avoir en ces estimations.

La régression bayésienne adopte une philosophie radicalement différente. Au lieu d’estimer un seul point pour les poids, on estime une distribution complète sur l’espace des poids possibles. Cette distribution capture tout ce que le modèle « sait » et tout ce qu’il « ignore ».

Imaginez que vous essayez de deviner la position d’un trésor enfoui. L’approche classique vous donnerait une coordonnée unique : « Le trésor est ici. » L’approche bayésienne dirait plutôt : « Il y a 68 % de chances qu’il soit dans ce cercle de 5 mètres, 95 % dans ce cercle de 10 mètres. » La seconde réponse est infiniment plus utile pour prendre des décisions éclairées.

Plus on dispose de données, plus la distribution a posteriori se resserre autour de la vraie valeur des poids. Avec un nombre infini de données, la posterior converge vers une masse de Dirac centrée sur les vrais paramètres. C’est cette propriété de convergence qui garantit la cohérence de l’approche bayésienne.

L’intuition géométrique est aussi éclairante : dans l’espace des poids, le prior gaussien pénalise les poids de grande norme (comme Ridge), mais la posterior explore un volume autour du maximum, plutôt qu’un simple pic. Les régions où les données sont peu informatives conservent une variance élevée, signalant clairement au praticien : « Attention, je ne suis pas sûr ici. »

Implémentation Python complète

Installation des dépendances

pip install scikit-learn numpy matplotlib seaborn

Comparaison BayesianRidge vs Ridge

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import make_regression
from sklearn.linear_model import BayesianRidge, Ridge
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score

# Configuration du style
sns.set_style("whitegrid")
sns.set_palette("deep")

# 1. Génération des données
random_state = 42
X, y, coef_true = make_regression(
    n_samples=150,
    n_features=1,
    n_informative=1,
    noise=15.0,
    coef=True,
    random_state=random_state
)

# Ajout de bruit hétéroscédastique pour illustrer l'incertitude
np.random.seed(random_state)
y += np.random.normal(0, np.abs(X).flatten() * 2, size=y.shape)

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=random_state
)

# 2. Entraînement — Régression bayésienne
bayesian_ridge = BayesianRidge(
    n_iter=300,
    tol=1e-3,
    alpha_1=1e-6,
    alpha_2=1e-6,
    lambda_1=1e-6,
    lambda_2=1e-6,
    compute_score=True,
    fit_intercept=True
)
bayesian_ridge.fit(X_train, y_train)

# 3. Entraînement — Ridge (pour comparaison)
ridge = Ridge(alpha=1.0)
ridge.fit(X_train, y_train)

# 4. Prédictions avec incertitude (BayesianRidge uniquement)
X_plot = np.linspace(X.min(), X.max(), 500).reshape(-1, 1)

y_pred_bayes, sigma_bayes = bayesian_ridge.predict(X_plot, return_std=True)
y_pred_ridge = ridge.predict(X_plot)

# Bandes d'incertitude à 95 % (μ ± 2σ)
y_lower = y_pred_bayes - 2 * sigma_bayes
y_upper = y_pred_bayes + 2 * sigma_bayes

# 5. Évaluation
y_pred_test_bayes, sigma_test = bayesian_ridge.predict(X_test, return_std=True)
y_pred_test_ridge = ridge.predict(X_test)

print("=" * 60)
print("RÉGRESSION BAYÉSIENNE VS RIDGE")
print("=" * 60)
print(f"BayesianRidge — MSE : {mean_squared_error(y_test, y_pred_test_bayes):.4f}")
print(f"BayesianRidge — R²  : {r2_score(y_test, y_pred_test_bayes):.4f}")
print(f"Ridge         — MSE : {mean_squared_error(y_test, y_pred_test_ridge):.4f}")
print(f"Ridge         — R²  : {r2_score(y_test, y_pred_test_ridge):.4f}")
print(f"\nPoids estimé (Bayes) : {bayesian_ridge.coef_}")
print(f"Poids estimé (Ridge) : {ridge.coef_}")
print(f"Poids réel           : {coef_true}")
print(f"\nAlpha estimé (précision du bruit) : {bayesian_ridge.alpha_:.4f}")
print(f"Lambda estimé (précision des poids) : {bayesian_ridge.lambda_:.4f}")
print(f"Score de l'évidence (log marginal) : {bayesian_ridge.scores_[-1]:.4f}")

# 6. Visualisation
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# Graphique 1 : Prédictions avec bandes d'incertitude
axes[0].scatter(X_train, y_train, color="steelblue", alpha=0.6, label="Données d'entraînement", zorder=5)
axes[0].scatter(X_test, y_test, color="coral", alpha=0.6, label="Données de test", zorder=5)
axes[0].plot(X_plot, y_pred_bayes, color="navy", linewidth=2, label="Prédiction moy. (Bayes)")
axes[0].fill_between(
    X_plot.flatten(), y_lower, y_upper,
    color="navy", alpha=0.15, label="Intervalle 95 % (μ ± 2σ)"
)
axes[0].set_title("Régression bayésienne — Incertitude", fontweight="bold")
axes[0].set_xlabel("X")
axes[0].set_ylabel("y")
axes[0].legend(fontsize=8)

# Graphique 2 : Comparaison Bayes vs Ridge
axes[1].scatter(X_train, y_train, color="steelblue", alpha=0.4, zorder=5)
axes[1].plot(X_plot, y_pred_bayes, color="navy", linewidth=2, label="BayesianRidge")
axes[1].plot(X_plot, y_pred_ridge, color="crimson", linewidth=2, linestyle="--", label="Ridge")
axes[1].set_title("BayesianRidge vs Ridge", fontweight="bold")
axes[1].set_xlabel("X")
axes[1].set_ylabel("y")
axes[1].legend(fontsize=8)

# Graphique 3 : Évolution du score d'évidence
axes[2].plot(bayesian_ridge.scores_, color="navy", linewidth=2)
axes[2].set_title("Score de l'évidence (log P(y|X))", fontweight="bold")
axes[2].set_xlabel("Itération")
axes[2].set_ylabel("Log marginal likelihood")
axes[2].axhline(y=bayesian_ridge.scores_[-1], color="crimson", linestyle="--", alpha=0.7)

plt.tight_layout()
plt.savefig("bayesian_ridge_results.png", dpi=150, bbox_inches="tight")
plt.show()

# 7. Utilisation de ARDRegression (Automatic Relevance Determination)
from sklearn.linear_model import ARDRegression

ard = ARDRegression(
    n_iter=300,
    tol=1e-3,
    alpha_1=1e-6,
    alpha_2=1e-6,
    lambda_1=1e-6,
    lambda_2=1e-6,
    threshold_lambda=10_000,
    fit_intercept=True
)
ard.fit(X_train, y_train)
y_pred_ard = ard.predict(X_test)

print(f"\nARDRegression — MSE : {mean_squared_error(y_test, y_pred_ard):.4f}")
print(f"ARDRegression — R²  : {r2_score(y_test, y_pred_ard):.4f}")
print(f"Coefficients ARD : {ard.coef_}")

Points clés du code

• return_std=True : méthode exclusive à BayesianRidge qui retourne l’écart-type prédictif σ pour chaque point. C’est ce qui permet de tracer les bandes d’incertitude.
• compute_score=True : enregistre le score de l’évidence (logarithme de la vraisemblance marginale) à chaque itération. Utile pour diagnostiquer la convergence.
• alpha_1, alpha_2, lambda_1, lambda_2 : hyperparamètres des distributions Gamma a priori sur les précisions α et λ. Des valeurs proches de zéro correspondent à des priors non-informatifs.
• ARDRegression : variante qui attribue un λ_j par coefficient, permettant l’Automatic Relevance Determination — les variables non pertinentes voient leur coefficient tendre vers zéro.

Hyperparamètres de BayesianRidge

Comment choisir les hyperparamètres

Les paramètres alpha_1, alpha_2, lambda_1 et lambda_2 contrôlent les priors Gamma sur les précisions. Par défaut, avec des valeurs très petites (1e-6), on obtient des priors faiblement informatifs qui laissent les données parler. Si vous avez une connaissance a priori forte sur le niveau de bruit ou l’ampleur des coefficients, vous pouvez augmenter ces valeurs.

Pour n_iter, 300 itérations sont généralement suffisantes. Surveillez la courbe d’évidence : si elle stagne avant 300 itérations, le modèle a convergé. Vous pouvez réduire n_iter pour accélérer l’entraînement.

Avantages et Limites

Avantages

1. Quantification de l’incertitude — C’est l’avantage principal. Chaque prédiction vient avec un intervalle de confiance naturel. En sécurité, en médecine ou en finance, savoir « je ne suis pas sûr » est parfois aussi important que la prédiction elle-même.
2. Régularisation automatique — Les hyperparamètres α et λ sont appris automatiquement à partir des données via la maximisation de la vraisemblance marginale (Type-II Maximum Likelihood). Pas besoin de validation croisée coûteuse pour trouver le bon paramètre de régularisation.
3. Interprétation probabiliste rigoureuse — Contrairement aux méthodes fréquentistes qui s’appuient sur des approximations asymptotiques, le cadre bayésien fournit une interprétation exacte (sous les hypothèses du modèle).
4. Équivalence avec Ridge au MAP — On retrouve les garanties classiques de régularisation L2 tout en bénéficiant de l’incertitude complète.
5. Robustesse avec peu de données — Le prior régularise naturellement le modèle, empêchant le surapprentissage même lorsque le nombre d’échantillons est faible par rapport au nombre de caractéristiques.

Limites

1. Coût computationnel — Le calcul de la posterior nécessite l’inversion d’une matrice de taille (D × D) où D est le nombre de caractéristiques. Cela demande O(D³) opérations, ce qui devient prohibitif au-delà de quelques milliers de dimensions.
2. Hypothèse gaussienne — La forme analytique fermée repose sur l’hypothèse d’un prior gaussien et d’une vraisemblance gaussienne. Si les données présentent des queues lourdes, des outliers importants, ou des relations non linéaires, le modèle peut être mal calibré.
3. Modèle linéaire — Comme toutes les régressions linéaires, la régression bayésienne ne capture pas les interactions complexes non linéaires sans ingénierie de caractéristiques explicite (polynômes, noyaux, etc.).
4. Sensibilité au choix des hyperpriors — Bien que α et λ soient appris automatiquement, les hyperparamètres des priors Gamma peuvent influencer les résultats, surtout avec peu de données.
5. Pas de sélection de variables native — Contrairement au Lasso, la régression bayésienne standard ne force pas les coefficients à zéro. Pour obtenir de la parcimonie, il faut utiliser ARDRegression qui attribue un prior individuel à chaque coefficient.

Cas d’usage concrets

1. Modélisation d’incertitude en médecine

Dans le diagnostic médical ou la prédiction de risques, il est crucial de connaître la marge d’erreur d’une prédiction. La régression bayésienne peut être utilisée pour prédire un indicateur de santé (par exemple, le taux de glycémie) à partir de variables cliniques, tout en fournissant un intervalle de confiance. Un intervalle large signale au médecin que la prédiction est incertaine et qu’un examen complémentaire est nécessaire.

Exemple : Prédiction du temps de récupération post-opératoire. Un intervalle de confiance à 95 % large permet de planifier des ressources hospitalières avec une marge de sécurité appropriée.

2. Données avec peu d’échantillons (small data)

Dans de nombreux domaines — recherche pharmaceutique, ingénierie aérospatiale, études climatiques rares — le nombre d’observations est limité. La régression bayésienne excelle dans ce régime grâce à son prior qui agit comme un régularisateur naturel. Là où une régression classique surajuste rapidement, la régression bayésienne maintient des estimations raisonnables en s’appuyant sur les connaissances a priori encodées dans les priors.

Exemple : Prédiction des performances d’un nouveau matériau à partir de seulement 30 expériences en laboratoire, en incorporant un prior informatif issu de simulations physiques.

3. A/B testing et expérimentation

Dans un contexte d’A/B testing, on veut non seulement estimer l’effet d’un traitement (la différence entre les groupes A et B) mais aussi quantifier la certitude de cette estimation. La régression bayésienne fournit naturellement cette information via la distribution a posteriori du coefficient associé à la variable de traitement.

Exemple : Mesure de l’impact d’une nouvelle fonctionnalité sur le temps de session des utilisateurs. La distribution a posteriori du coefficient donne directement la probabilité que l’effet soit positif (par exemple, 97 % de chances d’un effet positif).

4. Exploration bayésienne et optimisation

La régression bayésienne est souvent utilisée comme modèle de substitution (surrogate model) dans l’optimisation bayésienne. La variance épistémique guide l’exploration vers les régions de l’espace de recherche où le modèle est le plus incertain, permettant un équilibre optimal entre exploration et exploitation.

Exemple : Optimisation d’hyperparamètres d’un réseau de neurones. Le modèle bayésien prédit la performance et identifie les combinaisons de paramètres les plus prometteuses à évaluer ensuite.

Voir aussi

• Calcul de la Somme des Produits en Python : Guide Complet et Astuces Pratiques
• Explorer les Nombres de $k$-Markov avec Python : Guide Complet pour Algorithmes et Programmation