Fondations : estimateur, MSE, CRLB

Cette page pose les outils pour juger un estimateur, sans encore dire comment en construire un. La section 1 définit ce qu'est un estimateur et montre, sur un exemple, qu'il s'agit d'une variable aléatoire. La section 2 en tire les trois quantités qui le caractérisent, biais, variance et MSE. La section 3 montre comment les estimer par simulation lorsque le calcul analytique n'aboutit pas. La section 4 reporte cette incertitude sur le résultat lui-même, sous la forme d'un intervalle de confiance. La section 5 lève la dernière hypothèse commode, celle d'une variance connue.

1. Problématique

⚠ Bascule de notation

Le chapitre 1 écrivait $X$ (majuscule) la variable aléatoire et $x$ (minuscule) sa réalisation, une distinction utile pour apprendre à voir le hasard. À partir de maintenant, nous suivons la convention de Kay : une seule casse, $x[n]$, pour la variable aléatoire et sa réalisation ; le contexte tranche. De même $\hat{\theta }$ désignera aussi bien l'estimateur (une fonction des données, donc une v.a.) que l'estimation qu'il produit sur un jeu de données particulier (un nombre).

Modèle d'estimateur

Considérons que nous disposons de données $\mathbf{x}=[x[0],\dots ,x[N-1]{]}^T$, modélisées par une loi $p(\mathbf{x};\theta )$ à un paramètre $\theta$ inconnu. Un estimateur de $\theta$ est une fonction des données censée approcher la valeur de $\theta$. Mathématiquement, un estimateur peut s'exprimer sous la forme d'une fonction :

$$\hat{\theta }=g(\mathbf{x})$$

où

• $g$ est une fonction des données, choisie par nous,
• $\hat{\theta }$ est lui-même une variable aléatoire, puisque fonction de données aléatoires : il possède une loi, une espérance et une variance.

Confusion

Attention à ne pas confondre la valeur vraie $\theta$ et son estimation $\hat{\theta }$.

Exemple : estimation de la borne supérieure

Une machine produit des pièces dont la cote varie uniformément entre $0$ et une valeur maximale $\theta$ qui lui est propre ; nous mesurons $N$ pièces et cherchons $\theta$. Le modèle est donc $N$ mesures i.i.d. $x[n]\sim \mathcal{U}(0,\theta )$, une loi uniforme continue (fiche uniforme) dont le seul paramètre inconnu est aussi la borne supérieure du support.

Répétons trois fois la même expérience, huit mesures à chaque fois, et plaçons les valeurs obtenues sur un axe.

Figure 2.1a, Trois répétitions indépendantes de la même expérience, $N=8$ mesures tirées selon $\mathcal{U}(0,\theta )$ avec $\theta =10$ (script figures/ch2/01-echantillon-uniforme.py).

Les points se répartissent sans privilégier aucune zone, c'est le propre de la loi uniforme, et aucun ne dépasse $\theta$, qui borne le support. Chaque série suggère donc une réponse : $\theta$ se trouve quelque part au-dessus de la plus grande valeur observée. Mais cette plus grande valeur change d'une série à l'autre, $9,1$, $9,7$, puis seulement $6,6$, et elle reste toujours en deçà de $\theta$.

Trois candidats se présentent pour estimer $\theta$ :

$$\begin{array}{rl}{\hat{\theta }}_1& =\underset{n}{max}x[n],\\ {\hat{\theta }}_2& =\underset{n}{max}x[n]-\underset{n}{min}x[n],\\ {\hat{\theta }}_3& =2\overline{x}=\frac{2}{N}\sum _{n}x[n].\end{array}$$

${\hat{\theta }}_1$ mise sur l'observation la plus informative : aucune mesure ne peut dépasser $\theta$, donc la plus grande mesure observée doit en être proche. ${\hat{\theta }}_2$ retranche la plus petite mesure, dans l'idée de corriger un possible décalage. ${\hat{\theta }}_3$ mobilise toutes les mesures à la fois via la moyenne empirique (chapitre 1) : $E(\overline{x})=\theta /2$, donc ${\hat{\theta }}_3$ vise juste en espérance. Trois fonctions légitimes des données, donc trois variables aléatoires. Laquelle choisir ? Il faut un critère de comparaison.

Illustration

Fixons $\theta =10$ et tirons une seule série de 200 mesures ; à chaque $N$, recalculons les trois estimateurs sur les $N$ premières mesures de cette même série.

