Digital resources in the Social Sciences and Humanities OpenEdition Our platforms OpenEdition Books OpenEdition Journals Hypotheses Calenda Libraries OpenEdition Freemium Follow us

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é….

Discussion on stress scenarios

Friday morning, I had the honor to discuss a presentation by Alexander McNeil, on Stress Testing and Reverse Stress Testing, at the Financial Risks International Forum on Risk Dependencies (here).

This was an opportunity to rediscover techniques I have studied briefly a few years ago, on outliers detection, namely the bagplot (I will probably upload a post on that topic soon, in French unfortunately). The slides of my discussion are available here.

http://freakonometrics.hypotheses.org/5338

Quand on a les boules…..

Vendredi dernier, le bureau des élèves avait organisé un tournoi de speed-pétanque dans le cloître de la fac de Sciences Eco (comme l’an dernier, ici). Nous avions monté une équipe, qui n’a pas eu la chance de franchir les qualifications*…. Pourtant, pour reprendre l’expression de Nathalie, la pétanque c’est facile, il suffit de minimiser la distance entre la boule ou le cochonnet. Ou plus précisément, on recherche à résoudre un programme de la forme suivante,

Pour cela, il y a deux méthodes principalement,

  • essayer de minimiser , i.e. minimiser la distance entre ma boule et le cochonnet. On appelle ça pointer,
  • essayer de maximiser , i.e. maximiser la distance entre la boule de mon adversaire et le cochonnet. On appelle ça tirer,

Malheureusement, c’est formellement un peu plus compliqué que ça, car en vrai il s’agit d’un jeu successif, où l’équipe la plus loin du cochonnet joue. Et en pratique, ça se complique encore, car entre la pluie, le terrain en pente, le vent, les gens qui fond des bruits autour, c’est pas simple de gagner (qui a dit qu’on était des mauvais perdants ?).

En fait, et c’est là que je voulais en venir, si modéliser une partie de pétanque est complexe – de part l’aspect dynamique – quelqu’un avait suggéré vendredi un jeu un peu plus simple en apparence: pourquoi pas lancer le cochonnet à la fin, une fois toutes les boules placées ? Le dual du jeu de pétanque en quelque sorte… Par exemple, ci-dessous, on a les boules rouges, et on se demande où il faudrait lancer le cochonnet pour gagner,

La réponse est simple (et les organisatrices du tournoi connaissent la réponse puisque je l’ai enseigné l’an dernier en master 1): il s’agit d’un problème de plus proches voisins.

On sait alors numériquement tracer les différentes régions. On peut alors même complexifier le problème, en se demandant où lancer le cochonnet pour gagner de deux points, au moins,

ou de trois points, ou plus,

Bref, l’année prochaine, je demande à ce que l’on change les règles, ça nous laissera une chance de ne pas être ridicules ! Reste à voir si on est capable de lancer le cochonnet dans la région rouge, mais promis, on va s’entraîner….
* comme en témoignent certains documents (ici par exemple) nous n’avons pas forcément été ridicule, mais comme il nous fallait généralement un quart d’heure pour se chauffer à chaque partie, ce n’était pas simple de briller…

Exchange rates and correlation matrices

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.euro_dollar_yen_s.jpg I wanted to upload here a small problem I started to work on…. unfortunatley, I could not find (yet) a proper answer. Any comments and suggestions are welcomed. The problem is simple: consider 3 random variables that can be interpreted as exchange rates. Thus X is the exchange rate USD versus EUR, Y is the exchange rate EUR versus YEN, while Z is the exchange rate USD versus YEN. We then have the following constraint XY=Z. Given that constraint, what are the implied constraints on the correlation matrix of random vector (X,Y,Z) ?

  • Bounds on correlation

There are classically standard constraints on correlation. For instance, from Cauchy-Schwarz’s inequality

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

but those bounds are not necessarily sharp, given distributions of https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-02.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-03.png. From Fréchet-Hoeffding bounds, we know that

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

where https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-05.png. Further, correlation matrices are necessarily positive-semidefinite matrices.
To go further to find possible bounds, recall that

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

and therefore

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

i.e.

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

