Machine Learning Python
Niveau Expert Data Scientist

1. Régression linéaire — fondements et limites

La régression OLS minimise la somme des carrés des résidus : β = (XᵀX)⁻¹Xᵀy. Elle est non biaisée (théorème Gauss-Markov) mais souffre en haute dimension ou en présence de multicollinéarité. Avant toute modélisation, diagnostiquez les hypothèses : linéarité, homoscédasticité, normalité des résidus, indépendance.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from scipy import stats

model = LinearRegression()
model.fit(X_train, y_train)
y_pred = model.predict(X_test)

# Métriques complètes
print(f"R² : {r2_score(y_test, y_pred):.4f}")
print(f"MAE : {mean_absolute_error(y_test, y_pred):.4f}")
print(f"RMSE : {np.sqrt(mean_squared_error(y_test, y_pred)):.4f}")

# Diagnostic des résidus — test de normalité Shapiro-Wilk
residuals = y_test - y_pred
_, p_value = stats.shapiro(residuals[:5000])
print(f"Shapiro p-value : {p_value:.4f} ({'normal' if p_value > 0.05 else 'non-normal'})")

# Plot résidus vs valeurs prédites
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].scatter(y_pred, residuals, alpha=0.4, s=15)
axes[0].axhline(0, color='red', linestyle='--')
axes[0].set_title('Résidus vs Prédits')
stats.probplot(residuals, dist='norm', plot=axes[1])
axes[1].set_title('Q-Q Plot des résidus')
plt.tight_layout()

2. Ridge, Lasso, ElasticNet — quand utiliser lequel

La régularisation ajoute une pénalité sur les coefficients pour réduire la variance. Le choix dépend de la structure des features :

• Ridge (L2) : contracte tous les coefficients sans les annuler. Idéal quand toutes les features contribuent, en présence de multicollinéarité.
• Lasso (L1) : annule les coefficients des features non pertinentes (sélection de features intégrée). Optimal en haute dimension avec sparsité.
• ElasticNet : combine L1 + L2. Recommandé quand des groupes de features corrélées doivent être conservés ou supprimés ensemble.

from sklearn.linear_model import RidgeCV, LassoCV, ElasticNetCV
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

# RidgeCV trouve automatiquement le meilleur alpha par CV
pipe_ridge = Pipeline([
    ('scaler', StandardScaler()),
    ('ridge', RidgeCV(alphas=np.logspace(-3, 3, 50), cv=5))
])
pipe_ridge.fit(X_train, y_train)
print(f"Ridge meilleur alpha : {pipe_ridge.named_steps['ridge'].alpha_:.4f}")

# LassoCV — coefficients annulés = features non pertinentes
pipe_lasso = Pipeline([
    ('scaler', StandardScaler()),
    ('lasso', LassoCV(cv=5, max_iter=10000, random_state=42))
])
pipe_lasso.fit(X_train, y_train)
n_zero = (pipe_lasso.named_steps['lasso'].coef_ == 0).sum()
print(f"{n_zero} features éliminées sur {X_train.shape[1]}")

3. Régression quantile — prédire un intervalle de confiance

Contrairement à l’OLS qui prédit la moyenne conditionnelle, la régression quantile prédit un percentile arbitraire — essentiel pour la prévision de stocks (percentile 90% pour éviter les ruptures) ou la tarification (percentile 10% pour des prix compétitifs).

from sklearn.linear_model import QuantileRegressor

# Prédiction d’intervalle de confiance à 80%
results = {}
for q in [0.1, 0.5, 0.9]:
    qr = QuantileRegressor(quantile=q, alpha=0.01, solver='highs')
    qr.fit(X_train, y_train)
    results[q] = qr.predict(X_test)

# Taux de couverture réel de l’intervalle [Q10, Q90]
coverage = ((y_test >= results[0.1]) & (y_test <= results[0.9])).mean()
print(f"Couverture IC 80% : {coverage:.1%}")

1. K-Means++ — initialisation intelligente et évaluation

K-Means minimise l’inertie intra-cluster. L’initialisation K-Means++ (initialisation k-means++) évite les minima locaux en espaçant les centroïdes initiaux. La méthode du coude + Silhouette score détermine le k optimal.

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score, silhouette_samples
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