Figure 2.1b, Les trois estimateurs recalculés sur les $N$ premières mesures d'un même tirage ($\theta =10$), comparés à la valeur vraie (script figures/ch2/01-fondations.py).

Les trois courbes se rapprochent de $\theta =10$ quand $N$ grandit, mais à des vitesses très différentes, ${\hat{\theta }}_1$ colle à la cible presque immédiatement, ${\hat{\theta }}_3$ tâtonne plus longtemps. Ce graphe n'est qu'une seule réalisation pour une meme valeur de $\theta$ : un autre tirage déplacerait les trois courbes, sans changer leur allure générale. Comparer des estimateurs sur un tirage unique est donc fragile. Il faut un critère qui résume leur comportement sur l'ensemble des tirages possibles.

2. Métriques de performance : Biais, variance et MSE

Comment résumer, en un seul nombre, la qualité d'un estimateur sur l'ensemble de ses tirages possibles ? La mesure la plus naturelle est une distance moyenne à la cible, l'erreur quadratique moyenne, l'espérance du carré de l'écart entre l'estimateur et la vraie valeur.

Définition

Définition, Erreur quadratique moyenne (MSE)

$$\mathrm{MSE}(\hat{\theta })=E((\hat{\theta }-\theta {)}^2).$$

où

• $\hat{\theta }=g(\mathbf{x})$ est l'estimateur évalué,
• $\theta$ est la valeur vraie, inconnue en pratique.

⚠ Statut de $\theta$

La définition de $\mathrm{MSE}(\hat{\theta })$ suppose $\theta$ déterministe mais inconnu, l'approche fréquentiste, celle que rappelle le point-virgule de $p(\mathbf{x};\theta )$. Le §4 des estimateurs usuels assignera au contraire une loi à $\theta$, et la notation basculera alors en $p(\mathbf{x}\mid \theta )$.

Cette quantité peut se décomposer en deux termes de natures différentes.

Définition, Biais et variance

Pour un estimateur $\hat{\theta }$ de $\theta$ :

$$\mathrm{MSE}(\hat{\theta })=\mathrm{var}(\hat{\theta })+b(\hat{\theta }{)}^2$$

où

• $b(\hat{\theta })=E(\hat{\theta })-\theta$ correspond au biais de l'estimateur,
• $\mathrm{var}(\hat{\theta })=E((\hat{\theta }-E(\hat{\theta }){)}^2)$ correspond à la variance de l'estimateur.

Preuve

Posons $m=E(\hat{\theta })$ et insérons-le :

$$\hat{\theta }-\theta =(\hat{\theta }-m)+(m-\theta ).$$

Le second terme, $m-\theta$, est une constante (il ne dépend pas des données). En développant le carré et en prenant l'espérance, le terme croisé s'annule car $E(\hat{\theta }-m)=0$ :

$$\mathrm{MSE}(\hat{\theta })=\underset{\text{dispersion autour de }m}{\underbrace{E((\hat{\theta }-m{)}^2)}}+2(m-\theta )\underset{=0}{\underbrace{E(\hat{\theta }-m)}}+\underset{\text{écart entre }m\text{ et }\theta }{\underbrace{(m-\theta {)}^2}}.◼$$

Le biais mesure une erreur systématique (nous visons à côté, toujours du même côté) ; la variance mesure une erreur de fluctuation (nous dispersons autour de notre propre cible). La MSE les réunit en un seul nombre, mais un nombre ne dit pas lequel des deux défauts pèse le plus.

Propriétés

En pratique, nous n'aurons pas de technique sur étagère pour choisir l'estimateur optimal $g(.)$. Néanmoins, dans la classe des estimateurs sans biais, il sera possible de comparer la variance de n'importe quel estimateur avec une borne présentant la variance optimale. Lorsqu'un estimateur sans biais atteint cette borne, il deviendra de facto optimal.

Théorème, Borne de Cramér-Rao (CRLB)

Pour tout estimateur sans biais $\hat{\theta }=g(\mathbf{x})$ de $\theta$ (sous conditions de régularité) :

$$\mathrm{var}(\hat{\theta })\ge \frac{1}{I(\theta )}.$$

où $I(\theta )$ correspond à l'information de Fisher et est définie par :

$$I(\theta )=-E\left[\frac{{\partial }^2\ell (\theta )}{\partial {\theta }^2}\right].$$

et $\ell (\theta )=\ln p(\mathbf{x};\theta )$ correspond à la log-vraisemblance.

Illustration

