---
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 computationnelles

Avant l'ère des ordinateurs, les statisticiens devaient se contenter de résultats analytiques exacts — formules closes, tables de distributions, approximations asymptotiques. Aujourd'hui, la puissance de calcul disponible permet d'**estimer par simulation** ce qui ne peut pas être calculé analytiquement. Ce chapitre couvre les grandes familles de méthodes computationnelles : Monte Carlo, tests de permutation, bootstrap, jackknife, et validation croisée statistique.

```{code-cell} python3
: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 norm, t as t_dist, chi2
from sklearn.model_selection import KFold, LeaveOneOut, cross_val_score
from sklearn.linear_model import LinearRegression, LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_iris
import warnings
warnings.filterwarnings('ignore')

plt.rcParams.update({
    'figure.dpi': 110,
    'axes.spines.top': False,
    'axes.spines.right': False,
    'font.size': 11,
})
rng = np.random.default_rng(42)
```

## Simulation Monte Carlo

### Principe général

La méthode de Monte Carlo repose sur la **loi des grands nombres** : la moyenne empirique d'une grande quantité de réalisations indépendantes d'une variable aléatoire converge vers son espérance. Pour estimer $\mu = \mathbb{E}[f(X)]$ :

$$\hat{\mu}_N = \frac{1}{N}\sum_{i=1}^N f(X_i), \quad X_i \sim F$$

L'erreur d'estimation (écart-type de $\hat{\mu}_N$) décroît en $1/\sqrt{N}$ — indépendamment de la dimension du problème.

### Estimation de π

L'exemple classique : on tire des points uniformément dans $[-1,1]^2$. La proportion de points dans le disque de rayon 1 converge vers $\pi/4$.

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

def estimer_pi(n_points, graine=None):
    """Estimation de π par Monte Carlo."""
    rng_local = np.random.default_rng(graine)
    x = rng_local.uniform(-1, 1, n_points)
    y = rng_local.uniform(-1, 1, n_points)
    dans_cercle = x**2 + y**2 <= 1
    return 4 * dans_cercle.mean(), dans_cercle

# Convergence
n_max = 100_000
tailles = np.logspace(2, 5, 50).astype(int)
estimations = [estimer_pi(n, graine=i)[0] for i, n in enumerate(tailles)]

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

# Visualisation du tirage
ax = axes[0]
n_visu = 3000
x_v = rng.uniform(-1, 1, n_visu)
y_v = rng.uniform(-1, 1, n_visu)
in_v = x_v**2 + y_v**2 <= 1
ax.scatter(x_v[in_v], y_v[in_v], s=2, color='steelblue', alpha=0.5, label='Dans le cercle')
ax.scatter(x_v[~in_v], y_v[~in_v], s=2, color='lightgray', alpha=0.5, label='Hors du cercle')
theta = np.linspace(0, 2*np.pi, 300)
ax.plot(np.cos(theta), np.sin(theta), 'orangered', linewidth=1.5)
ax.set_aspect('equal')
ax.set_title(f'Estimation de π = {4*in_v.mean():.4f}\n(n = {n_visu:,})')
ax.legend(fontsize=9, markerscale=5)

# Convergence
ax = axes[1]
ax.semilogx(tailles, estimations, 'steelblue', linewidth=1.5, alpha=0.8, label='Estimation MC')
ax.axhline(np.pi, color='orangered', linewidth=2, linestyle='--', label=f'π = {np.pi:.5f}')
# Bande de confiance théorique ±2σ
sigma_n = np.sqrt(np.pi/4 * (1 - np.pi/4) / tailles)  # variance de la proportion
ax.fill_between(tailles, np.pi - 2*sigma_n*4, np.pi + 2*sigma_n*4,
                alpha=0.2, color='steelblue', label='±2σ théorique')
ax.set_xlabel('Nombre de points N (échelle log)')
ax.set_ylabel('Estimation de π')
ax.set_title('Convergence Monte Carlo (erreur ∝ 1/√N)')
ax.legend(fontsize=9)

