Comparez des lois entières avec Kolmogorov-Smirnov et les noyaux

Jusqu’ici, nous avons comparé des positions : une médiane, un décalage, un rang moyen. Mais on veut parfois comparer des lois entières. Les temps entre deux requêtes sur un serveur suivent-ils bien une loi exponentielle, comme le suppose le modèle de dimensionnement ? Les latences mesurées en Europe et aux États-Unis ont-elles la même distribution, pas seulement la même médiane ?

Ce chapitre présente deux outils complémentaires. Le test de Kolmogorov-Smirnov compare des fonctions de répartition. L’estimation de densité par noyau dessine la forme d’une loi à partir des données, sans lui imposer de famille paramétrique.

Le test de Kolmogorov-Smirnov à un échantillon

Le fil rouge : des temps entre requêtes

Un générateur de charge est censé envoyer des requêtes selon un processus de Poisson d’intensité 1 par seconde : les temps entre deux requêtes devraient suivre une loi exponentielle de moyenne 1 s, de fonction de répartition $F_0(x) = 1 - e^{-x}$. On relève dix intervalles (en secondes), déjà triés :

0,12 ; 0,35 ; 0,48 ; 0,71 ; 0,95 ; 1,30 ; 1,62 ; 2,10 ; 2,85 ; 3,90

La statistique D

Au chapitre 1, nous avons défini la fonction de répartition empirique $F_n$, qui monte de $1/n$ à chaque observation. Si les données suivent $F_0$, $F_n$ devrait rester proche de $F_0$ partout. Le test de Kolmogorov-Smirnov mesure le plus grand écart vertical entre les deux courbes :

$$D_n = \sup_x \, \lvert F_n(x) - F_0(x) \rvert$$

où $\sup_x$ désigne la plus grande valeur atteinte quand $x$ parcourt tous les réels. Comme $F_n$ est en escalier, il suffit de regarder, en chaque observation $x_{(i)}$, les écarts juste avant et juste après le saut :

$$D_n = \max_i \max\left( \frac{i}{n} - F_0(x_{(i)}),\ F_0(x_{(i)}) - \frac{i-1}{n} \right)$$

Exemple

Le plus grand écart apparaît en $x_{(6)} = 1{,}30$ s. Juste avant le saut, $F_n$ vaut $5/10 = 0{,}5$, alors que $F_0(1{,}30) = 1 - e^{-1{,}30} \approx 0{,}727$. D’où

$$D_{10} \approx 0{,}727 - 0{,}5 = 0{,}227$$

Les données ont un peu moins de petites valeurs que prévu par la loi exponentielle.

La loi de D sous l’hypothèse nulle

Le résultat remarquable, dû à Kolmogorov, est que si $F_0$ est continue et entièrement spécifiée, la loi de $D_n$ sous $H_0$ ne dépend pas de $F_0$. Le test est donc libre de distribution. Pour $n$ petit, la loi exacte est tabulée ou calculée numériquement ; pour $n$ grand, $\sqrt{n}\,D_n$ suit approximativement la loi de Kolmogorov, et la valeur critique à 5 % vaut environ $1{,}36/\sqrt{n}$.

Ici, la valeur critique exacte à 5 % pour $n = 10$ est d’environ 0,41 (l’approximation $1{,}36/\sqrt{10}$ donne 0,43). Avec $D = 0{,}227$, on est loin du seuil : la p-valeur exacte vaut environ 0,60. Rien ne contredit la loi exponentielle.

Attention

« Ne pas rejeter » ne veut pas dire « la loi est exponentielle ». Avec dix observations, le test est très peu puissant : il faudrait un écart de plus de 0,41 entre les courbes pour le détecter. Un test d’ajustement ne prouve jamais qu’une loi est la bonne, il signale seulement les écarts flagrants.

Le piège des paramètres estimés

Et si l’on ne connaît pas la moyenne de la loi exponentielle ? La tentation est de l’estimer sur les données (ici, la moyenne vaut 1,438 s), puis de lancer le test avec cette loi ajustée. On obtient $D \approx 0{,}116$ et une p-valeur affichée de 0,997.

Piège

