Estimez un paramètre avec un intervalle de confiance

Un fabricant de drones teste 10 batteries de son nouveau modèle et mesure leur autonomie : 26,95 minutes en moyenne. Peut-il écrire « autonomie : 27 minutes » sur la boîte ? Avec 10 autres batteries, il aurait obtenu une autre moyenne. Ce qu’il lui faut, c’est une fourchette accompagnée d’un niveau de confiance : « entre 26,3 et 27,6 minutes, avec 95 % de confiance ».

Jusqu’ici, nous partions d’une loi connue pour calculer des probabilités. Ce chapitre inverse la démarche : on observe des données, et l’on cherche à remonter à la loi qui les a produites. C’est le passage des probabilités à la statistique. Le fil rouge sera l’autonomie de ces batteries. Le cours de Statistiques approfondit ensuite l’estimation (maximum de vraisemblance, propriétés asymptotiques).

Estimateur et estimation

Les 10 autonomies mesurées (en minutes) sont :

26,1	27,4	25,8	28,2	26,9	27,7	25,5	26,6	28,0	27,3

On les modélise comme les réalisations $x_1, \dots, x_n$ de variables $X_1, \dots, X_n$ i.i.d., de loi inconnue, d’espérance $\mu$ (l’autonomie moyenne de toutes les batteries produites) et d’écart-type $\sigma$. Le but est d’estimer le paramètre $\mu$.

Un estimateur $\hat{\theta}$ d’un paramètre $\theta$ est une variable aléatoire calculée à partir de l’échantillon, comme $\bar{X}_n$. Une estimation est la valeur qu’il prend sur les données observées, comme $\bar{x} = 26{,}95$. La distinction compte : on juge un estimateur sur son comportement sur tous les échantillons possibles, pas sur un seul.

Deux critères dominent :

• le biais $E[\hat{\theta}] - \theta$ : un estimateur est sans biais s’il tombe juste en moyenne ;
• la variance $\mathrm{Var}(\hat{\theta})$ : un bon estimateur fluctue peu d’un échantillon à l’autre.

On les combine dans l’erreur quadratique moyenne, $E[(\hat{\theta} - \theta)^2] = \text{biais}^2 + \text{variance}$. D’après le chapitre 7, $\bar{X}_n$ est un estimateur sans biais de $\mu$, de variance $\sigma^2/n$ : il s’améliore quand $n$ grandit.

Pourquoi diviser par n − 1 ?

Pour estimer la variance $\sigma^2$, l’idée naturelle est la moyenne des carrés des écarts, $S^2 = \frac{1}{n} \sum_i (X_i - \bar{X}_n)^2$. Or on démontre que :

$$E[S^2] = \frac{n - 1}{n} \, \sigma^2$$

$S^2$ sous-estime $\sigma^2$ en moyenne. La raison est intuitive : on mesure les écarts à $\bar{X}_n$, qui est par construction le point le plus proche de l’ensemble des données, plutôt qu’à la vraie moyenne $\mu$ ; les écarts sont donc un peu trop petits. On corrige en divisant par $n - 1$ :

$$S'^2 = \frac{1}{n - 1} \sum_{i=1}^{n} (X_i - \bar{X}_n)^2, \qquad E[S'^2] = \sigma^2$$

C’est la variance empirique corrigée, celle que calculent var et sd en R. Pour nos batteries, $s'^2 \approx 0{,}869$ et $s' \approx 0{,}932$ minute (la version biaisée donnerait $0{,}782$).

Sur des échantillons de 5 valeurs, diviser par n donne des estimations en moyenne trop petites (3,2 au lieu de 4). Diviser par n − 1 décale tout l'histogramme juste ce qu'il faut.

Attention

En Python, numpy.var et numpy.std divisent par $n$ par défaut ; il faut passer ddof=1 pour obtenir la version corrigée. pandas, lui, divise par $n - 1$ par défaut. Sur de petits échantillons, l’écart n’est pas négligeable.

Qu’est-ce qu’un intervalle de confiance ?

