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

# Estimation et intervalles de confiance

```{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
from scipy.stats import bootstrap as scipy_bootstrap

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

Estimer une quantité à partir d'un échantillon, c'est inévitablement accepter une incertitude. Le **point estimé** — la valeur calculée sur l'échantillon — n'est qu'une réalisation parmi toutes celles qu'on aurait pu obtenir avec d'autres données. L'**intervalle de confiance** quantifie cette incertitude de façon opérationnelle. Ce chapitre couvre les méthodes classiques et computationnelles pour construire ces intervalles.

## Rappel : estimateur, biais, variance, MSE

Un **estimateur** $\hat{\theta}$ d'un paramètre $\theta$ est une statistique — une fonction des données — qui approche $\theta$. Sa qualité se mesure par :

- **Biais** : $B(\hat{\theta}) = E[\hat{\theta}] - \theta$. Un estimateur sans biais satisfait $E[\hat{\theta}] = \theta$.
- **Variance** : $\text{Var}(\hat{\theta})$ — variabilité de l'estimateur d'un échantillon à l'autre.
- **MSE** (Mean Squared Error) : $\text{MSE}(\hat{\theta}) = \text{Var}(\hat{\theta}) + B(\hat{\theta})^2$ — combine les deux.

```{admonition} Biais vs variance
:class: note
Il existe un arbitrage fondamental entre biais et variance. La moyenne empirique est sans biais pour $\mu$ ; l'estimateur shrinkage de James-Stein est biaisé mais peut avoir une MSE plus faible. Pour les démonstrations rigoureuses de ces propriétés, voir le livre *Mathématiques*, chapitre sur les statistiques mathématiques.
```

```{code-cell} python
# Illustration : biais de l'estimateur de variance
# s² avec n-1 est sans biais ; avec n est biaisé
sigma_vrai = 5.0
mu_vrai = 10.0
n_obs = 15

n_sim = 10000
var_n = []
var_n1 = []

for _ in range(n_sim):
    x = rng.normal(mu_vrai, sigma_vrai, n_obs)
    var_n.append(np.var(x, ddof=0))   # diviseur n (biaisé)
    var_n1.append(np.var(x, ddof=1))  # diviseur n-1 (sans biais)

print(f"Variance vraie σ²    : {sigma_vrai**2:.4f}")
print(f"E[s²_n]   (ddof=0)  : {np.mean(var_n):.4f}  → biais = {np.mean(var_n)-sigma_vrai**2:.4f}")
print(f"E[s²_n-1] (ddof=1)  : {np.mean(var_n1):.4f}  → biais = {np.mean(var_n1)-sigma_vrai**2:.4f}")
```

## Bootstrap : principe du rééchantillonnage

Le **bootstrap** (Efron, 1979) est une méthode computationnelle d'estimation de la distribution d'un estimateur. L'idée centrale : l'échantillon observé est la meilleure approximation disponible de la population. On peut donc **rééchantillonner avec remise** pour simuler la variabilité.

**Algorithme :**

1. À partir de l'échantillon $(x_1, \ldots, x_n)$, tirer $B$ échantillons bootstrap $(x_1^*, \ldots, x_n^*)$ avec remise.
2. Calculer l'estimateur $\hat{\theta}^*_b$ sur chaque échantillon.
3. La distribution des $\{\hat{\theta}^*_b\}$ approxime la distribution de $\hat{\theta}$.

```{code-cell} python
# Bootstrap from scratch
def bootstrap(data, statistic, n_boot=2000, random_state=None):
    """
    Bootstrap non paramétrique.

    Parameters
    ----------
    data : array-like
    statistic : callable — fonction(array) → scalaire
    n_boot : int — nombre de rééchantillons
    random_state : int ou Generator

    Returns
    -------
    array des statistiques bootstrap
    """
    rng_b = np.random.default_rng(random_state)
    data = np.asarray(data)
    n = len(data)
    boots = np.array([
        statistic(rng_b.choice(data, size=n, replace=True))
        for _ in range(n_boot)
    ])
    return boots

