1 Historique et principe

1.1 Anderson-Darling

Le test d’Anderson-Darling a pour origine les deux papiers suivants :

co-écrits par Theodore Wilbur Anderson (1918-2016), auteur de “An Introduction to Multivariate Statistical Analysis”

Figure 1. T. W. Anderson

Figure 1. T. W. Anderson

et Donald A. Darling (1915-2014), tous les deux professeurs à la Columbia.

Figure 2. D. A. Darling

Figure 2. D. A. Darling

Il existe deux versions du test :

  • Un test qui permet de vérifier si un échantillon est distribué selon une loi connue

  • Un test qui permet de vérifier si \(K\) échantillons sont distribués selon une même loi (pas nécessairement connue).

1.2 Cramér-Von Mises

Il y a deux papiers distincts à l’origine de ce test :

  • Cramér, H. (1928). On the Composition of Elementary Errors. Scandinavian Actuarial Journal. 13–74. doi:10.1080/03461238.1928.10416862.

  • von Mises, R. E. (1928). Wahrscheinlichkeit, Statistik und Wahrheit. Julius Springer.

Harald Cramér (1893–1985), professeur à l’université de Stockholm (a reçu la médaille Guy décernée par la société royale de la statistique - https://en.wikipedia.org/wiki/Guy_Medal)

Figure 3. H. Cramér

Figure 3. H. Cramér

et Richard von Mises (professeur à l’université de Harvard) a parmi ses nombreux travaux travaillé sur le paradoxe des anniversaires (https://fr.wikipedia.org/wiki/Paradoxe_des_anniversaires)

Figure 4. Richard Von Mises

Figure 4. Richard Von Mises

A l’origine, le test permet de vérifier si un échantillon est distribué selon une loi connue. Puis, il a été généralisé à deux échantillons par Anderson

  • Anderson, T. W. (1962). On the Distribution of the Two-Sample Cramer–von Mises Criterion (PDF). Annals of Mathematical Statistics. Institute of Mathematical Statistics. 33 (3): 1148–1159. doi:10.1214/aoms

On va voir que le test de Cramér-Von Mises peut être vu comme un cas particulier du test d’Anderson-Darling. On va donc présenter le principe du test d’Anderson-Darling et on présentera alors le test de Cramér-Von Mises.

2 Test de comparaison d’un échantillon à une distribution

Bibliographie : les pages suivantes présentent bien le principe du test d’Anderson-Darling

NB : la page Wikipedia française ne parle que du cas particulier où on teste la normalité.

Principe : on dispose d’un échantillon et on souhaite savoir si cet échantillon pourraît être issu d’une loi de distribution connue. L’hypothèse nulle testée est donc

\[H_o:\{\text{L'échantillon observée est distribuée selon } F\}\] contre

\[H_1:\{\text{L'échantillon observé n'est pas distribuée selon } F\}\] où \(F\) est la fonction de répartition \(F_X(x)=P[X\leq x]=\int_{-\infty}^xf_x(x)dx\) d’une loi continue connue.

Pour le test d’Anderson-Darling, il semble que les lois pour lesquelle une distribution asymptotique de la statistique de test ait été calculée sont : normale, uniforme, log-normale, exponentielle, Weibull, extreme value type I, generalized Pareto, and logistic distributions.

Pour le test de Cramér-von Mises, la statistique de test est “distribution-free” et semble marcher pour le cas variable continue, comme variable discrète.

Exemple : on considère un échantillon de taille \(n=100\) simulé selon une loi uniforme \(\mathbb{U}_{[0,1]}\).

n <- 100
set.seed(123)
u <- runif(n, min = 0, max = 1)

Question : on souhaite vérifier si l’échantillon pourrait être issu d’une loi gaussienne \(\mathcal{N}(\mu=0.5, \sigma=0.2)\).

2.1 Analyse descriptive

L’outil statistique utilisé pour comparer des distributions est la fonction de répartition : la version empirique \(F_n\) obtenue sur les données et la version théorique \(F\) obtenue en choisissant une loi connue (ici la loi gaussienne) et ses paramètres (\(\mu\) et \(\sigma\)).

On peut représenter ces courbes sur R :

  • la fonction ecdf() permet de calculer \(F_n\) qui vaut : \(F_n(x)=\displaystyle \frac{1}{n}\sum_{i=1}^n1(x_i\leq x)\) pour un échantillon \(x_1,\ldots,x_n\)

  • la fonction pnorm() permet de calculer \(F_{gaussien}\) et donne pour une valeur de \(x\) la probabilité d’observer une valeur inférieure ou égale à \(x\) si \(X\) était distribuée selon un loi gaussienne de paramètres \(\mu\) et \(\sigma\).

fdr.u <- ecdf(u)
plot(fdr.u, main = "Distribution de u et la loi théorique N(0.5,0.2)", 
     xlim = c(-2, 2),
     lwd = 2, col = "magenta")
plot(function(x) pnorm(x, mean = 0.5, sd = 0.25), 
     from = -2, to = 2, add = T,
     col = "black", lwd = 2)
legend("topleft", c("Fn(u)",  "F(gaussienne)"),
       lwd = 2, col = c("magenta", "cyan"), cex =0.7)

2.2 Principe du test

Il s’agit de savoir si la fonction de répartition empirique est proche de la théorique que l’on souhaite tester. Pour cela, on va définir une mesure de distance entre \(F\) et \(F_n\) :

\[n\int_{-\infty}^{+\infty}(F_n(x) - F(x))^2w(x)dF(x)\]

Si on prend \(w(x)=1\) on obtient la distance de Cramér-von Mises (dans ce cas toutes les observations ont le même poids).

\[\omega^2=n\int_{-\infty}^{+\infty}(F_n(x) - F(x))^2dF(x)\]

Si on prend \(w(x)=(F(x)(1-F(x)))^{-1}\) (on donne plus de poids aux points en queues de distribution), on obtient la statistique de test d’Anderson-Darling :

\[A^2=n\int_{-\infty}^{+\infty}\frac{(F_n(x) - F(x))^2}{F(x)(1-F(x))}dF(x)\]

Remarque: la statistique de test fait intervenir \(F\) qui dépend de \(\theta\) (par exemple, dans le cas d’une loi gaussienne \(\theta=(\mu,\sigma)\). Dans la forme basique du test, \(\theta\) est suposé connu (c’est l’exemple qu’on a pris où \(\mu=0.5\) et \(\sigma=0.2\)), mais selon qu’on doit estimer \(\theta\) à partir de l’échantillon cela aura un impact sur le test et les valeurs critiques.

Application :

  • on programme la fonction \(F_n\) ainsi que la fonction \(F\) (pour la loi gaussienne) et on les représente (ainsi que l’échantillon u avec la fonction rug():
Fn <- function(x, my_samp) {
    n <- length(my_samp)
    res <- numeric(length(x))
    for (i in 1:length(x)) {
      res[i] <- 1 / n * sum(my_samp <= x[i])
    }
    return(res)
}
Fx <- function(x, mean, sd) {
  pnorm(x, mean = mean, sd = sd)
}
plot(function(x) Fx(x, mean = 0.5, sd = 0.2), -0.5, 1.5, col = "magenta", ylab = "F")
plot(function(x) Fn(x, u), -0.5, 1.5, add = T, col = "cyan") 
rug(u)
legend("topleft", c("F",  "Fn"),
       lwd = 2, col = c("magenta", "cyan"), cex =0.7)

On représente la fonction à intégrer :

  • Pour Cramér-von Mises, on a \((F_n(x) - F(x))^2dF(x)=(F_n(x) - F(x))^2\frac{dF(x)}{dx}dx=(F_n(x) - F(x))^2f(x)dx\)
f_1 <- function(x, my_samp) {
  Fx_eval <- Fx(x, mean = 0.5, sd = 0.2) 
  return((Fn(x, my_samp) - Fx_eval) ^ 2 * dnorm(x, mean = 0.5, sd = 0.2))
}
plot(function(x) f_1(x, u), 0, 1) 
rug(u)

Et il suffit d’intégrer :

n * integrate(function(x) f_1(x, u), lower = min(u), upper = max(u), 
          subdivisions = 1000)$value
## [1] 1.041762
  • Pour Anderson-Darling, on a \(\frac{(F_n(x) - F(x))^2}{F(x)(1-F(x))}dF(x)=\frac{(F_n(x) - F(x))^2}{F(x)(1-F(x))}\frac{dF(x)}{dx}dx=\frac{(F_n(x) - F(x))^2}{F(x)(1-F(x))}f(x)dx\)
f_2 <- function(x, my_samp) {
  Fx_eval <- Fx(x, mean = 0.5, sd = 0.2) 
  return((Fn(x, my_samp) - Fx_eval) ^ 2 / 
    (Fx_eval * (1 - Fx_eval)) * dnorm(x, mean = 0.5, sd = 0.2))
}
plot(function(x) f_2(x, u), 0, 1) 
rug(u)

Remarque : au-delà de l’intervalle \([0,1]\), il y a des conditions d’intégralité qui ne sont pas respectées à cause du terme \(\frac{1}{F(x)(1-F(x))}\)

On calcule néanmoins l’intrégale sur un domaine où la fonction semble intégrable :

n * integrate(function(x) f_2(x, u), lower = min(u), upper = max(u), 
          subdivisions = 1000)$value
## [1] 9.705243

On remarque que cette intégrale peut être décomposée par morceaux :

u_sort <- sort(u)
my_int <- 0
for(i in 1:(n - 1)) {
  my_int <- my_int + integrate(function(x) f_2(x, u), 
              lower = u_sort[i], upper = u_sort[i + 1], 
              subdivisions = 2)$value
}
n * my_int
## [1] 9.704962

En pratique, on remplace les intégrales par des sommes en utilisant les simplifications suivantes :

  • Pour Cramér-von Mises, une simplification est la suivante (\(Y_i\) doit d’abord être ordonnée):

\[T=\omega^2=\frac{1}{12n}+\sum_{i=1}^n(\frac{2i-1}{2n}-F(Y_i))^2\]

u_order <- sort(u)
(my_T <- 1 / (12 * length(u)) + sum(((2 * (1:n) - 1) / (2 * n) - Fx(u_order, mean = 0.5, sd = 0.2)) ^ 2))
## [1] 1.04036
  • Pour Anderson-Darling, une simplification de la formule est la suivante (\(Y_i\) doit d’abord être ordonnée):

\[A^2=-n-\displaystyle\sum_{i=1}^n\frac{2i-1}{n}[\log(F(Y_i))+\log(1-F(Y_{n+1-i})))]\]

u_order <- sort(u)
-n - sum((2 * (1:n) - 1) / n * (log(Fx(u_order, mean = 0.5, sd = 0.2)) + 
                                log(1 - Fx(rev(u_order), mean = 0.5, sd = 0.2)))) 
## [1] 9.709207

L’intrepétation est la suivante et la même pour les deux tests : si \(A^2\) ou si \(\omega^2\) est trop grande, on ne pourra pas accepter \(H_0\).

2.3 Calcul de la p-value

2.3.1 Anderson-Darling

La statistique de test \(A^2\) peut être comparée avec une table théorique (obtenue en obtenant la distribution asymptotique de \(A^2\) à partir de la fonction caractérisique de \(A^2\)) qui dépend :

  • la famille de distribution choisie \(F\) (normale, uniforme, exponentielle, etc)

  • la méthode utilisée pour l’estimation des paramètres dans \(F\). Par exemple, dans le cas où on teste la normalité, on a 4 possibilités :

    • cas 0 : \(\mu\) et \(\sigma\) sont connus
    • cas 1 : \(\sigma\) est connue, mais la moyenne est estimée par \(\bar{x}\)
    • cas 2 \(\mu\) est connue, mais \(\sigma\) est estimée par \(\frac{1}{n}\sum_{i=1}^n(x_i-\mu)^2\)
    • cas 3 : \(\mu\) est estimée par \(\bar{x}\) et \(\sigma\) est estimée par \(\frac{1}{n}\sum_{i=1}^n(x_i-\bar{x})^2\)

La table (pour la loi normale) est la suivante (pour les autres familles de distribution, il faudra utiliser d’autres tables qui ont été présentées dans le chapitre de livre suivant écrit par M. A. Stephens en 1986 : Tests Based on EDF Statistics, dans Goodness-of-Fit Techniques) :

Dans notre exemple, on est dans la situation 0 (\(\mu\) et \(\sigma\) connus, qui est la situation prise par Anderson et Darling). La valeur de \(A^2\) est 9.709207, on a donc une \(p\)-value inférieure à \(1\%\), ce qui montre que cette valeur est anormalement grande : on ne peut donc pas accepter l’hypothèse nulle et l’échantillon u n’est pas distribué selon une loi normale de paramètre \(\mu=0.5\) et \(\sigma = 0.2\)

2.3.2 Cramér-von Mises

Si on ne connaissait pas la distribution de \(\omega^2\), on pourrait estimer la distribution de \(\omega^2\) sous \(H_0\) en faisant du MCMC.

En effet, on peut simuler 1000 échantillons distribués selon une \(\mathcal{N}(\mu=0.5, \sigma=0.2)\) et calculer la statistique de test pour chaque échantillon; on obtient ainsi une distribution MCMC de \(\omega^2\).

T_mcmc <- numeric(1000)
for (k in 1:1000) {
  my_s <- rnorm(length(u), 0.5, 0.2)
  my_s_order <- sort(my_s)
  T_mcmc[k] <- 1 / (12 * length(my_s_order)) + sum(((2 * (1:n) - 1) / (2 * n) - Fx(my_s_order, mean = 0.5, sd = 0.2)) ^ 2) 
}

On représente ensuite l’histogramme et on regarde où se situe la vraie valeur de \(\omega^2\) :

hist(T_mcmc, col = "royalblue")
abline(v = my_T, col = "red", lty = 2)

\(\omega^2\) a l’air bien trop grande par rapport à ce que cela aurait dû être si l’échantillon était distribué selon une loi normale de paramètres \(\mu=0.5\) et \(\sigma=0.2\). On pourrait prendre comme \(p\)-value le pourcentage de fois où les statistiques de tests sur échantillon simulé sont supérieurs à \(\omega^2\) :

length(which(T_mcmc > T)) / 1000
## [1] 0.003

Pleusieurs articles se sont employés à trouver la distribution asymptotique de \(\omega^2\) dont Smirnov dans un premier temps et Anderson et Darling (1952) ensuite ou encore Durbin. On pourra par exemple citer la formule 6.6 de https://bookdown.org/egarpor/NP-UC3M/nptests-gof-dist.html pour voir une formule de l’ecdf de \(\omega^2\) (cette statistique ne dépend pas de la distribution de \(F\) et est dite “distribution free”) et vaut

\[\lim_{n->\infty}P[\omega^2\leq x]=1 - \frac{1}{\pi}\sum_j(-1)^{j-1}\int_{(2j-1)^2\pi^2}^{4j^2\pi^2}\sqrt{\frac{-\sqrt{y}}{\sin\sqrt{y}}}\frac{e^{-xy/2}}{y}dy\]

Lorsqu’on teste l’hypothèse de normalité, Stephens a proposé une formule analytique de la \(p\)-value (programmée dans le package nortest)

  • Stephens M.A. (1974). EDF statistics for goodness of fit and some comparisons. J Am Stat Assoc, 69.

2.4 Implémentation

2.4.1 Anderson-Darling

  • le package nortest (Gross et Ligges, 2015) permet de calculer le test d’Anderson Darling pour l’hypothèse de normalité dans le cas 3 présenté ci-dessus (\(\mu\) est estimée par \(\bar{x}\) et \(\sigma\) est estimée par \(\frac{1}{n}\sum_{i=1}^n(x_i-\bar{x})^2\)).
nortest::ad.test(u)
## 
##  Anderson-Darling normality test
## 
## data:  u
## A = 1.2686, p-value = 0.002543

Le test confirme que l’échantillon n’est pas distribué selon une loi gaussienne de paramètre \(\mu=\bar{x}\) et \(\sigma^2=\frac{1}{n}\sum_{i=1}^n(x_i-\bar{x})^2\)).

On vérifie qu’on aurait trouvé la même chose en faisant du MCMC :

u_order <- sort(u)
my_t_AD <- -n - sum((2 * (1:n) - 1) / n * (log(Fx(u_order, mean = mean(u), sd = sd(u))) + 
                                log(1 - Fx(rev(u_order), mean = mean(u), sd = sd(u))))) 
my_samp <- numeric(10000)
for (k in 1:n) {
  u_mcmc <- rnorm(length(u), mean(u), sd = sd(u))
  u_order <- sort(u_mcmc)
  my_samp[k] <- -n - sum((2 * (1:n) - 1) / n * (log(Fx(u_order, mean = mean(u), sd = sd(u))) + 
                                log(1 - Fx(rev(u_order), mean = mean(u), sd = sd(u))))) 
}
hist(my_samp)
abline(v = my_t_AD, lty = 2)

length(which(my_samp > my_t_AD)) / 10000
## [1] 0.0029
  • le package goftest (Faraway et al., 2019) permet de calculer le test d’Anderson Darling pour n’importe quelle hypothèse sur la distribution de \(F\) (argument null=) et la méthode utilisée pour l’estimation des paramètres dans \(F\). Ici il y a deux cas possible : cas 0 en utilisant l’option estimated=F :
goftest::ad.test(u, null = pnorm, mean = 0.5, sd = 0.2,
                 estimated = FALSE)
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: Normal distribution
##  with parameters mean = 0.5, sd = 0.2
##  Parameters assumed to be fixed
## 
## data:  u
## An = 9.7092, p-value = 1.738e-05

ou cas 3 en utilisant l’option estimated = TRUE; dans ce dernier cas la p-value est calculée en utilisant une méthode basé sur du re-échantillonage :

goftest::ad.test(u, null = pnorm, mean = mean(u), sd = sd(u),
                 estimated = TRUE)
## 
##  Anderson-Darling test of goodness-of-fit
##  Braun's adjustment using 10 groups
##  Null hypothesis: Normal distribution
##  with parameters mean = 0.498558994238265, sd = 0.284995269404707
##  Parameters assumed to have been estimated from data
## 
## data:  u
## Anmax = 1.7798, p-value = 0.7303

Toutefois, cela conduit à des résultats très étranges car très différent de ce que donne la fonction nortest::ad.test() (voir https://stats.stackexchange.com/questions/392563/why-do-two-implementations-of-the-anderson-darling-test-produce-such-different-p un début d’explication)

  • Un autre package ADGofTest propose une implémentation du test AD (cas 0) basé sur l’article suivant :

G. and J. Marsaglia (2004). Evaluating the Anderson-Darling Distribution, Journal of Statistical Software.

ADGofTest::ad.test(u, distr.fun = pnorm, mean = 0.5, sd = 0.2)
## 
##  Anderson-Darling GoF Test
## 
## data:  u  and  pnorm
## AD = 9.7092, p-value = 1.738e-05
## alternative hypothesis: NA

2.4.2 Cramér-von Mises

  • Dans le package nortest, il y a aussi la fonction cvm.test() qui permet de calculer le test de Cramér-von Mises pour l’hypothèse de normalité. La statistique de test est légèrement modifiée (\(\omega^2(1+0.5/n)\))
nortest::cvm.test(u)
## 
##  Cramer-von Mises normality test
## 
## data:  u
## W = 0.18036, p-value = 0.009333
  • Dans le package goftest, il y a la fonction cvm.test() qui permet d’utiliser plusieurs hypothèses sur la loi (argument null=) et la méthode utilisée pour l’estimation des paramètres dans \(F\) (soit \(\theta\) est connu, soit il est estimé)
goftest::cvm.test(u, null = pnorm, mean = 0.5, sd = 0.2,
                 estimated = FALSE)
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: Normal distribution
##  with parameters mean = 0.5, sd = 0.2
##  Parameters assumed to be fixed
## 
## data:  u
## omega2 = 1.0404, p-value = 0.001911

ou cas 3 en utilisant l’option estimated = TRUE; dans ce dernier cas la p-value est calculée en utilisant une méthode basé sur du re-échantillonage :

goftest::cvm.test(u, null = pnorm, mean = mean(u), sd = sd(u), 
                  estimated = TRUE)
## 
##  Cramer-von Mises test of goodness-of-fit
##  Braun's adjustment using 10 groups
##  Null hypothesis: Normal distribution
##  with parameters mean = 0.498558994238265, sd = 0.284995269404707
##  Parameters assumed to have been estimated from data
## 
## data:  u
## omega2max = 0.67651, p-value = 0.1203

2.5 Efficacité du test par rapport à d’autres tests

Le test d’Anderson-Darling ne semble pas aussi bon que le test de Shapiro-Wilk, mais il semble meilleur que les autres tests de normalité. Dans la page Wikipedia, les tests connexes à celui d’Anderson-Darling sont : Kolmogorov–Smirnov test, Kuiper’s test, Shapiro–Wilk test, Jarque–Bera test.

Voir Stephens (1974) pour une comparaison approfondie de ces tests d’ajustements.

3 Test de comparaison de \(K\) échantillons

3.1 Principe

3.1.1 Anderson-Darling

Les deux articles d’Anderson-Darling ne traitent pas ce test. L’article qui présente ce test est le suivant :

  • Scholz, F. W. et Stephens, M. A. (1987). K-sample Anderson–Darling Tests. Journal of the American Statistical Association. 82 (399): 918–924. doi:10.1080/01621459.1987.10478517.

Les auteurs généralisent d’abord la statistique de test \(A^2\) et la généralise pour comparer deux échantillons entre eux :

\[A^2_{mn}=\frac{mn}{N}\int_{-\infty}^{\infty}\frac{(F_m(x)-G_n(x))^2}{H_N(x)(1-H_N(x))}dH_N(x)\]

où :

  • \(F_m(x)\) est la fonction de répartition empirique d’un échantillon \(X_1,\ldots,X_m\)
  • \(G_n(x)\) est la fonction de répartition empirique d’un échantillon \(Y_1,\ldots,Y_n\)
  • \(H_N(x)\) est la fonction de réparition empirique des échantillons mélangés et peut s’exprimer ainsi \(\frac{1}{N}(mF_m(x)+nG_n(x))\), avec \(N=m+n\)

Remarque : ici, on n’utilise pas la distribution \(F\) ou \(G\) théorique et on n’a donc pas d’hypothèses à faire sur le paramètre \(\theta\) des lois.

L’hypothèse nulle est alors :

\[H_o:\{F=G\text{ (les deux échantillons ont la même distribution)}\}\] contre

\[H_1:\{\text{Les deux échantillons n'ont pas la même distribution}\}\]

Les auteurs généralisent ensuite ce test à \(k\) échantillons de taille \(n_1,\ldots,n_k\) et de distributions \(F_1, \ldots, F_k\). L’hypothèse nulle est alors

\[H_o:\{F_1=\ldots =F_k\}\] contre

\[H_1:\{\text{Il existe au moins un échantillon différent des autres}\}\] La statistique de test est alors :

\[A_{kN}^2=\displaystyle\sum_{i=1}^k\int \frac{(F_{i}(x)-H_N(x))^2}{H_N(x)(1-H_N(x))}dH_N(x)\]

avec \(N=\sum_{i=1}^k n_i\) pour laquelle les auteurs trouvent une formule simplifiée :

Soit \(Z_1,\ldots,Z_N\) l’échantillon ordonné (sans ex-aquo) qui regroupe les \(k\) échantillons, la formule devient :

\[A_{kN}^2=\frac{1}{N}\displaystyle\sum_{i=1}^k\frac{1}{n_i}\sum_{j=1}^{N-1} \frac{(NM_{ij} - jn_i)^2}{j(N-j)}\] où \(M_{ij}\) est le nombre d’observations de l’échantillon \(i\) qui sont inférieures ou égales à \(Z_j\).

Il y a également une deuxième version en cas de valeurs discrètes avec des répétitions (voir formule 7 dans https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.1087.8833&rep=rep1&type=pdf)

Exemple :

On considère deux échantillons k1 de taille \(n_1=6\) et \(k2\) de taille \(n_2=8\) distribués selon deux lois uniformes \(\mathbb{U}_{[0,1]}\) et \(\mathbb{U}_{[0.1,1.1]}\).

n1 <- 6
n2 <- 8
n_i <- list(n1, n2)
N <- n1 + n2
set.seed(123)
k1 <- runif(n1, min = 0, max = 1)
k2 <- runif(n2, min = 0.1, max = 1.1)
k_i <- list(k1, k2)
Z <- sort(c(k1, k2))

On représente les distributions de k1, k2 et Z :

plot(function(x) Fn(x, k1), 0, 1.1, col = "magenta", lwd = 2, ylab = "Fn")
plot(function(x) Fn(x, k2), 0, 1.1, col = "cyan", add = T, lwd = 2) 
plot(function(x) Fn(x, Z), 0, 1.1, col = "grey", add = T, lwd = 2) 
rug(k1, col = "magenta", lwd = 1)
rug(k2, col = "cyan", lwd = 1)
legend("topleft", c("Fn(k1)",  "Fn(k2)", "Fn(Z)"),
       lwd = 2, col = c("magenta", "cyan", "grey"), cex =0.7)

On calcule la version 1 (sans ex-aquo):

A_2 <- 0
for (k in 1:2) {
  for (j in 1:(N - 1)) {
    A_2 <- A_2 + (N * length(which(k_i[[k]] <= Z[j])) - j * n_i[[k]]) ^ 2 / 
      (j * (N - j)) / n_i[[k]]
  }
}
A_2 / N
## [1] 1.222883

On calcule la version 2 (avec ex-aquo):

A_2 <- 0
Z <- sort(unique(c(k1, k2)))
f_ij <- matrix(1, 2, length(Z))
f_ij[1, ] <- Z %in% k1
f_ij[2, ] <- Z %in% k2
l_j <- apply(f_ij, 2, sum)
M_ij <- f_ij
M_ij[1, ] <- cumsum(M_ij[1, ]) - f_ij[1, ] / 2
M_ij[2, ] <- cumsum(M_ij[2, ]) - f_ij[2, ] / 2
B_j <- cumsum(l_j) - l_j / 2
for (k in 1:2) {
  for (j in 1:(N)) {
    A_2 <- A_2 + (N * M_ij[k, j] - B_j[j] * n_i[[k]]) ^ 2 / 
      (B_j[j] * (N - B_j[j]) - N * l_j[j] / 4) / n_i[[k]]
  }
}
(N - 1) / N * A_2 / N
## [1] 1.268899

Interprétation : les auteurs ont calculé la distribution asymptotique de la statistique de test normalisée \(T_{kN}\frac{A^2_{kN}-(k-1)}{\sigma_N}\) où la valeur de \(\sigma_N\) peut être trouvée dans Scholz et Stephens (1987). Si la valeur est trop grande (dans la table ci-dessous \(m=k-1\)), on rejette l’hypothèse nulle.

3.1.2 Cramér-von Mises

C’est Anderson (1962) qui est à l’origine de ce test (comparaison de 2 échantillons). Les hypothèses \(H_0\) et \(H_1\) sont identiques à celles du test précédent :

\[H_o:\{F=G\text{ (les deux échantillons ont la même distribution)}\}\] contre

\[H_1:\{\text{Les deux échantillons n'ont pas la même distribution}\}\]

La statistique de test est :

\[W_{n,m}^2=\frac{nm}{(n+m)}\int(F_n(z)-G_m(z))^2dH_{n+m}(z)\]

où \(H_{n+m}\) est la fonction de répartition empirique de l’échantillon mélangé \(Z_1,\ldots,Z_{n+m}\).

En partant de l’égalité \(H_{n+m}(z)=\frac{n}{n+m}F_n(z)+\frac{m}{n+m}G_m(z)\), on obtient la formule suivante pour la statistique de test :

\[W_{n,m}^2=\frac{nm}{(n+m)^2}\sum_{k=1}^{n+m}(F_n(Z_k)-G_m(Z_k))^2\] Autrement dit, on s’intéresse à la différence quatratique entre les deux courbes, calculé sur chaque point de l’échantillon mélangé :

plot(function(x) Fn(x, k1), 0, 1.1, n = 1001, col = "magenta", lwd = 2, ylab = "Fn")
plot(function(x) Fn(x, k2), 0, 1.1, n = 1001, col = "cyan", add = T, lwd = 2) 
rug(k1, col = "magenta", lwd = 1)
rug(k2, col = "cyan", lwd = 1)

for(k in 1:N) {
  segments(Z[k], Fn(Z[k], k1), Z[k], Fn(Z[k], k2), col = "grey", lwd = 2, lty = 2) 
}

legend("topleft", c("Fn(k1)",  "Fn(k2)"),
       lwd = 2, col = c("magenta", "cyan"), cex =0.7)

On calcule la statistique de test en créant notre fonction :

CvM_2samples <- function(k1, k2) {
  n <- length(k1)
  m <- length(k2)
  Z <- c(k1, k2)
  N <- length(Z)
  my_T <- 0
  for(k in 1:N) {
    my_T <- my_T + (Fn(Z[k], k1) - Fn(Z[k], k2)) ^ 2
  }
 return(my_T <- (n1 * n2) / N ^ 2 * my_T)
}
(my_T <- CvM_2samples(k1, k2))
## [1] 0.1904762

Anderson (1962), puis Burr (1963) ont étudié la distribution de \(W_{n,m}^2\) et ont donné des tables pour de petits échantillons.

Une façon de savoir si cette statistique est anormalement grande serait de permuter aléatoirement les valeurs des deux échantillons (on crée dans ce cas là deux échantillons qui seraient issus de la même distribution \(Z\)) et on calcule la statistique de test sur chaque de ces échantillons permutés :

permutation_test <- function(k1, k2) {
  B <- 1000
  my_mcmc <- numeric(B)
  n <- length(k1)
  m <- length(k2)
  Z <- c(k1, k2)
  N <- length(Z)
  for(k in 1:B) {
    permus <- sample(Z)
    my_mcmc[k] <- CvM_2samples(permus[1:n], permus[(n+1):N])
  }
  return(my_mcmc)
}

On représente la distribution des statistiques de test permutés :

my_perm <- permutation_test(k1, k2) 
hist(my_perm, col = "royalblue")
abline(v = my_T, col = "red", lty = 2)

Et pour la p-value, on compte le nombre de fois où la statistique de test des échantillons permutés est plus grande que la vraie valeur :

length(which(my_perm > my_T)) / 1000
## [1] 0.315

Ainsi, le test de Cramér-von Mises indique que les deux échantillons pourraient appartenir à la même distribution.

3.2 Implémentation du test sous R

3.2.1 Anderson-Darling

Le package kSamples (Fritz Scholz and Angie Zhu, 2019) contient une fonction ad.test() qui permet de réaliser ce test. Il utilise la syntaxe à la ggplot2 ou des modèles linéaires :

library(kSamples)
## Loading required package: SuppDists
my_df <- data.frame(y = c(k1, k2),
                    group = c(rep("1", n1), rep("2", n2)))
kSamples::ad.test(y ~ group, data = my_df)
## 
## 
##  Anderson-Darling k-sample test.
## 
## Number of samples:  2
## Sample sizes:  6, 8
## Number of ties: 0
## 
## Mean of  Anderson-Darling  Criterion: 1
## Standard deviation of  Anderson-Darling  Criterion: 0.67069
## 
## T.AD = ( Anderson-Darling  Criterion - mean)/sigma
## 
## Null Hypothesis: All samples come from a common population.
## 
##               AD   T.AD  asympt. P-value
## version 1: 1.223 0.3323           0.2480
## version 2: 1.270 0.4009           0.2306

Ici, la statistique de test normalisée n’est pas anormale et on ne peut donc pas rejeter l’hypothèse nulle; les deux échantillons pourraient donc appartenir à la même distribution.

3.2.2 Cramér-von Mises

A partir de la statistique de test calculé ci-dessus, on peut utiliser la fonction goftest::pCvM() pour obtenir la significativité car la statistique de test \(W_{n,m}^2\) a la même distribution que \(\omega^2\) (voir https://bookdown.org/egarpor/NP-UC3M/nptests-comp.html)

1 - goftest::pCvM(q = my_T)
## [1] 0.286705

Une alternative est d’utiliser le package twosamples (qui contient plusieurs tests de comparaison de deux échantillons dont Kolmogorov-Smirnov, Kuiper, Cramer-von Mises, Anderson-Darling, Wasserstein, and DTS)

twosamples::cvm_test(k1, k2, nboots = 1000)
## Test Stat   P-Value 
## 0.7777778 0.1670000 
## attr(,"details")
##      n1      n2 n.boots 
##       6       8    1000

Sinon, il existe un code en C++ qui permet de faire le job (https://www.jstatsoft.org/article/view/v017i08)

3.3 Comparaison avec d’autres test

Un test dans le même registre est le test de de Kruskal-Wallis qui teste si la médiane d’échantillons (dont on connaît pas la distribution) est la même ou non. A noter que ce test n’a pas exactement la même vocation que A-D ou C-vM et peuvent donc aboutir à des résultats différents.