plt.suptitle('Méthode de Monte Carlo — Estimation de π', fontsize=13, y=1.01)
plt.tight_layout()
plt.show()
```

### Intégration par Monte Carlo

Pour estimer $I = \int_a^b f(x)\,dx$, on peut écrire $I = (b-a)\,\mathbb{E}[f(U)]$ où $U \sim \mathcal{U}(a,b)$ et estimer par la moyenne empirique.

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

# Intégration Monte Carlo vs intégration numérique classique
# Exemple : ∫₀¹ exp(-x²) dx (pas de forme close simple)
from scipy.integrate import quad

def f(x):
    return np.exp(-x**2)

# Valeur exacte (numérique)
I_exact, _ = quad(f, 0, 1)

# Monte Carlo
n_vals = np.logspace(2, 6, 40).astype(int)
mc_estimates = []
mc_erreurs = []

for n in n_vals:
    u = rng.uniform(0, 1, n)
    vals = f(u)
    mc_estimates.append(vals.mean())
    mc_erreurs.append(vals.std() / np.sqrt(n))

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

ax = axes[0]
x_plot = np.linspace(0, 1, 300)
ax.plot(x_plot, f(x_plot), 'steelblue', linewidth=2)
ax.fill_between(x_plot, f(x_plot), alpha=0.3, color='steelblue',
                label=f'I = {I_exact:.6f}')
ax.set_title(r'$\int_0^1 e^{-x^2}\, dx$')
ax.set_xlabel('x'); ax.set_ylabel('f(x)')
ax.legend()

ax = axes[1]
erreur_absolue = np.abs(np.array(mc_estimates) - I_exact)
ax.loglog(n_vals, erreur_absolue, 'steelblue', linewidth=2, label='Erreur MC')
ax.loglog(n_vals, np.array(mc_erreurs), 'orangered', linestyle='--',
          linewidth=1.5, label='±1σ théorique')
# Référence 1/√n
ref_n = n_vals[5]
ref_e = erreur_absolue[5]
ax.loglog(n_vals, ref_e * np.sqrt(ref_n / n_vals), 'gray', linestyle=':',
          linewidth=1, label='∝ 1/√N')
ax.set_xlabel('N (log)')
ax.set_ylabel('Erreur absolue (log)')
ax.set_title('Convergence de l\'intégration Monte Carlo')
ax.legend(fontsize=9)

plt.suptitle('Intégration par Monte Carlo', fontsize=13, y=1.01)
plt.tight_layout()
plt.show()
```

---

## Tests de permutation

### Principe

Un test de permutation est une alternative **non paramétrique** aux tests classiques. L'idée : sous H₀ (les deux groupes sont échangeables), toutes les permutations des labels sont également probables. La p-valeur est la proportion de permutations qui donnent une statistique aussi ou plus extrême que celle observée.

**Avantages :**

- Ne suppose aucune distribution particulière
- Valide à taille d'échantillon finie (exactitude)
- Fonctionne pour n'importe quelle statistique (même sans distribution analytique connue)

**Inconvénient :** coûteux pour $n$ grand (mais on peut approcher par un sous-ensemble aléatoire de permutations).

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

def test_permutation(groupe_a, groupe_b, statistic='mean_diff', n_permutations=10000, graine=42):
    """Test de permutation pour deux groupes indépendants."""
    rng_p = np.random.default_rng(graine)
    combined = np.concatenate([groupe_a, groupe_b])
    n_a = len(groupe_a)

    if statistic == 'mean_diff':
        stat_obs = groupe_b.mean() - groupe_a.mean()
        stats_perm = np.array([
            rng_p.permutation(combined)[:n_a].mean() - rng_p.permutation(combined)[n_a:].mean()
            for _ in range(n_permutations)
        ])
    elif statistic == 'median_diff':
        stat_obs = np.median(groupe_b) - np.median(groupe_a)
        stats_perm = np.array([
            np.median(rng_p.permutation(combined)[:n_a]) - np.median(rng_p.permutation(combined)[n_a:])
            for _ in range(n_permutations)
        ])

    p_value = np.mean(np.abs(stats_perm) >= np.abs(stat_obs))
    return stat_obs, stats_perm, p_value

# Données : temps de réaction (ms) de deux groupes
groupe_A = rng.normal(350, 50, 30)   # contrôle
groupe_B = rng.normal(320, 55, 28)   # traitement (réduction ~30ms)

