Tests d'adéquation

Les cinq détecteurs de la page précédente choisissaient entre deux lois données. Le chapitre 1 a posé la question qui vient avant, la loi retenue décrit-elle les données, et l'a tranchée à l'œil. Cette page en fait une décision chiffrée : la section 1 pose l'adéquation comme un problème de détection, la section 2 construit le test du ${\chi }^2$, la section 3 en tire la p-value, la section 4 délimite ce qu'un tel test permet de conclure.

1. Problématique

Un serveur enregistre le nombre de requêtes reçues chaque seconde pendant $N=200$ secondes. La loi de Poisson est le modèle naturel de ce comptage, mais rien ne garantit que ses hypothèses tiennent : un trafic par rafales, par exemple, produit des comptages bien plus dispersés que ce qu'elle prévoit. Comment décider, sur les données, si le modèle est acceptable ?

La question se coule dans le cadre de la page précédente, à ceci près que l'hypothèse porte cette fois sur le modèle lui-même :

$$H_0:x[n]\sim p(x;{\theta }_0),$$$$H_1:\text{les données ne suivent pas cette loi}.$$

Les deux hypothèses ne sont pas de même nature. Seule $H_0$ est spécifiée : $H_1$ rassemble toutes les autres lois possibles, sans en désigner aucune. La probabilité de fausse alarme reste donc calculable, rejeter un modèle pourtant vrai, mais la probabilité de détection ne l'est pas, puisqu'elle dépend de l'alternative, qu'il faut choisir pour l'évaluer.

Cette asymétrie a une conséquence qui gouverne toute la page : le test rejette ou ne rejette pas, il ne valide jamais.

2. Test du ${\chi }^2$ d'adéquation

Hypothèses

Trois hypothèses définissent ce cadre :

• Données groupées. Les observations sont réparties en $K$ classes disjointes couvrant tout le support du modèle. Le symbole $K$ compte ici des classes, non des hypothèses comme au §3 des fondations.
• Effectifs suffisants. Chaque classe reçoit un effectif attendu ${\hat{n}}_k\ge 5$, faute de quoi les classes voisines sont fusionnées.
• Paramètres estimés par maximum de vraisemblance. Si $r$ paramètres du modèle sont inconnus, ils sont estimés sur les mêmes données, par MLE.

Expression

L'idée est de comparer, classe par classe, l'effectif observé à celui que le modèle prédit.

Définition, Statistique du ${\chi }^2$ d'adéquation

La statistique de test est la somme des écarts quadratiques normalisés :

$$T(\mathbf{x})=\sum _{k=1}^K\frac{(n_k-{\hat{n}}_k{)}^2}{{\hat{n}}_k},$$

où

• $n_k$ est l'effectif observé de la classe $k$,
• ${\hat{n}}_k=NP(\text{classe }k;\hat{\theta })$ est l'effectif attendu sous le modèle, dont le chapeau rappelle qu'il dépend du paramètre estimé,
• $K$ est le nombre de classes.

Chaque terme est un écart au carré rapporté à l'effectif attendu, la variance d'un comptage étant elle-même de l'ordre de cet effectif : une différence de dix unités ne pèse pas le même poids sur une classe qui en attend douze et sur une classe qui en attend mille. La statistique est nulle pour un accord parfait et grandit dès qu'une classe s'écarte, dans un sens ou dans l'autre. La règle de décision est celle de toute la page précédente, décider $H_1$ si $T(\mathbf{x})>\gamma$.

Propriétés

Proposition, Loi asymptotique sous $H_0$

Sous $H_0$, quand $N\to \mathrm{\infty }$, la statistique $T$ suit une loi du ${\chi }^2$ à

$$\nu =K-1-r$$

degrés de liberté, où $r$ est le nombre de paramètres du modèle estimés sur les données.

Les deux soustractions ont chacune leur cause. Le $-1$ vient de ce que les effectifs somment à $N$ : la dernière classe se déduit des autres. Le $-r$ vient de ce que chaque paramètre estimé rapproche mécaniquement le modèle des données, et diminue donc l'écart mesuré. Cette statistique est par ailleurs une approximation du GLRT, dont la dérivation fait l'objet du projet Tests classiques : le ${\chi }^2$ n'est pas une recette isolée, c'est une instance du cadre de la page précédente.

Confusion classique

« Le nombre de degrés de liberté vaut $K-1$. » Non, pas dès qu'un paramètre est estimé sur les données. Dans l'exemple ci-dessous, retenir $\nu =6$ au lieu de $5$ porte le seuil de $11,07$ à $12,59$ : le test ne rejette plus qu'une fois sur trente-sept sous $H_0$, au lieu d'une fois sur vingt.

Implémentation

Les quatre étapes du §2 de la page précédente valent ici sans changement, l'étape de calibration comprise :

1. grouper les observations en classes et compter les effectifs observés ;
2. estimer les paramètres du modèle, puis en déduire les effectifs attendus ;
3. évaluer la statistique $T$ ;
4. comparer au seuil, lu sur la loi du ${\chi }^2$ ou calibré par simulation sous $H_0$.

