Distributions en pratique#

Hide code cell source

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import seaborn as sns
from scipy import stats
from scipy.special import boxcox1p

sns.set_theme(style="whitegrid", palette="muted", font_scale=1.1)
rng = np.random.default_rng(42)

La théorie des probabilités décrit des distributions idéales. Les données réelles n’obéissent jamais exactement à ces distributions — mais elles s’en approchent souvent suffisamment pour que les outils associés soient valides. Ce chapitre enseigne à naviguer entre ces deux mondes : reconnaître la forme d’une distribution réelle, la tester, l’ajuster, et transformer les données quand nécessaire.

Reconnaître une distribution sur données réelles#

La première étape est toujours visuelle. L’histogramme, la densité estimée par noyau (KDE) et le QQ-plot donnent des informations complémentaires sur la forme d’une distribution.

# Génération de jeux de données avec différentes distributions
n = 500

# Revenus : log-normale (la loi des revenus par excellence)
revenus = rng.lognormal(mean=3.6, sigma=0.5, size=n)  # en dizaines de milliers €

# Temps d'attente : exponentielle (files d'attente, temps entre événements)
temps_attente = rng.exponential(scale=5, size=n)  # minutes

# Notes d'examen : bêta (bornée entre 0 et 20)
notes = 20 * rng.beta(a=5, b=2, size=n)

# Températures : normale (phénomène symétrique)
temperatures = rng.normal(loc=15, scale=5, size=n)

datasets = {
    "Revenus (log-normale)": revenus,
    "Temps d'attente (exponentielle)": temps_attente,
    "Notes d'examen (bêta)": notes,
    "Températures (normale)": temperatures,
}

Hide code cell source

fig, axes = plt.subplots(2, 4, figsize=(16, 7))

for col, (nom, data) in enumerate(datasets.items()):
    # Histogramme + KDE
    ax = axes[0, col]
    ax.hist(data, bins=30, density=True, alpha=0.55, color="steelblue")
    kde_x = np.linspace(data.min(), data.max(), 300)
    ax.plot(kde_x, stats.gaussian_kde(data)(kde_x), color="tomato", lw=2)
    skew = stats.skew(data)
    ax.set_title(f"{nom}\nskew={skew:.2f}", fontsize=8.5)
    ax.set_xlabel("Valeur")

    # QQ-plot par rapport à la normale
    ax2 = axes[1, col]
    (osm, osr), (slope, intercept, r) = stats.probplot(data, dist="norm")
    ax2.scatter(osm, osr, alpha=0.4, s=10, color="steelblue")
    x_line = np.array([osm.min(), osm.max()])
    ax2.plot(x_line, slope * x_line + intercept, "r-", lw=1.5)
    ax2.set_xlabel("Quantiles théoriques (normale)")
    ax2.set_ylabel("Quantiles observés")
    ax2.set_title(f"QQ-plot (r²={r**2:.3f})")

plt.suptitle("Histogrammes et QQ-plots pour quatre distributions", y=1.01, fontsize=12)
plt.tight_layout()
plt.show()
_images/de06d8545ea241d79571f509e77eb2a51f24ae6cabefe0bd5fe7dcac95bce8b9.png

Le QQ-plot : lecture et interprétation#

Le QQ-plot (quantile-quantile plot) compare les quantiles empiriques d’un échantillon aux quantiles théoriques d’une loi de référence. Si l’échantillon suit cette loi, les points s’alignent sur la droite \(y = x\) (ou sur la droite de régression ajustée).

Les déviations caractéristiques ont chacune une signature visuelle reconnaissable.

Hide code cell source

fig, axes = plt.subplots(2, 3, figsize=(15, 9))

cas = [
    ("Normale parfaite", rng.normal(0, 1, 500), "norm",
     "Points sur la droite — distribution normale"),
    ("Queues lourdes (t)", rng.standard_t(df=3, size=500), "norm",
     "Points s'éloignent aux deux extrémités (S)"),
    ("Asymétrie droite", rng.lognormal(0, 0.8, 500), "norm",
     "Points au-dessus de la droite à droite"),
    ("Asymétrie gauche", -rng.lognormal(0, 0.8, 500), "norm",
     "Points en-dessous de la droite à gauche"),
    ("Bimodale", np.concatenate([rng.normal(-2, 0.8, 250), rng.normal(2, 0.8, 250)]), "norm",
     "Sigmoïde avec inflexion centrale"),
    ("Discrète/bornée", rng.integers(1, 11, 500).astype(float), "norm",
     "Escaliers visibles"),
]

