---
jupytext:
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.16.0
kernelspec:
  name: python3
  display_name: Python 3
---

# Statistiques descriptives

```{code-cell} python
:tags: [hide-input]

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

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

Les statistiques descriptives ont un objectif modeste en apparence — résumer un jeu de données en quelques chiffres ou graphiques — mais cet objectif est en réalité fondamental. Avant tout test, toute modélisation, toute inférence, il faut *comprendre* ses données : leur forme, leurs valeurs atypiques, leur dispersion, leur asymétrie. Un modèle appliqué à des données mal comprises produit des résultats mal compris. Ce chapitre couvre les outils essentiels pour explorer des données numériques, en partant d'un jeu de données concret et en progressant vers les diagnostics avancés.

## Un jeu de données réaliste : salaires dans une entreprise

Imaginons une entreprise de 500 employés. Les salaires annuels (en milliers d'euros) suivent une distribution réaliste : la grande majorité tourne autour de 35 à 60 k€, avec quelques cadres supérieurs et dirigeants qui tirent la queue droite vers le haut.

```{code-cell} python
# Génération de salaires réalistes : mélange de deux populations
n_employes = 400
n_cadres = 80
n_dirigeants = 20

salaires_employes = rng.normal(loc=38, scale=8, size=n_employes)
salaires_cadres = rng.normal(loc=65, scale=12, size=n_cadres)
salaires_dirigeants = rng.normal(loc=130, scale=30, size=n_dirigeants)

salaires = np.concatenate([salaires_employes, salaires_cadres, salaires_dirigeants])
salaires = np.clip(salaires, 22, 350)  # salaires plausibles

departements = (
    ["Technique"] * 150 + ["Ventes"] * 120 + ["RH"] * 80 + ["Finance"] * 50 +
    ["Technique"] * 30 + ["Ventes"] * 20 + ["Finance"] * 20 + ["RH"] * 10 +
    ["Direction"] * 20
)

df = pd.DataFrame({
    "salaire": salaires,
    "departement": departements,
    "anciennete": rng.integers(1, 25, size=500),
    "age": np.clip(rng.normal(38, 9, 500).astype(int), 22, 62),
})

print(df.shape)
df.head()
```

## Tendance centrale : quelle est la valeur « typique » ?

Trois mesures coexistent pour résumer le centre d'une distribution. Elles coïncident pour une distribution symétrique, mais divergent dès qu'apparaît de l'asymétrie ou des valeurs aberrantes.

**La moyenne arithmétique** $\bar{x} = \frac{1}{n}\sum_{i=1}^n x_i$ est sensible à tous les points, y compris les extrêmes. Elle minimise la somme des carrés des écarts.

**La médiane** est la valeur qui partage l'échantillon en deux moitiés égales. Elle est robuste aux outliers : on peut multiplier le salaire du PDG par dix sans déplacer la médiane d'un centime.

**Le mode** est la valeur la plus fréquente. Utile pour des données catégorielles ou discrètes ; pour des données continues, on parle plutôt de pic de la densité estimée.

```{code-cell} python
moyenne = df["salaire"].mean()
mediane = df["salaire"].median()
mode_val = float(stats.mode(df["salaire"].round(), keepdims=True).mode[0])

print(f"Moyenne  : {moyenne:.1f} k€")
print(f"Médiane  : {mediane:.1f} k€")
print(f"Mode (arrondi) : {mode_val:.0f} k€")
print(f"Écart moyenne-médiane : {moyenne - mediane:.1f} k€  → asymétrie droite")
```

```{code-cell} python
:tags: [hide-input]

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

# Histogramme avec KDE et indicateurs
ax = axes[0]
ax.hist(df["salaire"], bins=40, density=True, alpha=0.55, color="steelblue",
        label="Distribution observée")
kde_x = np.linspace(df["salaire"].min(), df["salaire"].max(), 300)
kde = stats.gaussian_kde(df["salaire"])
ax.plot(kde_x, kde(kde_x), color="steelblue", lw=2, label="Densité estimée (KDE)")
ax.axvline(moyenne, color="tomato", lw=2, linestyle="--", label=f"Moyenne = {moyenne:.0f}")
ax.axvline(mediane, color="seagreen", lw=2, linestyle="-.", label=f"Médiane = {mediane:.0f}")
ax.set_xlabel("Salaire (k€)")
ax.set_ylabel("Densité")
ax.set_title("Distribution des salaires — tendance centrale")
ax.legend(fontsize=9)