Cette décomposition se lit aussi sur un schéma classique : imaginez trois séries de tirs sur une cible, chaque tir étant une estimation $\hat{\theta }$ obtenue sur un tirage différent, et le centre marquant la valeur vraie $\theta$.

Figure 2.1c, Trois cibles, de gauche à droite : tirs groupés mais décalés du centre (biais fort, variance faible) ; tirs dispersés mais centrés en moyenne (sans biais, variance forte) ; tirs groupés et proches du centre malgré un léger décalage (biais faible, MSE la plus faible des trois) (script figures/ch2/01-fondations.py).

Le biais déplace le nuage de tirs loin du centre ; la variance l'étale autour de sa propre moyenne. Les cibles de gauche et du centre n'illustrent chacune qu'un seul défaut, et leur MSE reste élevée, celle de gauche à cause du biais, celle du centre à cause de la variance. La cible de droite cumule un peu des deux, un biais non nul et une variance non nulle, mais tous deux petits : sa MSE est la plus faible des trois. Un peu de biais, bien maîtrisé, peut valoir mieux qu'aucun biais du tout.

Exemple : estimation de la borne supérieure

En utilisant les propriétés de l'espérance mathématique, il est parfois possible d'obtenir analytiquement l'expression de la MSE, du biais et de la variance d'un estimateur.

Pour les trois candidats du §1 ($x[n]$ i.i.d. $\mathcal{U}(0,\theta )$), le tableau suivant présente les propriétés de chaque estimateur.

	${\hat{\theta }}_1=max$	${\hat{\theta }}_2=max-min$	${\hat{\theta }}_3=2\overline{x}$
Biais	$-{\displaystyle \frac{\theta }{N+1}}$	$-{\displaystyle \frac{2\theta }{N+1}}$	$0$
Variance	${\displaystyle \frac{N{\theta }^2}{(N+1{)}^2(N+2)}}$	${\displaystyle \frac{2(N-1){\theta }^2}{(N+1{)}^2(N+2)}}$	${\displaystyle \frac{{\theta }^2}{3N}}$
MSE	${\displaystyle \frac{2{\theta }^2}{(N+1)(N+2)}}$	${\displaystyle \frac{6{\theta }^2}{(N+1)(N+2)}}$	${\displaystyle \frac{{\theta }^2}{3N}}$

Calcul des moments de ${\hat{\theta }}_1$ et ${\hat{\theta }}_3$

Sous indépendance, $P({\hat{\theta }}_1\le x)=P(\text{tous les }x[n]\le x)=F_X(x{)}^N=(x/\theta {)}^N$ sur $[0,\theta ]$, la CDF du maximum s'obtient en élevant celle d'une seule mesure à la puissance $N$ (fonction de répartition). En dérivant, la densité est $f_{(N)}(x)=N{x}^{N-1}/{\theta }^N$, d'où

$$E({\hat{\theta }}_1)={\int }_0^{\theta }x\frac{N{x}^{N-1}}{{\theta }^N}\mathrm{d}x=\frac{N}{N+1}\theta ,E({\hat{\theta }}_1^2)=\frac{N}{N+2}{\theta }^2,$$

et $\mathrm{var}({\hat{\theta }}_1)=E({\hat{\theta }}_1^2)-E({\hat{\theta }}_1{)}^2=N{\theta }^2/((N+1{)}^2(N+2))$ par König-Huygens. Pour ${\hat{\theta }}_3$, le corollaire sur la loi de la moyenne empirique donne $\mathrm{var}(\overline{x})={\sigma }^2/N$ avec ${\sigma }^2={\theta }^2/12$ (fiche uniforme), d'où $\mathrm{var}(2\overline{x})=4{\sigma }^2/N={\theta }^2/(3N)$. Les moments de ${\hat{\theta }}_2$ s'obtiennent de façon analogue, à partir de la loi jointe du couple (min, max), nous les admettons ici.

Un seul de ces estimateurs est sans biais, ${\hat{\theta }}_3$, et pourtant ce n'est pas lui qui présente la MSE la plus faible. En effet, $\mathrm{MSE}({\hat{\theta }}_1)$ et $\mathrm{MSE}({\hat{\theta }}_2)$ décroissent en $1/N^2$, bien plus vite que $\mathrm{MSE}({\hat{\theta }}_3)$, qui décroît seulement en $1/N$. Dès $N=3$, le maximum biaisé a donc une MSE plus petite que la moyenne sans biais, et l'écart se creuse ensuite. Entre les deux estimateurs biaisés, ${\hat{\theta }}_1$ l'emporte partout sur ${\hat{\theta }}_2$, d'un facteur $3$ exactement.