# Données : temps de chargement d'une page web (ms)
data_web = rng.lognormal(mean=np.log(250), sigma=0.6, size=80)
print(f"n = {len(data_web)}, médiane = {np.median(data_web):.1f} ms")
print(f"Skewness = {stats.skew(data_web):.2f}  → distribution asymétrique")
```

### Bootstrap percentile

La méthode la plus simple : les bornes de l'IC sont les quantiles $\alpha/2$ et $1-\alpha/2$ de la distribution bootstrap.

### Bootstrap BCa (Bias-Corrected and Accelerated)

Le bootstrap percentile peut être biaisé et sous-optimal. Le **BCa** (Efron & Tibshirani, 1993) corrige :

1. **Correction de biais** $z_0$ : proportion des valeurs bootstrap inférieures à l'estimateur observé.
2. **Accélération** $a$ : mesure la vitesse de variation de l'écart-type de l'estimateur.

C'est l'IC bootstrap recommandé pour la plupart des applications.

```{code-cell} python
def bootstrap_ic(data, statistic, n_boot=2000, alpha=0.05, methode="bca"):
    """
    IC bootstrap : méthode percentile ou BCa.

    Returns
    -------
    ic_low, ic_high, boots
    """
    boots = bootstrap(data, statistic, n_boot=n_boot, random_state=0)
    theta_hat = statistic(data)

    if methode == "percentile":
        ic_low = np.percentile(boots, 100 * alpha / 2)
        ic_high = np.percentile(boots, 100 * (1 - alpha / 2))

    elif methode == "bca":
        # Correction de biais z0
        z0 = stats.norm.ppf(np.mean(boots < theta_hat) + 1e-10)

        # Accélération a (jackknife)
        n = len(data)
        jk_vals = np.array([statistic(np.delete(data, i)) for i in range(n)])
        jk_mean = jk_vals.mean()
        num = np.sum((jk_mean - jk_vals)**3)
        den = 6 * (np.sum((jk_mean - jk_vals)**2))**1.5
        a = num / den if den != 0 else 0

        # Quantiles ajustés
        z_alpha_lo = stats.norm.ppf(alpha / 2)
        z_alpha_hi = stats.norm.ppf(1 - alpha / 2)

        def adj_quantile(z_in):
            num_q = z0 + z_in
            return stats.norm.cdf(z0 + num_q / (1 - a * num_q))

        q_lo = adj_quantile(z_alpha_lo)
        q_hi = adj_quantile(z_alpha_hi)
        q_lo = np.clip(q_lo, 0.001, 0.999)
        q_hi = np.clip(q_hi, 0.001, 0.999)

        ic_low = np.percentile(boots, 100 * q_lo)
        ic_high = np.percentile(boots, 100 * q_hi)

    return ic_low, ic_high, boots

# Calcul des IC pour la médiane (temps de chargement)
ic_p_low, ic_p_high, boots_perc = bootstrap_ic(data_web, np.median, methode="percentile")
ic_b_low, ic_b_high, boots_bca = bootstrap_ic(data_web, np.median, methode="bca")

print("IC bootstrap pour la médiane des temps de chargement :")
print(f"  Percentile    : [{ic_p_low:.1f}, {ic_p_high:.1f}] ms")
print(f"  BCa           : [{ic_b_low:.1f}, {ic_b_high:.1f}] ms")
print(f"  Médiane obs.  : {np.median(data_web):.1f} ms")
```

```{code-cell} python
# Comparaison IC bootstrap pour plusieurs statistiques
statistiques = {
    "Moyenne": np.mean,
    "Médiane": np.median,
    "Écart-type": np.std,
    "P90": lambda x: np.percentile(x, 90),
    "Skewness": stats.skew,
}

print(f"{'Statistique':15s} {'Estimé':>10} {'IC perc. bas':>14} {'IC perc. haut':>14} {'IC BCa bas':>12} {'IC BCa haut':>12}")
print("-" * 82)
for nom, stat in statistiques.items():
    val = stat(data_web)
    ic_p_l, ic_p_h, _ = bootstrap_ic(data_web, stat, methode="percentile")
    ic_b_l, ic_b_h, _ = bootstrap_ic(data_web, stat, methode="bca")
    print(f"{nom:15s} {val:>10.2f} {ic_p_l:>14.2f} {ic_p_h:>14.2f} {ic_b_l:>12.2f} {ic_b_h:>12.2f}")