Si les paramètres de $F_0$ sont estimés sur les mêmes données, la loi de Kolmogorov ne s’applique plus : la loi ajustée colle artificiellement aux données, $D$ est trop petit et la p-valeur standard est beaucoup trop grande. On rejette trop rarement. Il faut alors une version corrigée : le test de Lilliefors (pour la loi normale ou exponentielle), un test dédié comme Shapiro-Wilk pour la normalité, ou un bootstrap paramétrique (chapitre 7). Dans scipy, stats.goodness_of_fit fait ce travail par simulation.

Ajuster la loi sur les données la rapproche artificiellement de Fₙ : D est divisé par deux et la p-valeur standard (0,997) n'a plus de sens.

Kolmogorov-Smirnov à deux échantillons

Le même principe compare deux échantillons, sans référence théorique. Avec $F_{n}$ et $G_{m}$ les fonctions de répartition empiriques des deux échantillons :

$$D_{n,m} = \sup_x \, \lvert F_n(x) - G_m(x) \rvert$$

Sous $H_0$ (les deux échantillons viennent de la même loi continue), la loi de $D_{n,m}$ ne dépend pas non plus de cette loi.

Exemple

Latences (ms) d’une API mesurées depuis l’Europe et depuis les États-Unis :

• Europe : 82, 88, 92, 95, 97, 101, 105, 110 ;
• États-Unis : 99, 108, 115, 118, 122, 125, 131, 140.

Parcourons quelques valeurs de $x$, en notant $F$ la répartition européenne et $G$ l’américaine :

$x$	97	101	105	110	115
$F(x)$	5/8	6/8	7/8	8/8	8/8
$G(x)$	0	1/8	1/8	2/8	3/8
Écart	0,625	0,625	0,75	0,75	0,625

Le maximum vaut $D = 0{,}75$ : 7 latences européennes sur 8 sont inférieures ou égales à 105 ms, contre une seule américaine. La p-valeur exacte, obtenue en énumérant les $\binom{16}{8} = 12\,870$ répartitions possibles, vaut environ 0,019. Les deux lois diffèrent.

La courbe européenne monte bien avant l'américaine. Juste après 105 ms, 7 latences européennes sur 8 sont passées contre 1 américaine : D = 0,75.

Quelle différence avec Mann-Whitney ? Mann-Whitney est surtout sensible à un décalage de position. Kolmogorov-Smirnov détecte n’importe quelle différence entre les lois : position, dispersion ou forme. En contrepartie, pour un pur décalage, il est moins puissant que Mann-Whitney. Il est aussi plus sensible au centre de la loi qu’aux queues, où les deux répartitions sont toutes deux proches de 0 ou de 1.

À retenir

Kolmogorov-Smirnov mesure le plus grand écart vertical entre deux fonctions de répartition (empirique contre théorique, ou empirique contre empirique). Il est libre de distribution si la loi de référence est continue et fixée à l’avance, et il détecte tout type de différence entre lois.

import numpy as np
from scipy import stats

x = np.array([0.12, 0.35, 0.48, 0.71, 0.95, 1.30, 1.62, 2.10, 2.85, 3.90])
print(stats.kstest(x, "expon"))            # D ≈ 0.227, p ≈ 0.60 (loi exponentielle de moyenne 1)

europe = [82, 88, 92, 95, 97, 101, 105, 110]
usa    = [99, 108, 115, 118, 122, 125, 131, 140]
print(stats.ks_2samp(europe, usa))         # D = 0.75, p ≈ 0.019

x <- c(0.12, 0.35, 0.48, 0.71, 0.95, 1.30, 1.62, 2.10, 2.85, 3.90)
ks.test(x, "pexp", rate = 1)               # D = 0.227, p ≈ 0.60

europe <- c(82, 88, 92, 95, 97, 101, 105, 110)
usa    <- c(99, 108, 115, 118, 122, 125, 131, 140)
ks.test(europe, usa)                       # D = 0.75, p ≈ 0.019

En data science, le test à deux échantillons sert souvent à détecter une dérive de données (data drift) : on compare, variable par variable, la distribution des données de production à celle des données d’entraînement.

Estimer une densité par noyau

Un histogramme donne une idée de la forme d’une loi, mais il dépend beaucoup du choix des classes et produit une courbe en escalier. L’estimation de densité par noyau (kernel density estimation, KDE) le remplace par une courbe lisse.

Le principe

On pose une petite « bosse » sur chaque observation, puis on additionne toutes les bosses :

$$\hat f_h(x) = \frac{1}{nh} \sum_{i=1}^{n} K\!\left(\frac{x - X_i}{h}\right)$$