Similarity, https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-09.png will be a function of https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-10.png.
Thus, given marginal distributions of https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-02.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-03.png, covariances (and also correlations https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-11.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-12.png) are functions of L'image “https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-13.png” ne peut être affichée car elle contient des erreurs.,https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-14.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-15.png.
Somehow, we have obtained 2 degrees of freedom: if there are no additional constraints on those two coefficients, there should be no other constraints than having positive-semidefinite matrices.
I guess that https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-14.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-15.png should not be connected. Those quantities are related to the co-skewness* and they should be different when dependence is asymmetric.
Given a 3 parameter copula for pair https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-16.png, respectively describing global dependence (that can be related to L'image “https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-13.png” ne peut être affichée car elle contient des erreurs. and tail asymmetry (that should be related either to https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-14.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-15.pngdepending on the corner) it should be possible to derive any kind of triplets.

  • Running simulations

I have tried to generate a lot of samples. The idea is to draw pairs from mixtures of Marshall-Olkins’ copula (wich are asymmetric), where parameters where generated randomly. I draw 250 values according to that copula. Then, the third variable is obtained taking the product of normalized pairs (or Studentized). And then calculate the associated correlation matrix.

  1.  draw https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-17.png from https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-18.png (but actually any positive distribution should work) and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-19.png
  2. do 250 times the following (index https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-20.png)
  • draw https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-21.png with a Bernoulli https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-22.png distribution
  • draw&nbs
    p;https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-23.png from copula https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-24.png if https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-25.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-27.png if https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-26.png
  • set https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-28.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-29.png, and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-30.png
  1. estimate the correlation of sample https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-31.png

In dimension 3, I have the following graph,

In order to get a better understanding, I look at slices, for instance when https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-32.png. On the graph on the right, to dark area are values that can not be obtained since correlation have to be positive-semidefinite. Dots are points that have been obtained in a particular scenario. I have also included the convex hull if the points obtained with all the scenarios (here 50,000 scenarios), since I assume that the set of positive correlation is necessarily convex. Graphs below have been obtained when https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-33.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-34.png

So far I have not been able to obtain all possible correlations, I guess…  If anyone has suggestions…
* Recall that the co-skewness matrix, in dimension 2, for a pair https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-35.png is simply

https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-36.png

where https://perso.univ-rennes1.fr/arthur.charpentier/latex/change-corr-37.png.

On the return period of the 2003 heat wave

The paper on the 2003 heat wave has been submitted a few months ago, and now, the preprint can be downloaded (here for the preprint version, or here on Hal). The paper has a long story, that I will probably summarize here, some day…

Extremal events are difficult to model since it is difficult to characterize formally those events. The 2003 heat wave in Europe was not characterized by very high temperatures, but mainly the fact that night temperature were no cool enough for a long period of time. Hence, simulation of several models (either with heavy tailed noise or long range dependence) yielddifferent estimations for the return period of that extremal event. 

The paper has been presented in several conferences, in started a couple of days after the birth of Romane, in March 2005, in talks on dependencies in environmental series….

Lillie versus Kolmogorov-Smirnov

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 Lillie1, il faut revenir sur Kolmogorov-Smirnov. L’idée est de tester

https://perso.univ-rennes1.fr/arthur.charpentier/tex/KS-01b.png

contre

https://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-02.png

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éorème de Glivenko-Cantelli. En fait, Andreï Kolmogorov et Vladimir Smirnov ont montré que l’on pouvait encadrer la vitesse de convergence,

https://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-03.png

ce qui se montre en notant que https://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-04.png converge en loi (si l’hypothèse H0 était la bonne) vers un pont brownien changé de temps par la fonction quantile https://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-05.png.
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 H0)
> 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 enrouge pour l’échantillon Gaussien, enbleu pour l’échantillon uniforme et envert 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 ), 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é https://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-encore.png 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

https://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-06.png

contre

https://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-07.png

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 loirouge, et maintenant on va tester l’ajustement de la loibleue,

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)
&nbsp
;     95%
0.2943206
> quantile(SNML,.95)
95%
0.1923295
Sous R, deux fonctions distinctes sont à utiliser. Par exemple l’ajustement à une loi https://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-08.png 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 https://perso.univ-rennes1.fr/arthur.charpentier/tex/ks-09.png 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 H0 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.
1 Lillie pour Hubert Lilliefors qui, le premier, a proposé ce test

Ecart absolue ou écart type ?

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

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

alors que l’erreur absolue moyenne est définie comme

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

(que je définie ici comme l’écart à la moyenne, certain considérant la distance à la médiane, qui minimise cette quantiié, 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 appellerait “dispersion” autour de la moyenne.

Si l’on prend deux observations, ces quantités sont rigoureusement identiques, car

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

compte tenu du fait que

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

En revanche, comme le montre le graphique ci-dessous,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/eqt-00-col.png

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,

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

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’accroît (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“).

"sendo l'intento mio scrivere cosa utile a chi la intende…"