Une estimation seule (26,95) ne dit rien de sa précision. Un intervalle de confiance (IC) de niveau $1 - \alpha$ est un intervalle aléatoire $[L\,;\,U]$, calculé à partir de l’échantillon, tel que :

$$P(L \leq \theta \leq U) = 1 - \alpha$$

Le niveau le plus courant est 95 % ($\alpha = 0{,}05$). La recette est toujours la même : trouver une quantité, appelée pivot, qui dépend de $\theta$ et dont la loi est connue et ne dépend pas de $\theta$ ; encadrer ce pivot entre deux quantiles ; puis « isoler » $\theta$.

On note $u_{1-\alpha/2}$ le quantile d’ordre $1 - \alpha/2$ de $\mathcal{N}(0, 1)$, c’est-à-dire la valeur telle que $P(Z \leq u_{1-\alpha/2}) = 1 - \alpha/2$. Pour $\alpha = 0{,}05$, $u_{0{,}975} \approx 1{,}96$ : on laisse 2,5 % de probabilité de chaque côté.

Piège

Une fois calculé, l’intervalle $[26{,}3\,;\,27{,}6]$ contient $\mu$ ou ne le contient pas : il n’y a plus rien d’aléatoire. Dire « $\mu$ a 95 % de chances d’être dans cet intervalle » est donc un abus de langage. Le « 95 % » porte sur la procédure : si l’on répétait l’expérience 100 fois, environ 95 des 100 intervalles obtenus contiendraient $\mu$.

Intervalle pour une moyenne

Quand l’écart-type est connu

Supposons d’abord, pour simplifier, que le fabricant connaisse l’écart-type de l’autonomie grâce aux modèles précédents : $\sigma = 1$ minute. Si les $X_i$ sont normales (ou si $n$ est grand, grâce au TCL), le pivot

$$Z = \frac{\bar{X}_n - \mu}{\sigma / \sqrt{n}}$$

suit $\mathcal{N}(0, 1)$. On a donc $P(-1{,}96 \leq Z \leq 1{,}96) = 0{,}95$, et en isolant $\mu$ :

$$\mu \in \left[\bar{x} - u_{1-\alpha/2} \frac{\sigma}{\sqrt{n}}\,;\, \bar{x} + u_{1-\alpha/2} \frac{\sigma}{\sqrt{n}}\right]$$

Pour les batteries : $26{,}95 \pm 1{,}96 \times 1/\sqrt{10} = 26{,}95 \pm 0{,}62$, soit $[26{,}33\,;\,27{,}57]$.

La demi-largeur $e = u_{1-\alpha/2} \, \sigma/\sqrt{n}$ mesure la précision. Elle diminue en $1/\sqrt{n}$, ce qui permet de calculer la taille d’échantillon nécessaire pour une précision donnée : $n \geq (u_{1-\alpha/2} \, \sigma / e)^2$. Pour une précision de $\pm 0{,}3$ minute, il faudrait $n \geq (1{,}96/0{,}3)^2 \approx 42{,}7$, donc 43 batteries.

Quand l’écart-type est inconnu : la loi de Student

En réalité, $\sigma$ est presque toujours inconnu. On le remplace par son estimation $s'$, mais cela ajoute de l’incertitude : le pivot $\frac{\bar{X}_n - \mu}{S'/\sqrt{n}}$ ne suit plus une normale. Pour des données normales, il suit une loi de Student à $n - 1$ degrés de liberté, notée $t_{n-1}$.

Deux lois nouvelles apparaissent ici, toutes deux construites à partir de la normale :

• la loi du khi-deux à $k$ degrés de liberté, $\chi^2_k$, est la loi de $Z_1^2 + \dots + Z_k^2$ où les $Z_i$ sont des $\mathcal{N}(0, 1)$ indépendantes ;
• la loi de Student à $k$ degrés de liberté, $t_k$, est la loi de $Z / \sqrt{U/k}$, avec $Z \sim \mathcal{N}(0, 1)$ et $U \sim \chi^2_k$ indépendantes.

