Construisez un estimateur par les moments ou le maximum de vraisemblance

Au chapitre précédent, les estimateurs tombaient un peu du ciel : « prenons deux fois la moyenne », « prenons le maximum ». Pour un modèle plus riche, avec plusieurs paramètres ou une loi peu familière, il faut une méthode qui fabrique l’estimateur à coup sûr.

Ce chapitre en présente deux : la méthode des moments, simple et intuitive, et le maximum de vraisemblance, qui est la méthode de référence et le moteur caché de très nombreux modèles de machine learning. Notre fil rouge : des temps d’attente à un arrêt de bus.

Le fil rouge : attendre le bus

Sur une ligne sans horaires fixes (un bus toutes les quelques minutes), vous arrivez à l’arrêt à des instants quelconques et notez votre temps d’attente, en minutes, sur 8 jours :

3,1	7,4	1,2	5,8	0,6	9,3	2,5	4,1

La somme vaut 34, donc $\bar x_8 = 4{,}25$ minutes. Un modèle classique pour ce genre d’attente est la loi exponentielle de paramètre $\theta > 0$ (le taux de passage, en bus par minute), de densité

$$f(x; \theta) = \theta \, e^{-\theta x} \quad \text{pour } x \ge 0$$

avec $E[X] = 1/\theta$ et $V(X) = 1/\theta^2$. Le paramètre $\theta$ est inconnu : c’est lui qu’il faut estimer.

La méthode des moments

L’idée tient en trois temps :

1. écrire les moments théoriques ($E[X]$, $E[X^2]$, …) en fonction des paramètres ;
2. inverser ce système pour exprimer les paramètres en fonction des moments ;
3. remplacer chaque moment théorique par le moment empirique correspondant ($\frac{1}{n}\sum X_i$, $\frac{1}{n}\sum X_i^2$, …).

On utilise autant de moments qu’il y a de paramètres. Pourquoi ça marche ? Parce que les moments empiriques sont des estimateurs consistants des moments théoriques (loi des grands nombres), et qu’une fonction continue d’estimateurs consistants reste consistante.

Exemple

Exponentielle. Un seul paramètre, un seul moment : $E[X] = 1/\theta$, donc $\theta = 1/E[X]$, d’où $\hat\theta = 1/\bar X_n$. Pour nos attentes, $\hat\theta = 1/4{,}25 \approx 0{,}235$ bus par minute.

Bernoulli. $E[X] = \theta$, donc $\hat\theta = \bar X_n$ : la proportion observée.

Normale $\mathcal{N}(\mu, \sigma^2)$. Deux paramètres, deux moments : $E[X] = \mu$ et $E[X^2] = \sigma^2 + \mu^2$. En inversant, $\mu = E[X]$ et $\sigma^2 = E[X^2] - E[X]^2$. D’où $\hat\mu = \bar X_n$ et $\hat\sigma^2 = \frac{1}{n}\sum X_i^2 - \bar X_n^2 = S_n^2$, la variance non corrigée.

Uniforme sur $[0, \theta]$. $E[X] = \theta/2$, donc $\hat\theta = 2\bar X_n$ : c’est le premier estimateur des vélos du chapitre 2.

La méthode est simple et fonctionne même quand la vraisemblance est difficile à écrire. Mais elle a deux défauts : le choix des moments n’est pas unique, et les estimateurs obtenus sont souvent moins précis que ceux du maximum de vraisemblance. L’exemple uniforme le montre bien : $2\bar X_n$ a une EQM plus de deux fois supérieure à celle du maximum corrigé.

La vraisemblance : quelle valeur rend les données les plus plausibles ?

Changeons de point de vue. Les données sont fixées ; faisons varier $\theta$ et demandons-nous : pour quelle valeur de $\theta$ ces données étaient-elles les plus probables ?

La vraisemblance d’une observation $x$ est la fonction de $\theta$ :

• $L_1(\theta; x) = P_\theta(X = x)$ si la loi est discrète ;
• $L_1(\theta; x) = f(x; \theta)$ si la loi a une densité.

Pour un $n$-échantillon i.i.d., l’indépendance transforme la probabilité jointe en produit :