stat_obs, stats_perm, p_perm = test_permutation(groupe_A, groupe_B)

# Test t pour comparaison
t_stat, p_t = stats.ttest_ind(groupe_A, groupe_B)
# Test de Mann-Whitney
u_stat, p_mw = stats.mannwhitneyu(groupe_A, groupe_B, alternative='two-sided')

print("Comparaison des méthodes de test")
print(f"  Différence observée (B - A) : {stat_obs:.2f} ms")
print()
print(f"{'Méthode':>25} | {'Statistique':>12} | {'p-valeur':>10}")
print("-" * 55)
print(f"{'Test de permutation':>25} | {'|Δ| observé':>12} | {p_perm:>10.4f}")
print(f"{'Test t de Student':>25} | {t_stat:>12.3f} | {p_t:>10.4f}")
print(f"{'Mann-Whitney':>25} | {u_stat:>12.0f} | {p_mw:>10.4f}")

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

# Distribution de permutation
ax = axes[0]
ax.hist(stats_perm, bins=60, color='steelblue', alpha=0.7, density=True,
        label='Distribution de permutation')
ax.axvline(stat_obs, color='orangered', linewidth=2.5, label=f'Observé = {stat_obs:.2f}')
ax.axvline(-stat_obs, color='orangered', linewidth=2.5, linestyle='--')

# Superposition distribution normale théorique (test t)
x_range = np.linspace(stats_perm.min(), stats_perm.max(), 300)
se_th = np.sqrt(groupe_A.var()/len(groupe_A) + groupe_B.var()/len(groupe_B))
ax.plot(x_range, norm.pdf(x_range, 0, se_th), 'seagreen', linewidth=2,
        label='Distribution normale (test t)')
ax.fill_betweenx([0, ax.get_ylim()[1] if ax.get_ylim()[1] > 0 else 0.05],
                  -abs(stat_obs)*5, -abs(stat_obs), alpha=0.15, color='orangered')
ax.fill_betweenx([0, ax.get_ylim()[1] if ax.get_ylim()[1] > 0 else 0.05],
                  abs(stat_obs), abs(stat_obs)*5, alpha=0.15, color='orangered')
ax.set_xlabel('Différence de moyennes')
ax.set_ylabel('Densité')
ax.set_title(f'Test de permutation\np-valeur = {p_perm:.4f} ({(np.abs(stats_perm) >= np.abs(stat_obs)).sum()} / {len(stats_perm)})')
ax.legend(fontsize=9)

# Données
ax = axes[1]
ax.boxplot([groupe_A, groupe_B], labels=['Groupe A (contrôle)', 'Groupe B (traitement)'],
           patch_artist=True,
           boxprops=dict(facecolor='steelblue', alpha=0.5))
ax.scatter([1]*len(groupe_A), groupe_A, alpha=0.3, s=20, color='steelblue')
ax.scatter([2]*len(groupe_B), groupe_B, alpha=0.3, s=20, color='orangered')
ax.set_ylabel('Temps de réaction (ms)')
ax.set_title('Données des deux groupes')