Illustration

Figure 2.1d, Biais, variance et MSE des trois estimateurs en fonction de $N$ ($\theta =10$), $N\ge 2$ ; variance et MSE en échelle log-log (script figures/ch2/01-fondations.py).

En échelle log-log, une décroissance en $1/N^k$ devient une droite de pente $-k$ : les MSE de ${\hat{\theta }}_1$ et ${\hat{\theta }}_2$ y sont deux fois plus pentues que celle de ${\hat{\theta }}_3$, l'écart en $1/N^2$ contre $1/N$ se lisant directement. Le biais reste en échelle linéaire, ses valeurs étant négatives et nulles pour ${\hat{\theta }}_3$.

3. Simulation de Monte-Carlo

Les trois définitions du §2 posent une difficulté pratique immédiate. Elles reposent toutes sur une espérance, rarement calculable analytiquement dès que l'estimateur s'écarte des quelques cas d'école. La simulation de Monte-Carlo contourne cet obstacle : elle génère $B$ réalisations de la variable aléatoire $\hat{\theta }$, puis remplace chaque espérance par une moyenne empirique sur ces $B$ valeurs.

Définition

Définition, Estimation Monte-Carlo du biais, de la variance et de la MSE

Soient ${\hat{\theta }}^{(1)},\dots ,{\hat{\theta }}^{(B)}$ les estimations obtenues sur $B$ répétitions indépendantes de l'expérience, et $\bar{\hat{\theta }}=\frac{1}{B}\sum _b{\hat{\theta }}^{(b)}$ leur moyenne empirique. Le biais, la variance et la MSE s'estiment par

$$\hat{b}=\bar{\hat{\theta }}-\theta ,$$$$\widehat{\mathrm{var}}=\frac{1}{B}\sum _{b=1}^B({\hat{\theta }}^{(b)}-\bar{\hat{\theta }}{)}^2,$$$$\widehat{\mathrm{MSE}}=\frac{1}{B}\sum _{b=1}^B({\hat{\theta }}^{(b)}-\theta {)}^2.$$

Ces trois expressions se justifient par la loi des grands nombres : une moyenne empirique se rapproche de l'espérance qu'elle estime à mesure que $B$ grandit. La qualité de l'estimation dépend donc directement de $B$. Une valeur trop faible laisse subsister des fluctuations, une valeur trop élevée allonge le temps de simulation.

Implémentation

La mise en œuvre suit toujours les trois mêmes étapes, quelle que soit la loi des données et quel que soit l'estimateur étudié :

1. simuler un jeu de $N$ mesures selon le modèle, avec une valeur de $\theta$ que nous fixons ;
2. estimer $\hat{\theta }$ sur ce jeu, en appliquant l'estimateur à évaluer ;
3. répéter $B$ fois, puis moyenner selon les trois formules ci-dessus.

Seules les deux premières étapes dépendent du problème traité. Elles s'isolent donc dans deux fonctions, simulate et g, que le reste du code ne fait qu'appeler :

import numpy as np
rng = np.random.default_rng(2026)

theta, N, B = 10.0, 8, 100_000

def simulate(theta, N):
    """Draw N measurements from the model, here x[n] ~ U(0, theta)."""
    return rng.uniform(0, theta, size=N)

def g(x):
    """The estimator under study: from data to an estimate of theta."""
    return x.max()

# steps 1 and 2, repeated B times: B estimates of theta
theta_hat = np.empty(B)
for b in range(B):
    x = simulate(theta, N)      # one repetition of the experiment
    theta_hat[b] = g(x)         # the corresponding estimate

# step 3: each expectation becomes an average over these B values
bias = theta_hat.mean() - theta
variance = theta_hat.var()
mse = ((theta_hat - theta) ** 2).mean()

Changer d'estimateur revient à changer le corps de g, changer de modèle à changer celui de simulate : le reste du code est inchangé. Lorsque l'estimateur s'exprime de façon vectorielle, comme ici avec max, il est plus rapide de remplacer la boucle par un unique tirage dans un tableau $B\times N$, puis d'appliquer l'estimateur ligne par ligne.

Exemple

Appliquons cette méthode à ${\hat{\theta }}_1=\underset{n}{max}x[n]$ pour $N$ variable, en confrontant sa MSE théorique à son estimation Monte-Carlo pour deux valeurs de $B$.

