Arthur Charpentier

Toujours pour répondre à une question posée par mail ("ça quoi ça sert le test lillie sous R ? "), un court billet. Avant de parler de Lillie¹, il faut revenir sur Kolmogorov-Smirnov. L’idée est de tester

contre
Je prends une forme paramétrique pour la fonction que l’on cherche à tester car c’est ce qui est souvent fait en pratique. Et surtout, ça permettra de mieux comprendre (je pense) l’idée du test lilllie. L’idée du test repose sur le théroème de Glivenko-Cantelli. En fait, Andreï Kolmogorov et Vladimir Smirnov ont montré que l’on pouvait encadrer la vitesse de convergence,
ce qui se montre en notant que  converge en loi (si l’hypothèse H₀ était la bonnne) vers un pont brownien changé de temps par la fonction quantile .
La loi limite est appelée loi de Kolmogorv-Smirnov, et ne dépend pas de la loi sous-jacente. Par contre elle n’est qu’asymptotique. Pour récupérer la loi à distance finie, regardons les simulations suivantes. On notera en particulier que la loi de la différence ne dépend pas de la vraie loi des observations (même si c’est cette loi que l’on teste, comme toujours lorsque l’on cherche la loi sous l’hypothèse H₀)
> n=20
> ns=100000
> SU=SN=SE=rep(NA,ns)
> for(i in 1:ns){
+ X=runif(n)
+ SU[i]=max(c(abs(rank(X)/n-X),abs((rank(X)-1)/n-X)))
+ X=rnorm(n,0,1)
+ SN[i]=max(c(abs(rank(X)/n-pnorm(X,0,1)),abs((rank(X)-1)/n-pnorm(X,0,1))))
+ X=rexp(n,1)
+ SE[i]=max(c(abs(rank(X)/n-pexp(X,1)),abs((rank(X)-1)/n-pexp(X,1))))
+ }
> plot(density(SN),col="red",lwd=2)
> lines(density(SU),col="blue",lwd=2)
> lines(density(SE),col="green",lwd=2)
avec la densité de la loi du maximum entre la fonction de répartition théorique et la fonction de répartition empirique en rouge pour l’échantillon Gaussien, en bleu pour l’échantillon uniforme et en vert pour l’échantillon exponentiel. Peu importe la loi sous-jacente, la statistique de test est la même.

On peut en particulier récupérer les quantiles de cette lois pour n=20 observations (la loi du test dépend simplement du nombre d’observations).

> quantile(SE,.95)
      95%
0.2937023
> quantile(SE,.90)
     90%
0.264661
On peut d’ailleurs confirmer ces valeurs en regardant dans les tables (ici ou là), par exemple 0.294 à 95%,

Mais ce test ne sert que pour juger de l’adéquation à une loi précise (que nous avions noté  ici). En pratique, on ne connaît pas cette loi, mais on peut l’estimer, en particulier pour les lois paramétriques. On peut alors chercher à tester autre chose, comme
contre
qui s’apparente plus à un test d’ajustement à une famille de lois. Ce test étant différent du précédant, il n’y a pas de raison d’avoir la même loi pour la statistique de test (qui néanmoins, à  peut rester celle que nous avions considérée, à savoir la norme infinie entre les fonctions de répartition). Graphiquement par exemple, on cherchait aurparavant à tester l’ajustement de la loi rouge, et maintenant on va tester l’ajustement de la loi bleue,

qui va s’adapter à l’échantillon, contrairement au test initial,

Le code ici, dans le cas Gaussien par exemple (la loi sous-jacente n’est plus neutre dans ce cas) donne
> n=20
> ns=100000
> SN=SNML=rep(NA,ns)
> for(i in 1:ns){
+ X=rnorm(n,0,1)
+ SN[i]  =max(c(abs(rank(X)/n-pnorm(X,0,1)),
+               abs((rank(X)-1)/n-pnorm(X,0,1))))
+ SNML[i]=max(c(abs(rank(X)/n-pnorm(X,mean(X),sd(X))),
+               abs((rank(X)-1)/n-pnorm(X,mean(X),sd(X)))))
+ }
> plot(density(SN),col="blue",lty=2)
> lines(density(SNML),col="red",lwd=2)
Les densités sont sensiblement différentes,