for ax, (titre, data, dist_ref, interpretation) in zip(axes.flat, cas):
    (osm, osr), (slope, intercept, r) = stats.probplot(data, dist=dist_ref)
    ax.scatter(osm, osr, alpha=0.4, s=12, color="steelblue")
    x_line = np.array([osm.min(), osm.max()])
    ax.plot(x_line, slope * x_line + intercept, "r-", lw=2)
    ax.set_title(titre, fontsize=9.5, fontweight="bold")
    ax.set_xlabel("Quantiles théoriques")
    ax.set_ylabel("Quantiles observés")
    ax.annotate(interpretation, xy=(0.03, 0.97), xycoords="axes fraction",
                fontsize=7.5, va="top", style="italic",
                bbox=dict(boxstyle="round,pad=0.3", fc="lightyellow", alpha=0.8))

plt.suptitle("Signatures visuelles des déviations à la normalité dans les QQ-plots",
             y=1.01, fontsize=12)
plt.tight_layout()
plt.show()
_images/02e4499f8d9b96f4700de7b20b3c31e9a51baa346076a9ff10cfb5f11fe092e4.png

Lecture pratique du QQ-plot

  • Points sur la droite : la distribution correspond bien à la loi de référence.

  • Courbe en S : queues plus lourdes que la référence (leptokurtique).

  • Courbe en S inversé : queues plus légères (platykurtique).

  • Inflexion vers le haut à droite : asymétrie positive (queue droite étalée).

  • Inflexion vers le bas à gauche : asymétrie négative.

  • Marches d’escalier : données discrètes ou nombreuses ex-aequo.

Tests de normalité#

Les tests de normalité répondent formellement à la question : « peut-on rejeter l’hypothèse que ces données suivent une loi normale ? »

Attention à l’usage des tests de normalité

Ces tests détectent la non-normalité, pas la normalité. Ne pas rejeter H₀ ne prouve pas que les données sont normales — cela signifie seulement que l’échantillon est trop petit pour détecter la déviation. Avec un grand n, le test rejette H₀ pour des déviations pratiquement insignifiantes.

Shapiro-Wilk#

Optimal pour les petits échantillons (\(n < 50\)). Il compare les statistiques d’ordre à leurs espérances sous la normale. Il est le plus puissant des tests courants pour détecter des déviations à la normalité.

Kolmogorov-Smirnov (avec correction Lilliefors)#

Basé sur la distance maximale entre la ECDF et la CDF théorique. Dans scipy, stats.kstest suppose les paramètres connus ; utiliser stats.kstest avec les paramètres estimés à partir des données invalide le test (biais conservateur). La correction de Lilliefors (statsmodels.stats.diagnostic.lilliefors) est préférable.

D’Agostino-Pearson#

Combine le skewness et le kurtosis en une statistique de test. Robuste, applicable pour \(n > 20\), disponible via scipy.stats.normaltest.

def tests_normalite(data, nom="données"):
    """Applique trois tests de normalité et retourne les résultats."""
    resultats = {}

    # Shapiro-Wilk (recommandé pour n < 5000)
    if len(data) <= 5000:
        stat_sw, p_sw = stats.shapiro(data)
        resultats["Shapiro-Wilk"] = {"statistique": stat_sw, "p-valeur": p_sw}

    # D'Agostino-Pearson
    stat_dp, p_dp = stats.normaltest(data)
    resultats["D'Agostino-Pearson"] = {"statistique": stat_dp, "p-valeur": p_dp}

    # Jarque-Bera
    stat_jb, p_jb = stats.jarque_bera(data)
    resultats["Jarque-Bera"] = {"statistique": stat_jb, "p-valeur": p_jb}

    df_res = pd.DataFrame(resultats).T
    df_res["rejet H0 (α=0.05)"] = df_res["p-valeur"] < 0.05
    print(f"\n=== Tests de normalité : {nom} ===")
    print(df_res.round(4).to_string())
    return df_res

_ = tests_normalite(temperatures, "Températures (normale)")
_ = tests_normalite(revenus[:500], "Revenus (log-normale)")
_ = tests_normalite(temps_attente[:500], "Temps d'attente (exponentielle)")
=== Tests de normalité : Températures (normale) ===
                    statistique  p-valeur  rejet H0 (α=0.05)