Figure 2.2, MSE de ${\hat{\theta }}_1=\underset{n}{max}x[n]$ en fonction de $N$ ($\theta =10$), valeur théorique $2{\theta }^2/((N+1)(N+2))$ et estimations Monte-Carlo pour $B=10$ et $B=100$, en échelle log-log (script figures/ch2/01-monte-carlo-convergence.py).

Les deux estimations suivent la décroissance théorique sans jamais s'en écarter systématiquement : elles oscillent autour de la courbe, tantôt au-dessus, tantôt en dessous. C'est la signature d'un estimateur sans biais de la MSE, et c'est ce qui rend la méthode fiable. Leur dispersion, en revanche, diffère nettement. À $B=10$, l'écart à la valeur théorique atteint $104\mathrm{\%}$ et vaut $46\mathrm{\%}$ en médiane : dix répétitions ne permettent pas de conclure. À $B=100$, il tombe à $17\mathrm{\%}$ en médiane. Multiplier $B$ par dix a donc divisé l'écart typique par $2,7$, proche du facteur $\sqrt{10}\approx 3,2$ attendu, la précision d'une moyenne empirique progressant en $1/\sqrt{B}$. Atteindre une décimale sûre exige de ce fait des $B$ de l'ordre de ${10}^4$, ce qui reste bon marché : le coût est du temps de calcul, jamais des mesures supplémentaires.

Cette technique dépasse largement l'exemple : dès qu'un estimateur est trop compliqué pour un calcul analytique, ce qui est la règle plutôt que l'exception, la simulation de Monte-Carlo reste disponible, et c'est alors le seul accès à son risque.

4. Intervalle de confiance

Le biais, la variance et la MSE jugent un estimateur sur l'ensemble de ses tirages possibles. Le résultat livré à l'utilisateur reste pourtant un nombre unique, une estimation dont rien n'indique la précision. Comment reporter, à côté de l'estimation, l'incertitude qui l'accompagne ? La réponse consiste à remplacer ce nombre par un intervalle, dont les deux bornes sont elles aussi calculées sur les données.

Définition

Définition, Intervalle de confiance

Un intervalle de confiance de niveau $1-\alpha$ pour $\theta$ est un intervalle aléatoire $[{\hat{\theta }}_{inf},{\hat{\theta }}_{sup}]$ vérifiant

$$P({\hat{\theta }}_{inf}\le \theta \le {\hat{\theta }}_{sup})=1-\alpha ,$$

où

• ${\hat{\theta }}_{inf}=g_{inf}(\mathbf{x})$ et ${\hat{\theta }}_{sup}=g_{sup}(\mathbf{x})$ sont deux fonctions des données, donc deux variables aléatoires,
• $\theta$ est la valeur vraie, déterministe et inconnue,
• $1-\alpha$ est le niveau de confiance, et $\alpha$ le risque toléré.

Dans cette probabilité, $\theta$ est une constante : l'aléa provient exclusivement des deux bornes, qui changent d'un jeu de données à l'autre. Le niveau $1-\alpha$ mesure donc une fréquence de couverture, la proportion des expériences répétées où l'intervalle construit contient la valeur vraie. C'est une propriété de la procédure, au même titre que le biais ou la MSE du §2.

Confusion classique

« Il y a $95\mathrm{\%}$ de chances que $\theta$ appartienne à $[2,23;2,35]$. » Non : une fois les données observées, les deux bornes sont des nombres, et $\theta$ en est un également. Soit $\theta$ appartient à cet intervalle, soit il ne lui appartient pas ; cette proposition est vraie ou fausse, et aucune probabilité intermédiaire ne lui est attachée.

Propriétés

Construire un intervalle de confiance exige de connaître la loi de l'estimateur. Le cas de la moyenne empirique d'un échantillon Gaussien se traite exactement, à partir de la loi de la moyenne empirique.

Proposition, Intervalle de confiance de la moyenne d'une Gaussienne

Soient $x[n]\sim \mathcal{N}(\theta ,{\sigma }^2)$ i.i.d. avec $\sigma$ connu, et $\hat{\theta }=\overline{x}$. L'intervalle

$$[\overline{x}-z_{1-\alpha /2}\frac{\sigma }{\sqrt{N}},\overline{x}+z_{1-\alpha /2}\frac{\sigma }{\sqrt{N}}]$$

est de niveau exactement $1-\alpha$, où

• $z_{1-\alpha /2}$ est le quantile d'ordre $1-\alpha /2$ de la loi $\mathcal{N}(0,1)$,
• $z_{\mathrm{0,975}}=1,96$ pour un niveau de confiance de $95\mathrm{\%}$.

