Jugez la qualité d'un estimateur

Pendant la Seconde Guerre mondiale, les Alliés estimaient la production de chars allemands à partir des numéros de série des chars capturés. Le problème est le même pour un service de vélos en libre-service dont les vélos sont numérotés de 1 à $\theta$ : en voyant passer quelques numéros, combien de vélos compte la flotte ?

On peut imaginer plusieurs formules. Laquelle est la meilleure ? Pour trancher, il faut des critères. C’est l’objet de ce chapitre, et ces vélos numérotés seront notre fil rouge.

Estimateur et estimation

Vous observez 5 vélos, de numéros 112, 47, 305, 198 et 260. On modélise ces numéros comme un 5-échantillon $X_1, \dots, X_5$ de la loi uniforme sur $[0, \theta]$ (une approximation continue de « un numéro au hasard entre 1 et $\theta$ »). Le nombre $\theta$ est le paramètre inconnu.

Un estimateur de $\theta$ est une statistique $T_n = h(X_1, \dots, X_n)$ dont l’expression ne dépend pas de $\theta$. On le note souvent $\hat\theta$. Une estimation est la valeur $t_n$ qu’il prend sur les données observées.

Voici deux candidats raisonnables.

• Par la moyenne. Pour la loi uniforme sur $[0, \theta]$, l’espérance vaut $\theta/2$. D’où l’idée $T_n = 2 \bar X_n$. Ici $\bar x_5 = 184{,}4$, donc $t_5 = 368{,}8$.
• Par le maximum. Le plus grand numéro vu, $X_{(n)} = \max X_i$, est forcément inférieur à $\theta$. On le gonfle un peu : $T'_n = \frac{n+1}{n} X_{(n)}$. Ici $305 \times 6/5 = 366$.

Les deux estimations sont proches, mais sur un autre échantillon elles pourraient différer beaucoup. Comme $T_n$ est une variable aléatoire, c’est sa loi qu’il faut juger.

Le biais : viser juste en moyenne

Le biais d’un estimateur $T_n$ de $\theta$ est

$$b(T_n) = E[T_n] - \theta$$

où $E[T_n]$ est l’espérance de l’estimateur, c’est-à-dire la valeur moyenne qu’il prendrait sur une infinité d’échantillons. L’estimateur est sans biais si $b(T_n) = 0$, et asymptotiquement sans biais si $b(T_n) \to 0$ quand $n \to +\infty$.

Exemple

La moyenne empirique est sans biais pour l’espérance. Si $E[X] = m$, la linéarité de l’espérance donne

$$E[\bar X_n] = \frac{1}{n} \sum_{i=1}^{n} E[X_i] = \frac{1}{n} \times n m = m$$

Pour les vélos, $E[2 \bar X_n] = 2 \times \theta/2 = \theta$ : $T_n$ est sans biais.

Le maximum seul est biaisé. On montre que $E[X_{(n)}] = \frac{n}{n+1} \theta$, toujours inférieure à $\theta$. Son biais vaut $-\theta/(n+1)$ : il sous-estime la flotte. Le facteur $\frac{n+1}{n}$ de $T'_n$ corrige exactement ce biais.

Le cas célèbre de la variance

Voici le résultat qui justifie le $n - 1$ du chapitre précédent. Notons $\sigma^2 = V(X)$ et $m = E[X]$. On part de l’écriture $S_n^2 = \frac{1}{n}\sum X_i^2 - \bar X_n^2$ et de l’identité $E[Z^2] = V(Z) + E[Z]^2$ :

$$E[S_n^2] = E[X^2] - E[\bar X_n^2] = (\sigma^2 + m^2) - \left( \frac{\sigma^2}{n} + m^2 \right) = \frac{n-1}{n} \, \sigma^2$$

On a utilisé $V(\bar X_n) = \sigma^2/n$, qui vient de l’indépendance : $V\left(\frac{1}{n}\sum X_i\right) = \frac{1}{n^2} \sum V(X_i) = \frac{\sigma^2}{n}$.

Le biais de $S_n^2$ vaut donc $-\sigma^2/n$ : il est négatif (on sous-estime la dispersion) mais tend vers 0. En multipliant par $\frac{n}{n-1}$, on obtient $\tilde S_n^2$, dont l’espérance vaut exactement $\sigma^2$.

Avec n = 5, Sₙ² vise en moyenne 0,8 σ² : elle sous-estime la dispersion. Multiplier par n/(n − 1) recentre la distribution sur σ².