La loi de Student ressemble à une normale centrée réduite, mais avec des queues plus épaisses : c’est le prix de l’incertitude sur $\sigma$. Elle s’en rapproche quand $k$ augmente : le quantile d’ordre 0,975 vaut 2,78 pour $k = 4$, 2,26 pour $k = 9$, 2,05 pour $k = 29$, et tend vers 1,96. L’intervalle devient :

Les lois de Student ont des queues plus épaisses que la normale : pour laisser 2,5 % de chaque côté, il faut aller plus loin, et d'autant plus que k est petit.

$$\mu \in \left[\bar{x} - t_{n-1}(1 - \alpha/2) \frac{s'}{\sqrt{n}}\,;\, \bar{x} + t_{n-1}(1 - \alpha/2) \frac{s'}{\sqrt{n}}\right]$$

Exemple

Pour les batteries, $n = 10$, $\bar{x} = 26{,}95$, $s' \approx 0{,}932$ et $t_9(0{,}975) \approx 2{,}262$. La demi-largeur vaut $2{,}262 \times 0{,}932/\sqrt{10} \approx 0{,}67$, d’où l’intervalle $[26{,}28\,;\,27{,}62]$. Il est un peu plus large que celui obtenu avec $\sigma$ connu : ne pas connaître $\sigma$ se paie en précision.

Les deux intervalles sont centrés sur x̄ = 26,95 ; celui de Student est un peu plus large, car il paie l'estimation de σ. Tous deux sont bien plus étroits que la dispersion des batteries elles-mêmes.

En R, t.test fait tout le calcul :

autonomie <- c(26.1, 27.4, 25.8, 28.2, 26.9, 27.7, 25.5, 26.6, 28.0, 27.3)
n <- length(autonomie)

mean(autonomie); sd(autonomie)                          # 26.95 et 0.932 (sd divise par n - 1)
mean(autonomie) + qt(c(0.025, 0.975), df = n - 1) * sd(autonomie) / sqrt(n)   # 26.28 27.62
t.test(autonomie, conf.level = 0.95)$conf.int           # même résultat

Intervalle pour une variance

Le fabricant veut aussi savoir si l’autonomie est régulière d’une batterie à l’autre : il faut un intervalle pour $\sigma^2$. Pour des données normales, le pivot $(n-1) S'^2 / \sigma^2$ suit la loi $\chi^2_{n-1}$. En l’encadrant entre ses quantiles d’ordre $\alpha/2$ et $1 - \alpha/2$, puis en isolant $\sigma^2$ :

$$\sigma^2 \in \left[\frac{(n-1)\,s'^2}{\chi^2_{n-1}(1 - \alpha/2)}\,;\, \frac{(n-1)\,s'^2}{\chi^2_{n-1}(\alpha/2)}\right]$$

Attention à l’inversion : le grand quantile donne la borne inférieure. Pour les batteries, $(n-1)s'^2 = 9 \times 0{,}869 = 7{,}825$, et les quantiles de $\chi^2_9$ valent 2,70 et 19,02. D’où $\sigma^2 \in [0{,}41\,;\,2{,}90]$, soit $\sigma \in [0{,}64\,;\,1{,}70]$ minute. L’intervalle est large et asymétrique : 10 mesures renseignent mal sur la dispersion.

La loi du khi-deux n'est pas symétrique, donc l'intervalle pour σ² ne l'est pas non plus : l'estimation 0,869 est bien plus proche de la borne basse que de la borne haute.

Attention

Les intervalles pour la moyenne sont robustes : grâce au TCL, ils restent approximativement valables sans normalité dès que $n$ est assez grand. Celui pour la variance, non : il repose sur la normalité des données, même pour de grands échantillons. Vérifiez-la (histogramme, droite de Henry) avant de l’utiliser.

Intervalle pour une proportion

Enfin, le fabricant sonde 400 clients : 92 d’entre eux déclarent utiliser leur drone plus de 25 minutes d’affilée, soit une proportion observée $\hat{p} = 0{,}23$. Le nombre de « oui » suit $\mathcal{B}(n, p)$, que le TCL approche par une normale (chapitre 8) dès que $n \geq 30$, $np \geq 5$ et $n(1-p) \geq 5$. Le pivot $\frac{\hat{p} - p}{\sqrt{p(1-p)/n}}$ suit donc approximativement $\mathcal{N}(0, 1)$, et :