et les quantiles aussi,
> quantile(SN,.95)
      95%
0.2943206
> quantile(SNML,.95)
      95%
0.1923295
Sous R, deux fonctions distinctes sont à utiliser. Par exemple l’ajustement à une loi  donne
>  X=rnorm(20)
> mean(X)
[1] -0.553641
> sd(X)
[1] 0.8293167
>  ks.test(X,"pnorm", mean=0,sd=1)

        One-sample Kolmogorov-Smirnov test

data:  X
D = 0.3096, p-value = 0.03362
alternative hypothesis: two-sided
(on rejette ici l’ajustement à une loi normale centrée réduite) et l’ajustement à une loi  donne
> lillie.test(X)

        Lilliefors (Kolmogorov-Smirnov) normality test

data:  X
D = 0.1519, p-value = 0.2629
avec une p-value de l’ordre de 26%. Notons que cette p-value est interprétable, contrairement au test suivant, qui n’est pas valide,
>  ks.test(X,"pnorm", mean=mean(X),sd=sd(X))

        One-sample Kolmogorov-Smirnov test

data:  X
D = 0.1519, p-value = 0.69
alternative hypothesis: two-sided
Encore une fois, la statistique de test est la même, sauf que la loi de la statistique de test sous H₀ n’est pas juste dans le dernier cas. En particulier, n sera beaucoup plus exigent dans le test car on a choisi la meilleure loi normale possible (à l’aide d’un critère de maximisation de vraisemblance). Sinon je renvoie ici pour un document sur le même sujet.
¹ Lillie pour Hubert Lilliefors qui, le premier, a proposé ce test

Lors d’une discussion avec quelques anciens collègues samedi midi, nous évoquions des prénoms, et du fait que les prénoms rares (pour ne pas dire bizarres) étaient de plus en plus fréquents. Je me suis demandé si les gens donnaient vraiment plus des prénoms originaux qu’avant. Était-ce une idée en l’air (du buzz médiatique, genre ici ou là), ou au contraire est-elle (statistiquement) fondée ?

Pour mesurer la concentration des prénoms (ou de n’importe quoi d’ailleurs), le plus simple est peut être d’utiliser la courbe de Lorenz. En statistique descriptive, on construit cette courbe de la manière suivante: à partir d’un échantillon que l’on supposera ordonné, i.e. , la courbe de Lorenz est le nuage de points où et

La relecture probabiliste est que l’on trace où

Je renvoie ici ou là pour plus d’information sur la courbe de Lorenz.
Si je reviens à mon histoire de prénoms, on pourra noter qu’en 2003, par exemple¹, il y a eu 795493 naissances, donnant lieu à 59481 prénoms (j’entends par à des prénoms d’orthographe différente),
> head(base)
    sexe preusuel annais nombre rang
32     1    AADIL   2003      3 8775
71     1    AARON   2003    135  764
100    1    ABASS   2003      3 8301
133    1    ABBAS   2003      5 6897
171    1    ABBES   2003      3 8816
228    1      ABD   2003     11 3982
Si on trie, on obtient le classement suivant,
> head(baser)
    sexe preusuel annais nombre rang
300101    1      LEA   2003   8988    1
110747    1    LUCAS   2003   8292   2
168477    1     THEO   2003   7855    3
77273    1     HUGO   2003   7386    4
318115    1    MANON   2003   6917    5
Autrement dit, le top 5 des prénoms représente 5% des naissances, et les 100 premiers représentent 43% des naissances. En 1900, les 5 premiers représentaient 22% de l’ensemble, et 80% pour les 100 premiers. Afin de résumer l’information, on notera, par exemple, que 2,5% des prénoms (les prénoms les plus donnés) représentaient 50% des naissances en 2003. C’est ceci que l’on récupère via la courbe de Lorenz.
Le graphique ci-dessous montre combien représentaient 10% des prénoms (les prénoms le plus donnés, en rouge) au cours du temps, et 2,5% des prénoms (en bleu). Autrefois, si l’on prenait les 2,5% des prénoms les plus donnés, ils étaient portés par 9O% de la population en 1900 (100-10), contre 50% de la population en 2004 (100-50) [des compléments sur la lecture des graphiques ont été apportés ici]