```

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

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

# Distribution bootstrap de la médiane
ax = axes[0]
boots_med = bootstrap(data_web, np.median, n_boot=300, random_state=42)
theta_obs = np.median(data_web)

ax.hist(boots_med, bins=50, density=True, alpha=0.55, color="steelblue",
        label="Distribution bootstrap (médiane)")
ax.axvline(theta_obs, color="tomato", lw=2.5, label=f"Médiane obs. = {theta_obs:.1f}")
ax.axvline(ic_p_low, color="seagreen", lw=2, linestyle="--", label=f"IC perc. [{ic_p_low:.1f}, {ic_p_high:.1f}]")
ax.axvline(ic_p_high, color="seagreen", lw=2, linestyle="--")
ax.axvline(ic_b_low, color="darkorange", lw=2, linestyle=":", label=f"IC BCa [{ic_b_low:.1f}, {ic_b_high:.1f}]")
ax.axvline(ic_b_high, color="darkorange", lw=2, linestyle=":")

kde_x = np.linspace(boots_med.min(), boots_med.max(), 300)
ax.plot(kde_x, stats.gaussian_kde(boots_med)(kde_x), "k-", lw=1.5)
ax.set_xlabel("Médiane bootstrap (ms)")
ax.set_ylabel("Densité")
ax.set_title("Distribution bootstrap de la médiane\n(temps de chargement)")
ax.legend(fontsize=8)

# Comparaison IC bootstrap vs IC analytique pour la moyenne
ax2 = axes[1]
boots_moy = bootstrap(data_web, np.mean, n_boot=300, random_state=42)
mu_obs = np.mean(data_web)
se_obs = np.std(data_web, ddof=1) / np.sqrt(len(data_web))
t_crit = stats.t.ppf(0.975, df=len(data_web)-1)
ic_analytique = (mu_obs - t_crit * se_obs, mu_obs + t_crit * se_obs)
ic_boot_perc = np.percentile(boots_moy, [2.5, 97.5])

ax2.hist(boots_moy, bins=50, density=True, alpha=0.55, color="steelblue")
ax2.axvline(mu_obs, color="tomato", lw=2.5, label=f"Moyenne = {mu_obs:.1f}")
ax2.axvspan(*ic_analytique, alpha=0.15, color="tomato",
            label=f"IC t de Student [{ic_analytique[0]:.1f}, {ic_analytique[1]:.1f}]")
ax2.axvspan(*ic_boot_perc, alpha=0.1, color="steelblue",
            label=f"IC bootstrap [{ic_boot_perc[0]:.1f}, {ic_boot_perc[1]:.1f}]")
# Gaussienne de l'IC analytique
x_norm = np.linspace(mu_obs - 4*se_obs, mu_obs + 4*se_obs, 300)
ax2.plot(x_norm, stats.norm.pdf(x_norm, mu_obs, se_obs), "r--", lw=2,
         label="IC analytique (normale)")
ax2.set_xlabel("Moyenne bootstrap (ms)")
ax2.set_ylabel("Densité")
ax2.set_title("IC bootstrap vs IC analytique pour la moyenne")
ax2.legend(fontsize=8)

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

## IC pour les proportions

Estimer une proportion $p$ à partir de $k$ succès sur $n$ essais. Trois méthodes existent, avec des propriétés très différentes.

### Approximation normale (Wald)

$\hat{p} \pm z_{\alpha/2} \sqrt{\hat{p}(1-\hat{p})/n}$

Simple mais sous-couvre pour les $p$ proches de 0 ou 1, ou pour les petits $n$.

### Intervalle de Wilson

$\frac{\hat{p} + z^2/(2n) \pm z\sqrt{\hat{p}(1-\hat{p})/n + z^2/(4n^2)}}{1 + z^2/n}$

Bonne couverture même pour $p$ extrême ou petit $n$. C'est la méthode recommandée par défaut.

### Intervalle de Clopper-Pearson (exact)

Basé sur la distribution bêta : $[\text{Beta}(\alpha/2; k, n-k+1), \text{Beta}(1-\alpha/2; k+1, n-k)]$. Garantit une couverture au moins égale à $1-\alpha$ (conservatif).

```{code-cell} python
def ic_proportion(k, n, alpha=0.05):
    """Trois méthodes d'IC pour une proportion."""
    phat = k / n
    z = stats.norm.ppf(1 - alpha / 2)

    # Wald (normale)
    se_wald = np.sqrt(phat * (1 - phat) / n)
    wald = (phat - z * se_wald, phat + z * se_wald)

    # Wilson
    denom = 1 + z**2 / n
    centre = (phat + z**2 / (2*n)) / denom
    rayon = z * np.sqrt(phat*(1-phat)/n + z**2/(4*n**2)) / denom
    wilson = (centre - rayon, centre + rayon)

    # Clopper-Pearson (exact)
    if k == 0:
        cp_low = 0.0
    else:
        cp_low = stats.beta.ppf(alpha/2, k, n - k + 1)
    if k == n:
        cp_high = 1.0
    else:
        cp_high = stats.beta.ppf(1 - alpha/2, k + 1, n - k)

    return {
        "p_hat": phat,
        "Wald": wald,
        "Wilson": wilson,
        "Clopper-Pearson": (cp_low, cp_high),
    }