$$p \in \left[\hat{p} - u_{1-\alpha/2} \sqrt{\frac{p(1-p)}{n}}\,;\, \hat{p} + u_{1-\alpha/2} \sqrt{\frac{p(1-p)}{n}}\right]$$

Problème : $p$, l’inconnue, apparaît dans les bornes. Deux solutions :

• remplacer $p$ par $\hat{p}$ dans la racine (intervalle de Wald, très répandu mais approximatif, surtout quand $\hat{p}$ est proche de 0 ou de 1) : $0{,}23 \pm 1{,}96\sqrt{0{,}23 \times 0{,}77 / 400} = 0{,}23 \pm 0{,}041$, soit $[0{,}189\,;\,0{,}271]$ ;
• utiliser la borne $p(1 - p) \leq 1/4$, valable pour tout $p$, qui donne un intervalle prudent : $\hat{p} \pm u_{1-\alpha/2} / (2\sqrt{n}) = 0{,}23 \pm 0{,}049$, soit $[0{,}181\,;\,0{,}279]$.

La version prudente est un peu plus large, mais son niveau de confiance est garanti quel que soit $p$. Elle donne aussi la formule de taille des sondages : pour une précision de $\pm 2$ points à 95 %, il faut $n \geq (1{,}96/(2 \times 0{,}02))^2 \approx 2\,401$ personnes.