Mais pourquoi sous-estime-t-on ? Parce que les écarts sont mesurés par rapport à $\bar X_n$, qui est calculée sur les mêmes données et « colle » donc aux points plus que la vraie moyenne $m$. On dit qu’on a perdu un degré de liberté.

Attention

Un bon estimateur de $\theta$ ne donne pas forcément un estimateur sans biais de $\varphi(\theta)$. Par exemple, $\tilde S_n^2$ est sans biais pour $\sigma^2$, mais $\tilde S_n = \sqrt{\tilde S_n^2}$ est biaisé pour $\sigma$ : l’espérance d’une racine n’est pas la racine de l’espérance.

L’erreur quadratique moyenne : juste et précis

Un estimateur sans biais peut quand même être mauvais : s’il vise juste en moyenne mais s’éparpille énormément, une estimation isolée peut tomber très loin. On mesure donc l’erreur globale par l’erreur quadratique moyenne (EQM, ou risque quadratique, en anglais mean squared error) :

$$EQM(T_n) = E\left[(T_n - \theta)^2\right] = V(T_n) + b(T_n)^2$$

La décomposition vient du développement de $(T_n - E[T_n] + E[T_n] - \theta)^2$, dont le terme croisé a une espérance nulle. Elle dit que l’erreur se compose d’une partie variance (l’éparpillement) et d’une partie biais au carré (le décalage systématique).

Chaque point est l'estimation tirée d'un échantillon. Le biais décale le nuage, la variance l'étale : l'EQM pénalise les deux.

Entre deux estimateurs sans biais, on préfère celui de plus petite variance. Le meilleur possible, s’il existe, s’appelle estimateur sans biais de variance minimale (ESBVM). Mais un estimateur légèrement biaisé peut avoir une EQM plus faible qu’un estimateur sans biais. Les vélos en donnent un exemple frappant.

Exemple

Pour la loi uniforme sur $[0, \theta]$, $V(X) = \theta^2/12$. On calcule :

Estimateur	Biais	Variance	EQM pour $n = 5$
$2\bar X_n$	0	$\theta^2/(3n)$	$0{,}067\,\theta^2$
$X_{(n)}$	$-\theta/(n+1)$	$\frac{n\,\theta^2}{(n+1)^2(n+2)}$	$0{,}048\,\theta^2$
$\frac{n+1}{n} X_{(n)}$	0	$\theta^2/(n(n+2))$	$0{,}029\,\theta^2$

Le maximum brut, pourtant biaisé, fait mieux que $2\bar X_n$ sans biais. Et le maximum corrigé fait mieux que les deux. Une simulation de 200 000 échantillons de taille 5 retrouve ces trois valeurs à 0,001 près.

Le maximum brut est biaisé (sa moyenne tombe vers 250) mais bien moins dispersé que 2 X̄ₙ ; une fois corrigé, il vise juste et reste serré : la plus petite EQM.

import numpy as np

rng = np.random.default_rng(0)
theta, n, R = 300, 5, 200_000
X = rng.uniform(0, theta, size=(R, n))       # R échantillons de taille n

estimateurs = {
    "2 * moyenne": 2 * X.mean(axis=1),
    "maximum": X.max(axis=1),
    "maximum corrigé": (n + 1) / n * X.max(axis=1),
}
for nom, T in estimateurs.items():
    eqm = np.mean((T - theta) ** 2) / theta**2
    print(f"{nom:>16} : biais {T.mean() - theta:+6.1f}, EQM / theta^2 = {eqm:.3f}")

À retenir

$EQM = \text{variance} + \text{biais}^2$. Un estimateur sans biais n’est pas automatiquement le meilleur : on compare les estimateurs par leur EQM. C’est exactement le compromis biais-variance de l’apprentissage automatique, où la régularisation accepte un peu de biais pour réduire beaucoup la variance.

Que se passe-t-il quand l’échantillon grandit ?

On attend d’un estimateur raisonnable qu’il s’améliore avec plus de données. Il y a plusieurs façons de le formaliser.

Convergence en moyenne quadratique

$T_n$ converge en moyenne quadratique vers $\theta$ si $EQM(T_n) \to 0$. Pour un estimateur sans biais, il suffit que sa variance tende vers 0. C’est le cas de $\bar X_n$, dont la variance $\sigma^2/n$ tend vers 0 dès que $\sigma^2$ est finie.

Consistance (convergence en probabilité)