Shapiro-Wilk             0.9973    0.5917              False
D'Agostino-Pearson       0.5649    0.7539              False
Jarque-Bera              0.5131    0.7737              False

=== Tests de normalité : Revenus (log-normale) ===
                    statistique  p-valeur  rejet H0 (α=0.05)
Shapiro-Wilk             0.8742       0.0               True
D'Agostino-Pearson     189.4560       0.0               True
Jarque-Bera            711.3370       0.0               True

=== Tests de normalité : Temps d'attente (exponentielle) ===
                    statistique  p-valeur  rejet H0 (α=0.05)
Shapiro-Wilk             0.8013       0.0               True
D'Agostino-Pearson     237.9953       0.0               True
Jarque-Bera           1220.2718       0.0               True

Hide code cell source

# Illustration de l'effet de la taille d'échantillon sur les tests
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Pour des données légèrement non-normales (t avec 10 degrés de liberté)
tailles = [20, 50, 100, 200, 500, 1000, 2000]
n_sim = 200
puissances = []

for n_t in tailles:
    rejets = []
    for _ in range(n_sim):
        # t(10) est légèrement non-normale (kurtosis = 1)
        d = rng.standard_t(df=10, size=n_t)
        _, p = stats.shapiro(d) if n_t <= 5000 else stats.normaltest(d)
        rejets.append(p < 0.05)
    puissances.append(np.mean(rejets))

axes[0].plot(tailles, puissances, "o-", color="steelblue", lw=2)
axes[0].axhline(0.05, color="tomato", linestyle="--", label="Niveau α=0.05")
axes[0].axhline(1.0, color="seagreen", linestyle=":", alpha=0.5)
axes[0].set_xscale("log")
axes[0].set_xlabel("Taille d'échantillon n")
axes[0].set_ylabel("Fréquence de rejet de H₀")
axes[0].set_title("Puissance du test de Shapiro-Wilk\n(données t(10), légèrement non-normales)")
axes[0].legend()

# QQ-plot : effet de la taille d'échantillon
for n_t, alpha, color in [(20, 0.9, "tomato"), (100, 0.7, "steelblue"), (1000, 0.5, "seagreen")]:
    d = rng.standard_t(df=10, size=n_t)
    (osm, osr), (slope, intercept, _) = stats.probplot(d, dist="norm")
    axes[1].scatter(osm, osr, alpha=alpha, s=8, color=color, label=f"n={n_t}")

x_line = np.linspace(-4, 4, 100)
axes[1].plot(x_line, x_line, "k-", lw=1.5, label="Diagonale")
axes[1].set_xlabel("Quantiles théoriques")
axes[1].set_ylabel("Quantiles observés")
axes[1].set_title("QQ-plot : t(10) pour différents n\n(queue légèrement lourde)")
axes[1].legend(fontsize=9)

plt.tight_layout()
plt.show()
_images/8bf2aa4a347b9d14d6902d93e5a4f9211ef4d816b858fb4fc9b26ab94bd6e939.png

Ajustement d’une loi : scipy.stats.fit()#

Lorsqu’on souspecte une distribution particulière, on peut ajuster ses paramètres sur les données observées par maximum de vraisemblance (MLE).

# Ajustement de plusieurs lois sur les revenus (log-normaux)
lois_candidates = ["lognorm", "gamma", "expon", "norm", "weibull_min"]
data_revenus = revenus.copy()

resultats_fit = []
for nom_loi in lois_candidates:
    loi = getattr(stats, nom_loi)
    try:
        params = loi.fit(data_revenus)
        # Log-vraisemblance
        loglik = loi.logpdf(data_revenus, *params).sum()
        k = len(params)
        n = len(data_revenus)
        aic = 2 * k - 2 * loglik
        bic = k * np.log(n) - 2 * loglik
        # Test KS
        stat_ks, p_ks = stats.kstest(data_revenus, nom_loi, args=params)
        resultats_fit.append({
            "Loi": nom_loi,
            "Paramètres": str([round(p, 3) for p in params]),
            "Log-vraisemblance": round(loglik, 1),
            "AIC": round(aic, 1),
            "BIC": round(bic, 1),
            "KS p-valeur": round(p_ks, 4),
        })
    except Exception as e:
        pass