# Impact des outliers sur la moyenne
ax2 = axes[1]
n_sim = 200
moyennes_sans = []
moyennes_avec = []
for _ in range(n_sim):
    echantillon = rng.choice(salaires_employes, size=50, replace=False)
    moyennes_sans.append(echantillon.mean())
    avec_outlier = np.append(echantillon, 350)  # un dirigeant exceptionnel
    moyennes_avec.append(avec_outlier.mean())

ax2.hist(moyennes_sans, bins=25, alpha=0.6, color="seagreen", density=True,
         label="Sans outlier")
ax2.hist(moyennes_avec, bins=25, alpha=0.6, color="tomato", density=True,
         label="Avec outlier à 350 k€")
ax2.axvline(np.median(salaires_employes), color="gray", lw=2, linestyle=":",
            label=f"Médiane stable ≈ {np.median(salaires_employes):.0f}")
ax2.set_xlabel("Moyenne estimée (k€)")
ax2.set_ylabel("Densité")
ax2.set_title("Robustesse : médiane vs moyenne face aux outliers")
ax2.legend(fontsize=9)

plt.tight_layout()
plt.show()
```

La figure de droite illustre un fait fondamental : la moyenne fluctue énormément en présence d'un seul outlier, alors que la médiane reste stable. Pour des données salariales, il est donc préférable de rapporter la **médiane** comme mesure de tendance centrale.

```{admonition} Règle pratique
:class: tip
Choisissez la **médiane** dès que la distribution est asymétrique ou contient des valeurs aberrantes probables : revenus, prix immobiliers, durées de vie, tailles de fichiers. Réservez la **moyenne** aux distributions approximativement symétriques et bornées.
```

## Dispersion : à quel point les valeurs s'étalent-elles ?

La dispersion mesure la variabilité. Connaître la tendance centrale sans la dispersion ne suffit pas : deux entreprises peuvent avoir la même masse salariale médiane avec des inégalités très différentes.

**Variance et écart-type.** La variance $s^2 = \frac{1}{n-1}\sum (x_i - \bar{x})^2$ est la moyenne des carrés des écarts à la moyenne (divisée par $n-1$ pour l'estimateur sans biais). L'écart-type $s$ est sa racine carrée, dans la même unité que les données.

**IQR (interquartile range).** L'IQR = Q3 − Q1 est la longueur de la boîte du boxplot. Il mesure l'étendue des 50 % centraux et est insensible aux valeurs extrêmes.

**MAD (median absolute deviation).** $\text{MAD} = \text{médiane}(|x_i - \text{médiane}(x)|)$. Encore plus robuste que l'IQR, c'est l'estimateur de dispersion de référence pour les données contaminées. On le normalise parfois par 1,4826 pour l'aligner sur l'écart-type d'une gaussienne.

```{code-cell} python
s = df["salaire"]

variance = s.var()          # ddof=1 par défaut dans pandas
ecart_type = s.std()
iqr = s.quantile(0.75) - s.quantile(0.25)
mad = stats.median_abs_deviation(s)
mad_normalise = stats.median_abs_deviation(s, scale="normal")  # ≈ σ pour gaussienne