Démonstration

La moyenne empirique de $N$ mesures Gaussiennes i.i.d. suit $\overline{x}\sim \mathcal{N}(\theta ,{\sigma }^2/N)$, donc la variable centrée réduite

$$z=\frac{\overline{x}-\theta }{\sigma /\sqrt{N}}$$

suit $\mathcal{N}(0,1)$, quel que soit $N$. Par symétrie de cette loi, $P(|z|\le z_{1-\alpha /2})=1-\alpha$. Isoler $\theta$ dans cette double inégalité donne l'intervalle annoncé : le paramètre n'a pas changé de statut, seule l'écriture a déplacé l'aléa des données vers les bornes. $◼$

La demi-largeur $z_{1-\alpha /2}\sigma /\sqrt{N}$ vaut $z_{1-\alpha /2}\sqrt{\mathrm{var}(\hat{\theta })}$ : l'intervalle n'est rien d'autre que l'écart-type du §2, rendu lisible sur l'axe des valeurs du paramètre. Sa décroissance en $1/\sqrt{N}$ fixe le prix de la précision, quadrupler le nombre de mesures pour diviser la largeur par deux. Deux restrictions accompagnent cette construction. La loi de $\hat{\theta }$ est rarement Gaussienne exactement, et $\sigma$ est rarement connu ; l'un et l'autre se remplacent alors par une approximation valable pour $N$ grand, comme le fera l'intervalle asymptotique du maximum de vraisemblance.

Confusion classique

« L'intervalle contient $95\mathrm{\%}$ des mesures. » Non : il porte sur le paramètre $\theta$, pas sur les données. Sa largeur décroît d'ailleurs quand $N$ augmente, alors que la dispersion des mesures, elle, ne dépend pas de leur nombre.

Implémentation

La couverture annoncée se vérifie par simulation, en reprenant les trois étapes du §3 :

1. simuler un jeu de $N$ mesures selon le modèle, avec une valeur de $\theta$ que nous fixons ;
2. construire l'intervalle sur ce jeu, puis tester s'il contient $\theta$ ;
3. répéter $B$ fois, puis compter la proportion d'intervalles qui le contiennent.

import numpy as np
rng = np.random.default_rng(2026)

theta, sigma, N, B = 2.30, 0.10, 10, 100_000
z = 1.96

def simulate(theta, sigma, N):
    """Draw N measurements from the model, here x[n] ~ N(theta, sigma^2)."""
    return rng.normal(theta, sigma, size=N)

def confidence_interval(x, sigma, z):
    """Two-sided interval centred on the sample mean, half-width z*sigma/sqrt(N)."""
    half_width = z * sigma / np.sqrt(len(x))
    return x.mean() - half_width, x.mean() + half_width

# steps 1 and 2, repeated B times: count the intervals that catch theta
covered = 0
for b in range(B):
    x = simulate(theta, sigma, N)               # one repetition of the experiment
    lower, upper = confidence_interval(x, sigma, z)
    if lower <= theta <= upper:
        covered += 1

# step 3: the coverage is a frequency over repetitions, not a probability on theta
coverage = covered / B                          # >>> 0.9495  (predicted: 0.95)

Ce calcul n'est possible qu'en simulation, où $\theta$ est fixé par nous. Face à des données réelles, $\theta$ reste inconnu et la couverture n'est jamais observable.

Exemple : mesure d'une tension

Un voltmètre bruité mesure $N=10$ fois la même tension continue, avec un écart-type de bruit $\sigma =0,10$ V connu par la notice de l'appareil. Les dix mesures donnent $\overline{x}=\mathrm{2,288}$ V, et la demi-largeur vaut $1,96\times 0,10/\sqrt{10}=\mathrm{0,062}$ V, d'où l'intervalle de confiance à $95\mathrm{\%}$

$$[\mathrm{2,226};\mathrm{2,350}]\text{V}.$$

Cet intervalle contient-il la vraie tension ? Sur des données réelles, la question reste sans réponse : le résultat obtenu est un intervalle numérique, et la vraie tension y est ou n'y est pas.

Illustration

En simulation, en revanche, la vraie tension est connue, et la couverture se lit directement. Si le niveau de $95\mathrm{\%}$ décrit bien la procédure, une centaine de répétitions doit produire une poignée d'intervalles qui manquent la cible.