df_fit = pd.DataFrame(resultats_fit).set_index("Loi").sort_values("AIC")
print("Ajustement de lois sur les revenus — comparaison par AIC/BIC :")
print(df_fit.to_string())
Ajustement de lois sur les revenus — comparaison par AIC/BIC :
                                                             Paramètres  Log-vraisemblance     AIC     BIC  KS p-valeur
Loi                                                                                                                    
lognorm                  [np.float64(0.518), 2.432, np.float64(33.608)]            -2138.1  4282.3  4294.9       0.8426
gamma         [np.float64(2.398), np.float64(9.31), np.float64(13.156)]            -2139.5  4285.0  4297.7       0.3237
weibull_min  [np.float64(1.544), np.float64(10.024), np.float64(34.38)]            -2147.1  4300.2  4312.8       0.0548
expon                                                  [10.142, 30.713]            -2212.3  4428.7  4437.1       0.0000
norm                           [np.float64(40.854), np.float64(21.208)]            -2236.7  4477.3  4485.7       0.0000

Hide code cell source

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

x_plot = np.linspace(data_revenus.min(), data_revenus.max(), 400)
colors = ["tomato", "seagreen", "darkorange", "purple", "steelblue"]

# Histogramme + ajustements
ax = axes[0]
ax.hist(data_revenus, bins=40, density=True, alpha=0.4, color="gray", label="Données")
ax.plot(x_plot, stats.gaussian_kde(data_revenus)(x_plot), "k-", lw=2, label="KDE empirique")

for (nom_loi, row), color in zip(df_fit.iterrows(), colors[:len(df_fit)]):
    loi = getattr(stats, nom_loi)
    params_str = row["Paramètres"]
    params = loi.fit(data_revenus)
    ax.plot(x_plot, loi.pdf(x_plot, *params), lw=1.8, color=color,
            label=f"{nom_loi} (AIC={row['AIC']:.0f})", alpha=0.9)

ax.set_xlabel("Revenu (k€)")
ax.set_ylabel("Densité")
ax.set_title("Ajustement de plusieurs lois sur les revenus")
ax.legend(fontsize=8)
ax.set_xlim(data_revenus.min(), np.percentile(data_revenus, 99))

# QQ-plot pour les deux meilleures lois
ax2 = axes[1]
meilleures = df_fit.index[:2]
for nom_loi, color in zip(meilleures, ["tomato", "steelblue"]):
    loi = getattr(stats, nom_loi)
    params = loi.fit(data_revenus)
    (osm, osr), (slope, intercept, r) = stats.probplot(data_revenus, sparams=params,
                                                         dist=loi)
    ax2.scatter(osm, osr, alpha=0.3, s=10, color=color)
    x_line = np.array([osm.min(), osm.max()])
    ax2.plot(x_line, slope * x_line + intercept, "-", color=color, lw=2,
             label=f"{nom_loi} (r²={r**2:.4f})")

ax2.set_xlabel("Quantiles théoriques")
ax2.set_ylabel("Quantiles observés")
ax2.set_title("QQ-plots pour les meilleures lois")
ax2.legend(fontsize=9)

plt.tight_layout()
plt.show()
_images/8a8bd378a4cc3044db3744bf07747babefcd6688cf0de33a756a6bc4bc6467c3.png

Lois utiles en pratique#

Normale \(\mathcal{N}(\mu, \sigma^2)\)#

Omniprésente par le théorème central limite. Convient aux mesures physiques, tailles humaines, erreurs de mesure. À utiliser par défaut quand rien n’indique le contraire.

Log-normale \(\text{LogNorm}(\mu, \sigma)\)#

Si \(\log X \sim \mathcal{N}(\mu, \sigma^2)\), alors \(X\) est log-normale. Modèle naturel des revenus, prix immobiliers, durées de vie de composants, tailles de fichiers. La queue droite est étalée, la distribution est toujours positive.

Exponentielle \(\text{Exp}(\lambda)\)#

Modèle des temps entre événements dans un processus de Poisson : temps entre arrivées dans une file d’attente, durée entre appels téléphoniques. Sa propriété remarquable est l’absence de mémoire : \(P(X > s+t \mid X > s) = P(X > t)\).

Pareto \(\text{Pareto}(\alpha, x_m)\)#

Formalise la « loi de puissance » : 20 % des fichiers représentent 80 % de l’espace. Convient aux revenus très élevés, tailles de fichiers, fréquence des mots (loi de Zipf).

Weibull \(\text{Weibull}(k, \lambda)\)#