# Normalisation indispensable pour K-Means
X_scaled = StandardScaler().fit_transform(X)

# Méthode du coude + Silhouette
K_range = range(2, 11)
inertias, silhouettes = [], []

for k in K_range:
    km = KMeans(n_clusters=k, init='k-means++', n_init=10, random_state=42)
    labels = km.fit_predict(X_scaled)
    inertias.append(km.inertia_)
    silhouettes.append(silhouette_score(X_scaled, labels))

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
ax1.plot(list(K_range), inertias, 'bo-')
ax1.set_title('Méthode du coude (inertie)')
ax2.plot(list(K_range), silhouettes, 'ro-')
ax2.set_title('Silhouette score (plus haut = meilleur)')
plt.tight_layout()

2. DBSCAN — clustering basé sur la densité

DBSCAN ne requiert pas de spécifier k et détecte les clusters de forme arbitraire. Il identifie les points de bruit (label -1) — idéal pour la détection d’anomalies. Les deux paramètres clés : eps (rayon de voisinage) et min_samples (densité minimale).

from sklearn.cluster import DBSCAN
from sklearn.neighbors import NearestNeighbors

# Trouver eps optimal — coude de la courbe k-distances
nn = NearestNeighbors(n_neighbors=5)
nn.fit(X_scaled)
distances, _ = nn.kneighbors(X_scaled)
k_distances = np.sort(distances[:, -1])[::-1]

plt.plot(k_distances)
plt.xlabel('Points triés')
plt.ylabel('Distance au 5e voisin')
plt.title('Courbe k-distances — le coude donne eps optimal')
plt.show()

# Application DBSCAN
dbscan = DBSCAN(eps=0.5, min_samples=5)
labels = dbscan.fit_predict(X_scaled)
n_clusters = len(set(labels)) - (1 if -1 in labels else 0)
n_noise = (labels == -1).sum()
print(f"Clusters : {n_clusters} | Points de bruit : {n_noise} ({n_noise/len(labels):.1%})")

3. Gaussian Mixture Models — clustering probabiliste

Les GMM modélisent les données comme un mélange de gaussiennes. Contrairement à K-Means, chaque point appartient à chaque cluster avec une probabilité — utile pour la segmentation client où les frontières sont floues.

from sklearn.mixture import GaussianMixture

# Sélection du nombre de composantes par BIC
bic_scores = []
for n in range(1, 11):
    gmm = GaussianMixture(n_components=n, covariance_type='full', random_state=42)
    gmm.fit(X_scaled)
    bic_scores.append(gmm.bic(X_scaled))

best_n = np.argmin(bic_scores) + 1
print(f"Nombre optimal de composantes (BIC) : {best_n}")

# Probabilités d’appartenance par cluster
best_gmm = GaussianMixture(n_components=best_n, random_state=42)
best_gmm.fit(X_scaled)
proba = best_gmm.predict_proba(X_scaled)  # shape (n, k)

1. PCA — variance expliquée et choix du nombre de composantes

La PCA projette les données sur les axes de variance maximale (vecteurs propres de la matrice de covariance). Elle réduit le bruit, accélère l’entraînement et visualise des données haute dimension. La variance expliquée cumulée guide le choix de n_components.

from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

# PCA complète pour analyser la variance expliquée
pca_full = PCA().fit(X_scaled)
cumvar = np.cumsum(pca_full.explained_variance_ratio_)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
ax1.bar(range(1, 21), pca_full.explained_variance_ratio_[:20], color='#7c3aed')
ax1.set_title('Variance expliquée par composante')
ax2.plot(range(1, len(cumvar)+1), cumvar, 'o-', color='#7c3aed')
ax2.axhline(0.95, color='red', linestyle='--', label='95% variance')
ax2.set_title('Variance cumulée')
ax2.legend()
plt.tight_layout()

# Nombre de composantes pour 95% de variance
n95 = (cumvar >= 0.95).argmax() + 1
print(f"{n95} composantes pour 95% de la variance")

2. t-SNE vs UMAP — visualisation 2D de données haute dimension

