All posts by Arthur Charpentier

Arthur Charpentier, professor in Montréal, in Actuarial Science. Former professor-assistant at ENSAE Paristech, associate professor at Ecole Polytechnique and assistant professor in Economics at Université de Rennes 1.  Graduated from ENSAE, Master in Mathematical Economics (Paris Dauphine), PhD in Mathematics (KU Leuven), and Fellow of the French Institute of Actuaries.

De la vitesse de convergence, partie 1

Samedi après midi, à l’heure de la sieste pour certains, nous discutions avec Magali, Jean Baptiste et Florent des machines qui pouvaient trier les cartes, de là, la discussion a dévié sur les algorithmes de tri, et sur la difficulté d’expliquer l’importance de la vitesse d’un algorithme, ou plus généralement, la vitesse de convergence. Bref, avant de parler de vitesse de convergence en statistique – ou peut être avant sur les méthodes de calcul d’intégrale (ou d’espérance) – je voulais faire un petit billet rapide sur les algorithmes de tri. Je commence par des choses simples. De mémoire c’est des choses que je me souviens avoir vu au lycée dans les années 90, dans les cours d’informatique… Mais l’exemple peut être intéressant pour comprendre la différence, par exemple, entre https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-02.pnghttps://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-03.png et https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-04.png.
Le premier algorithme est appelé tri à bulles. C’est probablement l’algorithme le plus simple, et sa vitesse est en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-01.png. On va balayer autant de fois la base que l’on a d’observations, et à chaque passage, on récupère le maximum (ou le minimum, peu importe). L’algorithme ressemble à ça,

> n=25; s=1
> X=sample(1:n,size=n,replace=FALSE)
>          for(i in 1:(n-1)){
+          for(j in 1:(n-1)){
+          s=s+1
+          if(X[j] > X[j + 1]){X[c(j,j+1)]=X[c(j+1,j)]}
+          }}

Si on l’implémente, voilà ce que ça donne (tous les graphiques partent de la même permutation initiale, et tous contiennent autant d’images, de l’ordre de https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-01.png)

Bref, un algorithme en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-05.png.Bon, j’ai été un peu brutal dans mon algorithme. En fait, on voit que les derniers points sont rapidement bien classés, donc inutile que la seconde boucle aille jusqu’au bout… on peut donc faire du https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-06.png, en faisant

 