On peut alors effectivement parler de "course à l’originalité" depuis 1970. En 1900, seuls 5% des naissances donnaient lieu à des prénoms "originaux", contre 50% aujourd’hui. Étrangement, on notera qu’il y a eu des "pics d’originalité" pendant les guerres. Si quelqu’un a une explication ?... On peut aussi visualiser ci-dessous les courbes de Lorenz,

Les courbes montrent que les prénoms sont de moins en moins "concentrés" (pour reprendre une terminologie classique sur les courbes de Lorenz). On le retrouve également sur les prénoms à la mode qui furent très stables par le passé, et qui depuis les années 70 n’arrêtent pas de changer... Bref, il y a de plus en plus de prénoms originaux, certes, mais les prénoms classiques ou à la mode restent donnés par beaucoup de monde, l’indice de Gini restant très élevé. Et je ne suis pas assez calé en statistique textuelle pour voir si ce qui change, c’est le nombre de prénoms (avec un total de 6556 prénoms en 1900 contre 10 fois plus en 2004) ou simplement l’orthographe,
       sexe       preusuel annais nombre
 313765    1       MAEL   2004     37
313861    1      MAELE   2004     27
314056    1      MAELL   2004      3
314127    1     MAELLE   2004   1577
315525    1      MAHEL   2004      3
etc....
La prochaine fois, j’essayerais de reprendre une étude intéressante sur les liens entre les prénoms et les professions (ici) pour étudier les prénoms des actuaires.... oui, j’ai toute la base ! à suivre donc...

¹ Je passe sur les soucis techniques de programmation. En effet, il existe un libellé étrange pour les "prénoms rares",
       sexe       preusuel annais nombre 
408961    1 _PRENOMS_RARES   2003  26253    
185715    1 _PRENOMS_RARES   2003  23718    
mais comme on a le tail on peut construire la courbe de Lorenz proprement... enfin, au moins par la partie de droite car j’ai postulé car les prénoms rares étaient donnés une unique fois (mais a priori c’est une ou deux, et je n’ai pas moyen de le savoir.... je n’ai que le nombre de naissances avec un prénom rare). Mais ça n’influence que la partie gauche qui peut être un peu différente....

Un petit billet très court pour répondre à une question d’élèves de licence. L’erreur quadratique moyen est définie comme

alors que l’erreur absolue moyenne est définie comme
(que je définie ici comme l’écart à la moyenne, certain considérant la distance à la médiane, qui minimise cette quantitié, alors que la moyenne minimise l’écart quadratique, mais je reviendrais là dessus bientôt, sur l’histoire de la loi normale).
Souvent ces quantités sont présentées comme des mesures de la même quantité, que l’on appelerait "dispersion" autour de la moyenne.
Si l’on prend deux observations, ces quantités sont rigoureusement identiques, car
compte tenu du fait que
En revanche, comme le montre le graphique ci-dessous,
Aussi, avec davantage de termes, la différence peut devenir significativement différente. En particulier, comme le notent Goldstein et Taleb (2007), en ligne ici, dans le cas Gaussien centré réduit,
Autrement dit, ces quantités peuvent parfaitement être différentes. Sur l’exemple ci-dessous, on considère un groupe de 4 individus, trois ayant le même revenu, et le dernier ayant un revenu plus important. On s’arrange pour que le revenu moyen soit constant (100 ici) et on regarde comment se comportent ces deux quantités lorsque le revenu du dernier individu s’accroit (de celui des autres).

(j’ai autorisé les revenus négatifs pour la simplicité des calculs et pour avoir un joli dessin). On en conclue qu’un écart significativement différent entre l’écart-type (L2) et l’écart absolu (L1) signifie simplement qu’il y a probablement des points aberrants dans la base (aberrants au sens "sensiblement différents", des "outliers").