$$L_n(\theta; x_1, \dots, x_n) = \prod_{i=1}^{n} L_1(\theta; x_i)$$

L’estimateur du maximum de vraisemblance (EMV, en anglais maximum likelihood estimator) est la valeur qui rend $L_n$ maximale :

$$\hat\theta = \arg\max_{\theta} L_n(\theta; X_1, \dots, X_n)$$

En pratique, on maximise la log-vraisemblance $\ln L_n$, qui transforme le produit en somme sans changer la position du maximum (le logarithme est croissant). On annule sa dérivée, puis on vérifie qu’il s’agit bien d’un maximum (dérivée seconde négative, ou étude des variations).

Les données (points ambre) sont fixées, on fait varier θ. La vraisemblance multiplie les hauteurs de la densité aux points observés (traits verts) : θ = 0,235 donne le produit le plus grand.

Exemple

Exponentielle. La vraisemblance du 8-échantillon est

$$L_n(\theta) = \prod_{i=1}^{n} \theta e^{-\theta x_i} = \theta^n e^{-\theta \sum x_i}$$

d’où $\ln L_n(\theta) = n \ln \theta - \theta \sum x_i$. La dérivée vaut $\frac{n}{\theta} - \sum x_i$, qui s’annule en $\hat\theta = n / \sum x_i = 1/\bar x_n$. La dérivée seconde, $-n/\theta^2$, est négative : c’est bien un maximum.

Avec nos données, $\ln L_8(\theta) = 8 \ln \theta - 34\,\theta$. Quelques valeurs :

$\theta$	0,15	0,20	0,235	0,30	0,35
$\ln L_8(\theta)$	$-20{,}28$	$-19{,}68$	$-19{,}58$	$-19{,}83$	$-20{,}30$

Le maximum est atteint en $\hat\theta = 1/4{,}25 \approx 0{,}235$. Ici, l’EMV coïncide avec l’estimateur des moments.

La log-vraisemblance est une fonction de θ : l'EMV est l'abscisse de son sommet, là où la dérivée s'annule.

D’autres exemples classiques