print(f"Variance          : {variance:.1f} k€²")
print(f"Écart-type        : {ecart_type:.1f} k€")
print(f"IQR               : {iqr:.1f} k€")
print(f"MAD               : {mad:.1f} k€")
print(f"MAD normalisé     : {mad_normalise:.1f} k€  (comparable à l'écart-type)")
print(f"\nRatio σ/MAD_norm  : {ecart_type/mad_normalise:.2f}  (> 1 → queue droite lourde)")
```

Le ratio écart-type / MAD normalisé supérieur à 1 indique que l'écart-type est gonflé par les valeurs extrêmes. C'est une signature d'asymétrie ou de queues lourdes.

## Forme de la distribution : asymétrie et aplatissement

Deux distributions peuvent avoir la même moyenne et le même écart-type tout en différant radicalement de forme. Le **skewness** (asymétrie) et le **kurtosis** (aplatissement) capturent ces différences.

**Skewness.** Le coefficient d'asymétrie $\gamma_1 = \frac{\mu_3}{\sigma^3}$ où $\mu_3 = E[(X-\mu)^3]$ est le moment centré d'ordre 3. Un skewness positif signifie une queue étalée à droite (revenus, temps d'attente) ; négatif, à gauche (âge au décès dans un pays développé).

**Kurtosis.** Le coefficient d'aplatissement $\gamma_2 = \frac{\mu_4}{\sigma^4} - 3$ (excès de kurtosis, nul pour la gaussienne). Un kurtosis positif (leptokurtique) indique des queues plus épaisses et un pic plus prononcé que la normale. Attention : `scipy.stats.kurtosis` retourne l'excès de kurtosis par défaut.

```{code-cell} python
skewness = stats.skew(df["salaire"])
kurt = stats.kurtosis(df["salaire"])  # excès de kurtosis

print(f"Skewness  : {skewness:.3f}")
print(f"Kurtosis (excès) : {kurt:.3f}")
print()
print("Interprétation :")
print(f"  → Skewness {skewness:.2f} > 0 : queue droite étalée (salaires élevés rares mais présents)")
print(f"  → Kurtosis {kurt:.2f} > 0 : queues plus lourdes qu'une gaussienne")
```

```{code-cell} python
:tags: [hide-input]

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# Comparaison de skewness
ax = axes[0]
x = np.linspace(-4, 8, 300)
for skew_param, label, color in [
    (0, "Symétrique (skew=0)", "steelblue"),
    (4, "Asymétrie droite (skew>0)", "tomato"),
    (-4, "Asymétrie gauche (skew<0)", "seagreen"),
]:
    a_param = skew_param
    dist = stats.skewnorm(a=a_param, loc=2, scale=1.2)
    ax.plot(x, dist.pdf(x), lw=2, label=label, color=color)
ax.set_title("Effets du skewness")
ax.set_xlabel("Valeur")
ax.set_ylabel("Densité")
ax.legend(fontsize=8)

# Comparaison de kurtosis
ax2 = axes[1]
x2 = np.linspace(-6, 6, 300)
for df_t, label, color in [
    (2.1, "Queues lourdes (t, ν=2)", "tomato"),
    (10, "Queues légères (t, ν=10)", "steelblue"),
    (None, "Gaussienne (kurtosis=0)", "seagreen"),
]:
    if df_t is None:
        ax2.plot(x2, stats.norm.pdf(x2), lw=2, label=label, color=color)
    else:
        ax2.plot(x2, stats.t.pdf(x2, df=df_t), lw=2, label=label, color=color)
ax2.set_title("Effets du kurtosis (excès)")
ax2.set_xlabel("Valeur")
ax2.legend(fontsize=8)

# Salaires : skewness et kurtosis réels
ax3 = axes[2]
ax3.hist(df["salaire"], bins=35, density=True, alpha=0.5, color="steelblue")
kde_x = np.linspace(df["salaire"].min(), df["salaire"].max(), 300)
ax3.plot(kde_x, stats.gaussian_kde(df["salaire"])(kde_x), color="steelblue", lw=2)
# Gaussienne de référence
x_ref = np.linspace(10, 200, 300)
ax3.plot(x_ref,
         stats.norm.pdf(x_ref, df["salaire"].mean(), df["salaire"].std()),
         "r--", lw=1.5, label="Gaussienne de référence")
ax3.set_xlabel("Salaire (k€)")
ax3.set_title(f"Salaires : skew={skewness:.2f}, kurt={kurt:.2f}")
ax3.legend(fontsize=8)