• $K$ est le noyau (kernel), une densité symétrique : le plus souvent la densité normale centrée réduite $K(u) = e^{-u^2/2}/\sqrt{2\pi}$, ou le noyau d’Epanechnikov $K(u) = \frac{3}{4}(1 - u^2)$ pour $|u| \le 1$ ;
• $h > 0$ est la largeur de bande (bandwidth), qui règle la largeur de chaque bosse ;
• le facteur $1/(nh)$ garantit que $\hat f_h$ est d’intégrale 1.

Exemple

Six temps de trajet (en minutes) : 2,1 ; 2,5 ; 3,0 ; 3,2 ; 4,8 ; 5,5. Estimons la densité en $x = 3$ avec un noyau gaussien et $h = 0{,}5$. Chaque observation contribue $\frac{1}{nh} K\left(\frac{3 - X_i}{h}\right)$ :

$X_i$	2,1	2,5	3,0	3,2	4,8	5,5
$(3 - X_i)/h$	1,8	1,0	0	−0,4	−3,6	−5,0
Contribution	0,026	0,081	0,133	0,123	0,000	0,000

Total : $\hat f_{0,5}(3) \approx 0{,}363$. Les observations lointaines (4,8 et 5,5) ne contribuent presque rien : l’estimation est locale.

Chaque observation apporte une bosse d'aire 1/n ; l'estimation est leur somme. En x = 3, seules les bosses des observations proches comptent.

Le choix de h compte plus que celui du noyau

Le noyau a peu d’influence sur le résultat ; la largeur de bande, énormément. C’est un compromis biais-variance :

• $h$ trop petit : chaque observation produit un pic étroit, la courbe est hérissée (sous-lissage, forte variance) ;
• $h$ trop grand : toutes les bosses se fondent, la courbe est trop lisse et gomme les modes réels (sur-lissage, fort biais).

Sur nos six temps de trajet, avec $h = 0{,}3$, la densité estimée vaut 0,46 en $x = 3$ mais seulement 0,01 en $x = 4$ : deux pics séparés. Avec $h = 1{,}5$, elle vaut 0,20 en $x = 3$ et 0,19 en $x = 4$ : une seule bosse aplatie.

Des règles automatiques proposent un $h$ raisonnable quand la loi n’est pas trop éloignée d’une normale :

• la règle de Silverman : $h = 0{,}9 \min\left(\hat\sigma, \frac{\mathrm{IQR}}{1{,}34}\right) n^{-1/5}$, où $\hat\sigma$ est l’écart-type empirique et IQR l’écart interquartile. C’est le défaut de R (bw.nrd0). Ici, $\hat\sigma \approx 1{,}34$, IQR $\approx 1{,}78$ (quartiles calculés à la manière de R), d’où $h \approx 0{,}83$ ;
• la règle de Scott : $h = \hat\sigma\, n^{-1/5}$, le défaut de scipy.stats.gaussian_kde, ici $h \approx 0{,}94$.

Pour une loi multimodale ou très asymétrique, ces règles lissent trop ; on préfère alors la validation croisée ou la méthode de Sheather-Jones (bw = "SJ" en R).

La largeur de bande décide de ce que l'on voit : h = 0,3 hérisse la courbe (0,46 en 3 mais 0,01 en 4), h = 1,5 gomme toute structure (0,20 et 0,19).

import numpy as np
from scipy import stats

t = np.array([2.1, 2.5, 3.0, 3.2, 4.8, 5.5])
kde = stats.gaussian_kde(t)                  # règle de Scott par défaut
print(kde([3.0, 4.0]))                       # environ [0.26 0.19]
print(kde.factor * t.std(ddof=1))            # largeur de bande effective, environ 0.94

kde_fin = stats.gaussian_kde(t, bw_method=0.3 / t.std(ddof=1))   # h = 0.3
print(kde_fin([3.0, 4.0]))                   # environ [0.46 0.01]

Attention à la convention de scipy : bw_method est un facteur multiplié par l’écart-type des données, pas la largeur de bande elle-même. En R, density prend directement $h$ :

t <- c(2.1, 2.5, 3.0, 3.2, 4.8, 5.5)
d <- density(t)             # noyau gaussien, h de Silverman (bw.nrd0)
d$bw                        # environ 0.83
plot(d); rug(t)
plot(density(t, bw = 0.3))  # sous-lissage : deux pics