plt.suptitle('Test de permutation vs test paramétrique', fontsize=13, y=1.01)
plt.tight_layout()
plt.show()
```

---

## Bootstrap

### Bootstrap non paramétrique

Le **bootstrap** de Efron (1979) est une méthode de rééchantillonnage pour estimer la distribution d'une statistique quelconque. L'idée : l'échantillon observé **est** la meilleure approximation de la population. On peut donc simuler ce qui se passerait si on rééchantillonnait la population en tirant des échantillons **avec remise** de l'échantillon observé.

**Algorithme :**

1. Calculer la statistique $\hat{\theta}$ sur l'échantillon original
2. Pour $b = 1, \ldots, B$ : tirer un échantillon bootstrap $X^{*b}$ de taille $n$ avec remise depuis $X$
3. Calculer $\hat{\theta}^{*b}$ sur chaque échantillon bootstrap
4. Utiliser la distribution de $\{\hat{\theta}^{*b}\}$ pour estimer biais, variance, et intervalles de confiance

**Intervalle de confiance percentile :**
$$\text{IC}_{1-\alpha} = \left[\hat{\theta}^{*(\alpha/2 \cdot B)}, \hat{\theta}^{*(1-\alpha/2 \cdot B)}\right]$$

**Intervalle BCa** (*Bias-Corrected and accelerated*) : corrige le biais et l'accélération de la distribution bootstrap — méthode recommandée.

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

def bootstrap(data, statistic, n_boot=500, ci_level=0.95, graine=42):
    """Bootstrap non paramétrique."""
    rng_b = np.random.default_rng(graine)
    n = len(data)
    stat_obs = statistic(data)
    boot_stats = np.array([
        statistic(rng_b.choice(data, n, replace=True))
        for _ in range(n_boot)
    ])
    alpha = 1 - ci_level
    ic_low  = np.percentile(boot_stats, alpha/2 * 100)
    ic_high = np.percentile(boot_stats, (1 - alpha/2) * 100)
    biais   = boot_stats.mean() - stat_obs
    return stat_obs, boot_stats, ic_low, ic_high, biais

# Données : revenus mensuels (distribution log-normale)
revenus = rng.lognormal(mean=8.5, sigma=0.5, size=80)

# Statistiques d'intérêt
def mediane(x): return np.median(x)
def q25(x): return np.percentile(x, 25)
def kurtosis(x): return stats.kurtosis(x)
def cv(x): return x.std() / x.mean()  # coefficient de variation

statistiques = {
    'Médiane'                : mediane,
    '1er quartile'           : q25,
    'Coeff. de variation'    : cv,
    'Kurtosis'               : kurtosis,
}

fig, axes = plt.subplots(2, 2, figsize=(13, 9))
for ax, (nom_stat, func) in zip(axes.flatten(), statistiques.items()):
    stat_obs, boot_stats, ic_low, ic_high, biais = bootstrap(revenus, func)

    ax.hist(boot_stats, bins=50, color='steelblue', alpha=0.7, density=True,
            label=f'Distribution bootstrap\n(B=5000)')
    ax.axvline(stat_obs, color='orangered', linewidth=2.5, label=f'Observé = {stat_obs:.4f}')
    ax.axvline(ic_low,  color='seagreen', linewidth=2, linestyle='--',
               label=f'IC 95%: [{ic_low:.4f}, {ic_high:.4f}]')
    ax.axvline(ic_high, color='seagreen', linewidth=2, linestyle='--')
    ax.set_title(f'{nom_stat}\n(biais bootstrap = {biais:.4f})')
    ax.legend(fontsize=8)
    ax.set_xlabel(f'Valeur de la statistique')

plt.suptitle('Bootstrap non paramétrique — Distribution de statistiques sur revenus (n=80)',
             fontsize=13, y=1.01)
plt.tight_layout()
plt.show()
```

### Bootstrap paramétrique

Le **bootstrap paramétrique** suppose un modèle paramétrique ajusté aux données. Les échantillons bootstrap sont générés depuis ce modèle ajusté, et non par rééchantillonnage direct.

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

# Bootstrap paramétrique : inférence sur un paramètre de loi log-normale
# On suppose X ~ LogNormal(μ, σ²)
revenus_log = np.log(revenus)
mu_hat = revenus_log.mean()
sigma_hat = revenus_log.std()

n = len(revenus)
B = 500
mu_boot_param    = np.zeros(B)
sigma_boot_param = np.zeros(B)

for b in range(B):
    sample_b = rng.normal(mu_hat, sigma_hat, n)
    mu_boot_param[b]    = sample_b.mean()
    sigma_boot_param[b] = sample_b.std()

# IC bootstrap paramétrique pour μ (en espace log)
ic_mu_low   = np.percentile(mu_boot_param, 2.5)
ic_mu_high  = np.percentile(mu_boot_param, 97.5)
# IC analytique (référence)
ic_mu_ana_low  = mu_hat - 1.96 * sigma_hat / np.sqrt(n)
ic_mu_ana_high = mu_hat + 1.96 * sigma_hat / np.sqrt(n)