plt.tight_layout()
plt.show()
```

## Quantiles, percentiles et boxplot

Les **quantiles** découpent la distribution en tranches d'égale probabilité. Le quantile d'ordre $p$ est la valeur $Q_p$ telle que $P(X \leq Q_p) = p$. Les quartiles ($Q_{0.25}$, $Q_{0.5}$, $Q_{0.75}$) sont les plus utilisés.

Le **boxplot** (boîte à moustaches) résume une distribution en cinq chiffres : minimum, Q1, médiane, Q3, maximum — avec un traitement visuel des outliers. La moustache s'étend jusqu'à $\text{Q1} - 1{,}5 \times \text{IQR}$ (en bas) et $\text{Q3} + 1{,}5 \times \text{IQR}$ (en haut). Les points au-delà sont tracés individuellement.

```{code-cell} python
quantiles = df["salaire"].quantile([0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99])
print("Quantiles des salaires :")
for q, v in quantiles.items():
    print(f"  P{int(q*100):2d} : {v:.1f} k€")
```

```{code-cell} python
:tags: [hide-input]

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

# Boxplot comparatif par département
ax = axes[0]
order = df.groupby("departement")["salaire"].median().sort_values(ascending=False).index
sns.boxplot(data=df, y="departement", x="salaire", order=order, hue="departement", palette="muted", legend=False, ax=ax)
ax.set_xlabel("Salaire (k€)")
ax.set_ylabel("")
ax.set_title("Distribution des salaires par département")
# Annoter IQR
for i, dept in enumerate(order):
    sub = df[df["departement"] == dept]["salaire"]
    q1, q3 = sub.quantile(0.25), sub.quantile(0.75)
    ax.annotate(f"IQR={q3-q1:.0f}", xy=(q3, i), xytext=(q3+5, i),
                fontsize=7.5, va="center", color="gray")

# Violin plot : densité complète + boxplot intégré
ax2 = axes[1]
sns.violinplot(data=df, y="departement", x="salaire", order=order,
               inner="box", hue="departement", palette="pastel", legend=False, ax=ax2)
ax2.set_xlabel("Salaire (k€)")
ax2.set_ylabel("")
ax2.set_title("Violin plot : densité + quartiles")

plt.tight_layout()
plt.show()
```

```{admonition} Boxplot vs violin plot
:class: note
Le **boxplot** est lisible et compact : il fonctionne bien pour comparer beaucoup de groupes. Le **violin plot** montre la forme de la distribution (multimodalité, asymétrie) mais est plus difficile à lire avec de nombreux groupes. Pour une distribution bimodale, le boxplot n'en révèle rien là où le violin plot la révèle immédiatement.
```

## Détection d'outliers

Un outlier n'est pas forcément une erreur de mesure — c'est une valeur qui se distingue du comportement typique. Sa détection dépend du contexte et de la méthode choisie.

### Règle de l'IQR (Tukey)

Un point $x$ est déclaré outlier si $x < Q_1 - 1{,}5 \times \text{IQR}$ ou $x > Q_3 + 1{,}5 \times \text{IQR}$. Simple et non paramétrique, cette règle est celle utilisée par défaut dans les boxplots.

### Z-score classique

$z = \frac{x - \bar{x}}{s}$. Convention : $|z| > 3$ signale un outlier. Cette méthode suppose une distribution approximativement gaussienne et est sensible aux outliers eux-mêmes (le masquage).

### Z-score modifié (Iglewicz & Hoaglin)

$z_{\text{mod}} = \frac{0{,}6745 \times (x - \text{médiane})}{\text{MAD}}$. Le seuil habituel est $|z_{\text{mod}}| > 3{,}5$. Beaucoup plus robuste : ni la médiane ni la MAD ne sont affectées par les outliers.

```{code-cell} python
def detecter_outliers(x):
    """Trois méthodes de détection d'outliers."""
    x = np.asarray(x)

    # IQR
    q1, q3 = np.percentile(x, 25), np.percentile(x, 75)
    iqr_val = q3 - q1
    outliers_iqr = (x < q1 - 1.5 * iqr_val) | (x > q3 + 1.5 * iqr_val)

    # Z-score classique
    z = np.abs(stats.zscore(x))
    outliers_z = z > 3

    # Z-score modifié
    med = np.median(x)
    mad_val = stats.median_abs_deviation(x)
    z_mod = np.abs(0.6745 * (x - med) / mad_val) if mad_val > 0 else np.zeros_like(x)
    outliers_zmod = z_mod > 3.5

    return outliers_iqr, outliers_z, outliers_zmod