Généralisation de l’exponentielle pour la fiabilité et l”analyse de survie. Le paramètre de forme \(k\) contrôle si le taux de défaillance est croissant (\(k>1\)), constant (\(k=1\), exponentielle) ou décroissant (\(k<1\)).

Hide code cell source

fig, axes = plt.subplots(2, 3, figsize=(15, 9))
x = np.linspace(0.001, 10, 400)

# Normale
ax = axes[0, 0]
for mu, sigma, label in [(0, 1, "μ=0, σ=1"), (1, 0.5, "μ=1, σ=0.5"), (-1, 2, "μ=-1, σ=2")]:
    ax.plot(np.linspace(-8, 8, 400),
            stats.norm.pdf(np.linspace(-8, 8, 400), mu, sigma), lw=2, label=label)
ax.set_title("Loi Normale")
ax.legend(fontsize=8)

# Log-normale
ax = axes[0, 1]
for sigma, label in [(0.3, "σ=0.3"), (0.7, "σ=0.7"), (1.2, "σ=1.2")]:
    ax.plot(x, stats.lognorm.pdf(x, s=sigma), lw=2, label=label)
ax.set_title("Loi Log-normale (μ=0)")
ax.legend(fontsize=8)
ax.set_xlim(0, 8)

# Exponentielle
ax = axes[0, 2]
for rate, label in [(0.5, "λ=0.5"), (1.0, "λ=1.0"), (2.0, "λ=2.0")]:
    ax.plot(x, stats.expon.pdf(x, scale=1/rate), lw=2, label=label)
ax.set_title("Loi Exponentielle")
ax.legend(fontsize=8)

# Pareto
ax = axes[1, 0]
x_p = np.linspace(1, 6, 400)
for alpha, label in [(1.0, "α=1.0"), (2.0, "α=2.0"), (3.5, "α=3.5")]:
    ax.plot(x_p, stats.pareto.pdf(x_p, b=alpha), lw=2, label=label)
ax.set_title("Loi de Pareto (x_m=1)")
ax.legend(fontsize=8)

# Weibull
ax = axes[1, 1]
for k, label, color in [(0.5, "k=0.5 (DFR)", "tomato"),
                         (1.0, "k=1.0 (exponentielle)", "seagreen"),
                         (2.5, "k=2.5 (IFR)", "steelblue"),
                         (5.0, "k=5.0 (pic étroit)", "darkorange")]:
    ax.plot(x, stats.weibull_min.pdf(x, c=k), lw=2, label=label)
ax.set_xlim(0, 4)
ax.set_title("Loi de Weibull (λ=1)")
ax.legend(fontsize=8)

# Données log-normales avec ajustement
ax = axes[1, 2]
data_ln = rng.lognormal(1.5, 0.7, 300)
ax.hist(data_ln, bins=30, density=True, alpha=0.5, color="steelblue", label="Données")
params_ln = stats.lognorm.fit(data_ln, floc=0)
x_fit = np.linspace(0.01, data_ln.max(), 300)
ax.plot(x_fit, stats.lognorm.pdf(x_fit, *params_ln), "r-", lw=2, label="Ajustement log-normale")
ax.set_title("Ajustement log-normale : durées de projets")
ax.legend(fontsize=8)

plt.suptitle("Les lois statistiques utiles en pratique", y=1.01, fontsize=12)
plt.tight_layout()
plt.show()
_images/0febd559d4306fad8e00bed3be4a4066dbf60325fdd6aaf184c568ade9e0c95b.png

Transformations de données#

Quand les données ne sont pas normales alors qu’une méthode le requiert, on peut les transformer pour les rapprocher d’une gaussienne.

Transformation logarithmique#

\(y = \log(x + c)\) (souvent \(c=0\) si \(x > 0\)). Comprime les grandes valeurs, étire les petites. Idéale pour les distributions log-normales et les données à une ordre de grandeur de variation.

Transformation de Box-Cox#

\(y = \frac{x^\lambda - 1}{\lambda}\) pour \(\lambda \neq 0\), \(y = \log(x)\) pour \(\lambda = 0\). Le paramètre \(\lambda\) est estimé par MLE. Requiert \(x > 0\).

Transformation de Yeo-Johnson#

Extension de Box-Cox qui accepte les valeurs négatives : \(y = \begin{cases} \frac{(x+1)^\lambda - 1}{\lambda} & \text{si } x \geq 0, \lambda \neq 0 \\ \log(x+1) & \text{si } x \geq 0, \lambda = 0 \\ \frac{-((-x+1)^{2-\lambda} - 1)}{2-\lambda} & \text{si } x < 0, \lambda \neq 2 \end{cases}\)