print("Bootstrap paramétrique vs IC analytique (loi log-normale)")
print(f"  μ_chapeau = {mu_hat:.4f},  σ_chapeau = {sigma_hat:.4f}")
print(f"  IC 95% bootstrap paramétrique : [{ic_mu_low:.4f}, {ic_mu_high:.4f}]")
print(f"  IC 95% analytique (normale)   : [{ic_mu_ana_low:.4f}, {ic_mu_ana_high:.4f}]")

fig, ax = plt.subplots(figsize=(9, 4))
ax.hist(mu_boot_param, bins=50, color='steelblue', alpha=0.7, density=True,
        label='Bootstrap paramétrique (μ log-revenu)')
ax.axvline(mu_hat, color='orangered', linewidth=2.5, label=f'Estimé = {mu_hat:.4f}')
ax.axvline(ic_mu_low, color='seagreen', linewidth=2, linestyle='--',
           label=f'IC 95% bootstrap')
ax.axvline(ic_mu_high, color='seagreen', linewidth=2, linestyle='--')

x_range = np.linspace(ic_mu_low - 0.1, ic_mu_high + 0.1, 300)
ax.plot(x_range, norm.pdf(x_range, mu_hat, sigma_hat/np.sqrt(n)),
        'gray', linestyle=':', linewidth=2, label='IC analytique')
ax.set_xlabel('μ (log-revenu)')
ax.set_title('Bootstrap paramétrique pour le paramètre de localisation d\'une loi log-normale')
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()
```

### Erreur bootstrap vs taille d'échantillon

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

# Largeur de l'IC bootstrap en fonction de n
n_vals = [20, 30, 50, 80, 120, 200, 500]
largeurs_median = []
largeurs_cv     = []

for n_val in n_vals:
    largeurs_m, largeurs_c = [], []
    for graine in range(20):
        echantillon = rng.lognormal(8.5, 0.5, n_val)
        _, _, il_m, ih_m, _ = bootstrap(echantillon, mediane, n_boot=500, graine=graine)
        _, _, il_c, ih_c, _ = bootstrap(echantillon, cv,      n_boot=500, graine=graine)
        largeurs_m.append(ih_m - il_m)
        largeurs_c.append(ih_c - il_c)
    largeurs_median.append(np.median(largeurs_m))
    largeurs_cv.append(np.median(largeurs_c))

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
for ax, largeurs, titre, ref in zip(
    axes,
    [largeurs_median, largeurs_cv],
    ['Largeur IC médiane', 'Largeur IC coeff. variation'],
    [largeurs_median[0], largeurs_cv[0]]
):
    ax.loglog(n_vals, largeurs, 'o-', color='steelblue', linewidth=2, label='Bootstrap (médiane 50 rep.)')
    # Référence 1/√n
    ax.loglog(n_vals, ref * np.sqrt(n_vals[0] / np.array(n_vals)), '--',
              color='orangered', linewidth=1.5, label='∝ 1/√n')
    ax.set_xlabel('Taille d\'échantillon n')
    ax.set_ylabel('Largeur IC 95%')
    ax.set_title(titre)
    ax.legend(fontsize=9)

plt.suptitle('Précision du bootstrap en fonction de la taille d\'échantillon', fontsize=13, y=1.01)
plt.tight_layout()
plt.show()
```

---

## Jackknife

Le **jackknife** de Quenouille (1949) est une autre méthode de rééchantillonnage, plus ancienne que le bootstrap. Chaque échantillon jackknife omet un individu : $X_{(-i)} = (X_1, \ldots, X_{i-1}, X_{i+1}, \ldots, X_n)$.

**Estimation du biais :**
$$\widehat{\text{biais}} = (n-1)(\bar{\theta}_{(-\cdot)} - \hat{\theta})$$

**Estimation de la variance :**
$$\widehat{\text{Var}}(\hat{\theta}) = \frac{n-1}{n}\sum_{i=1}^n (\hat{\theta}_{(-i)} - \bar{\theta}_{(-\cdot)})^2$$

Le jackknife est **plus conservateur** que le bootstrap (il tend à sous-estimer la variance pour des statistiques non lisses comme la médiane). Il reste utile pour l'estimation du biais.