Figure 2.3, Cent répétitions indépendantes de la mesure ($\theta =2,30$ V, $\sigma =0,10$ V, $N=10$) et les intervalles de confiance à $95\mathrm{\%}$ correspondants ; en rouge, ceux qui ne contiennent pas $\theta$ (script figures/ch2/01-intervalle-confiance.py).

Chaque segment est un intervalle obtenu sur un jeu de données différent. Tous ont la même largeur, $2\times \mathrm{0,062}$ V, puisque celle-ci ne dépend que de $\sigma$ et de $N$ ; seule leur position change, entraînée par $\overline{x}$. Sept segments sur cent ne coupent pas la verticale : pris isolément, chacun est simplement faux, sans que rien dans les données ne le signale. La couverture observée, $93\mathrm{\%}$, est elle-même une fréquence sur $B=100$ répétitions, et fluctue autour de $95\mathrm{\%}$ comme toute estimation Monte-Carlo.

5. Variance inconnue

L'intervalle du §4 suppose $\sigma$ connu, ce que la notice d'un appareil fournit parfois. Le cas courant est l'autre : $\sigma$ doit être estimé sur les mesures elles-mêmes, celles qui servent déjà à estimer $\theta$. Cette section lève cette hypothèse, d'abord en estimant la variance, puis en corrigeant l'intervalle qui en découle.

Définition

Définition, Variance empirique

Pour $N$ mesures i.i.d. de moyenne inconnue, la variance empirique est

$${\hat{\sigma }}^2=\frac{1}{N-1}\sum _{n=0}^{N-1}(x[n]-\overline{x}{)}^2,$$

où

• $\overline{x}$ est la moyenne empirique, estimée sur les mêmes données,
• $N-1$ est le nombre de degrés de liberté restants une fois $\overline{x}$ fixée.

Le diviseur $N-1$ n'est pas une convention d'usage. Les écarts sont mesurés à $\overline{x}$ et non à $\theta$, qui est inconnu ; or $\overline{x}$ est précisément la valeur qui minimise cette somme de carrés. La somme obtenue est donc systématiquement trop petite, et retirer un degré de liberté compense exactement ce rétrécissement.

Propriétés

Proposition, Biais des deux diviseurs

Pour $N$ mesures i.i.d. de variance ${\sigma }^2$ :

$$E[\sum _{n=0}^{N-1}(x[n]-\overline{x}{)}^2]=(N-1){\sigma }^2.$$

L'estimateur à diviseur $N-1$ est donc sans biais, tandis que celui à diviseur $N$ a pour biais $-{\sigma }^2/N$.

Démonstration

En insérant $\theta$ dans chaque écart et en développant, le terme croisé se regroupe :

$$\sum _n(x[n]-\overline{x}{)}^2=\sum _n(x[n]-\theta {)}^2-N(\overline{x}-\theta {)}^2.$$

Le premier terme a pour espérance $N{\sigma }^2$, chaque écart à $\theta$ étant de variance ${\sigma }^2$. Le second vaut $N\mathrm{var}(\overline{x})={\sigma }^2$ en espérance, par la loi de la moyenne empirique. La différence vaut $(N-1){\sigma }^2$. $◼$

Le biais du diviseur $N$ décroît en $1/N$ et devient négligeable pour $N$ grand, mais il atteint $-10\mathrm{\%}$ à $N=10$. Cette différence est celle du paramètre ddof de NumPy : np.var(x, ddof=1) calcule l'estimateur sans biais, et le défaut ddof=0 l'estimateur biaisé.

Confusion classique

« $\hat{\sigma }$ est un estimateur sans biais de $\sigma$. » Non : le diviseur $N-1$ ne corrige que ${\hat{\sigma }}^2$. La racine d'une espérance n'étant pas l'espérance d'une racine, $E(\hat{\sigma })<\sigma$, et l'écart-type empirique reste biaisé quel que soit le diviseur retenu.

Remplacer $\sigma$ par $\hat{\sigma }$ dans l'intervalle du §4 introduit un second aléa, celui de $\hat{\sigma }$, que le quantile Gaussien ne prend pas en compte. La correction porte un nom.

Proposition, Intervalle de confiance à variance inconnue

Soient $x[n]\sim \mathcal{N}(\theta ,{\sigma }^2)$ i.i.d. avec $\sigma$ inconnu. La variable

$$t=\frac{\overline{x}-\theta }{\hat{\sigma }/\sqrt{N}}$$

suit une loi de Student à $\nu =N-1$ degrés de liberté, et l'intervalle