# Données log-normales : revenu brut
data_brut = rng.lognormal(mean=3.5, sigma=0.8, size=500)

# Log
data_log = np.log(data_brut)

# Box-Cox
data_bc, lambda_bc = stats.boxcox(data_brut)

# Yeo-Johnson
data_yj, lambda_yj = stats.yeojohnson(data_brut)

print(f"Lambda Box-Cox  : {lambda_bc:.3f}  (proche de 0 → log convient)")
print(f"Lambda Yeo-Johnson : {lambda_yj:.3f}")
print()
for nom, d in [("Brut", data_brut), ("Log", data_log), ("Box-Cox", data_bc), ("Yeo-Johnson", data_yj)]:
    sw_stat, sw_p = stats.shapiro(d[:500])
    sk = stats.skew(d)
    print(f"{nom:12s} : skewness={sk:.3f}, Shapiro p={sw_p:.4f}")
Lambda Box-Cox  : -0.098  (proche de 0 → log convient)
Lambda Yeo-Johnson : -0.136

Brut         : skewness=3.792, Shapiro p=0.0000
Log          : skewness=0.226, Shapiro p=0.1325
Box-Cox      : skewness=0.004, Shapiro p=0.8416
Yeo-Johnson  : skewness=0.008, Shapiro p=0.7388

Hide code cell source

fig, axes = plt.subplots(2, 4, figsize=(16, 7))
transformations = [
    ("Données brutes", data_brut),
    ("Log", data_log),
    ("Box-Cox", data_bc),
    ("Yeo-Johnson", data_yj),
]

for col, (nom, d) in enumerate(transformations):
    # Histogramme
    axes[0, col].hist(d, bins=30, density=True, alpha=0.55, color="steelblue")
    kde_x = np.linspace(d.min(), d.max(), 300)
    axes[0, col].plot(kde_x, stats.gaussian_kde(d)(kde_x), "r-", lw=2)
    sk = stats.skew(d)
    axes[0, col].set_title(f"{nom}\nskew={sk:.3f}")

    # QQ-plot
    (osm, osr), (slope, intercept, r) = stats.probplot(d, dist="norm")
    axes[1, col].scatter(osm, osr, alpha=0.3, s=10, color="steelblue")
    x_line = np.array([osm.min(), osm.max()])
    axes[1, col].plot(x_line, slope * x_line + intercept, "r-", lw=2)
    axes[1, col].set_title(f"QQ-plot (r²={r**2:.4f})")
    axes[1, col].set_xlabel("Théoriques (N)")

plt.suptitle("Effets des transformations sur des revenus log-normaux", y=1.01)
plt.tight_layout()
plt.show()
_images/6d97e22ff89ac8d613cc6f050489d6b59a854ab838dc866d1a820e00b07f2abc.png

Quand transformer ?

  • Log : données strictement positives avec queue droite lourde (revenus, prix, concentrations chimiques).

  • Box-Cox : données strictement positives, λ estimé par MLE — optimal pour normaliser.

  • Yeo-Johnson : idem mais accepte les valeurs nulles ou négatives.

  • Racine carrée : données de comptage (Poisson) — variance proportionnelle à la moyenne.

  • Arcsine : proportions entre 0 et 1.

Interprétez toujours les résultats dans l’espace original (dé-transformer). Les coefficients d’une régression log-linéaire s’interprètent comme des effets multiplicatifs.

Mélange de gaussiennes (GMM)#

Certaines distributions semblent bimodales parce qu’elles sont un mélange de deux (ou plus) populations gaussiennes. Le Gaussian Mixture Model (GMM) décompose cette distribution mixte.

# Exemple : temps de réponse d'un serveur web (bimodal)
# Requêtes normales + requêtes lentes (cache miss)
t_normales = rng.normal(loc=50, scale=10, size=350)   # ms
t_lentes = rng.normal(loc=200, scale=40, size=150)    # ms
temps_reponse = np.concatenate([t_normales, t_lentes])

# Ajustement GMM "à la main" via scipy
from scipy.stats import norm

def gmm_pdf(x, pi1, mu1, sigma1, mu2, sigma2):
    return pi1 * norm.pdf(x, mu1, sigma1) + (1 - pi1) * norm.pdf(x, mu2, sigma2)