Samedi soir, Djalil et Florent se demandaient ce que Karl Pearson avait inventé en premier, la loi du khi (ou du chi), ou la loi du khi-deux (pensant que la loi du chi-deux était relativement naturelle, et que la loi du khi découlait de la loi du khi-deux). La loi du khi admet pour densité
et est aussi parfois définie à partir de la fonction Gamma incomplète, i.e. sa fonction de répartition est
où
avec
Pour la loi du khi-deux, sa densité est
dont la fonction de répartition est
Bref, on peut définir ces deux lois complètement indépendemment de leurs applications: si on considère des variables indépendantes , alors
suit une loi du khi-deux, alors que
suit une loi du khi. Bref, intuitivement, on définit d’abord la loi du chi-deux avant d’introduire la loi du khi, si l’on part de l’interprétation. Et si l’on reprend le texte fondamental de Fisher, ici, datant de 1922, on ne retrouve que le khi-deux

Mais dans "The lady drinking tea", David Salsburg suggère que Karl Pearson avait introduit le khi avant le khi-deux... Si l’on y regarde de plus près, Karl Pearson a proposé le "test du chi-deux" en 1900 (ici), avec l’idée que l’exposant qui apparait dans la loi normale pouvait être noté, comme on le retrouve ci-dessous 

Effectivement, si l’on regarde ici par exemple, la loi du chi-deux est reliée à l’étude du volume d’un ellipsoïde. Mais ce n’est pas la première apparition de cette variation dite du khi-deux. Si on y regarde de plus près, par exemple ici, ou là, on notera que Karl Pearson avait défini, dès 1896 aussi bien la loi du khi que la loi du khi-deux, dans le cadre de ce qui s’appelle les distribution de Pearson de type III,

Autrement dit, Karl Pearson a commencé par définir les deux lois simultanément, via les fonction gamma incomplète, et autres équations d’analyse fonctionnelle, et c’est seulement plus tard que l’utilisation de la loi du khi-deux l’a emporté, de part son interprétation relativement élégante (somme de carrés de lois normales centrées réduites) et le test d’indépendance (ou d’adéquation) du khi-deux.
Mais en 1931, Karl Pearson renvoit à un papier plus ancien où l’on trouve pour la première fois la loi du khi-deux, suggérant d’ailleurs d’appeler cette loi la loi de Helmert, en mémoire à Friedrich Robert Helmert, qui avait également beaucoup travailler sur les moindres carrés (dont il est naturel d’arriver à des sommes de carrés de lois normales, ici), qui aurait trouvé cette loi en 1875. Mais on peut également remonter aux travaux d’Ernst Karl Abbe en 1863 qui avait également mentionné la loi.

Mais c’est de Karl Pearson que l’on a hérité la notation et le nom du khi-deux.

I will start here a short post on extreme values, with some historical perspective. In a recent paper (in French), I mentioned the use of the Pareto distribution as a standard model for extremes, but if reinsurers have been using the Pareto distribution for a long time (see here e.g.), the oldest mathematical models when dealing with extreme value should be related to work on maximum values in finite samples.

• The work of Ronald Fisher and Leonard Tippett

Leonard Henry Tippett, a former student of Karl Pearson published in Biometrika a note on extremes, in 1925. The goal was "the determination of the distribution of the range and the extremes for a large number of samples". In 1925, everyone was looking for the Gaussian distribution everywhere, and Leonard Tippett observed that the distribution of the largest value did not have a Gaussian distribution.
A few years after, a joint work with Ronald Fisher was presented to the Cambridge Philosophical Society. The starting point was the idea of "stability" (even if the term did not appear explicitely in their work): the limiting distribution the maximum should be of the "same type" as the underlying distribution. Thus, if  stands for the cumulative distribution function, it should satisfy functional equation
Solutions of that functional equation will give all possible limiting distributions. Thus, Fisher and Tippett obtained three possible limits,

• solutions of , i.e. 
• solutions of , i.e.  with  (i.e. finite lower bound for the support), i.e. 
• solutions of , i.e.  if  (i.e. finite upper bound for the support), i.e. 

Based on those possible limiting distributions, Fisher and Tippett wanted to derive what has been called later on the "domain of attraction" of those distributions.

• The work of Maurice Fréchet, at the same time

