Statistiques computationnelles

Contenu

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.

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$.

_images/70810a660c0b1d4171d8538ddb750086faa8661b7353069105427232f5cec01c.png

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.

_images/cb36f7fc21883933dcb11d364b7a3d963433de022db3d584915d88cac5907eca.png

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).

Show code cell source

Hide code cell source

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()

Comparaison des méthodes de test
  Différence observée (B - A) : -27.57 ms

                  Méthode |  Statistique |   p-valeur
-------------------------------------------------------
      Test de permutation |  |Δ| observé |     0.0011
        Test t de Student |        2.457 |     0.0171
             Mann-Whitney |          567 |     0.0226

_images/d08b6bbec7b9fcda0a93c6edba334266c5b7af7e396494de8158460c1368814c.png

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 :

Calculer la statistique $\hat{\theta}$ sur l’échantillon original
Pour $b = 1, \ldots, B$ : tirer un échantillon bootstrap $X^{*b}$ de taille $n$ avec remise depuis $X$
Calculer $\hat{\theta}^{*b}$ sur chaque échantillon bootstrap
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.

Show code cell source

Hide code cell source

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()

_images/83bfb672d0bc86947e64fd49744ebc9fefd74d142d3983af414e7466e1bcdf53.png

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.

Bootstrap paramétrique vs IC analytique (loi log-normale)
  μ_chapeau = 8.4239,  σ_chapeau = 0.4922
  IC 95% bootstrap paramétrique : [8.3252, 8.5389]
  IC 95% analytique (normale)   : [8.3161, 8.5318]

_images/7c42e9955ffddc9043ba8637a95946a15af4152cd9071af0dd3064fa2add0a6b.png

Erreur bootstrap vs taille d’échantillon#

_images/30077d3e9c0dcd5b9ec9e176976465e890df590f47e1f5a134ac6ba8574dfe30.png

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.

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}")

Comparaison Bootstrap vs Jackknife (n=50, loi exponentielle)

         Statistique |    Observé |  Biais Jack. |   SE Jack. |   SE Boot.
------------------------------------------------------------------------
             Moyenne |     2.3522 |     +0.00000 |    0.28805 |    0.27426
             Médiane |     1.9723 |     -0.00000 |    0.74851 |    0.38221
          Écart-type |     2.0163 |     -0.03726 |    0.25936 |    0.23917

            Skewness |     1.0850 |     -0.11739 |    0.40985 |    0.32628

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).

Show code cell source

Hide code cell source

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()

Validation croisée — Régression linéaire (n=100, p=5)

   Stratégie |   RMSE moyen |         SE |  N folds
--------------------------------------------------
         LOO |      15.4868 |    12.3368 |      100
      5-fold |      19.3773 |     1.2249 |        5
     10-fold |      19.4105 |     2.3846 |       10
     20-fold |      18.9748 |     5.5047 |       20

_images/41bce143e6757cc960b3890512d18b7bd0fb2103feb613e8e8f6b2157e9d7fab.png

Simulation de puissance#

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

Show code cell source

Hide code cell source

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()

Simulation de puissance (n=30 par groupe, α=5%, 2000 simulations)

 Cohen's d |  t (norm.) |   Wilcoxon |     t (t₅)
------------------------------------------------

       0.1 |      0.048 |      0.048 |      0.058

       0.2 |      0.100 |      0.108 |      0.092

       0.3 |      0.172 |      0.160 |      0.166

       0.4 |      0.288 |      0.260 |      0.246

       0.5 |      0.434 |      0.398 |      0.352

       0.6 |      0.590 |      0.558 |      0.476

       0.7 |      0.730 |      0.706 |      0.568

       0.8 |      0.842 |      0.836 |      0.676

       0.9 |      0.908 |      0.892 |      0.750

       1.0 |      0.962 |      0.950 |      0.848

       1.1 |      0.990 |      0.990 |      0.912

_images/b396fb181b30d4dcdda90cd5d3d2eb740619aff62d67fd8df4fd984649033102.png

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.

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 :

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

_images/1ef24981b935ae0e6af69b332b30367c695906769f676cd2c9192db7dfc9ffad.png

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)