# Exemples
print(f"{'Cas':35s} {'Wald':>25} {'Wilson':>25} {'Clopper-Pearson':>25}")
print("-" * 115)
for desc, k, n in [
    ("50/100 (p=0.50, grand n)", 50, 100),
    ("5/10 (p=0.50, petit n)", 5, 10),
    ("95/100 (p=0.95, proche bord)", 95, 100),
    ("1/10 (p=0.10, rare)", 1, 10),
    ("0/20 (zéro succès)", 0, 20),
]:
    res = ic_proportion(k, n)
    wald_s = f"[{res['Wald'][0]:.3f}, {res['Wald'][1]:.3f}]"
    wils_s = f"[{res['Wilson'][0]:.3f}, {res['Wilson'][1]:.3f}]"
    cp_s = f"[{res['Clopper-Pearson'][0]:.3f}, {res['Clopper-Pearson'][1]:.3f}]"
    print(f"{desc:35s} {wald_s:>25} {wils_s:>25} {cp_s:>25}")
```

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

# Couverture simulée : quelle méthode atteint vraiment 95% ?
fig, axes = plt.subplots(1, 2, figsize=(13, 5))

valeurs_p = np.linspace(0.02, 0.98, 30)
n_test = 20
n_sim = 1000
methodes = ["Wald", "Wilson", "Clopper-Pearson"]
couleurs = {"Wald": "tomato", "Wilson": "seagreen", "Clopper-Pearson": "steelblue"}

couvertures = {m: [] for m in methodes}

for p_vrai in valeurs_p:
    cov_par_methode = {m: 0 for m in methodes}
    for _ in range(n_sim):
        k = rng.binomial(n_test, p_vrai)
        res = ic_proportion(k, n_test)
        for m in methodes:
            lo, hi = res[m]
            if lo <= p_vrai <= hi:
                cov_par_methode[m] += 1
    for m in methodes:
        couvertures[m].append(cov_par_methode[m] / n_sim)

ax = axes[0]
for m in methodes:
    ax.plot(valeurs_p, couvertures[m], lw=2, color=couleurs[m], label=m)
ax.axhline(0.95, color="black", lw=1.5, linestyle="--", label="Niveau nominal 95%")
ax.set_xlabel("Vraie proportion p")
ax.set_ylabel("Taux de couverture simulé")
ax.set_title(f"Couverture des IC à 95%\n(n={n_test}, {n_sim} simulations)")
ax.legend(fontsize=8)
ax.set_ylim(0.75, 1.02)

# Largeur des IC
ax2 = axes[1]
for m in methodes:
    largeurs = []
    for p_vrai in valeurs_p:
        k = round(p_vrai * n_test)
        res = ic_proportion(k, n_test)
        lo, hi = res[m]
        largeurs.append(hi - lo)
    ax2.plot(valeurs_p, largeurs, lw=2, color=couleurs[m], label=m)
ax2.set_xlabel("Proportion p")
ax2.set_ylabel("Largeur de l'IC")
ax2.set_title(f"Largeur des IC selon p\n(n={n_test})")
ax2.legend(fontsize=8)

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

## Taille d'échantillon et précision de l'IC

La largeur d'un IC pour la moyenne est $2 \times t_{\alpha/2} \times \sigma/\sqrt{n}$. Pour réduire la marge d'erreur de moitié, il faut **quadrupler** $n$.

```{code-cell} python
def n_necessaire_moyenne(sigma, marge, alpha=0.05, approximation=True):
    """n pour obtenir un IC de ±marge avec écart-type σ connu."""
    z = stats.norm.ppf(1 - alpha/2)
    return int(np.ceil((z * sigma / marge)**2))