+          for(j in 1:(n-i)){

Bref, la constante dans les https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-05.png peut être relativement importante comme on le voit ici.

Mais on peut faire  mieux, voire beaucoup mieux comme on le verra ensuite.
Le second algorithme est le tri par insertion. C’est aussi un algorithme très simple:  on balaye l’ensemble en plaçant les valeurs rencontrées à leur place, un peu comme quand on tri ses cartes… L’algorithme est là aussi en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-05.png. Ici l’algorithme converge un peu plus rapidement que le précédant, mais cela vient du fait que les données étaient partiellement ordonnées sur la fin. Autrement dit, un autre échantillon aurait pu donner autre chose… Mais généralement, c’est le cas.

>          for(i in 1:n){
+           x = X[i]
+          for(j in (i+1):2){
+          s=s+1
+         if(X[j - 1] < x){X[j] = X[j - 1];  j = j - 1;  X[j] = x}
+         }
+         }

(le tri est fait à l’envers ici…. mais c’est un détail).

Enfin, le tri fusion est un autre algorithme en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-07.png. L’idée est de couper, de trier, puis de fusionner. L’idée étant que trier un tableau dont la première moitié des éléments sont ordonnés, et où la dernière moitié est ordonné doit être plus rapide que quand tout est dans le désordre. A chaque étape, l’algorithme est simple, on compare le minimum des deux listes, puis on le vire à l’étape suivante. Sous R, ça donne, quand on commence

 

> for(i in 0:(trunc(n/2)-1)){
+ Y=X[c(1,2)+2*i]
+ if(Y[1]>Y[2]){X[c(1,2)+2*i]=X[c(2,1)+2*i]}
+ }

Bref, en commençant par classer les paires (1,2) puis (3,4), etc, puis les listes (1,2) versus (3,4). Le code est alors de la forme

> k=4
> for(i in 0:(trunc(n/4)-1)){
+ Y=X[1:k+k*i]
+ Y1=c(Y[1:(k/2)],rep(999999,k/2));Y2=c(Y[(k/2+1):k],rep(999999,k/2));
+ for(j in 1:k+k*i){
+     X[j]=min(c(Y1[1],Y2[1]),na.rm=TRUE)
+    if(X[j]==Y1[1]){Y1=Y1[
-1]}
+    if(X[j]==Y2[1]){Y2=Y2[-1]}
+ }}

(avec une petite correction pour le dernier bloc, qui ne compare pas forcément un bloc de taille 2 contre un autre de taille 2 (mais l’adaptation du code est infime). Puis on continue avec (1,2,3,4) versus (5,6,7,8), c’est relativement rapide de tout classer, comme le montre l’animation ci-dessous.

Il existe de très nombreux autres algorithmes de tri, le but étant de monter qu’un algorithme en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-07.png va effectivement beaucoup plus vite qu’un algorithme en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-05.png. J’ai tout fait ici en R, mais pour les adeptes d’autres langages, il existe plein de code en ligne, ici par exemple. Pour conclure, le graphique ci-dessous montre les différences entre les vitesses de convergence en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-08.png (en noir), en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-09.png (en rouge) et en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tri-10..png (en bleu),

La suite sur la convergence des méthodes de monte carlo, puis des théorèmes limites en stats….

Foot, probabilité, et assurance, partie 1

Depuis la fenêtre de mon bureau, je pourrais presque voir la grande banderolle sur le fronton du centre commercial de la Visitation,

Un vendeur d’électroménager proposer de rembourser certains achats sir la France gagne la coupe du monde de foot ! Cela me fait penser au court billet que j’avais fait l’an dernier, ici, comme quoi certains camping pouvaient rembourser les vacances des vacanciers en cas de mauvais temps…. Oui, on pourrait s’en douter, il y a un contrat d’assurance derrière – afin d’éviter que ce revendeur (il n’est pas le seul d’ailleurs) ne fasse faillite ! Des amis qui suivent le foot me faisaient remarqué que la prime d’assurance ne devaient pas être trop élevée….
Pour être un peu plus précis, Berkshire Hathaway a annoncé avoir assuré un de ses souscripteurs en cas de victoire de l’équipe de France lors de la prochaine coupe du monde. La perte maximale ne semble pas collossale, comme le notait Warren Duffet: « je pense que nous allons perdre 30 millions de dollars, ou quelque chose comme ça, si la France gagne ».

Quant aux probabilités dans le cas des vacances (ici), Meteo France propose des estimations objectives des probabilités. Dans le cas du foot, on peut faire comme en finance, et aller extraire les probabilités des marchés, en l’occurence des sites de pari en ligne, comme ici. En l’occurence, on apprend que la cote de la France est de l’ordre de 21… mais je reviendrais là dessus dans un autre billet !

Discussions géométriques du matin (partie 2)

Suite rapide de mon billet sur la construction de triangles avec des bouts de bois (ici). Adrian me faisait remarquer () qu’il y avait plus simple pour calculer la probabilité de construire un triangle (je suis d’ailleurs affreusement jaloux de ne pas y avoir pensé tout seul). Mon principal soucis que j’ai en lisant le commentaire est “ça semble juste, mais est-ce que c’est juste ?” (mon expérience en TD de proba est que l’on peut souvent trouver un paquet de solutions qui “semblent” justes). Bref, faisons preuve de rigueur quelques minutes (même si France Inter continue à nous dire que le mot est toujours tabou ces temps-ci).
Rappelons que si on dispose d’un échantillon L'image “https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-01.png” ne peut être affichée car elle contient des erreurs. i.i.d. de loi https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-02.png, et de densité https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-03.png, et que l’on définie les statistiques d’ordre

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-04.png

la loi du couple https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-05.png, avec https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-06.png,est alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-07.png

à condition que https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-08.png (sinon la densité sera nulle). C’est jsute des histoires de compter les classements possibles, et puis de calculer quelques probas simples. Un cas particulier connu de tous est la loi du couple https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-09.png, i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-10.png

Dans le cas qui nous intéresse, supposons que l’on tire deux bâtons, alors la loi du couple s’écrit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-11.png

alors que la densité du maximum était

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-12.png

donc la loi du minimum conditionnellement au maximum doit d’écrire, par la formule de Bayes

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-13.png

qui est effectivement la densité de la loi uniforme sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-14.png. Bref, ça marche pour deux bâtons.
Dans le cas de trois bâtons, la loi du triplet s’écrit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-15.png

alors que la loi du maximum devient

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-30.png

aussi la loi conditionnelle du couple des deux premiers s’écrit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-31.png

Or si l’on avait deux lois uniformes sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-16.png, compte tenu de ce que l’on a écrit auparavant, la loi du couple des deux statistiques d’ordre seraient

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-17.png

soit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-18.png

Bref, il semble qu’effectivement, le résultat soit correct: conditionnellement au maximum, les autres variables sont indépendantes, et suivent un loi uniforme sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-19.png.
Moralité, effectivement, les non-maximum sont indépendants, et leur loi est uniforme sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-16.png. On peut alors sommer ces deux lois uniformes et indépendantes, ce qui donne une loi triangulaire, sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-20.png, avec le sommet en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-21.png. Et par des raisons de symmétrie, la probabilité de dépasser https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-21.png est alors https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-22.png.Merci Adrian. C’est effectivement plus élégant que mon calcul moche…

De là à en faire un papier…

mais mon blog est l’endroit rêvé pour exposer mes interrogations récurentes. Je ne sais pas si vous avez remarqué, mais au toilettes, au bureau, il y a (généralement) du papier toilette dans un distributeur, et souvent, un rouleau qui traîne à côté….  et régulièrement, je me demande quelle peut être la stratégie optimale: prendre la papier du distributeur ou le papier du rouleau qui traîne…? car mine de rien, la personne chargée du réapprovisionement ne semble intervenir que si l’un des deux est vide.. autrement dit, si les deux diminuent trop vite ensemble, on risque d’avoir un problème un jour. Bon, ce qui me rassure, c’est que Donald Knuth avait remarqué des choses similaires, et avait eu la prétention de publier ses observations dans le American Mathematical Monthly. En fait, pour formaliser le problème, on note https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ02.png la probabilité que les gens utilisent du papier dans le rouleau le plus plein (on suppose connaître la quantité restant dans le distributeur qui dispose d’une fenêtre transparente), et https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ04.png la probabilité de piocher dans le plus petit rouleau. Et quand les deux font la même taille, on choisit au hasard (avec probabilités 1/2). On suppose que les rouleaux sont valables pour https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ60.png utilisateurs (je n’introduis pas de d’aléa, voire de dépendance en fonction du menu à la cantine). On parlera de portions. En gardant les notation de Knuth, on note https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ01.png l’espérance du nombre de portions restant sur le rouleau restant dès que l’un est terminé. Le but du papier de Knuth était d’étudier le comportement de https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ01.png quand https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ61.png. Mais comme ça n’a que peu d’intérêt, je trouve, on va regarder ce qui se passe à https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ60.png fixé.
Knuth note que https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ05.png et que https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ06.png. En fait, histoire d’épater les filles, Donald prétend que on peut obtenir facilement les résultats suivants

Pour le premier, c’est trivial, et pour le second, on peut faire un dessin (oui, je suis très dessin en ce moment),

A partir de la diagonale (points verts), il y a égalité entre les paquets, et donc on choisit avec probabilité 1/2. Pour calculer https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ62.png, on note qu’il y a deux issues possibles, soit il reste 1 portion sur le rouleau non vide (les points rouges) soit il reste 2 portions sur le rouleau non vide (les points bleus). Les probabilités de passer d’un point à l’autre sont indiquées sur le dessin. Bref, on a https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ03.png chances d’atteindre un point bleu, et https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ02.png chances d’atteindre un point rouge, dont l’espérance vaut ici

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ07.png

En continuant un peu, on obtient que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ08.png

En effet, sur notre petit dession, on rajoute juste les points autours, les valeurs à l’intérieur restant inchangées, et pour faire simple, je m’éloigne de la diagonale avec probabilité https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ03.png, et je m’en rapproche avec probabilité https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ02.png.

Comme tout à l’heure, on compte les probas d’arriver sur les différents points. On a https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ09.png chances d’atteindre un point violet, https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ10.png chances d’atteindre un point bleu, et enfin https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ11.png chances d’atteindre un point rouge (il suffit d’étudier toutes les trajectoires possibles). Au total, on arrive à

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ12.png

qui se simplifie en

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ13.png

Quant aux conditions de bord, on notera que https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ14.png et que https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ1.png. En fait, on peut étudier un problème plus général, mais plus réaliste: les paquets dans les distributeurs sont souvent beaucoup plus gros que celui qui traîne. On peut donc noter https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ16.png la même quantité, avec m le nombre de portions dans le distributeur et https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ60.png dans le rouleau qui traîne. Si on a compris le dessin ci-dessus, on se doute que ça ne change pas grande, chose, ormis que l’on ne commence pas sur la diagonale…
En fait, on note que https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ17.png, avec les conditions de bord de la forme https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ18.png.
On peut aussi noter (il suffit de regarde sur le dessin en haut) que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ19.png

La preuve plus propre est que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ20.png

or comme https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ21.png par une symmétrie du problème, on a l’égalité. En utilisant l’idée de la preuve du résultat précédant (et le dessin, comme toujours), si https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ22.png, alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ23.png

Bref, on est capable de faire tous les calculs de manière récursive et relativement simple (pour ceux qui compte le faire la prochaine fois qu’ils vont aux toilettes, il faudra peut être emmener quand même leur ordinateur). On note que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ24.png

or https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ25.png (comme on l’a montré) et

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ26.png

On a tout cette fois ci pour vérifier le calcul “facile” évoqué par Donald,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ27.png

Bref, on retombe sur le résultat précédant. On peut alors faire le calcul en fonction des valeurs de https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ02.png (que l’on fixera, puis de m et https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ60.png),

> p=0.4
> q=1-p
> m=50
> n=50
> M=matrix(NA,m+1,n+1)
> M[,1]=seq(0,m)
> M[1,]=seq(0,n)
> M[2,2]=2-p
> for(i in 3:(m+1)){
+ for(j in 2:(i-1)){
+ M[i,j]=p*M[i-1,j]+q*M[i,j-1]
+ }
+ M[i,i]=M[i,i-1]
+ }
> for(i in 2:(m)){
+ for(j in (i+1):(n+1)){
+ M[i,j]=M[j,i]
+ }}
> colnames(M)=seq(0,n)
> rownames(M)=seq(0,m)
> image(seq(0,n),seq(0,n),M,)
> contour(seq(0,n),seq(0,n),M,add=TRUE,
+ levels=seq(0,50,by=5))

Voilà ce que l’on obtient si https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ28.png. le nombre de portions qui reste dès que l’un des rouleaux est fini peut se visualiser sur le graphique suivant, la couleur rouge correspond aux cas un peu dangereux où il reste peu de feuilles dans le rouleau, les zones jaunes indiquant au contraire des niveaux plus élevés,

Avec https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ29.png, on obtient

Autrement dit, il reste toujours du papier sur l’autre dès lors que les gens pensent à alterner suffisement (ici 95% des utilisateurs) alors que pour https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ30.png,

autrement dit, on risque constament de se retrouver à court… Ce qui est marrant (?) c’est le cas où https://perso.univ-rennes1.fr/arthur.charpentier/latex/PQ31.png: s’il y autant de personnes qui piochent dans le gros paquet ou dans le petit, alors il y a de fortes chances pour que les deux rouleaux soient vides dès que l’un l’est…

[ce billet a été initialement posté sur freakonometrics.blog.free.fr/… en mai 2010]

Moyenne géométrique versus moyenne arithmétique…

L’autre jour, j’étais de surveillance d’examen. Pour passer le temps, j’ai survolé l’énoncé, et je me souviens avoir lu “montrez que la moyenne arithmétique est plus grande que la moyenne géométrique“. Ayant généralement tendance à préférer les copies où il y a peu à lire, mais avec de jolis dessins, j’ai eu envie de proposer une correction sans écrire un mot aux étudiants (qui semblaient souffrir dans leur démonstration). Il suffit de monter que

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/moyennes-latex.png

Ben voilà,

quod erat demonstrandum. Non ? On note en plus que l’on a égalité si et seulement si les nombres sont égaux…. Allez, pour convaincre les sceptiques, je rajoute une deuxième preuve… Bon, en fait je peux prouver que

mais c’est la même chose, non ?

Discussions géométriques du matin (partie 1)

Petit résumé du débat du matin avec mon fils …
– papa, papa j’ai trouvé trois bouts de bois,
– ben… fais une figure géométrique avec *
– ça marche pas !
– comment ça, “ça marche pas” ?!?
– j’arrive pas à faire un triangle !
– … hein ? attends, j’arrive
genre “mais qu’est-ce que tu racontes ? t’es pas réveillé ce matin ou quoi ?“. Bref, après avoir vu qu’effectivement ça n’était pas possible (un des bouts de bois était beaucoup plus long que les deux autres),  je me suis demandé si mon fils n’avait vraiment pas eu de chance… ou si au contraire, c’était plutôt “normal” de ne pas pouvoir faire un triangle à partir de trois bouts de bois….
Soient https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-01.pnghttps://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-02.png et https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-03.png les longueurs des trois morceaux de bois (oui, une fois le café avalé, je peux prendre un peu de temps pour formaliser). Ces trois bouts de bois forment  un triangle si et seulement si la somme des deux plus petits est plus grande que le plus grand. Autrement dit, si je m’intéresse aux statistiques d’ordre

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-04.png

il faut que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-05.png

Pour faire simple, on supposera que les https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-50.png suivent une loi uniforme sur [0,1].

  • Approche par simulation

Par simulation, ça ne devrait pas être trop compliqué d’avoir un ordre de grandeur de la probabilité

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-06.png

Le code est tout simplement

> n=1000000
> U=matrix(runif(3*n),n,3)
> U1=apply(U,1,max)
> U2=apply(U,1,sum)-U1
> U3=U1<U2
> mean(U3)
[1] 0.501067

Finalement la probabilité semble être proche de 1/2… ça doit pouvoir se montrer….

  • Approche par des calculs de proba

Il doit être possible de calculer directement

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-06.png

ou bien, on peut aussi noter que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-07.png

ce qui ne fait intervenir que deux variables (certes, toujours non indépendantes, mais on devrait s’en sortir malgré tout): le maximum des trois variables, et la somme des trois variables. On cherche alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-08.png

Je pense qu’il est plus simple de travailler sur la première expression, en notant qu’il existe 6 permutations de ces 3 valeurs, et cette probabilité s’écrit simplement

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-10.png

(qui  est un des 6 arrangements possibles, tous équiprobables). Cette probabilité est alors tout simplement, compte tenu de l’indépendance des tailles des bouts de bois, ce qui est peut être une hypothèse un peu forte, certes,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-11.png

qui doit pouvoir se réécrire

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-12.png

L’intégrale au centre se réécrit alors simplement, à condition de savoir comment https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-13.png se positionne par rapport à 1. Deux cas assez différents sont à distinguer: comme on veut https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-134.png, soit https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-15.png, soit https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-16.png. Autrement dit on écrira l’ensemble sur lequel on calcule la première intégrale comme la réunion de deux ensembles (disjoints),

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-18.png

et

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-19.png

Sur le premier ensemble,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-20.png

et sur le second

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-21.pngBref, l’intégrale se retrouve sous la forme de deux intégrales

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-22.png

et

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-23.png

La première intégrale correspond à un cas où https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-24.png est petit, ou plutôt https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-25.png. La seconde intégrale est un cas où https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-26.png est grand, ou plutôt https://perso.univ-rennes1.fr/arthur.charpentier/latex/b<br />
ois-27.png. Bref, la première intégrale peut s’écrire

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-28.png

et la seconde

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-29.png

(on aura en mémoire le facteur 6, que je ne mes pas pour alléger un peu). Bref, on a deux fois la même intégrale, qui s’écrire (dans le second cas à une translation des indices près)

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-30.png

Notre probabilité est alors 6 fois le double de cette intégrale, soit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-31.png

Bon, finalement c’est beaucoup plus improbable de pouvoir construire un triangle que je ne le croyais au départ…. Damned ! Moralité, mieux vaut ne pas trop réfléchir avant d’avoir pris un café…. Et comme “papa” est très fort en mauvaise foi, “papa” expliquera que si la longueur des bois de bois n’était pas distribuée suivant une loi uniforme, l’intuition initiale n’était pas complètement stupide… Non mais !
* pour ma défense, c’était ça ou la révision de la table de multiplication de 5… en gros, “papa a pas encore pris son café, c’est pas le moment“.

Histoire éthylique, et loi de l’arcsinus

Avant de revenir sur la marche aléatoire et le problème de la ruine d’un assureur comme problème de temps d’arrêt (je n’ai pas insisté en cours car les temps d’arrêt nécessitent un peu de recul sur les processus), je voulais prendre 5 minutes pour faire un exemple qui devrait intéresser les étudiants. On va se placer un peu avant les billets sur l’étude de l’ébriété rennaise, ici et .

  • De la probabilité de sortir un jour de la rue de la soif

Supposons que notre homme soit encore rue de la soif, un jeudi soit à 23 heures, dans un état d’ébriété assez avancé. Souhaitant quitter cette rue, il fait un pas (disons 1 mètre pour faire simple), puis hésite sur sa direction: il a 1 chance sur 2 de continuer dans la même direction, et 1 chance sur 2 de faire demi-tour.  Il se trouve à une distance L'image “https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-05.png” ne peut être affichée car elle contient des erreurs. de la place Sainte Anne et à une distance L'image “https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-06.png” ne peut être affichée car elle contient des erreurs.de la place des Lices. On se demande s’il arrivera un jour à quitter la rue de la soif, et surtout au bout de combien de temps (ou après avoir parcouru quel distance)… Au delà du temps moyen passé à errer dans la rue, on pourra se demander la probabilité pour qu’il soit sorti de la rue pour 8 heures, histoire d’avoir une chance d’arriver à l’heure en cours.

Il s’agit du problème de la marche aléatoire (qui, ma foi, est assez approprié au problème de l’ébriété) de type pile ou face. On note tau l’instant (exprimé en nombre de pas) de la sortie de la rue de la soif. Formellement, on note https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-01.png la position à la date https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-02.png. On suppose qu’à la date https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-02.png=0, il est à la position 0 (disons sa chaise à la terrasse du 1929). tau est l’instant où https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-01.png quitte l’intervalle https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-03.png. S’il sort par la partie inférieure (en https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-05.png) il sera place Sainte Anne, alors que s’il sort en https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-06.png, il sort par la place de Lices (la longueur de la rue étant https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-04.png).
Notons que https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-07.png où https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-08.png prend les valeurs https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-09.png avec probabilités https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-10.pnghttps://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-50.png s’il se dirige vers Sainte Anne, https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-11.png s’il part place des Lices.
On cherche

https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-12.png

C’est ce qu’on appelle un temps d’arrêt. Il s’agit d’une variable aléatoire, positive, qui prend valeurs dans https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-555.png.
Le problème est simple à formuler, mais pour le résoudre, c’est une autre histoire…
On peut d’abord montrer que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-13.png

ce qui est un peu rassurant (c’est le lemme de Fatou). On peut alors montrer que https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-14.png appartient forcément à https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-15.png, et donc

https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-16.png

aussi,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-17.png

ce qui est assez intuitif compte tenu de la symétrie naturel du problème. On a alors d’autant plus de sortir d’un coté que l’on part d’un endroit proche de ce côté.
On peut aussi monter que https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-18.pngest une martingale, nulle à la date initiale. Aussi,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-19.png

(je passe les difficultés techniques, mais on peut montrer que ce résultat est rigoureusement juste en travaillant sur la martingale arrêtée). Bon, on a pas mal d’information…. mais on peut aussi regarder la loi. On peut faire quelques simulations, dans un premier temps, afin de visualiser la densité (sous forme d’histogramme), avec ci-dessous la position de notre individu en ordonné, et en abscisse le temps depuis qu’il a quitté sa table.

https://perso.univ-rennes1.fr/arthur.charpentier/image-st-michel-v2.gif

L’histogramme ci-dessous
est le comptage des moments où l’individu a quitté la rue.

https://perso.univ-rennes1.fr/arthur.charpentier/compte-image-st-michel-v2.gif

Je me suis arrêté assez tôt sinon le fichier image serait trop lourd.
Pour étudier la distribution, mieux vaut éviter de travailler sur la fonction de répartition ou la densité du temps d’arrêt. En fait, le plus simple est de récupérer la fonction génératrice de tau. Pour cela, on va utiliser la fonction génératrice des incréments,
https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-20.png

en utilisant la notation du cosinus hyperbolique. On peut alors montrer

https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-21.png

Pour tout https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-22.png,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-23.png

En fait, on retrouve ici quelque chose qui est très proche de la probabilité qu’a le mouvement brownien de quitter un corridor (qui ont été beaucoup étudiée de part les applications en finance, non pas sur l’étude des gens bourrés, via les options parisiennes, ici ou , mais qui avait été étudié dans les revues de physique dans les années 60). Bon, comme on quitte le monde discret pour le monde continu, on ne travaille plus sur la fonction génératrice des moments https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-24.png, mais sur la transformée de Laplace https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-25.png dont on peut montrer que la transformée de Laplace du temps de sortie d’un corridor de la forme https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-521.pngest

https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-30.png

qui se rapproche de la formule précédente. Pour faire le lien, notons que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-31.png

et que pour https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-32.png,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-33.png

Bref, en bricolant, on peut relier assez clairement les deux notions (ce qui n’a rien d’une démonstration, j’en conviens). Mais on peut aussi se poser d’autres questions fondamentales……

  • De la probabilité de récupérer ses clés

J’avais évoqué ici les problèmes de clés de notre indigène. Mais j’avais oublié de dire qu’il avait oublié ses clés sur sa table en quittant le bar. On peut se demander quelle est la chance pour qu’il ne repasse jamais à sa table…. Autrement dit, on souhaite calculer

https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-522.png

ou éventuellement, si on sait que le bar ferme à l’heure https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-523.png,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-524.png

Dans le cas du mouvement Brownien, on peut montrer que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-525.png

On utilise la propriété de Markov, qui garantie que https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-526.png est indépendant de https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-527.png. On peut alors montrer que si https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-528.png est le temps d’arrêt défini comme le premier instant où le processus atteint https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-529.png, alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-530.png

Dans le cas du mouvement Brownien, cette propriété est une des lois de l’arcsinus, dont on trouvera des démonstrations ici ou . L’arcsinus apparaît ici dans un calcul d’intégral, en notant que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/alcool-531.png

Compte tenu des similitudes que nous avions observées dans la section précédente entre le cas continu et le cas discret, notre étudiant devrait pouvoir calculer assez simplement la probabilité qu’il a de récupérer ses clés…. enfin, dès qu’il aura dessaoulé un peu afin de faire les calculs….

Faut-il se fâcher avec tout le monde ?

(s’il le faut allons-y) Samedi matin, j’ai aperçu avec un grand plaisir le dernier numéro de Fakir chez mon buraliste préféré (i.e. le “journal fâché avec tout le monde“). Je n’ai pas vocation à me fâcher avec tout le monde, loin de là, mais je dois avoir un côté pitbull que j’ai du mal à cacher…

Et si quelques pages ont rapidement retenu mon attention, je dois avouer que j’en suis sorti très déçu (par ces quelques pages). Bon, autant le dire tout de suite, les pages en questions ne sont pas écrite par les journalistes de Fakir, mais correspondent à un “courrier des lecteur” (les “meilleurs” paraît-il… je me sens flaté, je suis un lecteur régulier). Le “papier” de deux pages part d’une attention louable (me semble-t-il), à savoir essayer de repérer des classes (au sens marxiste) en étudiant la distribution des revenus. On nous explique que l’auteur est “directeur de recherche au CNRS“, ce qui vise à donner a priori une légitimité scientifique à ses propos. N’ayant pas encore l’accord pour mettre en ligne les deux pages en question (je vais quand même demander l’accord à la rédaction), je me contente de ces deux copies de mauvaise qualité.

Pour résumer, l’auteur nous explique que l’INSEE ment, ou triche, dans sa façon de produire ses graphiques (comme je l’avais dit ici, je rejoins l’auteur sur ce point: j’avais déjà dit que la façon de présenter des informations chiffrées n’était pas neutre sur le message que l’on souhaite diffuser).

  • Densité ou fonction quantile ?

L’auteur commence par nous présenter deux graphiques. Le premier tiré de Max Weber montre la répartition de la taille de fourmis (géantes !), et observe que la distribution bimodale s’interprète comme la présence de deux classes de fourmis.

Ensuite, il nous présente la fonction quantile du revenu des français, et s’étonne de ne pas retrouver de distribution bimodale.

Il nous explique qu’avec ce genre de graphique, on ne trouverait pas les modes observés sur les tailles des fourmis,

plus précisément, l’auteur explique

Pour formaliser un peu (et revenir sur ce qu’est la “moulinette de l’insee“), rappelons le premier graphique est une densité https://perso.univ-rennes1.fr/arthur.charpentier/latex/fakir-01.png (je passe sur le fait que l’on travaille sur des échantillons, je vais faire une explication sur la modélisation sous-jacente, i.e. en travaillant sur la densité sous-jacente). Pour interpréter le second, rappelons que l’on construit la fonction de répartition en intégrant la densité,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fakir-02.png

La fonction quantile est alors l’inverse de cette fonction

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fakir-03.png

autrement dit, d’un point de vue géométrique, le symétrique par rapport à la première bissectrice de la fonction de répartition.  Avec des graphiques un peu plus clairs, je pense, voici une densité de mélange, en noir

Les deux populations sont visualisées avec les deux couleurs, rouge et bleu. Si je trace la fonction cumulative, j’obtiens

Malheureusement, on ne dispose pas de données précises sur les revenus (et je suis le premier à m’en désoler), et on observe une version discrète de cette fonction, à savoir les niveaux des déciles,

En inversant cette focntion, on obtient alors la fonction quantile, correspondant au graphique publié par l’insee,

Ma première observation est que si https://perso.univ-rennes1.fr/arthur.charpentier/latex/fakir-01.png est effectivement bimodale, cela signifie que la dérivée de https://perso.univ-rennes1.fr/arthur.charpentier/latex/fakir-04.png possède deux modes, ou plutôt deux maximum locaux. Si on se souvient des trucs que l’on apprend au lycée (et non pas à un niveau “bac+5” comme l’indique le haut du document), si on fait le tableau de signe (la première ligne étant le signe de la dérivée de la densité) on obtient

https://perso.univ-rennes1.fr/arthur.charpentier/latex/fakir-tableau-signe.png

Autrement dit, on ne voit rien. Formellement, cela se traduira uniquement sur la pente de https://perso.univ-rennes1.fr/arthur.charpentier/latex/fakir-04.png qui diminuera après le premier mode, qui augmentera juste avant d’atteindre le deuxième, avant de diminuer à nouveau après le second. Sur les courbes en escaliers, c’est vrai qu’il est difficile de visualiser une distribution bimodale, i.e. des accroissements sur la “pente” (ce n’était déjà pas forcément simple avec l’ensemble des valeurs possibles, mais sur une fonction en escalier….). Mais quoi qu’il en soit, je ne suis pas sûr que cette recherche de la bimodalité soit fondamental dans cette étude….

  • Bimodalité et mélange

En fait, ma seconde critique est un peu plus technique. Pour faire simple, la présence de “classes” de traduit par une présence de “mélange” (c’est la base de la tarification où l’on essaye d’identifier les classes tarifa
ires, en assurance, comme j’en ai parlé à de nombreuses reprises ici ou , par exemple). Mais ce “mélange” n’a aucune raison de se traduire par une distribution bimodale. On retrouve ici un argument fallacieux que l’on retrouve couramment, y compris dans des revues scientifiques, par exemple dans ce papier sur l’âge des tentatives de suicides, publié dans Journal of Psychiatric Research (ici). On y retrouve une tentative de retrouver des mélanges,

(je reviendrais un jour sur cette étude). Pour essayer de mieux comprendre en quoi l’argument est fallacieux, regardons quelques cas simples.
Le premier point est que la densité d’un mélange n’a aucune raison d’être bimodale. Même avec des lois sous-jacentes unimodales, comme la loi normale par exemple,

La courbe en noir est – comme auparavant – la distribution sur l’ensemble de la population. Elle est bimodale si les groupes sont bien discriminés, mais pas si les groupes sont proches. Comme je l’explique dans toutes mes interventions sur  les modèles linéaires généralisés, il est difficile de se faire une idée sur la loi par sous-groupe en regardant la distribution globale. Si la loi n’est pas unimodale (à en croire les travaux de Vilfredo Pareto – mentionnés ici ou – j’aurais tendance à croire que c’est le cas quand on parle des revenus), cette recherche des modes a encore moins de sens, me semble-t-il…. Le graphique ci-dessous correspond à mélanges de lois exponentielles (les traits en pointillés correspondant, comme auparavant, aux moyennes par sous groupe).

De plus, pour revenir aux exemples mentionnés dans le papier, à on avis, il y a une différence fondamentale entre le revenu et la taille: si la taille peut être supposé suivre une loi normale (par sous-population) je doute que ce soit le cas pour le revenu…Et encore, j’ai fait simple, car on suppose qu’il existe uniquement deux classes, à peu près équiréparties (ici 40-60 dans les deux animations précédentes).

Sur l’exemple ci-dessus, on suppose qu’il existe un petit groupe  hétérogène (i.e. avec une plus forte variance). Cette hypothèse ne me semble pas absurde si on raisonne sur les revenus. Bref, on ne trouve presque jamais de distribution bimodale… Et ce d’autant plus que le petit groupe est hétérogène. Si en plus on suppose qu’il existe non pas deux groupes, mais trois (voire plus), c’est encore pire: il devient difficile de visualiser les sous-groupes…

  • Faute de logique ou de statistique ?

J’ai l’impression que l’auteur part d’une affirmation qui me semble juste, à savoir “une distribution bimodale peut être modélisée par un mélange de deux distributions unimodales“. Sauf qu’il fait un raisonnement étrange en réécrivant cette affirmation sous la forme “une distribution qui n’est pas bimodale ne peut pas être un mélange” (je passe sur le premier point que j’avais mentionné comme quoi chercher un mélange sur une discrétisation de la fonction quantile, ou de la fonction de répartition, me semble étrange). Bref, je pense que l’ensemble de l’étude statistique repose sur des arguments fallacieux, ce qui a malheureusement tendance à décrébiliser l’ensemble de l’analyse. C’est dommage… Pourquoi ne pas avoir fait appel à un chercheur connaissant statistiques pour faire l’étude ? Mais beaucoup sont passés par l’INSEE, ce qui semblerait ôter tout légitimité à être critique (à en croire le petit encadré en bas)… Je crois qu’il serait temps d’arrêter de stigmatiser ainsi l’insee, beaucoup de gens y travaillant ayant gardé un esprit ouvert et critique ! Il serait dommage de ne pas s’appuyer sur leur expertise. L’autre regret est que la caution intellectuelle viennent d’un ” directeur de recherche du CNRS” (éventuellement en linguistique ou en anthropologie, on n’a malheureusement pas plus d’information).  Je trouve regrettable qu’un titre académique permette de se dédouaner de toute relecture critique… par des statisticiens….

Histoire éthylique, et problèmes de clés

Pour poursuivre le billet commencé ici, avant d’expliquer comment notre personne imbibée d’alcool est allée du bar à son appartement, je voulais raconter une petite mésaventure qu’il a eu.

Lors de ses multiples essais de clés, des agents de la maréchaussée ont trouvé étrange qu’il ait besoin d’autant d’essais pour trouver la bonne clé.

On part d’une distribution a priori d’une chance sur deux,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ivre-01.png

On sait que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ivre-02.png

alors que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ivre-03.png

et 0 sinon.
On peut alors utiliser la formule de Bayes pour déterminer la probabilité qu’une personne soit ivre sachant qu’il lui a fallu https://perso.univ-rennes1.fr/arthur.charpentier/latex/ivre-04.pngessais,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ivre-06.png

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ivre-07.png

Je laisse les plus courageux écrire la forme de cette probabilité, mais ça n’est pas sorcier….
Si on regarde ce qui se passe avec https://perso.univ-rennes1.fr/arthur.charpentier/latex/ivre-05.png=15 clés, par exemple, on obtient la probabilité suivante, en fonction du nombre d’essais nécessaire pour ouvrir la porte,

Pour la partie de droite, on a un résultat assez intuitif: s’il faut plus d’essais que le nombre de clés dans le trousseau, c’est qu’on est ivre. En revanche, tant qu’on n’atteint pas la valeur fatidique de https://perso.univ-rennes1.fr/arthur.charpentier/latex/ivre-08.png, alors plus le nombre de clés essayé est grand, plus cela “prouve” que l’on a de forte chances d’être à jeun. On retrouve ça sur le graphique ci-dessous,

Moralité, si on se fait arrêté par la police après avoir eu besoin de 12 essais pour trouver la bonne clé, ce petit laïus permet de prouver que l’on a peu de chances d’être ivre (malgré les apparences, j’en convient).
Mais pour l’instant, nous avons supposé que la maréchaussée intervenait lors de l’ouverture. Que se passe-t-il si elle intervient avantl’ouverture de la porte ? Autrement dit, on cherche ici

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ivre-09.png

dont l’évolution pour notre trousseau de https://perso.univ-rennes1.fr/arthur.charpentier/latex/ivre-05.png=15 clés donne

ou, en fonction de https://perso.univ-rennes1.fr/arthur.charpentier/latex/ivre-05.png et https://perso.univ-rennes1.fr/arthur.charpentier/latex/ivre-04.png,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso/desproges.jpgOn retrouve alors quelque chose qui pourrait paraître plus intuitif… Quoi que. Si au 12ème essai, la porte n’était pas ouverte, la police a tout intérêt à arrêter notre jeune héros, car la probabilité qu’il soit ivre est supérieure à la probabilité a priori. En revanche, s’il réussit à l’ouvrir au 13ème essai, elle a intérêt à ne pas l’arrêter, car il y a de très fortes chances pour qu’il soit à jeun. Étonnant non ?

Histoire éthylique, et de clés

En tant que statisticien, j’aime beaucoup les alcooliques… ils nous permettent de voir certains problèmes de modélisation sous un autre angle.
Considérons une personne ayant un trousseau de 10 clés, ou disons n pour généraliser. Il souhaite ouvrir la porte de chez lui, mais ne sait plus laquelle est la bonne. Combien d’essais – en moyenne – faudra-t-il avant de tomber sur la bonne clé ?
Pour une personne à jeun (et saine d’esprit), le plus simple est de tenter une clé au hasard, puis, si ça ne marche pas, de prendre celle à sa droite (ou à sa gauche, peut importe, mais il faut une stratégie et s’y tenir, à moins de pouvoir noter quelle clé a été utilisée). Bref, au bout de n essais, on est certain de tomber sur la bonne clé. Si le fait de manière un peu formelle (histoire de justifier une intuition), et que l’on note https://perso.univ-rennes1.fr/arthur.charpentier/latex/etyl-01.png l’évènement “la bonne clé est tirée au bout de https://perso.univ-rennes1.fr/arthur.charpentier/latex/etyl-02.png essais“, alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/etyl-03.png

Raisonnons alors simplement,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/etyl-04.png

puis

https://perso.univ-rennes1.fr/arthur.charpentier/latex/etyl-05.png

puis

https://perso.univ-rennes1.fr/arthur.charpentier/latex/etyl-06.png

etc. Aussi, on peut légitimement penser que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/etyl-07.png

Moralité (il s’agit de la moyenne d’une loi uniforme)

https://perso.univ-rennes1.fr/arthur.charpentier/latex/etyl-08.png

Pour une personne ivre (ou un peu moins à jeun), à chaque essai, elle prend une clé au hasard dans le trousseau (on pourrait appeler ça la stratégie du poisson rouge). Le nombre d’essais peut être potentiellement infini. Et en reprenant les notations précédentes, on note que si

https://perso.univ-rennes1.fr/arthur.charpentier/latex/etyl-09.png

en revanche

https://perso.univ-rennes1.fr/arthur.charpentier/latex/etyl-10.png

puisqu’à l’étape 2, on n’a pas écarté la clé qui ne marche pas… Aussi,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/etyl-11.png

qui correspond à une loi géométrique, ou loi de Pascal, de paramètre https://perso.univ-rennes1.fr/arthur.charpentier/latex/etyl-12.png. Pour ceux qui ne savent plus comment calculer l’espérance d’une telle loi, on peut prendre trois lignes pour refaire le calcul

https://perso.univ-rennes1.fr/arthur.charpentier/latex/etyl-13.png
soit tout simplement

https://perso.univ-rennes1.fr/arthur.charpentier/latex/etyl-14.png

i.e.
https://perso.univ-rennes1.fr/arthur.charpentier/latex/etyl-15.png

Aussi https://perso.univ-rennes1.fr/arthur.charpentier/latex/etyl-16.png.
Grosso modo, il faut deux fois plus de temps pour trouver une clé quand on est rond comme une queue d’pelle, que quand on est à jeun (sans compter le temps où l’on parle à son trousseau de clé, en espérant qu’il trouvera un moyen de désigner la bonne clé). La prochaine partie du billet parlera de “marche aléatoire“, et de difficultés pour quelqu’un d’ivre mort à quitter a rue Saint Michel….

Mais quand va-t-on enfin se décider à fermer les universités en France ?

J’ai été un peu titillé ce matin en voyant la une s’étaler sur la devanture du buraliste en bas de la fac,

Comme le dit Wikipedia, “l’expression grande école désigne en premier lieu les établissements d’enseignement supérieur français recrutant majoritairement sur concours parmi les élèves des classes préparatoires aux grandes écoles. Par extension, le terme est également associé à un petit nombre d’établissements de formation des hauts fonctionnaires français recrutant sur concours parmi les titulaires de la maîtrise ou d’un diplôme considéré comme equivalent. Enfin, de manière plus récente, le terme peut être associé à certains établissements d’enseignement supérieurs recrutant sur dossier ou concours directement, ou indirectement via des classes dites “préparatoires intégrées”, au niveau du baccalauréat et délivrant un diplôme après généralement cinq années d’études postérieures au baccalauréat“. Bref, au vu de cette dernière définition, j’ai du mal à voir clairement le lien entre la notion de “grande école” et celle d’élite, évoqué dans le premier paragraphe de de l’article, ou dans l’édito de la première page,

ÉLITE, subst. fém.

I.− [Avec ou sans déterminant] Ce qu’il y a de meilleur dans un ensemble composé d’êtres ou de choses; produit d’une élection qui, d’un ensemble d’êtres ou de choses, ne retient que les meilleurs sujets. A.− Usuel. [À propos d’êtres hum.] 1. Au sing. Minorité d’individus auxquels s’attache, dans une société donnée, à un moment donné, un prestige dû à des qualités naturelles (race, sang) ou à des qualités acquises (culture, mérites).
2. Au plur. avec art. déf. (néol. début xxe s.). Classe minoritaire composée de gens qui, du fait de leur naissance et de leurs mérites, de leur culture et de leur capacité sont reconnus (ou se reconnaissent) comme les plus aptes soit à occuper les premières places de la société à laquelle ils appartiennent, soit à donner le ton à leur milieu
Moi qui ai enseigné dans de nombreuses grandes écoles, et dans plusieurs universités, je suis un peu agacé par cette ambiance actuelle de dévalorisation de l’université, car pour être honnête, j’ai eu à l’université plusieurs très bons élèves (meilleurs que des bons élèves d’école), et réciproquement, j’ai vu des élèves en école que je n’aurais jamais pensé avoir une licence à l’université… Bref, je trouve dangereux de mentir ainsi aux parents en dénigrant autant l’université (car ces journaux sont plus lus par les parents que par les lycéens qui s’interrogent – éventuellement – sur leur avenir). Ce qui me gêne également, c’est lorsque les décideurs politiques lisent ce genre d’article (voire certains responsable universitaires) et décident d’affecter des moyens colossaux dans des “écoles” (au détriment de l’université, car comme toujours depuis quelques années, le jeu est à somme nulle). Si vraiment l’université était si mauvaise que ça, je me demande ce qu’on attend pour toutes les fermer….
Fort heureusement, les recruteurs sont – pour l’instant – un peu moins bornés que certains…. Espérons que ça dure…

De la qualité d’un score de classification

Un petit mot sur les courbes dites ROC, pour Receiver Operating Characteristic. Pour cela, on suppose que l’on dispose d’un prédicteur d’un variable prenant des valeurs 0 et 1 (pour simplifier), ou mieux encore “positif” et “négatif”. Peu importe le prédicteur, on peut considérer une régression logistique, une analyse discriminante, un classificateur nonparamétrique…. Bref, pour l’ensemble de nos observations, on a une valeur observée https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-01.png et une valeur prédite https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-02.png. En fait, comme je l’expliquais ici, on dispose plus précisément d’un score https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-03.png. La règle d’affectation est alors simple: on se fixe un seuil https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-04.png, et

  • si https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-05.png, alors  https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-02.png est “positif”
  • si https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-06.png, alors  https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-02.png est “négatif”

On peut alors construire une matrice dite de confusion, qui est simplement un table de contingence,

                valeur observée https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-01.png
valeur prédite
 https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-02.png
“positif” “négatif”
“positif” TP FP
“négatif” FN TN

où TP désigne les vrais positifs (true positive), TN  les vrais négatifs (true negative),FP désigne les faux positifs (false positive) ou erreur de type I (dans une terminologie de théorie de la décision, ou de théorie des tests), et FN désigne les faux négatifs (false negative) ou erreur de type II.
On peut alors définir toute une batterie d’indicateurs permettant de juger de la qualité de notre prédicteur (ou plutôt de notre score),

  • TPR = TP / P = TP / (TP + FN) appelé sensibilité, correspondant au taux de vrais positifs (true positive rate)
  • FPR = FP / N = FP / (FP + TN) correspondant au taux de faux positifs (false positive rate)
  • ACC = (TP + TN) / (P + N) appelé précision ou accuracy,
  • SPC = TN / N = TN / (FP + TN) = 1 − FPR appelé spécificité ou taux de vrais négatifs (True Negative Rate)
  • PPV = TP / (TP + FP) le taux de positifs prédits (positive predictive value)
  • NPV = TN / (TN + FN) le taux de négatifs prédits (negative predictive value)
  • FDR = FP / (FP + TP) correspondant au false discovery rate

Bref, on convertit cette matrice en probabilités conditionnelles, et beaucoup de notions peuvent être définie en changeant le conditionnement.
On peut aussi essayer de visualiser ces quantités. La représentation graphique la plus connue est probablement la courbe ROC, Receiver Operating Characteristic. L’idée est simple: il s’agit de juger du modèle, indépendmment du seuil s choisi. Ou plutôt de se donner un outils permettant de choisir le seuil s. On définie alors la fonction de sensibilité

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-07.png

et la fonction de spécificité

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-08.png

La courbe ROC est alors la courbe

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-09.pngSi cette courbe coïncide avec la diagonale, c’est que le modèle n’est pas plus performant qu’un modèle aléatoire (où on attribue la classe au hasard). Plus la courbe ROC s’approche du coin supérieur gauche, meilleur est le modèle, car il permet de capturer le plus possible de vrais positifs avec le moins possible de faux positifs. Notons que de part sa construction, la courbe ROC est invariante par toute transformation monotone croissante de la fonction de score (ce qui peut être particulièrement intéressant si on s’amuse à “normaliser” la fonction de score).
De plus, on notera que l’aire sous la courbe ROC doit pouvoir être vu comme une mesure de la qualité de l’ajustement. Le dessin ci-dessous montre la construction de la courbe ROC pour un échantillon discriminé de manière assez simple, par analyse discriminante.

Si on regarde une relecture probabiliste, on obtient la courbe suivante,

Il existe une autre courbe relativement classique, appelée courbe de lift. Cette dernière correspond à la courbe de Lorenz (que j’avais évoquée ici). On pose alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-10.png

et la courbe de lift est alors la courbe https://perso.univ-rennes1.fr/arthur.charpentier/latex/ROC-11.png. Là aussi, un calcul d’aire donne un indicateur de la qualité, et on retrouve l’indice de Gini.

There are known knowns. There are things we know that we know….

Après la publication de l’opuscule ce qui vient à nous, i.e. le premier volet du catalogue de la Biennale d’Art Contemporain (ici), le second opuscule ce qui devient a été mis en ligne voilà 15 jours (ici), et finalement, le troisième – ce qui survient – vient d’être mis en ligne ce week-end (). J’en parle parce que Julien Prévieux y a repris une interview que nous avions enregistré il y a quelques semaines, où nous discutions de la gestion des risques, des catastrophes improbables, et de la prédiction du futur en environnement incertain… Une expérience géniale, d’autant plus que Julien a une culture fabuleuse sur ces sujets, et surtout un œil très différents des chercheurs académiques, des étudiants, ou des praticiens de l’assurance ! Les discussions qu’on a pu avoir ont beaucoup nourri ma réflexion sur la prédiction des catastrophes, ceux qui survolent régulièrement le blog ont du s’en rendre compte !

La suite à l’ouverture de la biennale, avec en particulier le travail de Goldin & Senneby….

PS pour les amateurs, la vidéo dont est extrait le titre est en ligne ici,

Acheter des pochettes Pokemon ? et puis quoi encore !

Il y a quelque temps, Olivier Gergaud était venu à Rennes présenter son papier Stardom and Perceived Talent.Nous en avions profité pour discuter un peu d’un ancien papier qu’il avait fait sur les Pokemon (ici). Il faut avouer que je suis un peu jaloux car il essaye de répondre à une question que je me posais: le prix des cartes Pokemon a-t-il un lien avec les “pouvoirs” de la carte ou sa rareté.

Mais aujourd’hui je voulais répondre à une question simple que je me pose depuis que mon fils me dit (sans grand succès) “Papa, achète moi des cartes Pokemon” (ou “Naruto“, ou je ne sais quoi encore). Mon argument, en tant que père rationnel, est que ces cartes couchent très chers. Aujourd’hui, je vais le prouver. Supposons qu’il existe 200 cartes distinctes, et que ces cartes sont vendus 5 euros le paquet de 10. Pour simplifier, supposons que ces cartes soient vendues individuellement dans des cachés hermétiques à 50 centimes pièce. On va supposer que les cartes soient équi-réparties. Supposons que je puisse acheter une vignette spécifique (par exemple en ligne), dont je connais le personnage, pour 1 euro (soit le double du prix si le personnage est inconnu à l’achat). Quelle stratégie doit-je adopter ? et combien ça va me coûter ?
Je note https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-17.png le nombre de carte total (ici 200 car je suppose qu’il existe 200 personnages). Je note https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-01.png le coût d’une carte cachée (50 centimes ici) et https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-02.png le coût d’une carte non-cachée (1 euro). Supposons que j’ai à en ma possession les cartes https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-03.png (et que je n’ai pas les suivantes). Si je veux récupérer une nouvelle carte (que je noterais https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-05.png), en achetant au hasard i.e. dans des sachets hermétiques, le nombre de cartes à acheter est https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-06.png, qui suit une loi géométrique, de paramètre https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-07.png, où

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-08.png

Par exemple, si j’ai https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-09.pngcartes, et qu’il m’en manque simplement une, alors la probabilité de la tirer est

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-10.png

Pour revenir au tirage d’une nouvelle carte que je n’avais pas (alors que j’en https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-04.png), le coût moyen de ma stratégie est

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-11.png

J’ai donc intérêt à acheter cette carte directement dès lors que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-12.png

ce qui devient intéressant dès lors que

 https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-13.png

Dans mon exemple numérique, https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-02.png étant le double de https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-01.png, je commence à acheter directement dès lors que j’ai en ma possession 100 cartes distinctes. Le coût moyen est alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-14.png

Aussi, dans le cas que j’évoquais, on obtient 171€,
> .5*sum(200/(200-(1:101)-1))+1*(200-101+1)
[1] 171.0933

Je peux comparer à deux stratégies que l’on connaît. La première est de tout acheter à un collectionneur (i.e. à découvert, ce qui limite l’intérêt du “jeu“, j’en conviens). Dans ce cas, le coût total sera

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-21.png

soit 200€. L’autre stratégie est d’acheter continuellement des paquets clos chez mon buraliste. Dans ce cas, le coût moyen serait

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-16.png

Le calcul exact donne ici 587,80€, comme  le montre le calcul ci-dessous,
> 0.5*200*sum(1/1:200)
[1] 587.8031

mais on peut aussi utiliser une approximation de la forme

https://perso.univ-rennes1.fr/arthur.charpentier/latex/poke-22.png

(ce qui donnerait 529€ environ). Pour ceux qui ont besoin de voir pour croire, j’ai représenté ci-dessous le prix en fonction du moment où je décide de passer des cartes fermées aux cartes visibles.

Ça y est, j’ai enfin un argument rationnel à faire valoir à mon fils…. même si je ne suis pas certain de réussir à le convaincre aussi facilement…

Le bagplot, comme extension du boxplot

Lundi 15, Michel Grun Rehomme avait présenté les travaux qu’il avait obtenu avec Noureddine Ben Lagha sur des problèmes d’écrêtement de coûts de sinistres, que l’on peut rapprocher de la détection de points atypiques, ou outliers. Ce matin, Alexander McNeil nous a rappelé que les bagplot pouvaient être très intéressants en gestion des risques (mentionné ici). Je ne peux pas m’empêcher de revenir un petit peu sur cette notion, en m’écartant des applications financières proposées par Alexander pour proposer des applications en démographie et en climatologie (inspiré par des travaux de Rob Hyndman).

  • Le boxplot et la régression quantile

Un outil classique (en dimension 1) est le boxplot, proposé par John Tukey en 1977. Par exemple, sur les coûts de sinistres en assurance auto, on obtient

On peut aussi faire un boxplot par classes (ou tranches) d’âge,

qui correspond, en passant à la limite (i.e. en considérant l’âge comme une variable continue), à un ensemble de régressions quantiles.

On peut aussi regarder l’impact de l’ancienneté du véhicule (à la souscription du contrat) sur le coût individuel des sinistres,

avec en rouge la régression (linéaire) du log du coût sur l’âge du véhicule. Là aussi, on peut faire des régressions quantiles,

Le code pour obtenir les régressions quantiles lissées est assez simple, une fois installé les library(quantreg) et library(splines),

> library(quantreg);  library(splines)
> fit = rq(log(cout)~bs(agecond,df=8),data=BASECOUT,tau=.5)
> xu=seq(10,110)
> yu=predict(fit,newdata=data.frame(agecond=xu))
> lines(xu,yu,col="red",lwd=2)

Cette notion a été généralisée en multivariée sous le nom de bagplot. Par exemple, pour déterminer la “médiane” en dimension 2, on va recherche un ensemble (appelé bag) contenant 50% des observations.

  • La profondeur

Cela se fait à l’aide du concept de profondeur statistique, introduite par John Tukey en 1975. En dimension 1, si

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bag-01.png

désigne la fonction de répartition, la “profondeur de demi-espace” s’écrit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bag-02.png

qui est maximal pour la médiane (on obtient alors https://perso.univ-rennes1.fr/arthur.charpentier/latex/bag-03.png). Plus généralement, pour un quantile d’ordre https://perso.univ-rennes1.fr/arthur.charpentier/latex/bag-04.png,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bag-05.png

Cette notion se généralise en dimension supérieure en posant

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bag-06.png

On notera que Fraiman et Muñiz ont proposé une extension en 2001. S’il n’y a pas d’ex-aequos, en dimension 1,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bag-07.png

qui appartient à https://perso.univ-rennes1.fr/arthur.charpentier/latex/bag-08.png, et dont la médiane est la plus profonde. Liu a proposé une profondeur dite simpliciale en 1984,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bag-09.png

D’un point de vue statistique, la fonction de répartition est inconnue, mais on dispose de https://perso.univ-rennes1.fr/arthur.charpentier/latex/bag-10.png. On note alors https://perso.univ-rennes1.fr/arthur.charpentier/latex/bag-11.png l’équivalent empirique de la profondeur. En dimension 2, ça sera https://perso.univ-rennes1.fr/arthur.charpentier/latex/bag-12.png où https://perso.univ-rennes1.fr/arthur.charpentier/latex/bag-13.png. La méthode de construction est détaillée ci-dessous,

  • Construction d’un bagplot

Pour cela, on note https://perso.univ-rennes1.fr/arthur.charpentier/latex/bag-14.png le plus petit nombre d’observations contenu dans un demi plan passant par https://perso.univ-rennes1.fr/arthur.charpentier/latex/bag-15.png. On note https://perso.univ-rennes1.fr/arthur.charpentier/latex/bag-16.png l’ensemble des points https://perso.univ-rennes1.fr/arthur.charpentier/latex/bag-15.png tels que https://perso.univ-rennes1.fr/arthur.charpentier/latex/bag-17.png contienne au moins https://perso.univ-rennes1.fr/arthur.charpentier/latex/bag-18.png points. On notera qu’on peut montrer que https://perso.univ-rennes1.fr/arthur.charpentier/latex/b ag-20.png, et que ces régions sont convexes. Rousseeuw et Ruts ont proposé un algorithme permettant de construire un tel ensemble de région en 1998. Sur l’exemple suivant, on considère 20 points en dimension 2,

Notons qu’en changeant de méthode, on obtient des résultats qui peuvent être très changeant.

> library(mnormt)
> Z=rmnorm(20,mean=c(0,0),varcov=matrix(c(1,.6,.6,1),2,2))
> source("http://www.wiwi.uni-bielefeld.de/~wolf/software/R-wtools/bagplot/bagplot.R")
> bagplot(Z,factor=1,dkmethod=1)
> bagplot(Z,factor=1,dkmethod=2)

 

La méthodologie a été proposée par Rousseeuw, Ruts et Tukey en 1999 (ici).

  • Copules Gaussiennes, versus copules de Clayton

On peut essayer de voir si les bags changent de forme lorsque la copule sous-jacente change. L’exemple suivant montre le bagplot d’un vecteur gaussien,

alors que l’exemple suivant montre le bagplot pour une copule duale de Clayton, avec des marges Gaussiennes.

Dans le premier cas, on retrouve que les bagplots ressemblent des ellipses (ce qui correspond à ce que j’avais raconté ici sur les ellipses de confiance). On retrouve cette ellipse sur le graphique dessous, avec 5000 points (simulés) au lieu de 200,

pour la copule Gaussienne, et

pour la copule de Clayton. La petite animation ci-dessous montre l’évolution des bag en fonction de la corrélation du vecteur Gaussien,

  • Illustration sur les données Loss-ALAE

Considérons la base (célèbre) des coûts de sinistres, et des frais associés.

 

>  library(evd); data(lossalae)

La figure ci-contre montre l’échantillon de 1500 points, avec le logarithme des coûts, et le logarithme des frais. Classiquement, on s’intéresse à la loi conditionnelle des frais sachant le coût du sinistre. On peut alors tracer des boxplot conditionnels (comme pour le coût des sinistres conditionnels à l’âge du véhicule),

Mais l’idée des bagplot est sensemblement différentes. On s’intéresse véritablement au couple de variables.

  • Du bagplot fonctionnel

Rob Hyndman a également utilisé ces méthode de bagplot quand la taille devient infiniment grande, i.e. dans un cadre fonctionnel. Par exemple, on peut utiliser cette technique pour détecter des années atypiques sur des évolutions de température, ou dans un ajustement de modèle de type Lee-Carter en mortalité. Ces deux cas sont abordées par Rob Hyndman dans la library(rainbow). Pour faire simple et rapide, on peut définir une profondeur fonctionnelle,

Pour commencer, on peut regarder les température moyenne mensuelles de la température à la surface de la mer, entre 1951 et 2007, obtenue ici.

On note que 4 années se détachent des autres, 1982/1983 et 1997/1998, qui sont particulièrement élevées, mais pas seulement. Certaines semblent décallées de quelques mois. Le code est ici simplement

 

> library(rainbow)
> fboxplot(data = ElNino, plot.type = "functional", type = "bag")
> fboxplot(data = ElNino, plot.type = "bivariate", type = "bag")

Dans le cas des données démographiques, il faut bricoler un peu. Repartons des données français que nous avons étudiées plus d’une fois,

> DEATH=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/R/FR-Deaths_1x1.txt",header=TRUE)
> EXPOSURE=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/R/FR-Exposures_1x1.txt",header=TRUE)
> DEATH$Age=as.numeric(as.character(DEATH$Age))
> DEATH$Age[is.na(DEATH$Age)]=110
> EXPOSURE$Age=as.numeric(as.character(EXPOSURE$Age))
> EXPOSURE$Age[is.na(EXPOSURE$Age)]=110
> YEAR=unique(DEATH$Year);nC=length(YEAR)
> AGE =unique(DEATH$Age);nL=length(AGE)
> rownames(MUH)=AGE
> colnames(MUH)=YEAR
> rownames(MUF)=AGE
> colnames(MUF)=YEAR
> MUH=MUH[1:90,]
> MUF=MUF[1:90,]
> MUHF=fts(x = AGE[1:90], y = log(MUH), xname = "Age", 
+          yname = "Log Mortality Rate")
> MUFF=fts(x = AGE[1:90], y = log(MUF), xname = "Age", 
+          yname = "Log Mortality Rate")
> fboxplot(data = MUHF, plot.type = "functional", type = "bag")
> fboxplot(data = MUHF, plot.type = "bivariate", type = "bag")


avec un code similaire pour les femmes. Le graphique suivant montre l’évolution des taux de mortalité des hommes en fonctions de leur âge, pour différentes années,

On observe que certaines années se détachent assez nettement, en particulier sur les âges allant de 15 à 40 ans. Comme le montre le graphique à gauche, ce sont les années de guerre qui ressortent fortement (surtout la première guerre mondiale). Pour les femmes, l’évolution des taux de mortalités en fonction des âges semble en revanche beaucoup plus régulière, puisqu’aucune année atypique ne semble ressortir,

Ceci se confirme sur le graphique de droite, où, contrairement aux hommes, aucune année ne semble réellement se détacher. Notons qu’au lieu de chercher des années atypiques, on peut chercher faire une lecture duale, en cherchant les âges qui se ont eut un changement en terme d’évolution du taux de mortalité sur un siècle différent. Par exemple pour les femmes, les taux de mortalité ont l’évolution suivante au cours du siècle,

ce qui permet de repérer quelques âges un peu atypiques, i.e. les bébés, qui ont eu un gain en terme de taux de mortalité considérable, relativement aux autres âges. En revanche, pour les hommes, il n’y a pas particulièrement d’âge qui soit statistiquement différent des autres,

  • Applications en terme de stress scenarios

Nous avons vu qu’il y avait presque toujours (pour les hommes comme pour les femmes) des outliers, i.e. des courbes atypiques. Il serait donc normal de s’attendre à en trouver dans les scénarios de prédiction. Les bagplot sont alors un outils intéressant. Si on regarde les prédictions obtenues par la méthode de Lee & Carter par exemple, ou plutôt des scénarios générés par la méthode de Lee & Carter,

> LCH <- fdm(BASEH,order=10)
> LCHf=forecast(LCH,100,method="ets")
> LCHg=simulate(LCHf,1)
> LCHg1=as.matrix(LCHg[,,1])
> MUHf=LCHg1[1:90,]
> MUHF=fts(x = AGE[1:90], y = log(MUHf), xname = "Age",
+         yname = "Log Mortality Rate")
> library(rainbow)
> fboxplot(data = MUHF, plot.type = "functional", type = "bag")
> fboxplot(data = MUHF, plot.type = "bivariate", type = "bag")

On peut alors identifier clairement des années atypiques – pour un scenario donné. Bref, on a un outils pour vérifier que l’on génère effectivement des points atypiques, un peu extrêmes, mais surtout, on peut les identifier. On peut alors creuser un peu pour mieux comprendre ce qui se passe….

Pour aller plus loin sur ces histoires de bagplot et aux régressions quantiles, je renvoie ici ou (Thierry Roncalli avait remis la régression quantile au goût du jour il y a quelques années, à l’époque où je programmais encore en Gauss, je renvoie aussi ici pour des applications financières).

  • Les quantiles bivariés

En fait, le point de départ du billet était très lié d’une notion qui méritera un autre billet: comment définir des quantiles (et plus généralement des mesures de risques) dans un cadre multivarié….