Paramètre	Hypothèse	Pivot et loi	Intervalle à $1 - \alpha$
$\mu$	$\sigma$ connu	$\frac{\bar{X} - \mu}{\sigma/\sqrt{n}} \sim \mathcal{N}(0,1)$	$\bar{x} \pm u_{1-\alpha/2}\,\sigma/\sqrt{n}$
$\mu$	$\sigma$ inconnu	$\frac{\bar{X} - \mu}{S'/\sqrt{n}} \sim t_{n-1}$	$\bar{x} \pm t_{n-1}(1-\alpha/2)\,s'/\sqrt{n}$
$\sigma^2$	données normales	$\frac{(n-1)S'^2}{\sigma^2} \sim \chi^2_{n-1}$	voir la formule ci-dessus
$p$	$n$ grand	$\frac{\hat{p} - p}{\sqrt{p(1-p)/n}} \approx \mathcal{N}(0,1)$	$\hat{p} \pm u_{1-\alpha/2}/(2\sqrt{n})$ (prudent)

La simulation suivante illustre la bonne interprétation du niveau de confiance : on tire 100 échantillons d’une loi dont on connaît la moyenne, et l’on compte combien d’intervalles la contiennent.

set.seed(15)
contient <- replicate(100, {
  x <- rnorm(10, mean = 27, sd = 1)       # un échantillon de 10 batteries, vraie moyenne 27
  ic <- t.test(x)$conf.int
  ic[1] <= 27 && 27 <= ic[2]
})
sum(contient)                             # environ 95 intervalles sur 100 contiennent 27

Chaque échantillon donne un intervalle différent. La procédure capture la vraie moyenne 95 fois sur 100 ; les 5 intervalles roses sont des échantillons malchanceux, et rien ne permet de les reconnaître sans connaître μ.

En entretien

« Un intervalle de confiance à 95 % est produit par une procédure qui capture le vrai paramètre dans 95 % des échantillons. Quand l’écart-type est inconnu, on utilise la loi de Student, dont les queues plus épaisses élargissent l’intervalle ; sa largeur décroît en $1/\sqrt{n}$. »

Pour aller plus loin

Le cours de Statistiques reprend ces notions en détail : qualité d’un estimateur, méthode des moments et maximum de vraisemblance et intervalles de confiance.

En résumé

• Un estimateur est une variable aléatoire fonction de l’échantillon ; on le juge par son biais et sa variance.
• $\bar{X}_n$ estime $\mu$ sans biais ; $S'^2$ (division par $n - 1$) estime $\sigma^2$ sans biais, alors que la division par $n$ sous-estime.
• Un intervalle de confiance se construit à partir d’un pivot de loi connue ; son niveau porte sur la procédure, pas sur l’intervalle calculé.
• Moyenne : loi normale si $\sigma$ est connu, loi de Student $t_{n-1}$ sinon ; la demi-largeur décroît en $1/\sqrt{n}$.
• Variance : loi du khi-deux $\chi^2_{n-1}$, intervalle asymétrique et sensible à la non-normalité.
• Proportion : approximation normale, avec la borne prudente $\hat{p} \pm u_{1-\alpha/2}/(2\sqrt{n})$.

Exercices

Exercice · Un intervalle avec écart-type connu

facile

Une machine remplit des paquets avec un écart-type connu de 5 g. Sur 25 paquets, le poids moyen est de 512,3 g. Donnez un intervalle de confiance à 95 % pour le poids moyen, puis dites combien de paquets il faudrait peser pour une précision de $\pm 1$ g.

Voir la solution

Demi-largeur : $1{,}96 \times 5/\sqrt{25} = 1{,}96$ g. L’intervalle est $[510{,}34\,;\,514{,}26]$ g.

Pour une demi-largeur de 1 g : $n \geq (1{,}96 \times 5/1)^2 = 96{,}04$, donc 97 paquets. Diviser la demi-largeur par deux demande quatre fois plus de mesures.

En code

Python

import numpy as np
from scipy import stats

sigma, n, moyenne = 5, 25, 512.3
u = stats.norm.ppf(0.975)                      # 1.96
demi_largeur = u * sigma / np.sqrt(n)
print(demi_largeur)                            # 1.96
print(moyenne - demi_largeur, moyenne + demi_largeur)   # 510.34 514.26

# Taille d'échantillon pour une demi-largeur de 1 g
print((u * sigma / 1) ** 2, np.ceil((u * sigma / 1) ** 2))   # 96.04 97.0

R

sigma <- 5; n <- 25; moyenne <- 512.3
u <- qnorm(0.975)                              # 1.96
demi_largeur <- u * sigma / sqrt(n)
demi_largeur                                   # 1.96
moyenne + c(-1, 1) * demi_largeur              # 510.34 514.26

# Taille d'échantillon pour une demi-largeur de 1 g
(u * sigma / 1)^2                              # 96.04
ceiling((u * sigma / 1)^2)                     # 97

Exercice · Un sondage serré

moyen

Avant un vote d’association, on interroge 1 000 membres tirés au hasard : 540 sont favorables à la proposition.

1. Donnez l’intervalle de Wald et l’intervalle prudent à 95 %.
2. Peut-on affirmer, avec 95 % de confiance, que la proposition est majoritaire ?

Voir la solution

1. $\hat{p} = 0{,}54$. Wald : $0{,}54 \pm 1{,}96\sqrt{0{,}54 \times 0{,}46/1\,000} = 0{,}54 \pm 0{,}031$, soit $[0{,}509\,;\,0{,}571]$. Prudent : $0{,}54 \pm 1{,}96/(2\sqrt{1\,000}) = 0{,}54 \pm 0{,}031$, soit $[0{,}509\,;\,0{,}571]$ aussi. Les deux coïncident presque, car $\hat{p}$ est proche de $1/2$, où $p(1-p)$ atteint justement son maximum $1/4$.
2. Les deux intervalles sont entièrement au-dessus de 0,5 : oui, au niveau 95 %, les données sont compatibles avec une proportion favorable supérieure à 50 %, et incompatibles avec 50 % ou moins. Mais la marge est mince (la borne basse vaut 50,9 %), et le raisonnement suppose un tirage réellement aléatoire et des réponses sincères. Le chapitre 10 formalise ce type de conclusion sous forme de test.

En code

Python

import numpy as np
from scipy import stats

n, oui = 1000, 540
p_chapeau = oui / n
u = stats.norm.ppf(0.975)

wald = u * np.sqrt(p_chapeau * (1 - p_chapeau) / n)
prudent = u / (2 * np.sqrt(n))
print(p_chapeau - wald, p_chapeau + wald)          # 0.5091 0.5709
print(p_chapeau - prudent, p_chapeau + prudent)    # 0.5090 0.5710
# Les deux bornes basses dépassent 0.5 : majorité au niveau 95 %

R

n <- 1000; oui <- 540
p_chapeau <- oui / n
u <- qnorm(0.975)

wald <- u * sqrt(p_chapeau * (1 - p_chapeau) / n)
prudent <- u / (2 * sqrt(n))
p_chapeau + c(-1, 1) * wald       # 0.5091 0.5709
p_chapeau + c(-1, 1) * prudent    # 0.5090 0.5710
# Les deux bornes basses dépassent 0.5 : majorité au niveau 95 %

Exercice · Le biais de la variance empirique

difficile

Soient $X_1, \dots, X_n$ i.i.d. d’espérance $\mu$ et de variance $\sigma^2$.

1. Montrez que $\sum_i (X_i - \bar{X})^2 = \sum_i (X_i - \mu)^2 - n(\bar{X} - \mu)^2$.
2. Déduisez-en $E\left[\sum_i (X_i - \bar{X})^2\right] = (n-1)\sigma^2$, puis le biais de $S^2$.
3. Interprétez le cas extrême d’une seule observation ($n = 1$) : que vaut alors $S^2$ ?

Voir la solution

1. On écrit $X_i - \bar{X} = (X_i - \mu) - (\bar{X} - \mu)$ et on développe : $\sum_i (X_i - \bar{X})^2 = \sum_i (X_i - \mu)^2 - 2(\bar{X} - \mu)\sum_i (X_i - \mu) + n(\bar{X} - \mu)^2$. Or $\sum_i (X_i - \mu) = n(\bar{X} - \mu)$, donc le terme du milieu vaut $-2n(\bar{X} - \mu)^2$, d’où le résultat.
2. En prenant l’espérance : $E[(X_i - \mu)^2] = \sigma^2$ et $E[(\bar{X} - \mu)^2] = \mathrm{Var}(\bar{X}) = \sigma^2/n$. Donc $E\left[\sum_i (X_i - \bar{X})^2\right] = n\sigma^2 - n \times \sigma^2/n = (n - 1)\sigma^2$. On en déduit $E[S^2] = \frac{n-1}{n}\sigma^2$ : le biais vaut $-\sigma^2/n$, et il disparaît quand $n$ grandit.
3. Avec une seule observation, $\bar{X} = X_1$ et $S^2 = 0$ quelle que soit la dispersion réelle : une seule mesure ne dit rien de la variabilité. La version corrigée divise par $n - 1 = 0$ et n’est pas définie, ce qui est plus honnête. On dit qu’on a « consommé un degré de liberté » en estimant $\mu$ par $\bar{X}$ : c’est l’origine des $n - 1$ degrés de liberté de la loi de Student.

En code

Python

import numpy as np

# 200 000 échantillons de taille n = 5, de variance sigma^2 = 4
rng = np.random.default_rng(0)
n, sigma2 = 5, 4
x = rng.normal(10, np.sqrt(sigma2), size=(200_000, n))

s2_biaise = x.var(axis=1, ddof=0)          # division par n
s2_corrige = x.var(axis=1, ddof=1)         # division par n - 1
print(s2_biaise.mean(), (n - 1) / n * sigma2)   # 3.21  3.2
print(s2_corrige.mean())                        # 4.01, sans biais

# Cas n = 1 : la version biaisée vaut toujours 0
print(np.var([7.3]), np.var([7.3], ddof=1))     # 0.0 nan (avec un avertissement)

R

# 200 000 échantillons de taille n = 5, de variance sigma^2 = 4
set.seed(0)
n <- 5; sigma2 <- 4
x <- matrix(rnorm(200000 * n, mean = 10, sd = sqrt(sigma2)), ncol = n)

s2_corrige <- apply(x, 1, var)             # var divise par n - 1
s2_biaise <- s2_corrige * (n - 1) / n      # division par n
c(mean(s2_biaise), (n - 1) / n * sigma2)   # 3.20  3.2
mean(s2_corrige)                           # 4.00, sans biais

# Cas n = 1 : la version corrigée n'est pas définie
var(7.3)                                   # NA