5.2.8 Test de Kolmogorov-Smirnov

Ce test est plus naturel que celui du $\chi ^2$ lorsque la loi attendue a une fonction de répartition continue. Soit X une v.a. réelle dont la loi a une fonction de répartition continue F et X₁, X₂,..., X_m, m réalisations indépendantes d'un processus dont on suppose qu'il suit la même loi que X (on cherche à tester cette hypothèse). Le test consiste à mesurer l'écart entre la fonction de répartition exacte et la fonction de répartition empirique :

$$K_m = \sqrt{m} \sup _{-\infty < x < +\infty} \vert F(x) - F_m(x) \vert$$

et à le comparer à une valeur << admissible >>. À la limite, la loi de la variable aléatoire K_m a la fonction de répartition suivante :

$$\lim_{m \rightarrow +\infty} P(K_m \le x) = H(x) = 1 - 2 \sum_{j=1}^{+\infty} (1)^{j-1}e^{-2j^2x^2}$$

et de même que pour le test du $\chi ^2$ on rejette l'hypothèse lorsque :

$$K_m \gt H^{-1}(1 - \alpha)$$

avec $\alpha = 0.05$ par exemple. Si on utilise l'approximation $H(x) \simeq 1 - e^{-2x^2}$. alors la valeur seuil à ne pas dépasser est :

$$K_{seuil} = \sqrt{\frac{1}{2} \mbox{ln}(\frac{1}{\alpha})}$$

On dispose cependant d'une expression asymptotique pour ce seuil correspondant à la loi exacte suivie par K_m (suffisamment précise pour m>30) :

$$K_{seuil} = \sqrt{\frac{1}{2} \mbox{ln}(\frac{1}{\alpha})} - \frac{1}{6 \sqrt{m}} + O(\frac{1}{m})$$

Le calcul de K_m ne pose pas de problème si on trie le vecteur $X^r = (X_1,\; X_2,\; \dots,\; X_m)$.Supposons ce tri effectué, en remarquant que :

$$\sup _{x \in [X_i,X_{i+1}[} F_m(x) - F(x) = \frac{i}{m} - F(... ..._{x \in [X_i,X_{i+1}[} F(x) - F_m(x) = F(X_{i+1}) - \frac{i}{m}$$

les deux quantités suivantes se calculent facilement :

$$K_m^+ = \sqrt{m} \sup _{-\infty < x < +\infty} ( F_m(x) - F(... ... = \sqrt{m} \max _{1 \le j \le m} ( \frac{j}{m} - F(X_j) ) \\$$

$$K_m^- = \sqrt{m} \sup _{-\infty < x < +\infty} ( F(x) - F_m(x) ) = \sqrt{m} \max _{1 \le j \le m} (F(X_j) - \frac{j-1}{m} )$$

et l'on obtient alors $K_m = \max(K_m^+,\; K_m^-)$.

Nous allons illustrer ce test sur le processus stochatisque suivant (où U désigne la loi uniforme sur [0,1]) :

$$X_n(t) = \frac{1}{\sqrt{n}} \sum_{i=1}^n (1_{\{U_i \le t\}} - t)$$

qui est tel que pour $t \in ]0,1[$ fixé, on a :

$$\lim _{n \rightarrow +\infty} X_n(t) = {\cal N} (0, t(1-t))$$

Pour cela on se propose de prendre n << assez grand >> et de réaliser m simulations à n fixé et finalement d'observer la convergence en loi, d'abord graphiquement (en superposant la fonction de répartition empirique avec la fonction de répartition exacte) puis en réalisant ce test de Kolmogorov-Smirnov. Avec les raccourcis d'écriture permis en Scilab, la fonction principale pour obtenir une réalisation de X_n(t) peut s'écrire de la façon suivante :

function [X] = pont_brownien(t,n)
  //
  X = sum(bool2s(grand(n,1,"def") <= t) - t)/sqrt(n)

Voici maintenant un script possible, qui, après m simulations du processus, affiche le graphe de la fonction de répartition empirique puis lui superpose celui de la fonction de répartition attendue et finalement effectue le test statistique :

  // Le pont Brownien
  // Petite simulation pour illustrer que : Xn(t) --> N(0,t(1-t))  pour n -> infini
  //
  // où  Xn(t) = somme_{i=1}^n ( 1_(Ui <= t) - t ) /sqrt(n)
  //
  t = 0.3;
  sigma = sqrt(t*(1-t));  // l'écart type attendu
  n = 1000;               // n "grand"
  m = 4000;               // le nb de simulations
  X = zeros(m,1);         // initialisation du vecteur des realisations
  for k=1:m
    X(k) = pont_brownien(t,n);  // la boucle pour calculer les realisations
  end
  repartition_empirique(X)      // le dessin de la fonction de repartition empirique
  x = linspace(min(X),max(X),60)';               // les abscisses et
  [P,Q]=cdfnor("PQ",x,0*ones(x),sigma*ones(x));  // les ordonnees pour la fonction exacte
  plot2d(x,P,2,"000")                            // on l'ajoute sur le dessin initial 

  // mise en place du test KS
  alpha = 0.05 
  X = - sort(-X); // tri
  FX = cdfnor("PQ",X,0*ones(X),sigma*ones(X));
  Dplus = max( (1:m)'/m - FX );
  Dmoins = max( FX - (0:m-1)'/m );
  Km = sqrt(m)*max([Dplus ; Dmoins]);
  K_seuil = sqrt(log(1/alpha)/2) - 1/(6*sqrt(m)) ;

  // affichage des resultats
  //
  write(%io(2)," Test KS : ")
  write(%io(2)," -------- ")
  write(%io(2)," valeur obtenue par le test : "+string(Km))
  write(%io(2)," valeur seuil a ne pas depasser : "+string(K_seuil))
  if (Km > K_seuil) then
    write(%io(2)," Conclusion provisoire : Hypothese rejetee !")
  else
    write(%io(2)," Conclusion provisoire : Hypothese non rejetee !")
  end

Voici les résultats obtenus avec n=1000 et m = 4000 (cf figure(5.6)) :

 Test KS : 
 -------- 
 valeur obtenue par le test : 1.1204036
 valeur seuil a ne pas depasser : 1.2212382
 Conclusion provisoire : Hypothese non rejetee !

  

Figure: fonctions de répartition empirique et exacte

Pour finir nous allons, pour n fixé, dessiner la courbe $X_n(t), t \in [0,1]$. Il est assez simple de voir, qu'après le tri croissant des U_i (en supposant aussi les U_i tous distincts et différents de et 1 et en posant U₀ = 0 et U_n+1 = 1) que :

$$X_n(t) = \frac{i - n t}{\sqrt{n}}, \; \mbox{ pour } U_i \le t < U_{i+1}$$

On peut alors calculer cette courbe (en reliant les discontinuités par des segments verticaux) avec la fonction suivante :

 function [] = dessin_pont_brownien(n)
   //
   U = -sort(-rand(1,n))
   Xbas = ((0:n-1) - n*U)/sqrt(n)
   Xhaut = ((1:n) - n*U)/sqrt(n)
   absc = [0 ; matrix([U ; U], 2*n, 1) ; 1]
   ord  = [0 ; matrix([Xbas ; Xhaut], 2*n, 1) ; 0]
   xbasc()
   plot2d(absc,ord)

La figure (5.7) montre une courbe obtenue avec n=10000 avec 3 zooms successifs qui peuvent laisser entrevoir le caractère fractal de la << courbe limite >>.

  

Figure 5.7: Un pont Brownien...