In 1926, Maurice Fréchet wrote a paper on "la loi de probabilité de l’écart maximum". That paper, as well as the one by Fisher and Tippett (wrote at the same time), investigated asymptotic limits. Both obtained functional equations, but only Maurice Fréchet understood the importance of the stability concept, pointed out by Paul Levy in the context of sums. Thus, Maurice Fréchet introduced the concept of what is called now "max-stability". But Fréchet solve only functional equation . The point is that Fréchet studied absolute values of errors, i.e. strictly positive random variables. Thus, Maurice Fréchet considered distribution
where is an arbitrary positive constant. The "2" comes from the fact that Fréchet considered errors with respect to the median. But he did not introduced that new distribution function, he also proved that the distribution appears as a limit when the underlying distribution of the ’s has an algebraic behavior at infinity, i.e. equivalent to , for some . I.e. he proved that Pareto-type tailed distibutions where in the domain of attraction of the Fréchet distribution.

•  Later on, the work of Emil Gumbel

In 1932, Emil Gumbel gave a talk in France on the "âge limite". But as he wrote it "on peut donc supposer que la distribution de l’âge limite - c’est à dire la probabilité que la probabilité de cet âge ait une valeur donnée - soit Gaussienne". But a few years after, he read about Fisher’s work, and observed also that "la distribution d’une valeur extrêmes peut être représentée pour un nombre suffisant d’observations par la formule doublement exponentielle, pourvu que la distribution initiale se comporte asymptotiquement comme une exponentielle. La formule devient rigoureuse si la distribution initiale est exponentielle", as he wrote in 1935. Thus, as Fréchet proved that Pareto type distribution were in the max-domain of attraction of Fréchet’s distribution, Gumbel obtained that exponential type distributions were in the max-domain of attraction of Gumbel’s distribution. He also introduced the term "distribution de type exponentiel"
For Emil Gumbel, it was natural to study the logarithmic derivative of the distribution, since it is the mortality rate in demography (area that Emil Gumbel studied previously). As he mentioned "d’un point de vue théorique, il est intéressant de noter que M. Fréchet a construit une distribution initiale d’’une variable aléatoire pour laquelle la valeur absolue de la dérivée logarithmique diminue sans limite". But since it was not a valuable property for practical applications, he decided that "nous nous bornerons au traitement des données de type exponentiel". Emil Gumbel always tried to relate his work on extremes and what he did on demograpy.
For instance in 1937, he wrote a paper on "les centennaires" that can also be related to the work of Bortkiewicz on rare events. He also applied his work on radioactivity, and hydrology.
In the 30’s, hydrographs as Hazen or Graszberger introduced the concept of "yearly maximum" of a river level. They actually proposed to look for actuarial models to study decennial or centennial floods.  But they only used the lognormal distribution to model yearly maxima. In 1936, French hydrologist Aimé Coutagne met Emil Gumbel (who was teaching at the ISFA, in Lyon). At that time, Emil Gumbel was looking for possible applications (outside demography) for his doubly exponential distribution. As as pointed out by Aimé, "sa formule devait être applicable au cas des crues; c’est à dire des plus grands débits, problème analogue à celui des plus grands âges". Not only Gumbel’s distribution gave better empirical results, but also it came with a theoritical justification.

• Gumbel’s distribution properties

Consider the Gumbel distribution, with location and scale parameters \alpha and \beta respectively, i.e.
Note that the associated quantile function is
with mean
and variance

• The work of Waloddi Weibull

Waloddi Weibull, a Swedish physict proposed a distribution in 1939, to represent the distribution of breaking strength of materials. He used it in the 50’s in reliability concept. Actually, Weibull appeared late in the story of extremes, since Fréchet, Fisher and Tippett mentioned it already in the mid-20’s.

• From the central limit theorem (on the average) to Fisher-Tippett theorem (on the maxima)

In order to visualize those two theorem, consider the following animation, where samples of 20 exponential variables are generated. From those 20 values, we plot the maximum in blue, and the average in red, on top. Just below, be rescale those points by considering , and below again, }. When then look at the position of  and the one of the mean of . We then build an histogram to visualize the distribution of the rescaled maximum (in blue) and the rescale average (in red).