def n_necessaire_proportion(marge, p_estime=0.5, alpha=0.05):
    """n pour IC de Wilson de ±marge pour une proportion."""
    z = stats.norm.ppf(1 - alpha/2)
    return int(np.ceil(z**2 * p_estime * (1 - p_estime) / marge**2))

print("Taille d'échantillon nécessaire (IC 95%) :\n")
print("=== Pour une moyenne (σ = 10) ===")
for marge in [5, 2, 1, 0.5]:
    n = n_necessaire_moyenne(sigma=10, marge=marge)
    print(f"  Marge d'erreur ±{marge:4.1f} : n = {n:6d}")

print("\n=== Pour une proportion (p ≈ 0.5, pire cas) ===")
for marge in [0.10, 0.05, 0.03, 0.01]:
    n = n_necessaire_proportion(marge=marge, p_estime=0.5)
    print(f"  Marge d'erreur ±{marge:.2f} : n = {n:6d}")
```

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

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

# IC 95% pour la moyenne en fonction de n
sigma_fixe = 15
tailles_n = np.arange(10, 501, 5)
largeurs_ic = [2 * stats.t.ppf(0.975, df=n-1) * sigma_fixe / np.sqrt(n)
               for n in tailles_n]

ax = axes[0]
ax.plot(tailles_n, largeurs_ic, color="steelblue", lw=2)
ax.fill_between(tailles_n, 0, largeurs_ic, alpha=0.15, color="steelblue")
for n_mark, color in [(30, "tomato"), (100, "seagreen"), (200, "darkorange")]:
    larg = 2 * stats.t.ppf(0.975, df=n_mark-1) * sigma_fixe / np.sqrt(n_mark)
    ax.scatter(n_mark, larg, s=80, color=color, zorder=5)
    ax.annotate(f"n={n_mark}\n±{larg:.1f}", xy=(n_mark, larg),
                xytext=(n_mark+20, larg+1), fontsize=8, color=color,
                arrowprops=dict(arrowstyle="->", color=color))
ax.set_xlabel("Taille d'échantillon n")
ax.set_ylabel("Largeur totale de l'IC (unités de σ)")
ax.set_title(f"Largeur de l'IC 95% pour la moyenne (σ={sigma_fixe})")
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)

# Simulation de couverture des IC à 95% pour la moyenne
n_list = [10, 20, 50, 100]
n_exp = 1000
mu_vrai = 50
sigma_vrai = 10

ax2 = axes[1]
for j, n_t in enumerate(n_list):
    couvertures_n = []
    for _ in range(n_exp):
        x = rng.normal(mu_vrai, sigma_vrai, n_t)
        mu_hat = x.mean()
        se = x.std(ddof=1) / np.sqrt(n_t)
        t_crit = stats.t.ppf(0.975, df=n_t-1)
        ic_lo = mu_hat - t_crit * se
        ic_hi = mu_hat + t_crit * se
        couvertures_n.append(ic_lo <= mu_vrai <= ic_hi)
    couv = np.mean(couvertures_n)
    ax2.bar(j, couv, color="steelblue" if abs(couv - 0.95) < 0.02 else "tomato", width=0.6)
    ax2.text(j, couv + 0.002, f"{couv:.3f}", ha="center", fontsize=9)

ax2.axhline(0.95, color="tomato", lw=2, linestyle="--", label="Niveau nominal 95%")
ax2.set_xticks(range(len(n_list)))
ax2.set_xticklabels([f"n={n}" for n in n_list])
ax2.set_ylabel("Taux de couverture")
ax2.set_ylim(0.90, 1.00)
ax2.set_title(f"Couverture de l'IC t de Student\n({n_exp} expériences simulées)")
ax2.legend()
ax2.spines["top"].set_visible(False)
ax2.spines["right"].set_visible(False)

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

## IC pour la différence de moyennes : Welch vs Student

Comparer deux groupes est l'une des opérations les plus fréquentes. Le **test de Student** suppose des variances égales ; le **test de Welch** (recommandé par défaut) ne le suppose pas.

```{code-cell} python
# Deux groupes avec variances différentes
groupe_A = rng.normal(52, 8, 40)
groupe_B = rng.normal(58, 15, 35)