oi, oz, ozm = detecter_outliers(df["salaire"])

print(f"Méthode IQR          : {oi.sum()} outliers ({oi.mean()*100:.1f}%)")
print(f"Z-score classique    : {oz.sum()} outliers ({oz.mean()*100:.1f}%)")
print(f"Z-score modifié (MAD): {ozm.sum()} outliers ({ozm.mean()*100:.1f}%)")
print()
print("Outliers IQR (salaires en k€) :")
print(np.sort(df["salaire"][oi].values)[:10], "...")
```

```{code-cell} python
:tags: [hide-input]

fig, axes = plt.subplots(1, 3, figsize=(15, 4))
methodes = [("IQR", oi), ("Z-score", oz), ("Z-score modifié", ozm)]

for ax, (nom, masque) in zip(axes, methodes):
    ax.scatter(range(len(df)), np.sort(df["salaire"]),
               c=["tomato" if m else "steelblue" for m in masque[np.argsort(df["salaire"])]],
               alpha=0.5, s=15)
    ax.set_title(f"{nom} — {masque.sum()} outliers")
    ax.set_xlabel("Rang")
    ax.set_ylabel("Salaire (k€)")

plt.suptitle("Comparaison des méthodes de détection d'outliers", y=1.02)
plt.tight_layout()
plt.show()
```

```{admonition} Isolation Forest et méthodes non paramétriques
:class: note
Pour des données multivariées, les méthodes univariées sont insuffisantes : un point peut être normal sur chaque dimension séparément mais aberrant par leur combinaison. **Isolation Forest** (`sklearn.ensemble.IsolationForest`) isole les outliers en les coupant aléatoirement dans l'espace des features ; ils nécessitent moins de coupes que les points normaux. **Local Outlier Factor (LOF)** compare la densité locale d'un point à celle de ses voisins. Ces méthodes seront couvertes dans le livre *Data Science*.
```

## Résumé complet avec `describe()` enrichi

La méthode `describe()` de pandas donne un résumé rapide, mais elle peut être enrichie pour inclure skewness, kurtosis et MAD.

```{code-cell} python
def describe_enrichi(series, nom="variable"):
    """Résumé statistique enrichi."""
    s = series.dropna()
    q = s.quantile
    mad_v = stats.median_abs_deviation(s)

    resume = {
        "n": len(s),
        "n_manquants": series.isna().sum(),
        "min": s.min(),
        "P5": q(0.05),
        "Q1 (P25)": q(0.25),
        "médiane": q(0.50),
        "Q3 (P75)": q(0.75),
        "P95": q(0.95),
        "max": s.max(),
        "moyenne": s.mean(),
        "écart-type": s.std(),
        "IQR": q(0.75) - q(0.25),
        "MAD": mad_v,
        "skewness": stats.skew(s),
        "kurtosis (excès)": stats.kurtosis(s),
    }
    return pd.Series(resume, name=nom)

desc = describe_enrichi(df["salaire"], "salaire (k€)")
print(desc.to_string())
```

## Statistiques groupées

L'analyse par groupe est une des opérations les plus fréquentes en data analyse. `groupby` de pandas permet de calculer n'importe quelle statistique par sous-groupe.

```{code-cell} python
# Statistiques par département
stats_groupe = df.groupby("departement")["salaire"].agg([
    ("n", "count"),
    ("médiane", "median"),
    ("moyenne", "mean"),
    ("écart-type", "std"),
    ("IQR", lambda x: x.quantile(0.75) - x.quantile(0.25)),
    ("skewness", stats.skew),
    ("P90", lambda x: x.quantile(0.90)),
]).round(1)

print(stats_groupe.sort_values("médiane", ascending=False).to_string())
```

```{code-cell} python
:tags: [hide-input]

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