$T_n$ est consistant (ou convergent) si, pour toute marge $a > 0$,

$$P\left(\lvert T_n - \theta \rvert > a\right) \xrightarrow[n \to +\infty]{} 0$$

Autrement dit, la probabilité de se tromper de plus de $a$ devient négligeable. L’outil de base pour le montrer est l’inégalité de Bienaymé-Tchebychev : pour toute variable $X$ de variance finie et tout $a > 0$,

$$P\left(\lvert X - E[X] \rvert \ge a\right) \le \frac{V(X)}{a^2}$$

On en déduit que la convergence en moyenne quadratique entraîne la consistance.

Quand n grandit, la loi de l'estimateur se resserre autour de θ : la probabilité de se tromper de plus de 20 vélos tend vers 0. C'est la consistance.

Exemple

Pour estimer la proportion $\theta$ d’usagers satisfaits (loi de Bernoulli), on prend $\bar X_n$. Bienaymé-Tchebychev donne

$$P\left(\lvert \bar X_n - \theta \rvert > a\right) \le \frac{\theta(1-\theta)}{n a^2} \le \frac{1}{4 n a^2}$$

car $\theta(1-\theta) \le 1/4$. Avec $n = 1\,000$ et $a = 0{,}05$, la probabilité de se tromper de plus de 5 points est au plus $0{,}1$. La borne tend vers 0 : $\bar X_n$ est consistant.

Une propriété très pratique : si $T_n$ est consistant pour $\theta$ et $f$ est continue, alors $f(T_n)$ est consistant pour $f(\theta)$. C’est ainsi qu’on voit directement que $2\bar X_n$ est consistant pour $\theta$ dans le modèle uniforme.

Convergence en loi et théorème central limite

La consistance dit que $T_n$ se rapproche de $\theta$, mais pas comment il fluctue autour. Pour cela, on utilise la convergence en loi : $Z_n$ converge en loi vers $Z$ si $P(Z_n \le t) \to P(Z \le t)$ en tout point $t$ où la fonction de répartition de $Z$ est continue.

Le résultat phare est le théorème central limite (TCL) : si $X_1, \dots, X_n$ sont i.i.d. d’espérance $m$ et de variance $\sigma^2$ finie,

$$\frac{\sqrt{n}\,(\bar X_n - m)}{\sigma} \xrightarrow[n \to +\infty]{\mathcal{L}} \mathcal{N}(0, 1)$$

Quelle que soit la loi de départ, la moyenne fluctue autour de $m$ selon une loi à peu près normale, d’écart-type $\sigma/\sqrt{n}$. C’est ce résultat qui rendra possibles les intervalles de confiance et les tests asymptotiques.

Partie d'une loi très dissymétrique, la moyenne standardisée prend déjà une forme de cloche à n = 30 : c'est le théorème central limite à l'œuvre.

Plus généralement, on dit que $T_n$ est asymptotiquement normal si $\sqrt{n}\,(T_n - \theta)$ converge en loi vers $\mathcal{N}(0, \Sigma(\theta))$ ; $\Sigma(\theta)$ est la variance asymptotique.

La delta-méthode

Et si l’on s’intéresse à une fonction du paramètre ? La delta-méthode répond : si $\sqrt{n}\,(T_n - \theta) \to \mathcal{N}(0, \Sigma(\theta))$ en loi et si $g$ est dérivable avec $g'(\theta) \ne 0$, alors