```{code-cell} python3
def jackknife(data, statistic):
    """Estimation jackknife du biais et de la variance."""
    n = len(data)
    stat_obs = statistic(data)
    jack_stats = np.array([
        statistic(np.delete(data, i)) for i in range(n)
    ])
    mean_jack  = jack_stats.mean()
    biais_jack = (n-1) * (mean_jack - stat_obs)
    var_jack   = (n-1)/n * np.sum((jack_stats - mean_jack)**2)
    return stat_obs, jack_stats, biais_jack, np.sqrt(var_jack)

# Comparaison bootstrap vs jackknife pour plusieurs statistiques
n_data = 50
data_test = rng.exponential(scale=2.0, size=n_data)

print(f"Comparaison Bootstrap vs Jackknife (n={n_data}, loi exponentielle)\n")
print(f"{'Statistique':>20} | {'Observé':>10} | {'Biais Jack.':>12} | {'SE Jack.':>10} | {'SE Boot.':>10}")
print("-" * 72)

for nom, func in [('Moyenne', np.mean), ('Médiane', np.median),
                  ('Écart-type', np.std), ('Skewness', stats.skew)]:
    stat_obs, _, biais_j, se_j = jackknife(data_test, func)
    _, boot_stats, _, _, _ = bootstrap(data_test, func, n_boot=500)
    se_b = boot_stats.std()
    print(f"{nom:>20} | {stat_obs:>10.4f} | {biais_j:>+12.5f} | {se_j:>10.5f} | {se_b:>10.5f}")
```

---

## Validation croisée statistique

La validation croisée évalue la **capacité de généralisation** d'un modèle statistique en estimant son erreur sur des données non utilisées pour l'ajustement.

### Stratégies

**Leave-One-Out (LOO)** : chaque observation est laissée de côté une fois. Peu biaisé mais très variable (et coûteux pour $n$ grand).

**k-fold** : partition aléatoire en $k$ groupes ; on entraîne sur $k-1$ groupes et évalue sur le groupe restant. $k = 5$ ou $k = 10$ offrent le meilleur compromis biais-variance.

**Stratifiée** : la partition respecte la proportion des classes (recommandée pour la classification).

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

from sklearn.datasets import make_regression
from sklearn.metrics import mean_squared_error

# Comparaison LOO, 5-fold, 10-fold sur une régression linéaire
X_reg, y_reg = make_regression(n_samples=100, n_features=5, noise=20, random_state=42)
model_lr = LinearRegression()

def rmse(y_true, y_pred): return np.sqrt(mean_squared_error(y_true, y_pred))

strategies = {
    'LOO'     : LeaveOneOut(),
    '5-fold'  : KFold(n_splits=5,  shuffle=True, random_state=42),
    '10-fold' : KFold(n_splits=10, shuffle=True, random_state=42),
    '20-fold' : KFold(n_splits=20, shuffle=True, random_state=42),
}

resultats_cv = {}
for nom, cv_strat in strategies.items():
    scores = cross_val_score(model_lr, X_reg, y_reg,
                             cv=cv_strat, scoring='neg_root_mean_squared_error')
    resultats_cv[nom] = -scores

print("Validation croisée — Régression linéaire (n=100, p=5)\n")
print(f"{'Stratégie':>12} | {'RMSE moyen':>12} | {'SE':>10} | {'N folds':>8}")
print("-" * 50)
for nom, scores in resultats_cv.items():
    print(f"{nom:>12} | {scores.mean():>12.4f} | {scores.std():>10.4f} | {len(scores):>8}")

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

# Distribution des erreurs par fold
ax = axes[0]
positions = range(len(resultats_cv))
ax.boxplot([v for v in resultats_cv.values()],
           labels=list(resultats_cv.keys()), patch_artist=True,
           boxprops=dict(facecolor='steelblue', alpha=0.5))
ax.set_ylabel('RMSE par fold')
ax.set_title('Distribution des erreurs par fold')

# Biais-variance selon k
ax = axes[1]
ks = [2, 5, 10, 20, 50]
rmse_means = []
rmse_stds  = []
for k in ks:
    cv_k = KFold(n_splits=k, shuffle=True, random_state=42)
    sc = -cross_val_score(model_lr, X_reg, y_reg, cv=cv_k,
                          scoring='neg_root_mean_squared_error')
    rmse_means.append(sc.mean())
    rmse_stds.append(sc.std())