Attention

Un noyau gaussien place de la masse partout, y compris là où la variable ne peut pas aller : des durées ou des prix négatifs, par exemple. Près d’une borne naturelle, on estime la densité du logarithme de la variable, ou on utilise une correction de bord.

En entretien

« Une KDE pose un noyau sur chaque point et fait la moyenne. Le vrai réglage, c’est la largeur de bande : trop petite, on voit du bruit ; trop grande, on efface la structure. Je pars de la règle de Silverman ou de Scott, puis je vérifie visuellement ou par validation croisée. »

En résumé

• Le test de Kolmogorov-Smirnov mesure le plus grand écart vertical entre fonctions de répartition : $D_n = \sup_x |F_n(x) - F_0(x)|$.
• Sa loi sous $H_0$ ne dépend pas de $F_0$ si celle-ci est continue et fixée à l’avance ; avec des paramètres estimés sur les données, il faut Lilliefors ou un bootstrap.
• La version à deux échantillons détecte toute différence de lois (position, dispersion, forme), mais est moins puissante que Mann-Whitney pour un simple décalage.
• Ne pas rejeter une loi ne prouve pas qu’elle est la bonne, surtout avec peu de données.
• L’estimation par noyau $\hat f_h(x) = \frac{1}{nh}\sum K((x - X_i)/h)$ lisse les données ; la largeur de bande $h$ compte bien plus que le noyau.
• scipy.stats.kstest, ks_2samp, gaussian_kde ; en R, ks.test et density.

Exercices

Exercice · Calculer D à la main

facile

On teste si quatre valeurs, 0,1 ; 0,4 ; 0,5 ; 0,9, suivent la loi uniforme sur $[0, 1]$, pour laquelle $F_0(x) = x$. Calculez $D_4$.

Voir la solution

En chaque observation $x_{(i)}$, on calcule $i/n - F_0(x_{(i)})$ et $F_0(x_{(i)}) - (i-1)/n$ :

$i$	$x_{(i)}$	$i/4 - x_{(i)}$	$x_{(i)} - (i-1)/4$
1	0,1	0,15	0,10
2	0,4	0,10	0,15
3	0,5	0,25	0,00
4	0,9	0,10	0,15

$D_4 = 0{,}25$, atteint juste après $x = 0{,}5$ (trois valeurs sur quatre sont inférieures ou égales à 0,5, alors qu’on en attendrait la moitié). C’est très loin de la valeur critique à 5 % pour $n = 4$ (environ 0,62) : on ne rejette pas l’uniformité.

En code

Python

import numpy as np
from scipy import stats

x = np.sort([0.1, 0.4, 0.5, 0.9])
n = len(x)
i = np.arange(1, n + 1)
print(i / n - x)                  # [0.15 0.1  0.25 0.1 ]
print(x - (i - 1) / n)            # [0.1  0.15 0.   0.15]

res = stats.kstest(x, "uniform")  # loi uniforme sur [0, 1]
print(res.statistic, res.pvalue)  # 0.25 0.906
print(stats.kstwo.ppf(0.95, n))   # valeur critique à 5 % : 0.624

R

x <- sort(c(0.1, 0.4, 0.5, 0.9))
n <- length(x); i <- 1:n
i / n - x                  # 0.15 0.10 0.25 0.10
x - (i - 1) / n            # 0.10 0.15 0.00 0.15

ks.test(x, "punif")        # D = 0.25, p-value = 0.9062

Exercice · Une densité par noyau à la main

moyen

Avec les six temps de trajet de l’exemple (2,1 ; 2,5 ; 3,0 ; 3,2 ; 4,8 ; 5,5), estimez la densité en $x = 3$ avec le noyau d’Epanechnikov $K(u) = \frac{3}{4}(1 - u^2)$ pour $|u| \le 1$ (0 sinon) et $h = 1$. Comparez au résultat gaussien avec $h = 0{,}5$.

Voir la solution

Avec $h = 1$, $u_i = 3 - X_i$ : 0,9 ; 0,5 ; 0 ; −0,2 ; −1,8 ; −2,5. Seules les quatre premières valeurs vérifient $|u| \le 1$ :