# Student (variance poolée)
diff_obs = groupe_A.mean() - groupe_B.mean()
n_A, n_B = len(groupe_A), len(groupe_B)
s_A, s_B = groupe_A.std(ddof=1), groupe_B.std(ddof=1)

# Variance poolée (Student)
sp2 = ((n_A - 1) * s_A**2 + (n_B - 1) * s_B**2) / (n_A + n_B - 2)
se_student = np.sqrt(sp2 * (1/n_A + 1/n_B))
t_crit_s = stats.t.ppf(0.975, df=n_A + n_B - 2)
ic_student = (diff_obs - t_crit_s * se_student, diff_obs + t_crit_s * se_student)

# Welch (variances distinctes)
se_welch = np.sqrt(s_A**2/n_A + s_B**2/n_B)
df_welch = (s_A**2/n_A + s_B**2/n_B)**2 / (
    (s_A**2/n_A)**2 / (n_A-1) + (s_B**2/n_B)**2 / (n_B-1)
)
t_crit_w = stats.t.ppf(0.975, df=df_welch)
ic_welch = (diff_obs - t_crit_w * se_welch, diff_obs + t_crit_w * se_welch)

# scipy (Welch est la valeur par défaut)
t_stat, p_val = stats.ttest_ind(groupe_A, groupe_B, equal_var=False)

print(f"Groupe A : n={n_A}, μ={groupe_A.mean():.1f}, σ={s_A:.1f}")
print(f"Groupe B : n={n_B}, μ={groupe_B.mean():.1f}, σ={s_B:.1f}")
print(f"Ratio des variances : {(s_B/s_A)**2:.2f}  → variances clairement différentes")
print()
print(f"IC Student 95% (variances égales) : {ic_student[0]:.2f} à {ic_student[1]:.2f}")
print(f"IC Welch 95%   (variances libres) : {ic_welch[0]:.2f} à {ic_welch[1]:.2f}")
print(f"Degrés de liberté Welch : {df_welch:.1f}  (vs {n_A+n_B-2} pour Student)")
```

```{admonition} Welch ou Student ?
:class: tip
Par défaut, utilisez **Welch** (`scipy.stats.ttest_ind(equal_var=False)`). Il est robuste lorsque les variances sont égales (très peu de perte de puissance) et correct lorsqu'elles sont différentes. Le test de Student avec variances poolées peut être sévèrement anticonservatif si les variances et les tailles d'échantillon sont toutes deux différentes.
```

## Le Jackknife

Le **jackknife** (Quenouille, 1949) est l'ancêtre du bootstrap. Il consiste à recalculer l'estimateur en supprimant tour à tour chaque observation.

$$\hat{\theta}_{-i} = \text{statistique}(x_1, \ldots, x_{i-1}, x_{i+1}, \ldots, x_n)$$

**Estimation du biais :**

$$\hat{B}_{\text{jk}} = (n-1)(\bar{\hat{\theta}}_{-\cdot} - \hat{\theta})$$

**Estimation de la variance :**

$$\hat{V}_{\text{jk}} = \frac{n-1}{n} \sum_{i=1}^n (\hat{\theta}_{-i} - \bar{\hat{\theta}}_{-\cdot})^2$$

```{code-cell} python
def jackknife(data, statistic):
    """
    Estimation par jackknife du biais et de la variance.

    Returns
    -------
    dict avec theta_hat, biais_jk, var_jk, se_jk, vals_jk
    """
    data = np.asarray(data)
    n = len(data)
    theta_hat = statistic(data)
    vals_jk = np.array([statistic(np.delete(data, i)) for i in range(n)])
    theta_jk_mean = vals_jk.mean()
    biais = (n - 1) * (theta_jk_mean - theta_hat)
    var_jk = (n-1) / n * np.sum((vals_jk - theta_jk_mean)**2)
    se_jk = np.sqrt(var_jk)
    return {
        "theta_hat": theta_hat,
        "biais_jk": biais,
        "var_jk": var_jk,
        "se_jk": se_jk,
        "theta_corrige": theta_hat - biais,
        "vals_jk": vals_jk,
    }