$$[\overline{x}-t_{1-\alpha /2,\nu }\frac{\hat{\sigma }}{\sqrt{N}},\overline{x}+t_{1-\alpha /2,\nu }\frac{\hat{\sigma }}{\sqrt{N}}]$$

est de niveau exactement $1-\alpha$, où $t_{1-\alpha /2,\nu }$ est le quantile d'ordre $1-\alpha /2$ de cette loi.

Illustration

Si l'ignorance de $\sigma$ se paie, elle doit se lire sur la loi de Student : des queues plus épaisses que celles de la Gaussienne, donc un quantile plus éloigné de zéro à niveau de confiance égal.

Figure 2.4, Loi de Student à $\nu =9$ degrés de liberté et loi $\mathcal{N}(0,1)$, avec leurs quantiles d'ordre $\mathrm{0,975}$ ; en bleu, les queues de Student au-delà du quantile Gaussien (script figures/ch2/01-student-gaussienne.py).

Les deux densités se ressemblent au centre et diffèrent dans les queues, là où se lisent les quantiles. Les aires colorées représentent $8,2\mathrm{\%}$ de la masse de Student, alors que la Gaussienne n'en met que $5\mathrm{\%}$ au-delà des mêmes bornes : employer $z_{\mathrm{0,975}}=1,96$ avec $\hat{\sigma }$ ne couvre donc que $91,8\mathrm{\%}$ des expériences, trois points sous le niveau annoncé. Le quantile de Student, $2,26$, restaure les $95\mathrm{\%}$ en élargissant l'intervalle de $15\mathrm{\%}$. Cet élargissement s'atténue avec $N$ : $t_{\mathrm{0,975};29}=\mathrm{2,045}$ à $N=30$, soit $4\mathrm{\%}$ seulement, et la loi de Student tend vers la Gaussienne.

Implémentation

La couverture des deux intervalles se compare par la simulation du §4, en construisant les deux à chaque répétition :

1. simuler un jeu de $N$ mesures, avec des $\theta$ et $\sigma$ que nous fixons ;
2. estimer $\overline{x}$ et $\hat{\sigma }$ sur ce jeu, puis construire les deux intervalles ;
3. répéter $B$ fois, puis compter la proportion de chacun qui contient $\theta$.

import numpy as np
import scipy.stats as sta
rng = np.random.default_rng(2026)

theta, sigma, N, B = 2.30, 0.10, 10, 100_000
z = sta.norm.ppf(0.975)
t = sta.t.ppf(0.975, N - 1)              # >>> 2.2622  (Gaussian quantile: 1.96)

covered_z, covered_t = 0, 0
for b in range(B):
    x = rng.normal(theta, sigma, N)       # one repetition of the experiment
    mean, sigma_hat = x.mean(), x.std(ddof=1)
    half_width = sigma_hat / np.sqrt(N)
    if abs(mean - theta) <= z * half_width:
        covered_z += 1                    # naive interval, Gaussian quantile
    if abs(mean - theta) <= t * half_width:
        covered_t += 1                    # Student interval

coverage_z = covered_z / B                # >>> 0.9179  (predicted: 0.9184)
coverage_t = covered_t / B                # >>> 0.9493  (predicted: 0.95)

L'intervalle naïf n'est pas seulement imprécis, il est trop étroit : il annonce $95\mathrm{\%}$ et en tient $92\mathrm{\%}$. Une couverture inférieure au niveau affiché est le défaut le plus grave d'un intervalle, puisqu'il conduit à surestimer la confiance accordée au résultat.

Exemple : mesure d'une tension

Reprenons les dix mesures du §4, en supposant cette fois la notice du voltmètre indisponible. Les mêmes données donnent $\overline{x}=\mathrm{2,288}$ V et $\hat{\sigma }=\mathrm{0,087}5$ V, d'où une demi-largeur $\mathrm{2,262}\times \mathrm{0,087}5/\sqrt{10}=\mathrm{0,063}$ V et l'intervalle de confiance à $95\mathrm{\%}$

$$[\mathrm{2,225};\mathrm{2,351}]\text{V}.$$

Il est presque identique à celui du §4, $[\mathrm{2,226};\mathrm{2,350}]$ V, car deux effets se compensent : $\hat{\sigma }=\mathrm{0,087}5$ V sous-estime ici $\sigma =0,10$ V, et le quantile de Student, plus grand, rattrape cet écart. Cette coïncidence ne vaut que pour ce tirage : sur l'ensemble des répétitions, seule la version de Student tient son niveau.