$$\sqrt{n}\,\big(g(T_n) - g(\theta)\big) \xrightarrow{\mathcal{L}} \mathcal{N}\big(0,\ \Sigma(\theta)\, g'(\theta)^2\big)$$

L’idée est un développement limité : près de $\theta$, $g(T_n) - g(\theta) \approx g'(\theta)(T_n - \theta)$.

Exemple

Le nombre de vélos déposés en une minute à une station suit une loi de Poisson de paramètre $\theta$ (espérance et variance égales à $\theta$). Par le TCL, $\sqrt{n}\,(\bar X_n - \theta) \to \mathcal{N}(0, \theta)$.

On veut estimer la probabilité qu’aucun vélo n’arrive en une minute, $P(X = 0) = e^{-\theta}$. On prend $g(x) = e^{-x}$, de dérivée $-e^{-x}$ : l’estimateur $e^{-\bar X_n}$ a pour variance asymptotique $\theta \, e^{-2\theta}$.

Sur $n = 100$ minutes avec $\bar x = 2{,}3$, on estime $P(X = 0) \approx e^{-2{,}3} \approx 0{,}100$, avec un écart-type approché de $\sqrt{2{,}3}\, e^{-2{,}3}/\sqrt{100} \approx 0{,}015$.

Près de θ, la courbe se confond avec sa tangente : l'écart-type de g(X̄ₙ) est celui de X̄ₙ multiplié par |g′(θ)| = exp(−2,3) ≈ 0,1.

En entretien

« Un estimateur se juge par son erreur quadratique moyenne, qui vaut variance plus biais au carré. On lui demande aussi d’être consistant, c’est-à-dire de se rapprocher de la vraie valeur quand les données augmentent, et le théorème central limite dit comment il fluctue autour : en loi normale, avec un écart-type en $1/\sqrt{n}$. »

En résumé

• Un estimateur est une statistique qui ne dépend pas du paramètre ; une estimation est sa valeur sur les données.
• Biais : $E[T_n] - \theta$. $\bar X_n$ est sans biais pour l’espérance ; $S_n^2$ est biaisée ($-\sigma^2/n$), d’où la variance corrigée $\tilde S_n^2$.
• EQM $= V(T_n) + b(T_n)^2$ : un estimateur biaisé peut battre un estimateur sans biais.
• Convergence en moyenne quadratique ($EQM \to 0$) entraîne la consistance, que l’on montre souvent avec Bienaymé-Tchebychev.
• Le TCL donne la loi approchée de $\bar X_n$ : normale, d’écart-type $\sigma/\sqrt{n}$ ; la delta-méthode l’étend à $g(T_n)$.

Exercices

Exercice · Deux estimateurs de la moyenne

facile

Avec un 3-échantillon $X_1, X_2, X_3$ d’espérance $m$ et de variance $\sigma^2$, on propose $T = \bar X_3$ et $T' = \frac{1}{2} X_1 + \frac{1}{4} X_2 + \frac{1}{4} X_3$. Sont-ils sans biais ? Lequel préférer ?

Voir la solution

Par linéarité, $E[T] = m$ et $E[T'] = \frac{1}{2}m + \frac{1}{4}m + \frac{1}{4}m = m$ : les deux sont sans biais. On compare donc leurs variances, grâce à l’indépendance :

• $V(T) = \sigma^2/3 \approx 0{,}333\,\sigma^2$ ;
• $V(T') = \frac{1}{4}\sigma^2 + \frac{1}{16}\sigma^2 + \frac{1}{16}\sigma^2 = \frac{3}{8}\sigma^2 = 0{,}375\,\sigma^2$.

On préfère $T$ : à biais égal (nul), sa variance est plus faible. Plus généralement, parmi les moyennes pondérées de poids positifs de somme 1, les poids égaux minimisent la variance.

En code

Python

import numpy as np

# Variances exactes (sigma^2 = 1) : somme des poids au carré
poids_T = np.array([1/3, 1/3, 1/3])
poids_Tp = np.array([1/2, 1/4, 1/4])
print(np.sum(poids_T**2), np.sum(poids_Tp**2))     # 0.333  0.375

# Vérification par simulation : m = 10, sigma = 1
rng = np.random.default_rng(0)
X = rng.normal(10, 1, size=(100_000, 3))
T, Tp = X @ poids_T, X @ poids_Tp
print(T.mean(), Tp.mean())                         # 10.000  10.001 (sans biais)
print(T.var(), Tp.var())                           # 0.333  0.376

R

# Variances exactes (sigma^2 = 1) : somme des poids au carré
poids_T  <- c(1/3, 1/3, 1/3)
poids_Tp <- c(1/2, 1/4, 1/4)
sum(poids_T^2); sum(poids_Tp^2)                 # 0.333  0.375

# Vérification par simulation : m = 10, sigma = 1
set.seed(0)
X <- matrix(rnorm(3 * 100000, mean = 10, sd = 1), ncol = 3)
T  <- drop(X %*% poids_T)
Tp <- drop(X %*% poids_Tp)
mean(T); mean(Tp)                               # 9.999  9.999 (sans biais)
var(T); var(Tp)                                 # 0.333  0.376

Exercice · Consistance d'un estimateur uniforme

moyen

Dans le modèle uniforme sur $[0, \theta]$, montrez que $T_n = 2\bar X_n$ converge en moyenne quadratique vers $\theta$, puis qu’il est consistant.

Voir la solution

$T_n$ est sans biais, donc $EQM(T_n) = V(T_n) = 4\,V(\bar X_n) = 4 \times \frac{\theta^2}{12 n} = \frac{\theta^2}{3n}$, qui tend vers 0 : il y a convergence en moyenne quadratique.

Par Bienaymé-Tchebychev, pour tout $a > 0$ : $P(\lvert T_n - \theta \rvert \ge a) \le \frac{\theta^2}{3 n a^2} \to 0$. L’estimateur est donc consistant.

On pouvait aussi dire : $\bar X_n$ est consistant pour $\theta/2$ et $f(x) = 2x$ est continue, donc $2\bar X_n$ est consistant pour $\theta$.

En code

Python

import numpy as np

rng = np.random.default_rng(0)
theta, a, R = 300, 20, 20_000
for n in [10, 100, 1000]:
    T = 2 * rng.uniform(0, theta, size=(R, n)).mean(axis=1)
    eqm = np.mean((T - theta)**2)
    p_erreur = np.mean(np.abs(T - theta) >= a)
    borne = theta**2 / (3 * n * a**2)            # Bienaymé-Tchebychev
    print(n, round(eqm), round(theta**2 / (3 * n)), p_erreur, round(borne, 3))
# n     EQM simulée  theta^2/(3n)  P(|T - theta| >= 20)  borne
# 10    2994         3000          0.712                 7.5   (borne inutile ici)
# 100    299          300          0.248                 0.75
# 1000    30           30          0.00015               0.075 (tout tend vers 0)

R

set.seed(0)
theta <- 300; a <- 20; R <- 20000
for (n in c(10, 100, 1000)) {
  T <- 2 * rowMeans(matrix(runif(R * n, 0, theta), nrow = R))
  eqm <- mean((T - theta)^2)
  p_erreur <- mean(abs(T - theta) >= a)
  borne <- theta^2 / (3 * n * a^2)             # Bienaymé-Tchebychev
  cat(n, round(eqm), theta^2 / (3 * n), p_erreur, borne, "\n")
}
## n     EQM simulée  theta^2/(3n)  P(|T - theta| >= 20)  borne
## 10    2983         3000          0.717                 7.5   (borne inutile ici)
## 100    306          300          0.257                 0.75
## 1000    30           30          0.0004                0.075 (tout tend vers 0)

Exercice · Delta-méthode pour une loi exponentielle

difficile

Les durées d’appel à une hotline suivent une loi exponentielle de paramètre $\lambda$ : $E[X] = 1/\lambda$ et $V(X) = 1/\lambda^2$. On estime $\lambda$ par $1/\bar X_n$. Donnez sa variance asymptotique.

Voir la solution

Par le TCL, $\sqrt{n}\,(\bar X_n - 1/\lambda) \to \mathcal{N}(0, 1/\lambda^2)$. On applique la delta-méthode avec $g(x) = 1/x$, donc $g'(x) = -1/x^2$, évaluée en $1/\lambda$ : $g'(1/\lambda) = -\lambda^2$.

La variance asymptotique vaut $\frac{1}{\lambda^2} \times (\lambda^2)^2 = \lambda^2$. Ainsi $\sqrt{n}\,(1/\bar X_n - \lambda) \to \mathcal{N}(0, \lambda^2)$ : l’erreur-type approchée de $\hat\lambda$ est $\hat\lambda/\sqrt{n}$. Nous retrouverons ce résultat au chapitre suivant par le maximum de vraisemblance.

En code

Python

import numpy as np

rng = np.random.default_rng(0)
lam, n, R = 2.0, 500, 20_000
X = rng.exponential(scale=1 / lam, size=(R, n))
lam_hat = 1 / X.mean(axis=1)                    # estimateur 1 / moyenne
print(np.var(np.sqrt(n) * (lam_hat - lam)))     # 4.02, proche de lambda^2 = 4
print(lam_hat.std(), lam / np.sqrt(n))          # 0.0896  0.0894 : simulée puis lambda / sqrt(n)

R

set.seed(0)
lam <- 2; n <- 500; R <- 20000
X <- matrix(rexp(R * n, rate = lam), nrow = R)
lam_hat <- 1 / rowMeans(X)                      # estimateur 1 / moyenne
var(sqrt(n) * (lam_hat - lam))                  # 3.98, proche de lambda^2 = 4
sd(lam_hat); lam / sqrt(n)                      # 0.0892  0.0894 : simulée puis lambda / sqrt(n)