# Application sur les temps de chargement
print("=== Jackknife sur les temps de chargement ===\n")
for nom, stat in [("Moyenne", np.mean), ("Médiane", np.median),
                   ("Écart-type", np.std), ("Skewness", stats.skew)]:
    res = jackknife(data_web, stat)
    print(f"{nom:12s} : estimé = {res['theta_hat']:8.3f},  "
          f"biais = {res['biais_jk']:8.4f},  "
          f"SE jackknife = {res['se_jk']:8.4f},  "
          f"estimé corrigé = {res['theta_corrige']:8.3f}")
```

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

# Comparaison bootstrap vs jackknife : simulation de couverture (1000 expériences)
fig, axes = plt.subplots(1, 2, figsize=(13, 5))

n_experience = 30
mu_vrai = 50
sigma_vrai = 12
n_obs = 30

couv_t = 0
couv_boot_perc = 0
couv_boot_bca = 0
couv_jk = 0

for _ in range(n_experience):
    x = rng.normal(mu_vrai, sigma_vrai, n_obs)
    mu_hat = x.mean()
    se = x.std(ddof=1) / np.sqrt(n_obs)
    t_c = stats.t.ppf(0.975, df=n_obs-1)

    # IC t
    if mu_hat - t_c*se <= mu_vrai <= mu_hat + t_c*se:
        couv_t += 1

    # IC bootstrap percentile
    b_p_l, b_p_h, _ = bootstrap_ic(x, np.mean, n_boot=200, methode="percentile")
    if b_p_l <= mu_vrai <= b_p_h:
        couv_boot_perc += 1

    # IC bootstrap BCa
    b_b_l, b_b_h, _ = bootstrap_ic(x, np.mean, n_boot=200, methode="bca")
    if b_b_l <= mu_vrai <= b_b_h:
        couv_boot_bca += 1

    # IC jackknife
    res_jk = jackknife(x, np.mean)
    z_jk = stats.norm.ppf(0.975)
    if mu_hat - z_jk*res_jk["se_jk"] <= mu_vrai <= mu_hat + z_jk*res_jk["se_jk"]:
        couv_jk += 1

couvertures_finales = {
    "t de Student": couv_t / n_experience,
    "Boot. percentile": couv_boot_perc / n_experience,
    "Boot. BCa": couv_boot_bca / n_experience,
    "Jackknife": couv_jk / n_experience,
}

# Graphique de couverture
ax = axes[0]
noms_m = list(couvertures_finales.keys())
covs = list(couvertures_finales.values())
colors = ["seagreen" if abs(c - 0.95) < 0.015 else "tomato" for c in covs]
bars = ax.bar(noms_m, covs, color=colors, width=0.55)
ax.axhline(0.95, color="black", lw=2, linestyle="--", label="Niveau nominal 95%")
ax.set_ylim(0.85, 1.00)
for bar, c in zip(bars, covs):
    ax.text(bar.get_x() + bar.get_width()/2, c + 0.003, f"{c:.3f}",
            ha="center", va="bottom", fontsize=9.5)
ax.set_ylabel("Taux de couverture")
ax.set_title(f"Couverture des IC pour la moyenne\n(n={n_obs}, {n_experience} expériences, N(50,12))")
ax.legend()
ax.tick_params(axis="x", rotation=20)
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)

# Illustration jackknife vs bootstrap pour une statistique difficile
ax2 = axes[1]
# Corrélation : statistique pour laquelle jackknife et bootstrap diffèrent
x_corr = rng.normal(0, 1, 40)
y_corr = 0.5 * x_corr + rng.normal(0, 0.8, 40)
data_corr = np.column_stack([x_corr, y_corr])

def pearson_r(data_2d):
    return stats.pearsonr(data_2d[:, 0], data_2d[:, 1])[0]

r_obs = pearson_r(data_corr)

# Bootstrap
boots_r = []
for _ in range(500):
    idx = rng.integers(0, len(data_corr), size=len(data_corr))
    boots_r.append(pearson_r(data_corr[idx]))
boots_r = np.array(boots_r)
ic_r_boot = np.percentile(boots_r, [2.5, 97.5])

# IC de Fisher (transformation de Fisher)
z_r = np.arctanh(r_obs)
se_z = 1 / np.sqrt(len(data_corr) - 3)
ic_r_fish = (
    np.tanh(z_r - 1.96 * se_z),
    np.tanh(z_r + 1.96 * se_z)
)

ax2.hist(boots_r, bins=40, density=True, alpha=0.55, color="steelblue",
         label="Bootstrap (corrélation)")
ax2.axvline(r_obs, color="tomato", lw=2.5, label=f"r observé = {r_obs:.3f}")
ax2.axvspan(*ic_r_boot, alpha=0.2, color="steelblue",
            label=f"IC boot. [{ic_r_boot[0]:.3f}, {ic_r_boot[1]:.3f}]")
ax2.axvspan(*ic_r_fish, alpha=0.15, color="tomato",
            label=f"IC Fisher [{ic_r_fish[0]:.3f}, {ic_r_fish[1]:.3f}]")
x_kde = np.linspace(boots_r.min(), boots_r.max(), 300)
ax2.plot(x_kde, stats.gaussian_kde(boots_r)(x_kde), "k-", lw=1.5)
ax2.set_xlabel("r de Pearson")
ax2.set_ylabel("Densité")
ax2.set_title("Distribution bootstrap de la corrélation\n(n=40, 2000 rééchantillons)")
ax2.legend(fontsize=8)

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

## Récapitulatif des méthodes

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

# Tableau récapitulatif
recap = pd.DataFrame({
    "Méthode": [
        "IC t de Student (moyenne)",
        "IC Welch (diff. de moyennes)",
        "IC Wilson (proportion)",
        "IC Clopper-Pearson (proportion)",
        "Bootstrap percentile",
        "Bootstrap BCa",
        "Jackknife",
    ],
    "Hypothèses": [
        "Normalité ou grand n",
        "Indépendance, normalité approx.",
        "Indépendance des observations",
        "Indépendance des observations",
        "IID, n suffisant",
        "IID, n ≥ 20 recommandé",
        "IID, régularité de la statistique",
    ],
    "Quand l'utiliser": [
        "Moyenne, n ≥ 30 ou données normales",
        "Comparaison deux groupes (défaut)",
        "Proportions, usage général",
        "Proportions, garantie de couverture",
        "Toute statistique, données simples",
        "Toute statistique, IC précis",
        "Estimation du biais, petit n",
    ],
    "Limites": [
        "Sensible aux outliers pour petit n",
        "Idem Student",
        "Léger sous-couvrement aux extrêmes",
        "Conservatif (IC trop larges)",
        "Biaisé si distribution asymétrique",
        "Coûteux, jackknife pour a",
        "Inexact pour statistiques non-lisses (médiane)",
    ],
}).set_index("Méthode")

print(recap.to_string())
```

L'intervalle de confiance est une réponse calibrée à une question fondamentale : « quelle est notre incertitude sur ce que les données nous disent ? » Sa construction, qu'elle soit analytique ou computationnelle, repose toujours sur le même principe — simuler (analytiquement ou par rééchantillonnage) ce qui se passerait si on répétait l'expérience. Comprendre ce principe, c'est comprendre la logique entière de l'inférence statistique fréquentiste.