# Estimation simple par inspection visuelle / EM simplifié
# (sklearn.mixture.GaussianMixture fait ça proprement en multidimensionnel)
from scipy.optimize import minimize

def neg_log_vraisemblance(params):
    pi1, mu1, log_s1, mu2, log_s2 = params
    if not (0.01 < pi1 < 0.99):
        return 1e10
    s1, s2 = np.exp(log_s1), np.exp(log_s2)
    pdf = gmm_pdf(temps_reponse, pi1, mu1, s1, mu2, s2)
    pdf = np.clip(pdf, 1e-300, None)
    return -np.sum(np.log(pdf))

x0 = [0.7, 50, np.log(10), 200, np.log(40)]
res = minimize(neg_log_vraisemblance, x0, method="Nelder-Mead",
               options={"maxiter": 10000, "xatol": 1e-6})
pi1, mu1, log_s1, mu2, log_s2 = res.x
s1, s2 = np.exp(log_s1), np.exp(log_s2)
print(f"Composante 1 : π={pi1:.2f}, μ={mu1:.1f} ms, σ={s1:.1f} ms")
print(f"Composante 2 : π={1-pi1:.2f}, μ={mu2:.1f} ms, σ={s2:.1f} ms")
Composante 1 : π=0.70, μ=49.5 ms, σ=10.0 ms
Composante 2 : π=0.30, μ=194.2 ms, σ=40.7 ms

Hide code cell source

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

# Histogramme avec décomposition GMM
ax = axes[0]
x_plot = np.linspace(temps_reponse.min(), temps_reponse.max(), 400)
ax.hist(temps_reponse, bins=40, density=True, alpha=0.45, color="gray", label="Données")
ax.plot(x_plot, stats.gaussian_kde(temps_reponse)(x_plot), "k-", lw=2, label="KDE total")
ax.fill_between(x_plot, pi1 * norm.pdf(x_plot, mu1, s1),
                alpha=0.4, color="steelblue", label=f"C1: μ={mu1:.0f}ms, π={pi1:.0%}")
ax.fill_between(x_plot, (1-pi1) * norm.pdf(x_plot, mu2, s2),
                alpha=0.4, color="tomato", label=f"C2: μ={mu2:.0f}ms, π={1-pi1:.0%}")
ax.plot(x_plot, gmm_pdf(x_plot, pi1, mu1, s1, mu2, s2), "k--", lw=2, label="GMM total")
ax.set_xlabel("Temps de réponse (ms)")
ax.set_ylabel("Densité")
ax.set_title("Décomposition GMM : temps de réponse serveur")
ax.legend(fontsize=8)

# QQ-plot : données bimodales vs gaussienne
ax2 = axes[1]
(osm, osr), (slope, intercept, r) = stats.probplot(temps_reponse, dist="norm")
ax2.scatter(osm, osr, alpha=0.4, s=10, color="steelblue")
x_line = np.array([osm.min(), osm.max()])
ax2.plot(x_line, slope * x_line + intercept, "r-", lw=2)
ax2.annotate("Signature bimodale :\ncassure centrale", xy=(0, 150),
             xytext=(1, 80), fontsize=8.5,
             arrowprops=dict(arrowstyle="->", color="tomato"),
             color="tomato")
ax2.set_xlabel("Quantiles théoriques (normale)")
ax2.set_ylabel("Quantiles observés")
ax2.set_title("QQ-plot : distribution bimodale vs normale")

plt.tight_layout()
plt.show()
_images/c232a08a68514d6b954b2686db2f07b2fd5154946843ca89fd47f9f44dcede98.png

GMM en pratique

sklearn.mixture.GaussianMixture offre une implémentation robuste du GMM par algorithme EM (Expectation-Maximisation) pour des données multivariées. On choisit le nombre de composantes \(K\) par critère BIC. Le GMM est aussi utilisé comme méthode de clustering doux : chaque point appartient à chaque composante avec une probabilité. Il sera étudié en détail dans le livre Data Science.

Les outils de ce chapitre — QQ-plots, tests de normalité, ajustements de lois, transformations — constituent la boîte à outils de diagnostic indispensable avant toute modélisation. Connaître la loi (ou l’approximation raisonnable) d’une variable permet de choisir les bons tests, d’interpréter correctement les résultats, et d’éviter les erreurs classiques de l’analyse de données.