• Bernoulli. Sur 200 visiteurs d’une page, 37 cliquent sur un bouton. $L_n(\theta) = \theta^{\sum x_i} (1 - \theta)^{n - \sum x_i}$, dont la log-vraisemblance s’annule en $\hat\theta = \bar x_n = 37/200 = 0{,}185$.
• Normale. $\ln L_n(\mu, \sigma^2) = -\frac{n}{2}\ln(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum (x_i - \mu)^2$. En annulant les deux dérivées partielles, on obtient $\hat\mu = \bar x_n$ et $\hat\sigma^2 = \frac{1}{n}\sum (x_i - \bar x_n)^2 = s_n^2$. L’EMV de la variance est donc biaisé.
• Uniforme sur $[0, \theta]$. La densité vaut $1/\theta$ sur $[0, \theta]$ et 0 ailleurs, donc $L_n(\theta) = \theta^{-n}$ si $\theta \ge \max x_i$, et 0 sinon. Cette fonction décroît sur $[\max x_i, +\infty[$ : le maximum est au bord, en $\hat\theta = X_{(n)}$. Pour les vélos, 305. Ici, on ne trouve pas l’EMV en annulant une dérivée.

La vraisemblance saute de 0 à son maximum en θ = 305, puis décroît : le maximum est au bord, là où la dérivée ne s'annule pas.

Piège

Annuler la dérivée ne marche que si le maximum est intérieur et la vraisemblance dérivable. Quand le support de la loi dépend du paramètre (loi uniforme), le maximum est souvent au bord du domaine : il faut étudier la fonction, pas seulement sa dérivée.

Quand il n’y a pas de formule : maximiser numériquement

Pour beaucoup de lois (loi gamma, mélanges, régression logistique), l’équation « dérivée nulle » n’a pas de solution explicite. On maximise alors la log-vraisemblance numériquement.

import numpy as np
from scipy import optimize, stats

attentes = np.array([3.1, 7.4, 1.2, 5.8, 0.6, 9.3, 2.5, 4.1])

# Opposé de la log-vraisemblance exponentielle (on minimise)
def nll(theta):
    return -np.sum(stats.expon.logpdf(attentes, scale=1 / theta))

res = optimize.minimize_scalar(nll, bounds=(1e-6, 5), method="bounded")
print(res.x, 1 / attentes.mean())       # les deux valent environ 0.235

attentes <- c(3.1, 7.4, 1.2, 5.8, 0.6, 9.3, 2.5, 4.1)
loglik <- function(theta) sum(dexp(attentes, rate = theta, log = TRUE))
optimize(loglik, interval = c(0.001, 5), maximum = TRUE)$maximum

Les propriétés de l’EMV

L’invariance

Si $\hat\theta$ est l’EMV de $\theta$, alors pour une bijection $g$, l’EMV de $g(\theta)$ est $g(\hat\theta)$. C’est très pratique.

Exemple

Quelle est la probabilité d’attendre plus de 5 minutes ? Pour la loi exponentielle, $P(X > 5) = e^{-5\theta}$. Par invariance, son EMV est $e^{-5\hat\theta} = e^{-5/4{,}25} \approx 0{,}31$. Environ une attente sur trois dépasserait 5 minutes.

Le comportement pour de grands échantillons

Un EMV n’est pas forcément unique, ni sans biais, ni de variance minimale à $n$ fixé. Mais quand la loi est suffisamment régulière (en particulier, son support ne dépend pas de $\theta$), il a d’excellentes propriétés asymptotiques :

• il est consistant ;
• il est asymptotiquement sans biais et normal :

$$\sqrt{n}\,(\hat\theta_n - \theta) \xrightarrow{\mathcal{L}} \mathcal{N}\left(0,\ \frac{1}{I_1(\theta)}\right)$$

où $I_1(\theta)$ est l’information de Fisher d’une observation :

$$I_1(\theta) = V\left( \frac{\partial}{\partial\theta} \ln L_1(\theta; X) \right) = E\left[ \left( \frac{\partial}{\partial\theta} \ln L_1(\theta; X) \right)^2 \right]$$

La dérivée de la log-vraisemblance s’appelle le score. L’information de Fisher mesure à quel point la vraisemblance est « pointue » autour du vrai paramètre : plus elle est grande, plus une observation renseigne sur $\theta$, et plus l’EMV est précis.

Avec dix fois plus de données, la courbe est dix fois plus courbée au sommet : l'information totale n I₁(θ) augmente et l'intervalle des valeurs plausibles de θ se resserre.

Exemple

Pour l’exponentielle, $\ln L_1(\theta; X) = \ln\theta - \theta X$, donc le score vaut $\frac{1}{\theta} - X$. Sa variance est $V(X) = 1/\theta^2$, d’où $I_1(\theta) = 1/\theta^2$ et une variance asymptotique de $\theta^2$. C’est exactement le résultat obtenu par la delta-méthode au chapitre 2.

Pour nos 8 attentes, l’erreur-type approchée de $\hat\theta$ vaut $\hat\theta/\sqrt{8} \approx 0{,}235/2{,}83 \approx 0{,}083$. Avec si peu de données, l’estimation reste très imprécise.

Un complément : la borne de Cramér-Rao

Ce résultat classique, absent de nos supports de cours, éclaire la place de l’EMV. Pour tout estimateur sans biais $T_n$ d’un modèle régulier,

$$V(T_n) \ge \frac{1}{n \, I_1(\theta)}$$

Aucun estimateur sans biais ne peut descendre sous cette variance. Un estimateur qui l’atteint est dit efficace. L’EMV l’atteint asymptotiquement : pour de grands échantillons, on ne peut pas faire mieux. Par exemple, pour la loi de Poisson, $I_1(\theta) = 1/\theta$ et $V(\bar X_n) = \theta/n$ : la moyenne empirique est efficace.

À retenir

La méthode des moments égale moments théoriques et empiriques. Le maximum de vraisemblance choisit le paramètre qui rend les données les plus probables : on maximise $\ln L_n$. L’EMV est invariant par reparamétrage et, pour un modèle régulier, asymptotiquement normal de variance $1/(n I_1(\theta))$, la meilleure possible.

	Méthode des moments	Maximum de vraisemblance
Principe	égaler moments théoriques et empiriques	maximiser la probabilité des données
Calcul	souvent explicite et simple	parfois explicite, souvent numérique
Précision	souvent moindre	asymptotiquement optimale (modèle régulier)
Usage typique	point de départ, modèles complexes	méthode de référence, logiciels

En entretien

« Le maximum de vraisemblance choisit les paramètres qui rendent les données observées les plus probables. Minimiser l’entropie croisée d’une régression logistique ou les moindres carrés d’une régression gaussienne, c’est exactement maximiser une vraisemblance. »

En résumé

• La méthode des moments exprime les paramètres en fonction des moments, puis remplace ceux-ci par leurs versions empiriques.
• La vraisemblance $L_n(\theta) = \prod L_1(\theta; x_i)$ mesure la plausibilité des données pour chaque valeur de $\theta$ ; l’EMV la maximise, en général via la dérivée de $\ln L_n$.
• Exemples à connaître : Bernoulli ($\bar x$), exponentielle ($1/\bar x$), normale ($\bar x$ et $s_n^2$, biaisée), uniforme ($\max x_i$, maximum au bord).
• L’EMV est invariant : l’EMV de $g(\theta)$ est $g(\hat\theta)$.
• Pour un modèle régulier, l’EMV est consistant et asymptotiquement normal, de variance asymptotique $1/I_1(\theta)$ (information de Fisher) ; il atteint asymptotiquement la borne de Cramér-Rao.

Exercices

Exercice · Méthode des moments pour une loi de Poisson

facile

Le nombre de tickets de support reçus par heure suit une loi de Poisson de paramètre $\theta$. Sur 10 heures, on compte 2, 0, 3, 1, 4, 2, 2, 1, 3, 2 tickets. Estimez $\theta$ par la méthode des moments, puis la probabilité qu’une heure se passe sans ticket.

Voir la solution

Pour une loi de Poisson, $E[X] = \theta$, donc $\hat\theta = \bar X_n$. La somme vaut 20, d’où $\hat\theta = 2$ tickets par heure.

$P(X = 0) = e^{-\theta}$. En remplaçant $\theta$ par son estimation (c’est aussi l’EMV par invariance, voir l’exercice suivant), on obtient $e^{-2} \approx 0{,}135$ : environ une heure sur sept sans ticket.

En code

Python

import numpy as np
from scipy import stats

tickets = np.array([2, 0, 3, 1, 4, 2, 2, 1, 3, 2])
theta_hat = tickets.mean()                         # E[X] = theta
print(theta_hat)                                   # 2.0
print(np.exp(-theta_hat), stats.poisson.pmf(0, theta_hat))   # 0.135  0.135

R

tickets <- c(2, 0, 3, 1, 4, 2, 2, 1, 3, 2)
theta_hat <- mean(tickets)                  # E[X] = theta
theta_hat                                   # 2
exp(-theta_hat); dpois(0, theta_hat)        # 0.135  0.135

Exercice · Maximum de vraisemblance pour une loi de Poisson

moyen

Montrez que l’EMV du paramètre d’une loi de Poisson est $\bar X_n$, puis calculez l’information de Fisher $I_1(\theta)$.

Voir la solution

$P(X = x) = e^{-\theta}\theta^x / x!$, donc

$$\ln L_n(\theta) = -n\theta + \left(\sum x_i\right)\ln\theta - \sum \ln(x_i!)$$

La dérivée $-n + \frac{\sum x_i}{\theta}$ s’annule en $\hat\theta = \bar x_n$ ; la dérivée seconde $-\sum x_i/\theta^2$ est négative (dès qu’une observation est non nulle), c’est donc un maximum.

Pour une observation, le score vaut $-1 + X/\theta$, de variance $V(X)/\theta^2 = \theta/\theta^2 = 1/\theta$. Donc $I_1(\theta) = 1/\theta$ : la variance asymptotique de l’EMV est $\theta$, et $V(\bar X_n) = \theta/n$ atteint exactement la borne de Cramér-Rao.

En code

Python

import numpy as np
from scipy import optimize, stats

tickets = np.array([2, 0, 3, 1, 4, 2, 2, 1, 3, 2])

# L'EMV numérique retombe sur la moyenne
nll = lambda theta: -np.sum(stats.poisson.logpmf(tickets, theta))
res = optimize.minimize_scalar(nll, bounds=(0.01, 10), method="bounded")
print(res.x, tickets.mean())                       # 2.000  2.0

# Information de Fisher : variance du score -1 + X / theta
rng = np.random.default_rng(0)
theta = 2.0
X = rng.poisson(theta, size=200_000)
score = -1 + X / theta
print(score.var(), 1 / theta)                      # 0.500  0.5

R

tickets <- c(2, 0, 3, 1, 4, 2, 2, 1, 3, 2)

# L'EMV numérique retombe sur la moyenne
loglik <- function(theta) sum(dpois(tickets, theta, log = TRUE))
optimize(loglik, interval = c(0.01, 10), maximum = TRUE)$maximum   # 2.000

# Information de Fisher : variance du score -1 + X / theta
set.seed(0)
theta <- 2
X <- rpois(200000, theta)
score <- -1 + X / theta
var(score); 1 / theta                        # 0.500  0.5

Exercice · Une loi à densité puissance

difficile

Une proportion $X \in [0, 1]$ (le taux de remplissage d’un entrepôt) a pour densité $f(x; \theta) = \theta x^{\theta - 1}$ sur $[0, 1]$, avec $\theta > 0$. Déterminez l’estimateur des moments et l’EMV de $\theta$. Sont-ils égaux ?

Voir la solution

Moments. $E[X] = \int_0^1 \theta x^\theta \, dx = \frac{\theta}{\theta + 1}$. En posant $m = E[X]$, on inverse : $\theta = \frac{m}{1 - m}$. D’où $\hat\theta_{mom} = \frac{\bar X_n}{1 - \bar X_n}$.

Vraisemblance. $\ln L_n(\theta) = n\ln\theta + (\theta - 1)\sum \ln x_i$. La dérivée $\frac{n}{\theta} + \sum \ln x_i$ s’annule en

$$\hat\theta_{MV} = -\frac{n}{\sum_{i=1}^{n} \ln X_i}$$

(positif, car les $\ln X_i$ sont négatifs). La dérivée seconde $-n/\theta^2$ est négative : c’est un maximum.

Les deux estimateurs sont différents : l’un utilise la moyenne des $X_i$, l’autre la moyenne des $\ln X_i$. Ils sont tous deux consistants, mais l’EMV est asymptotiquement le plus précis.

En code

Python

import numpy as np

def theta_mom(x):                     # moyenne / (1 - moyenne)
    return x.mean(axis=-1) / (1 - x.mean(axis=-1))

def theta_mv(x):                      # -n / somme des log
    return -x.shape[-1] / np.log(x).sum(axis=-1)

# Simulation par inversion : F(x) = x^theta, donc X = U^(1 / theta)
rng = np.random.default_rng(0)
theta, n = 0.5, 50
X = rng.uniform(size=(20_000, n)) ** (1 / theta)
for nom, f in [("moments", theta_mom), ("vraisemblance", theta_mv)]:
    T = f(X)
    print(nom, round(T.mean(), 3), round(np.mean((T - theta)**2), 4))
# moments       0.507  0.0094
# vraisemblance 0.510  0.0055   (estimateurs différents, EMV plus précis)

R

theta_mom <- function(x) mean(x) / (1 - mean(x))   # moyenne / (1 - moyenne)
theta_mv  <- function(x) -length(x) / sum(log(x))  # -n / somme des log

# Simulation par inversion : F(x) = x^theta, donc X = U^(1 / theta)
set.seed(0)
theta <- 0.5; n <- 50
echantillons <- replicate(20000, runif(n)^(1 / theta))   # une colonne par échantillon
T_mom <- apply(echantillons, 2, theta_mom)
T_mv  <- apply(echantillons, 2, theta_mv)
c(mean(T_mom), mean((T_mom - theta)^2))   # moyenne et EQM
c(mean(T_mv),  mean((T_mv - theta)^2))
## moments       0.506  0.0093
## vraisemblance 0.510  0.0056   (estimateurs différents, EMV plus précis)