from scratch

import numpy as np
import scipy.stats as sta

edges = np.array([0, 2, 3, 4, 5, 6, 7])       # class lower bounds, last one open

def observed(x):
    """Counts per class; the last class gathers everything above its bound."""
    return np.histogram(x, bins=np.append(edges, 1e9))[0]

def expected(x):
    """Counts predicted by a Poisson model whose parameter is fitted on x."""
    lam = x.mean()                             # MLE of the Poisson parameter
    hi = sta.poisson.cdf(np.append(edges[1:] - 1, np.inf), lam)
    lo = sta.poisson.cdf(edges - 1, lam)
    return len(x) * (hi - lo)

def T(x):
    """Chi-square goodness-of-fit statistic."""
    o, e = observed(x), expected(x)
    return np.sum((o - e) ** 2 / e)

nu = len(edges) - 1 - 1                        # K - 1 - r, one parameter fitted
gamma = sta.chi2.ppf(0.95, nu)                 # >>> 11.0705  (nu = 5)

scipy

import scipy.stats as sta

# ddof counts the parameters fitted on the data; scipy then uses K - 1 - ddof
stat, p_value = sta.chisquare(observed(x), expected(x), ddof=1)

Le seuil peut aussi se calibrer sans invoquer la loi asymptotique, par la simulation de Monte-Carlo : engendrer $B$ échantillons sous $H_0$, réestimer les paramètres sur chacun, et lire le quantile empirique de $T$. Cette voie est plus lente mais plus sûre, puisqu'elle ne suppose ni $N$ grand ni effectifs suffisants.

Illustration

Si la proposition est exacte, l'histogramme de $T$ obtenu par simulation sous $H_0$ doit se superposer à la densité du ${\chi }^2$ à $5$ degrés de liberté, et le quantile empirique à $95\mathrm{\%}$ retomber sur le seuil théorique.

Figure 3.9, Loi de $T$ sous $H_0$ ($\lambda =4$, $N=200$, $50000$ répétitions) et densité ${\chi }^2(5)$ ; en rouge la zone de rejet au niveau $\alpha =5\mathrm{\%}$, au-delà de $\gamma =11,07$ (script figures/ch3/03-adequation-loi-h0.py).

L'histogramme épouse la densité théorique, et le quantile empirique à $95\mathrm{\%}$ vaut $11,02$ contre $11,07$ pour le ${\chi }^2(5)$ : la loi asymptotique est déjà acquise à $N=200$. La zone rouge représente exactement le budget de fausse alarme, la fraction des échantillons conformes au modèle que le test rejetterait pourtant.

Exemple : requêtes sur un serveur

Deux séries de $N=200$ comptages sont soumises au même test, au niveau $\alpha =5\mathrm{\%}$. La première est réellement issue d'une loi de Poisson ; la seconde provient d'un trafic par rafales, de même moyenne mais de variance double.

Série	Moyenne	Variance	$T$	Décision
Trafic régulier	$4,23$	$4,51$	$5,28$	non rejeté
Trafic par rafales	$4,06$	$7,26$	$49,1$	rejeté

La première série donne $T=5,28$, bien en deçà de $\gamma =11,07$ : rien dans ces données ne contredit le modèle de Poisson. La seconde donne $T=49,1$, plus de quatre fois le seuil. La colonne des variances dit pourquoi : la loi de Poisson impose $\mathrm{var}(X)=E(X)$, or la seconde série affiche une variance de $7,26$ pour une moyenne de $4,06$. Le test du ${\chi }^2$ chiffre ce que la comparaison de la moyenne et de la variance suggérait déjà en fiche Poisson.

3. La p-value

Le test précédent exige de fixer $\alpha$ avant de regarder les données. Rendre compte du résultat sans imposer ce choix au lecteur demande de publier, non pas la décision, mais le budget de fausse alarme auquel elle bascule.

Définition

Définition, p-value

La p-value d'une observation est la probabilité, sous $H_0$, d'obtenir une statistique au moins aussi grande que celle observée :

$$p=P(T\ge T_{\mathrm{obs}};H_0),$$

où

• $T_{\mathrm{obs}}$ est la valeur de la statistique sur les données reçues,
• la probabilité est calculée sous le modèle $H_0$, seul entièrement spécifié.

La p-value est donc la $P_{FA}$ du test dont le seuil vaudrait exactement $T_{\mathrm{obs}}$ : sur la figure 3.9, c'est l'aire sous la densité à droite du triangle bleu. Décider $H_1$ quand $p<\alpha$ est rigoureusement la règle $T>\gamma$, réécrite sur l'axe des probabilités. Pour les deux séries de l'exemple, $p=0,38$ et $p=2\cdot {10}^{-9}$.

Confusion classique