rmse_means = np.array(rmse_means)
rmse_stds  = np.array(rmse_stds)

ax.errorbar(ks, rmse_means, rmse_stds, fmt='o-', color='steelblue', capsize=5,
            linewidth=2, label='RMSE moyen ± SE')
ax.set_xlabel('Nombre de folds k')
ax.set_ylabel('RMSE')
ax.set_title('RMSE moyen et variabilité selon k')
ax.legend()

plt.suptitle('Validation croisée : impact du nombre de folds', fontsize=13, y=1.01)
plt.tight_layout()
plt.show()
```

---

## Simulation de puissance

Calculer la puissance d'un test complexe par simulation est souvent plus simple et plus flexible que les formules analytiques.

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

def puissance_par_simulation(n, effet, distribution='normale', n_simulations=500,
                              alpha=0.05, test='t', graine=42):
    """
    Estime la puissance d'un test par simulation.
    effet : différence des moyennes (Cohen's d si distribution normale et σ=1)
    """
    rng_sim = np.random.default_rng(graine)
    rejets = 0
    for _ in range(n_simulations):
        if distribution == 'normale':
            a = rng_sim.normal(0, 1, n)
            b = rng_sim.normal(effet, 1, n)
        elif distribution == 'exponential':
            a = rng_sim.exponential(1, n)
            b = rng_sim.exponential(1 / (1 + effet), n)  # décalage de la médiane approx
        elif distribution == 't5':  # t à 5 degrés de liberté (queues épaisses)
            a = rng_sim.standard_t(5, n)
            b = rng_sim.standard_t(5, n) + effet

        if test == 't':
            _, p = stats.ttest_ind(a, b)
        elif test == 'wilcoxon':
            _, p = stats.mannwhitneyu(a, b, alternative='two-sided')
        elif test == 'permutation':
            _, _, p = test_permutation(a, b, n_permutations=500, graine=rng_sim.integers(0, 10000))

        if p < alpha:
            rejets += 1
    return rejets / n_simulations

# Courbes de puissance simulées
effets = np.arange(0.1, 1.2, 0.1)
n_test = 30

print("Simulation de puissance (n=30 par groupe, α=5%, 2000 simulations)")
print(f"\n{'Cohen\'s d':>10} | {'t (norm.)':>10} | {'Wilcoxon':>10} | {'t (t₅)':>10}")
print("-" * 48)

puissances = {'t_normale': [], 'wilcoxon_normale': [], 't_t5': []}
for d in effets:
    p_t_n  = puissance_par_simulation(n_test, d, 'normale', test='t')
    p_mw_n = puissance_par_simulation(n_test, d, 'normale', test='wilcoxon')
    p_t_t5 = puissance_par_simulation(n_test, d, 't5', test='t')
    puissances['t_normale'].append(p_t_n)
    puissances['wilcoxon_normale'].append(p_mw_n)
    puissances['t_t5'].append(p_t_t5)
    print(f"{d:>10.1f} | {p_t_n:>10.3f} | {p_mw_n:>10.3f} | {p_t_t5:>10.3f}")

fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(effets, puissances['t_normale'],       'o-', color='steelblue', linewidth=2,
        label='Test t, loi normale')
ax.plot(effets, puissances['wilcoxon_normale'], 's-', color='seagreen', linewidth=2,
        label='Wilcoxon, loi normale')
ax.plot(effets, puissances['t_t5'],            '^-', color='orangered', linewidth=2,
        label='Test t, loi t₅ (queues épaisses)')
ax.axhline(0.80, color='gray', linestyle='--', linewidth=1, label='Puissance cible 80%')
ax.axhline(0.05, color='lightgray', linestyle=':', linewidth=1, label='Niveau α = 5%')
ax.set_xlabel('Taille d\'effet (Cohen\'s d)')
ax.set_ylabel('Puissance (probabilité de rejeter H₀)')
ax.set_title(f'Courbes de puissance simulées (n={n_test} par groupe, α=5%)')
ax.legend(fontsize=9)
ax.set_ylim(0, 1)
plt.tight_layout()
plt.show()
```