t-SNE et UMAP ne conservent pas les distances globales — ils ne sont utilisés que pour la visualisation, jamais comme preprocessing pour un modèle supervisé. UMAP est généralement préféré : plus rapide, déterministe et conserve mieux la structure globale.

from sklearn.manifold import TSNE
import umap

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# t-SNE — lent, stochastique, perplexity = voisinage local
tsne = TSNE(n_components=2, perplexity=30, learning_rate=200, random_state=42)
X_tsne = tsne.fit_transform(X_scaled[:5000])
axes[0].scatter(X_tsne[:, 0], X_tsne[:, 1], c=y[:5000], cmap='tab10', s=5, alpha=0.7)
axes[0].set_title('t-SNE')

# UMAP — rapide, déterministe, conserve la structure globale
reducer = umap.UMAP(n_neighbors=15, min_dist=0.1, random_state=42)
X_umap = reducer.fit_transform(X_scaled)
axes[1].scatter(X_umap[:, 0], X_umap[:, 1], c=y, cmap='tab10', s=5, alpha=0.7)
axes[1].set_title('UMAP')
plt.tight_layout()

1. Feature engineering temporel — extraire de l’information de la date

La majorité de la valeur dans les séries temporelles vient du feature engineering, pas du modèle. Un gradient boosting avec de bonnes features temporelles bat souvent ARIMA ou LSTM sur des données tabulaires.

import pandas as pd
import numpy as np

def create_temporal_features(df, date_col='date'):
    df = df.copy()
    df[date_col] = pd.to_datetime(df[date_col])

    # Features calendaires
    df['year'] = df[date_col].dt.year
    df['month'] = df[date_col].dt.month
    df['day_of_week'] = df[date_col].dt.dayofweek
    df['week_of_year'] = df[date_col].dt.isocalendar().week.astype(int)
    df['is_weekend'] = (df['day_of_week'] >= 5).astype(int)
    df['quarter'] = df[date_col].dt.quarter

    # Encodage cyclique (mois, jour) — préserve la continuité circulaire
    df['month_sin'] = np.sin(2 * np.pi * df['month'] / 12)
    df['month_cos'] = np.cos(2 * np.pi * df['month'] / 12)
    df['dow_sin'] = np.sin(2 * np.pi * df['day_of_week'] / 7)

    # Features de lag (valeurs passées)
    for lag in [1, 7, 14, 28]:
        df[f'lag_{lag}'] = df['target'].shift(lag)

    # Rolling statistics
    for w in [7, 30]:
        df[f'rolling_mean_{w}'] = df['target'].shift(1).rolling(w).mean()
        df[f'rolling_std_{w}'] = df['target'].shift(1).rolling(w).std()
    return df

2. Validation walk-forward — éviter le data leakage temporel

Le split aléatoire est interdit sur les séries temporelles — il révèle le futur au modèle. La validation walk-forward (expanding ou sliding window) simule le déploiement réel.

from sklearn.model_selection import TimeSeriesSplit
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.metrics import mean_absolute_percentage_error

tscv = TimeSeriesSplit(n_splits=5, gap=7)  # gap = 7 jours de séparation
scores = []

for fold, (train_idx, val_idx) in enumerate(tscv.split(X)):
    X_tr, X_val = X[train_idx], X[val_idx]
    y_tr, y_val = y[train_idx], y[val_idx]

    model = GradientBoostingRegressor(n_estimators=200, learning_rate=0.05,
                                         max_depth=4, random_state=42)
    model.fit(X_tr, y_tr)
    mape = mean_absolute_percentage_error(y_val, model.predict(X_val))
    scores.append(mape)
    print(f"Fold {fold+1} MAPE : {mape:.2%}")

print(f"MAPE moyen : {np.mean(scores):.2%} ± {np.std(scores):.2%}")

1. Pipeline complet — numérique + catégoriel + modèle

Le pipeline garantit que le scaler et l’encodeur sont fités uniquement sur les données d’entraînement. Sans pipeline, un fit_transform manuel sur l’ensemble du dataset introduit du data leakage — l’erreur classique qui fait paraître les modèles plus performants qu’ils ne le sont réellement.

from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.impute import SimpleImputer
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import cross_val_score

num_features = ['age', 'revenu', 'anciennete']
cat_features = ['region', 'segment', 'canal']