« $p=0,38$, donc le modèle a $38\mathrm{\%}$ de chances d'être vrai. » Non : la p-value est une probabilité sur les données, calculée en supposant $H_0$ vraie, jamais une probabilité sur l'hypothèse. C'est l'erreur déjà rencontrée sur l'intervalle de confiance, où les $95\mathrm{\%}$ portaient sur la procédure et non sur $\theta$.

4. Portée et limites

Illustration

Un test qui ne rejette pas ne dit rien tant que sa capacité à rejeter n'est pas connue. Confrontons donc le test au trafic par rafales, un modèle faux par construction, pour différentes tailles d'échantillon.

Figure 3.10, Fréquence de rejet du modèle de Poisson au niveau $\alpha =5\mathrm{\%}$ en fonction de $N$, sur des données réellement Poisson et sur des données surdispersées de même moyenne ($4000$ répétitions par point) (script figures/ch3/03-adequation-puissance.py).

La courbe orange reste collée au niveau visé quelle que soit la taille de l'échantillon : le test tient sa promesse de fausse alarme. La courbe bleue, elle, part de très bas : à $N=20$, un modèle pourtant faux n'est rejeté que dans $39\mathrm{\%}$ des cas, et il faut $N=100$ pour que le rejet devienne quasi certain. Sur un petit échantillon, un test qui ne rejette pas n'établit donc rien du tout, il constate seulement qu'il n'avait pas les moyens de conclure.

Confusion classique

« Le test ne rejette pas, donc le modèle est bon. » Non : ne pas rejeter, c'est manquer de preuves contre le modèle, jamais en avoir pour lui. La figure 3.10 le chiffre : à $N=20$, six fois sur dix, un modèle faux traverse le test sans encombre.

Le défaut symétrique guette à l'autre extrémité. Aucun modèle n'est exactement vrai, et la puissance croît avec $N$ : sur un échantillon assez grand, tout modèle finit par être rejeté, y compris un excellent modèle dont l'écart au réel est sans conséquence pratique. La décision d'un test d'adéquation se lit donc toujours avec la taille de l'échantillon à côté, et un écart statistiquement significatif n'est pas nécessairement un écart qui compte.

Extension : données continues (Kolmogorov-Smirnov)

Le ${\chi }^2$ exige un groupement en classes, naturel pour un comptage, arbitraire pour une variable continue : deux découpages différents peuvent conduire à deux conclusions différentes. Le test de Kolmogorov-Smirnov s'en dispense en comparant directement les deux fonctions de répartition :

$$D=\underset{x}{sup}|F_N(x)-F(x;{\theta }_0)|,$$

où $F(\cdot ;{\theta }_0)$ est la fonction de répartition du modèle et $F_N$ la fonction de répartition empirique, définie par

$$F_N(x)=\frac{\mathrm{\#}\{n:x[n]\le x\}}{N},$$

c'est-à-dire la proportion des mesures qui ne dépassent pas $x$. Comme toute fonction de répartition, elle croît de $0$ à $1$, ici par paliers de $1/N$. La statistique $D$ est l'écart vertical maximal entre les deux courbes, et sta.kstest la calcule.

Une réserve accompagne son usage courant. La loi de $D$ sous $H_0$ n'est tabulée que si ${\theta }_0$ est entièrement spécifié ; dès que des paramètres sont estimés sur les données, les seuils des tables deviennent trop permissifs et le seuil doit être calibré par simulation. Sur les chronos de marathon du chapitre 1, $D=\mathrm{0,060}$ pour $1998$ coureurs, à comparer au seuil tabulé de $\mathrm{0,030}$ et au seuil calibré par Monte-Carlo, $\mathrm{0,020}$ : la loi normale est rejetée par les deux, mais le seuil correct est une fois et demie plus sévère que celui de la table. Le diagramme quantile-quantile du chapitre 1 et cette valeur de $D$ disent la même chose, l'un en montrant où le modèle échoue, l'autre en chiffrant de combien.

5. Résumé

Deux tests ont été construits dans cette page, tous deux instances du cadre de la page précédente, appliqué cette fois au modèle lui-même.

Test	Statistique $T(\mathbf{x})$	Données	Loi sous $H_0$
${\chi }^2$	$\sum _k(n_k-{\hat{n}}_k{)}^2/{\hat{n}}_k$	groupées en classes	${\chi }^2(K-1-r)$, asymptotique
Kolmogorov-Smirnov	$\underset{x}{sup}|F_N(x)-F(x;{\theta }_0)|$	continues, non groupées	tabulée si ${\theta }_0$ connu, simulée sinon

Le choix entre les deux suit la nature des données, comptages contre mesures continues, et non la commodité. Dans les deux cas, le seuil se lit sur une loi de référence quand ses conditions sont réunies, et se calibre par simulation sous $H_0$ sinon, ce qui reste la voie la plus sûre dès que des paramètres sont estimés. Et dans les deux cas la conclusion s'énonce dans un seul sens : le modèle est rejeté, ou il résiste à ce test, sur cet échantillon, avec cette puissance.