For those who might be busy, after 1000 generations of samples, we obtain the following histograms (below), including the Gaussian distribution below (i.e. the average of exponential variables looks Gaussian, even with only 20 observations, actually the Gaussian distribution is only asymptotic, i.e. we should consider samples of size 2000), and the maximum over 20 observations of exponential variables (on top) looks like a Gumbel distribution (actually, here it is the exact distribution, and it is the asymptotic distribution for exponential type variables).

• The GEV distribution

The unified expression of those three distributions is call the GEV distribution. The generalized extreme value distribution has cumulative distribution function
for , where  is the location parameter,  the scale parameter and  the shape parameter. Note that the expected value is

Avec une petite semaine d’avance, un billet rapide pour préparer la Saint Valentin, i.e. récupérer des conseils d’experts. Il existe un site http://www.okcupid.com/ de Free Online Dating pour reprendre le titre. Mais son grand intérêt, c’est surtout le blog associé, qui propose des analyses statistiques sur les matching, http://blog.okcupid.com/.

Par exemple, les fondateurs (qui se sont amusés à faire beaucoup de statistiques) regardent les photos, et notent que pour espérer obtenir beaucoup de messages (ou tout du moins plus que la moyenne), les femmes doivent regarder l’objectif tout en souriant, en adoptant une moue un peu dragueuse ("flirty face") (cf Myth 1 ici). Les hommes, eux, doivent plutôt tirer la tronche et regarder ailleurs.... Mais le plus amusant, c’est - je m’adresse aux hommes qui lisent mon blog - il y a une étude sur ce qui doit apparaître sur la photo,

Bon, comme son nom l’indique, c’est un "blog perso" alors aujourd’hui, je vais parler de mon fils. Le 31 décembre est un jour un peu particulier, car c’est son anniversaire... Il est né avec plusieurs jours d’avance, et depuis qu’il va à l’école, on se dit que les choses auraient sûrement été plus simples s’il était né 3 heures plus tard... Bref, je me demandais si ce "retard" allait se combler un jour, ou bien s’il en "souffrirait" toujours¹ ...

Pour lui faire plaisir, j’ai voulu savoir si l’avenir dont il rêvait était compromis. Pour l’instant il rêve d’être footballeur (pour faire comme les copains, car entre nous, je le crois aussi peu doué que son père).

Loin de moi l’idée de lui ruiner totalement ses rêves, mais j’ai regardé un peu les joueurs de foot (tous les joueurs qui ont participé à la coupe du monde de 2006, i.e. ici, ou plutôt pour être honnête, tous ceux qui viennent d’un pays de l’hémisphère nord). Si on regarde le mois de naissance, on note qu’effectivement, peu de grands joueurs (au sens ’retenus dans leur équipe nationale, même s’ils ont passé toute la coupe du monde sur le banc de touche’) sont de la fin de l’année...

Bon, sans faire trop de psychologie de café du commerce, cela s’explique assez simplement: pour les jeunes enfants, disons de l’âge de mon fils, il y a un grand écart entre ceux qui sont nés en janvier et ceux qui sont né en décembre de la même année, et pourtant, ils sont dans la même catégorie sportive. Bref, dans une sélection (je parle des vrais clubs qui élèvent de la graîne de champion), on préfère les enfants du début d’année....
En regardant le jour et le mois de naissance, et en faisant une régression lowess, ainsi qu’une régression spline, on observe encore plus finement la tendance,

(pour être là aussi honnête, j’ai multiplié par 4 le nombre de naissances observé un 29 février). Bref, si mon fils était né trois heures plus tard, il aurait eu 60% de chances en plus de devenir un grand footballeur..! Bon, qu’à cela ne tienne, je vais lui trouver une autre vocation... En trainant toujours sur internet, on peut trouver des données similaires pour les basketteurs américains (ici). J’ai donc pris les joueurs NBA et ABA nés depuis 1946. Cette fois-ci, la tendance mensuelle semble disparaître...

et de même pour les tendances en regardant jour par jour,

Damned, il ne nous reste plus qu’à émigrer outre-atlantique et de convaincre mon fils qu’il a toutes ses chances au basket !

¹ Il existe un certain nombre d’études sur le sujet, avec ici ou là par exemple des liens entre le mois de naissance et la réussite scolaire.