• $K(0{,}9) = 0{,}75 \times 0{,}19 = 0{,}1425$ ;
• $K(0{,}5) = 0{,}75 \times 0{,}75 = 0{,}5625$ ;
• $K(0) = 0{,}75$ ;
• $K(-0{,}2) = 0{,}75 \times 0{,}96 = 0{,}72$.

Somme : 2,175. Donc $\hat f_1(3) = 2{,}175/(6 \times 1) \approx 0{,}363$. C’est presque exactement la valeur obtenue avec le noyau gaussien et $h = 0{,}5$ (0,363) : un noyau d’Epanechnikov de largeur 1 a un écart-type de $1/\sqrt 5 \approx 0{,}45$, comparable à celui du noyau gaussien de largeur 0,5. Ce qui compte, c’est l’échelle de lissage, pas la forme du noyau.

En code

Python

import numpy as np
from scipy import stats

t = np.array([2.1, 2.5, 3.0, 3.2, 4.8, 5.5])

def epanechnikov(u):
    return np.where(np.abs(u) <= 1, 0.75 * (1 - u ** 2), 0.0)

h = 1
contrib = epanechnikov((3 - t) / h)
print(contrib)                          # [0.1425 0.5625 0.75 0.72 0. 0.]
print(contrib.sum() / (len(t) * h))     # 0.3625

# Noyau gaussien avec h = 0.5, pour comparer
print(stats.norm.pdf((3 - t) / 0.5).sum() / (len(t) * 0.5))   # 0.3629

R

t <- c(2.1, 2.5, 3.0, 3.2, 4.8, 5.5)
epanechnikov <- function(u) ifelse(abs(u) <= 1, 0.75 * (1 - u^2), 0)

h <- 1
contrib <- epanechnikov((3 - t) / h)
contrib                              # 0.1425 0.5625 0.7500 0.7200 0 0
sum(contrib) / (length(t) * h)       # 0.3625

# Noyau gaussien avec h = 0.5, pour comparer
sum(dnorm((3 - t) / 0.5)) / (length(t) * 0.5)   # 0.3629

Exercice · Choisir entre KS et Mann-Whitney

difficile

Deux serveurs ont la même latence médiane, 100 ms. Sur le serveur 1, les latences sont très regroupées (entre 95 et 105 ms). Sur le serveur 2, elles sont très dispersées (entre 40 et 160 ms), mais symétriques autour de 100 ms.

1. Un test de Mann-Whitney sur deux grands échantillons a-t-il de bonnes chances de détecter une différence ? Et Kolmogorov-Smirnov ?
2. Lequel choisiriez-vous pour surveiller une dérive des latences en production ?

Voir la solution

1. Mann-Whitney porte sur $P(X > Y)$. Par symétrie autour de la même médiane, une latence du serveur 2 a à peu près une chance sur deux de dépasser une latence du serveur 1 : $\theta \approx 1/2$, et le test ne détectera pratiquement rien, même avec beaucoup de données. Kolmogorov-Smirnov compare les fonctions de répartition entières : celle du serveur 1 monte brutalement autour de 100 ms, celle du serveur 2 monte doucement de 40 à 160 ms. L’écart maximal est grand (juste avant 95 ms, par exemple, $F_1 \approx 0$ alors que $F_2$ vaut près de la moitié), et le test rejettera nettement.
2. Pour une surveillance de dérive, on veut détecter tout changement de distribution, y compris une hausse de variabilité sans changement de médiane (qui dégrade l’expérience des utilisateurs dans la queue). Kolmogorov-Smirnov est donc plus adapté, en complément d’un suivi direct des quantiles élevés (95e ou 99e centile).

En code

Python

import numpy as np
from scipy import stats

# Même médiane (100 ms), dispersions très différentes
rng = np.random.default_rng(0)
serveur1 = rng.uniform(95, 105, size=500)
serveur2 = rng.uniform(40, 160, size=500)

print(stats.mannwhitneyu(serveur1, serveur2).pvalue)   # 0.94 : rien détecté
res = stats.ks_2samp(serveur1, serveur2)
print(res.statistic, res.pvalue)                       # D = 0.466, p ≈ 2e-49

R

# Même médiane (100 ms), dispersions très différentes
set.seed(0)
serveur1 <- runif(500, 95, 105)
serveur2 <- runif(500, 40, 160)

wilcox.test(serveur1, serveur2)$p.value   # 0.96 : rien détecté
ks.test(serveur1, serveur2)               # D = 0.462, p-value < 2.2e-16