# Médiane et moyenne par département
ax = axes[0, 0]
x_pos = np.arange(len(stats_groupe))
dept_order = stats_groupe.sort_values("médiane", ascending=True).index
med_vals = stats_groupe.loc[dept_order, "médiane"]
moy_vals = stats_groupe.loc[dept_order, "moyenne"]
ax.barh(x_pos, med_vals.values, alpha=0.7, color="seagreen", label="Médiane")
ax.barh(x_pos, moy_vals.values, alpha=0.5, color="tomato", label="Moyenne")
ax.set_yticks(x_pos)
ax.set_yticklabels(dept_order)
ax.set_xlabel("Salaire (k€)")
ax.set_title("Médiane vs Moyenne par département")
ax.legend()

# Boxplot par département avec swarmplot
ax2 = axes[0, 1]
dept_order2 = df.groupby("departement")["salaire"].median().sort_values(ascending=False).index
sns.boxplot(data=df, x="departement", y="salaire", order=dept_order2,
            hue="departement", palette="Set2", legend=False, width=0.5, ax=ax2)
sns.stripplot(data=df, x="departement", y="salaire", order=dept_order2,
              color="black", alpha=0.2, size=3, jitter=True, ax=ax2)
ax2.tick_params(axis="x", rotation=20)
ax2.set_xlabel("")
ax2.set_ylabel("Salaire (k€)")
ax2.set_title("Distribution par département")

# Salaire vs ancienneté (corrélation brute)
ax3 = axes[1, 0]
scatter = ax3.scatter(df["anciennete"], df["salaire"],
                      c=df["age"], cmap="YlOrRd", alpha=0.5, s=20)
plt.colorbar(scatter, ax=ax3, label="Âge")
ax3.set_xlabel("Ancienneté (années)")
ax3.set_ylabel("Salaire (k€)")
ax3.set_title("Salaire vs ancienneté (couleur = âge)")
r, p = stats.pearsonr(df["anciennete"], df["salaire"])
ax3.annotate(f"r = {r:.2f} (p={p:.3f})", xy=(0.05, 0.92), xycoords="axes fraction",
             fontsize=9)

# Distribution cumulée (ECDF)
ax4 = axes[1, 1]
for dept in dept_order2[:4]:
    sub = df[df["departement"] == dept]["salaire"].sort_values()
    ecdf = np.arange(1, len(sub) + 1) / len(sub)
    ax4.step(sub, ecdf, where="post", lw=1.8, label=dept)
ax4.set_xlabel("Salaire (k€)")
ax4.set_ylabel("Fréquence cumulée")
ax4.set_title("ECDF par département")
ax4.legend(fontsize=9)
ax4.axhline(0.5, color="gray", lw=1, linestyle=":")
ax4.axhline(0.25, color="gray", lw=1, linestyle=":")
ax4.axhline(0.75, color="gray", lw=1, linestyle=":")

plt.tight_layout()
plt.show()
```

## Tableau récapitulatif des mesures

```{code-cell} python
:tags: [hide-input]

# Tableau synthétique des mesures et de leurs propriétés
synthese = pd.DataFrame({
    "Mesure": [
        "Moyenne", "Médiane", "Mode",
        "Variance", "Écart-type", "IQR", "MAD",
        "Skewness", "Kurtosis (excès)"
    ],
    "Catégorie": [
        "Tendance centrale", "Tendance centrale", "Tendance centrale",
        "Dispersion", "Dispersion", "Dispersion", "Dispersion",
        "Forme", "Forme"
    ],
    "Robustesse": [
        "Faible", "Forte", "Forte",
        "Faible", "Faible", "Forte", "Très forte",
        "Faible", "Faible"
    ],
    "Quand l'utiliser": [
        "Distribution symétrique, sans outliers extrêmes",
        "Distribution asymétrique, revenus, prix",
        "Variables discrètes, catégorielles",
        "Calculs théoriques, ANOVA",
        "Distribution ≈ normale",
        "Boxplots, détection d'outliers Tukey",
        "Données contaminées, outliers probables",
        "Diagnostic d'asymétrie",
        "Diagnostic de queues lourdes"
    ]
}).set_index("Mesure")

print(synthese.to_string())
```

Le choix de la bonne mesure statistique n'est pas anodin. Rapporter la moyenne des salaires d'une entreprise peut être politiquement trompeur si quelques dirigeants très bien payés la tirent vers le haut — la médiane révèle alors une réalité plus fidèle. Comprendre ces nuances, c'est déjà une grande partie du travail de statisticien.