---

## Nombres pseudo-aléatoires et reproductibilité

### Générateurs modernes

Python utilise par défaut le **Mersenne Twister** (MT19937) — une période de $2^{19937}-1$, adapté à la plupart des simulations. NumPy recommande depuis la version 1.17 le **PCG64** (Permuted Congruential Generator), plus rapide et avec de meilleures propriétés statistiques.

```python
rng = np.random.default_rng(seed=42)  # PCG64, recommandé
# Ne pas utiliser : np.random.seed(42)  # MT19937, interface legacy
```

### Reproductibilité

Pour garantir la reproductibilité d'une simulation :

1. **Fixer le seed** au début de la simulation
2. **Documenter le seed** dans le rapport / notebook
3. **Ne pas réutiliser** le même générateur global dans des fonctions appelées en parallèle

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

# Quasi-Monte Carlo : convergence plus rapide
# (mention uniquement — pas de bibliothèque standard)
from scipy.stats import qmc

# Halton vs uniforme : convergence de l'estimation de π
def estimer_pi_qmc(n, engine='halton'):
    if engine == 'halton':
        sampler = qmc.Halton(d=2, scramble=True, seed=42)
    else:
        sampler = qmc.LatinHypercube(d=2, seed=42)
    pts = sampler.random(n) * 2 - 1  # espace [-1,1]²
    return 4 * np.mean(pts[:, 0]**2 + pts[:, 1]**2 <= 1)

n_vals = np.logspace(2, 5, 40).astype(int)
errs_mc   = []
errs_qmc  = []

for n in n_vals:
    est_mc, _ = estimer_pi(n, graine=n)
    est_qmc = estimer_pi_qmc(n)
    errs_mc.append(abs(est_mc - np.pi))
    errs_qmc.append(abs(est_qmc - np.pi))

fig, ax = plt.subplots(figsize=(9, 4))
ax.loglog(n_vals, errs_mc, 'steelblue', linewidth=2, alpha=0.8, label='Monte Carlo standard')
ax.loglog(n_vals, errs_qmc, 'orangered', linewidth=2, alpha=0.8, label='Quasi-Monte Carlo (Halton)')

# Références de convergence
n0 = n_vals[10]
ax.loglog(n_vals, errs_mc[10] * np.sqrt(n0/np.array(n_vals)), 'gray', linestyle='--',
          linewidth=1, label='∝ 1/√N')
ax.loglog(n_vals, errs_qmc[10] * (n0/np.array(n_vals)), 'gray', linestyle=':',
          linewidth=1, label='∝ 1/N')
ax.set_xlabel('N (log)'); ax.set_ylabel('Erreur absolue (log)')
ax.set_title('Monte Carlo vs Quasi-Monte Carlo — convergence pour l\'estimation de π\n'
             '(QMC converge en O(1/N) vs O(1/√N) pour MC standard)')
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()
```

```{note}
Le Quasi-Monte Carlo utilise des **séquences à faible discrépance** (Halton, Sobol, Latin Hypercube) qui remplissent l'espace de façon plus uniforme que les nombres pseudo-aléatoires, d'où une convergence théorique en $O((\log N)^d / N)$ au lieu de $O(1/\sqrt{N})$. En pratique, le gain est substantiel en faible dimension ($d \leq 10$).
```

---

## Résumé des méthodes computationnelles

| Méthode | Objectif principal | Hypothèses | Complexité |
|---------|--------------------|-----------|------------|
| Monte Carlo | Estimer E[f(X)], intégrer | Aucune | O(N) |
| Test de permutation | Tester H₀ : groupes échangeables | Échangeabilité | O(n × B) |
| Bootstrap non-param. | IC, biais, variance de tout estimateur | n suffisant, i.i.d. | O(n × B) |
| Bootstrap paramétrique | IC sous modèle paramétrique | Modèle correctement spécifié | O(n × B) |
| Jackknife | Biais, variance (lissé) | Statistique lisse | O(n²) |
| Validation croisée | Erreur de généralisation | i.i.d. | O(k × coût_modèle) |
| Simulation de puissance | Puissance de tests complexes | Modèle de simulation valide | O(nsim × n) |