num_pipeline = Pipeline([
    ('imputer', SimpleImputer(strategy='median')),
    ('scaler', StandardScaler()),
])

cat_pipeline = Pipeline([
    ('imputer', SimpleImputer(strategy='most_frequent')),
    ('encoder', OneHotEncoder(handle_unknown='ignore', sparse_output=False)),
])

preprocessor = ColumnTransformer([
    ('num', num_pipeline, num_features),
    ('cat', cat_pipeline, cat_features),
])

full_pipeline = Pipeline([
    ('preprocessor', preprocessor),
    ('model', RandomForestRegressor(n_estimators=200, random_state=42)),
])

# CV sans leakage — le pipeline entier est fité sur chaque fold
scores = cross_val_score(full_pipeline, X_df, y, cv=5,
                         scoring='neg_root_mean_squared_error')
print(f"RMSE CV : {-scores.mean():.4f} ± {scores.std():.4f}")

2. Hyperparameter tuning avec Optuna — optimisation bayésienne

Optuna explore l’espace des hyperparamètres intelligemment via l’algorithme TPE (Tree-structured Parzen Estimator), bien plus efficace que GridSearch ou RandomizedSearch sur des espaces larges.

import optuna
from sklearn.model_selection import cross_val_score
from xgboost import XGBRegressor

def objective(trial):
    params = {
        'n_estimators': trial.suggest_int('n_estimators', 100, 1000),
        'learning_rate': trial.suggest_float('learning_rate', 1e-3, 0.3, log=True),
        'max_depth': trial.suggest_int('max_depth', 3, 10),
        'subsample': trial.suggest_float('subsample', 0.6, 1.0),
        'colsample_bytree': trial.suggest_float('colsample_bytree', 0.6, 1.0),
        'reg_alpha': trial.suggest_float('reg_alpha', 1e-8, 10.0, log=True),
    }
    model = XGBRegressor(**params, random_state=42, verbosity=0)
    score = cross_val_score(model, X_train, y_train, cv=5,
                            scoring='neg_root_mean_squared_error').mean()
    return -score

study = optuna.create_study(direction='minimize')
study.optimize(objective, n_trials=100, show_progress_bar=True)
print(f"Meilleurs params : {study.best_params}")

1. Sérialisation du pipeline avec joblib

import joblib
from datetime import datetime
import json

# Sauvegarder le pipeline complet + métadonnées
model_artifact = {
    'pipeline': full_pipeline,
    'feature_names': num_features + cat_features,
    'trained_at': datetime.now().isoformat(),
    'cv_rmse': float(-scores.mean()),
}
joblib.dump(model_artifact, 'model_v1.joblib', compress=3)

# Chargement et prédiction
artifact = joblib.load('model_v1.joblib')
pipeline = artifact['pipeline']
prediction = pipeline.predict(new_data_df)

2. API de prédiction avec FastAPI

from fastapi import FastAPI
from pydantic import BaseModel
import joblib, pandas as pd

app = FastAPI()
artifact = joblib.load('model_v1.joblib')

class PredictRequest(BaseModel):
    age: float
    revenu: float
    anciennete: float
    region: str
    segment: str
    canal: str

@app.post("/predict")
async def predict(req: PredictRequest):
    df = pd.DataFrame([req.model_dump()])
    pred = artifact['pipeline'].predict(df)[0]
    return {"prediction": round(float(pred), 4), "model_version": "v1"}

3. Détection de data drift — surveiller la dégradation en production

from scipy.stats import ks_2samp
import pandas as pd

def detect_drift(X_train, X_prod, threshold=0.05):
    """Test KS pour détecter un glissement de distribution."""
    results = {}
    for col in X_train.select_dtypes(include='number').columns:
        stat, p_value = ks_2samp(X_train[col].dropna(), X_prod[col].dropna())
        results[col] = {
            'ks_stat': round(stat, 4),
            'p_value': round(p_value, 4),
            'drift_detected': p_value < threshold
        }
    drift_df = pd.DataFrame(results).T
    drifted = drift_df[drift_df['drift_detected']]
    if len(drifted) > 0:
        print(f"⚠ Drift détecté sur {len(drifted)} features : {list(drifted.index)}")
    return drift_df