Les aventuriers du λ perdu

jeudi 11 avril 2013
par  Alain BUSSER

Sans aller jusqu’à tester l’adéquation à une loi exponentielle (hors programme), on peut chercher expérimentalement si les intervalles entre éruptions successives du Piton de la Fournaise suivent une loi exponentielle. La même question peut être vue pour la durée des éruptions. Il se pose alors deux questions :

  1. les durées suivent-elles, ou non, une loi exponentielle ? (statistiquement, c’est un test)
  2. si oui, quel est son paramètre ? (là c’est une estimation)

L’observatoire volcanologique de la Réunion a pour rôle de collecter des données sismographiques. Grâce à son directeur technique Philippe Kowalski, j’ai pu obtenir les dates des débuts et fins des éruptions du Piton de la Fournaise sur une vingtaine d’années. Ce qui permet de calculer par différences

  • la durée d’une éruption (intervalle entre le début de l’éruption et la fin de l’éruption)
  • l’intervalle entre deux éruptions successives (intervalle entre le début de l’éruption et le début de l’éruption suivante).

Le but du TP était de voir si les durées suivent bien une loi exponentielle, au moins approximativement, et le cas échéant, d’estimer la valeur du paramètre de cette loi [1].

La mission de l’observatoire volcanologique

Le Piton de la Fournaise est surveillé par tout un réseau de capteurs sismographiques, qui envoient en permanence des signaux à l’observatoire. Chaque secousse est donc captée, à des instants différents, par différents capteurs. Le temps mis par la secousse pour aller d’un capteur à l’autre, avec la distance entre capteurs, permet donc d’estimer la vitesse de propagation des ondes sismiques entre les capteurs. Ce sont les variations de cette vitesse, symptomatiques de variations de la structure du sous-sol, qui permettent de savoir s’il y a ou non risque d’éruption.

Coïncidence [2] : Le sujet du Rallye mathématique 2013 illustre ceci (exercice 3) !

PDF - 144.2 ko
le TP en version pdf
L’énoncé donné aux élèves
OpenDocument Spreadsheet - 32 ko
dates d’éruptions
Le fichier de travail, au format odt

Les élèves, habitant la même commune que le volcan qui est ici l’objet d’étude, s’intéressent très vite à ces données ; en particulier, comme il se trouve qu’en ce moment le volcan est particulièrement calme, le grand nombre d’éruptions répertoriées depuis les années 1980 les étonne immédiatement : Mais alors, le Piton de la Fournaise entrerait-il parfois en éruption en cachette, à l’insu de mon plein gré ?

L’allure de l’histogramme montre immédiatement que la loi suivie par les intervalles entre éruptions n’est ni uniforme, ni normale. Elle a l’air exponentielle (seule loi au programme, à part celles ci-dessus).

L’idée de départ était la suivante : En supposant que la loi est vraiment exponentielle, les effectifs théoriques avec des classes de largeur 30 (jours) sont de la forme 94∫30k30(k+1)λe-λxdx=94(1-e-30λ)e-30λk, donc une suite géométrique. Alors les logarithmes des effectifs devraient être fonction affine de k, et les points les représentant devraient être alignés. Cet alignement n’est pas évident :

on en a la confirmation

Le carré du coefficient de corrélation n’est pas si élevé :

On lit sur la « courbe de tendance » une valeur de λ proche de 0,3 en mois, donc 0,01 en jours (soit une espérance de 100 jours).

Ceci dit, comme un histogramme est une fonction en escalier, on est en train de chercher à approcher une fonction en escalier par une droite, ce qui est un choix quelque peu étrange... D’un autre côté, les escaliers, pour escalader un volcan, c’est plus pratique !

Cependant, les points d’abscisses 0, 90, 120, 150 et 210 sont mieux alignés, ce qui donne alors λ=1/60. La valeur donnée dans l’énoncé a été trouvée par tâtonnements sous tableur.

Mais un élève a fait la remarque que lorsqu’on affiche une « courbe de tendance », d’autres choix sont possibles, et qu’il est plus simple de faire une régression exponentielle sur les effectifs, plutôt qu’une régression linéaire sur les logarithmes des effectifs :

Ce qui suggère une valeur de 0,21/30=0,007 pour λ soit environ 147 jours pour l’espérance.


test du Χ2 sous tableur

On commence par vérifier qu’il n’y a pas de corrélation entre les durées et les dates :

On calcule ensuite les effectifs théoriques, qui sont les produits des fréquences théoriques par l’effectif total ; puis les résidus, qui sont les carrés des écarts entre les effectifs réels et théoriques, divisés par les effectifs constatés. Enfin, leur somme ; on trouve environ 6,9.

À moins d’être un expert en test de Χ2, il reste maintenant à estimer la probabilité d’avoir une telle valeur. Heureusement, le tableur de Libre Office sait calculer la probabilité critique, avec la loi du Χ2, mais également la valeur jusqu’à laquelle on accepte d’aller si les écarts à la loi exponentielle constatés sont dus à la seule fluctuation d’échantillonnage ; avec la loi du Χ2 inverse :

On a entré le risque 0,05 et le nombre de degrés de libertés, 15 parce qu’il y a 16 classes dans l’histogramme [3] ; la valeur du Χ2 admissible est 7,26 qui est supérieure à la valeur constatée, ce qui amène à accepter l’hypothèse selon laquelle la loi suivie est exponentielle, avec un risque inférieur à 5%.

Mais si la loi suivie par les intervalles entre éruptions était exponentielle, on s’attendrait à ce que sa moyenne et son écart-type soient tous deux proches de l’inverse du paramètre, et donc proches l’un de l’autre, ce qui n’est clairement pas le cas, la moyenne valant 97,8 jours, et l’écart-type 184 jours, soit presque le double ! D’ailleurs, en simulant une loi exponentielle de même paramètre avec l’algorithme du cours, on a des courbes visuellement identiques à celle ci-dessus, mais avec un écart-type voisin de l’inverse du paramètre. Ce qui montre que le regard peut tromper...

test du Χ2 avec GNU R

Faire un test de Χ2 avec le tableur, c’est un peu comme réinventer la roue : On ferait peut-être mieux de laisser le logiciel R faire le travail, vu que ce logiciel est spécialisé dans ce genre de tâches [4]. Voici comment Hubert Raymondaud s’y prend (le détail est lisible dans le pdf téléchargeable en bas d’article) :

  • Il exporte les données, depuis le tableur, au format comma-separated values, qu’il est facile d’importer dans toutes sortes de logiciels, en particulier R ;
  • Il choisit comme « working directory », dans R, le dossier où il a placé le fichier csv en question (téléchargeable en bas de l’article, avec le fichier R de la session).
  • Il lit le fichier de données au format csv, avec EruptA <- read.csv("EruptDuree_2.csv", sep=";", dec = ",") (en signalant à R que le séparateur de données est le point-virgule, et la virgule décimale une vraie virgule parce que cocorico, non mais)
  • il enlève les durées nulles (de moins d’un jour), celles-ci étant parfois la détection d’une même éruption par plusieurs capteurs ; le code est EruptA <- EruptionA[EruptionA$duree != 0, ]
  • il signale à R que la session se fera sur ces données, en les « attachant », avec attach(EruptA).
  • il construit les classes pour l’histogramme, en mettant à part (dans une seule classe large) la donnée la plus grande, isolée des autres. Cette liste de bornes d’intervalles se fait en concaténant (« c ») les valeurs construites par une séquence à la GeoGebra (ou calculatrice) : subdiv1 <- c(seq(0, 360, 30), 1560) ;
  • ce qui permet d’avoir l’histogramme avec la syntaxe suivante :
    histo <- hist(duree, breaks = subdiv1, freq = FALSE, right = FALSE)
    grid()
  • Comme estimation de λ, il choisit l’inverse de la moyenne de la série restante ; il s’en sert pour définir la loi de densité exponentielle ayant cette valeur pour paramètre, et la dessiner sur l’histogramme, avec ce code :
    lambda <- 1 / mean(duree)
    abscissesX <- seq(min(duree), max(duree), length = 1000)
      lines(abscissesX, dexp(abscissesX, lambda), col = "red", lwd = 2)

ce qui donne ce graphique [5] :

  • Il se sert de la fonction de répartition de cette loi exponentielle pour calculer les effectifs théoriques :
    OrdreSubdiv1 <- length(subdiv1)
    EffTheo <- NULL
    for(i in 1:(OrdreSubdiv1 - 1)){
      EffTheo <- c(EffTheo,
         (pexp(histo$breaks[i + 1], lambda) - pexp(histo$breaks[i], lambda)) * length(duree))
    }

 ; les effectifs sont alors les suivants [6] :

classes effectifs effectifs théoriques
moins de 30 jours 41 24.830499
de 30 à 60 jours 14 18.271417
de 60 à 90 jours 8 13.444945
de 90 à 120 jours 8 9.893406
de 120 à 150 jours 5 7.280021
de 150 à 180 jours 3 5.356972
de 180 à 210 jours 4 3.941905
de 210 à 240 jours 1 2.900634
de 240 à 270 jours 2 2.134419
de 270 à 300 jours 2 1.570603
de 300 à 330 jours 1 1.155722
de 330 à 360 jours 1 0.850433
plus de 360 jours 4 2.369012
  • avec ces données, il calcule le critère du Χ2 : (khi2Observe <- sum((histo$counts - EffTheo)^2 / EffTheo))
  • Et là toute la puissance de R entre en jeu, avec l’application d’un test de Χ2 à ce critère, pour avoir la probabilité critique (à comparer ensuite avec le risque) : pchisq(khi2Observe, length(histo$counts) - 2, lower.tail = FALSE)

Et comme la probabilité critique (10%) est supérieure aux 5% de risque admis, on admet que la loi est exponentielle. Mais les effectifs étant parfois trop petits (en tout cas, trop disparates), la fiabilité de ce test n’est pas extraordinaire.

Dans le document téléchargeable ci-dessous, il propose alors deux améliorations :

  1. Changer les classes pour que les effectifs soient mieux répartis, et là la probabilité critique est réduite d’un facteur 3, et passe donc en-dessous du seuil de risque de 5% : Cette fois-ci, le test de Χ2 échoue...
  2. Modifier imperceptiblement les effectifs pour qu’un test de Kolmogorov-Smirnov, basé sur la notion de fonction de répartition, puisse s’appliquer ; là encore, le test conclut à une non-exponentionnalité.

En fait, R est muni d’un outil d’estimation de paramètres pour une adéquation à une loi, ici exponentielle. En entrant

fitdistr(duree,"exponential")

On obtient à la fois une estimation de λ et la probabilité qu’elle ne vaille rien ; cette probabilité est bien trop élevée pour défendre un modèle exponentiel.

Pourquoi la loi n’est-elle pas exponentielle ?

En croyant expliquer pourquoi la loi suivie est exponentielle, un élève [7] ouvre la voie vers une explication à la non-exponentialité des intervalles entre éruptions :

Il faut dire qu’avec plus de 500000 ans d’âge, le Piton de la Fournaise pourrait être soumis à un vieillissement. Dans ce cas,

  • une loi autre qu’exponentielle (loi de Weibull par exemple, loi d’Erlang, loi log-logistique...) serait peut-être plus appropriée ici ;
  • l’intervention d’autres paramètres que λ seul : Le Piton de la Fournaise est peut-être juste un système trop complexe pour être décrit par une seule loi exponentielle.

Conclusion

Même en grimpant au sommet du volcan, les aventuriers n’ont pas réussi à retrouver le λ perdu...


Un autre exemple

Le geyser old faithful dans le parc national de Yellowstone, tire son nom de la régularité de ses éruptions. L’histogramme des intervalles entre éruptions en 2011 suggère plutôt une loi normale :

On peut charger ces données dans R en entrant

install.packages("alr3", dependencies=TRUE)

puis

library(alr3)

La variable oldfaith contient alors les données dont on a besoin, sous la forme durée/intervalle ; ce qui permet de voir si les durées des éruptions sont corrélées aux intervalles entre éruptions.

La loi suivie semble plus bimodale que normale. Un histogramme plus précis montre d’ailleurs que les faibles durées entre éruptions sont plus fréquentes que prévu (entre 52 minutes et 76 minutes)


[1on effectue donc à la fois un test et une estimation dans ce TP

[2ou pas : Après tout 2013 est à la fois l’année des maths et celle de la planète...

[3Remarques d’Hubert Raymondaud :

  1. en estimant la valeur de λ, on a perdu un degré de liberté supplémentaire et on aurait dû prendre 14 plutôt que 15 ;
  2. les dernières classes ont des effectifs trop faibles pour être fiables et auraient dû être écartées.

[4ce qui n’est absolument pas mon cas ; mais en réinventant la méthode avec le tableur, on apprend mieux comment elle fonctionne...

[5au format jpg pour affichage rapide dans le navigateur ; les exports svg, pdf et eps sont de bien meilleure qualité

[6les effectifs théoriques sont tellement théoriques qu’ils ne sont pas entiers ! En tout cas, on voit bien que plusieurs de ces effectifs sont inférieurs à 5 ce qui limite sérieusement la validité du test du Χ2...

[7Il sait lire sur Wikipedia, où se trouve la phrase suivante : Une loi exponentielle modèlise la durée de vie d’un phénomène sans mémoire, ou sans vieillissement, ou sans usure


Documents joints

le problème traité par Hubert Raymondaud avec (...)
le problème traité par Hubert Raymondaud avec (...)
R permet de faire les graphiques, mais aussi de faire les tests d’adéquation (khi-deux, (...)
les fichiers pour une séance R
les fichiers pour une séance R
le fichier des dates d’éruptions, au format csv pour import ; le fichier de la session (...)

Portfolio

FournaiseSecondesARVI

Commentaires

Annonces

Prochains rendez-vous de l’IREM

Séminaire EDIM-IREM

- Mercredi 11 octobre 2017, 14h-18h, campus du Tampon
- Mercredi 22 novembre 2017, 14h-18h, campus du Tampon
- Mercredi 7 février 2018, PTU, Saint-Denis, salle S23.6
- Mercredi 7 mars 2018, 14h-18h, campus du Tampon
- Mercredi 4 avril 2018, PTU, Saint-Denis, salle S23.6
- Mercredi 2 mai, 14h-18h, campus du Tampon
- Mardi 5 juin 2018, PTU, Saint-Denis, salle S23.6
- Mercredi 6 juin, 14h-18h, campus du Tampon

Fête de la science

Du 13 au 18 novembre 2017.
Thème : « La recherche à l’heure du numérique »

Semaine des mathématiques

Du 26 au 31 mars 2018.
Thème : « Mathématiques et mouvement »


Brèves

Comprendre la logique Shadok

dimanche 30 avril

Le meilleur cours de maths Shadok jamais réalisé...
Vidéo succulente ajoutée sur la canal des archives de l’INA le 26 avril 2017.

Décès de Kenneth Arrow

mercredi 15 mars

La vedette de la théorie du choix social, bien connue de nos lecteurs, est décédée récemment, à l’âge respectable de 95 ans.

CHAOS : une aventure mathématique

vendredi 8 mars 2013

CHAOS est un film mathématique constitué de neuf chapitres de treize minutes chacun. Il s’agit d’un film tout public autour des systèmes dynamiques, de l’effet papillon et de la théorie du chaos. Tout comme DIMENSIONS, ce film est diffusé sous une licence Creative Commons et a été produit par Jos Leys, Étienne Ghys et Aurélien Alvarez.

Sur le Web : CHAOS

Rencontres Mondiales Du Logiciel Libre Décentralisées à Saint-Joseph

mardi 28 juin 2011

C’est une manifestation qui aura lieu sur 3 jours, avec de nombreux
ateliers et conférences sur les logiciels libres.
C’est vendredi 1, samedi 2 et dimanche 3 juillet.

C’est une première dans l’île, petite soeur des Rencontres Mondiales du Logiciel Libre nationales qui se déroulent chaque année.

Le site des rencontres réunionnaises se trouve ici :
http://2011.d.rmll.info/

Yves Martin y donnera une conférence d’introduction à la géométrie hyperbolique avec CarMetal, Alain Busser parlera de sa contribution en tant que développeur à CarMetal et Nathalie Carrié présentera un logiciel d’élaboration de connaissances.
De nombreux ateliers vous y attendent : Ruby, Smalltalk, Stellarium, Audacity, Freeplane et d’autres encore...

Il y aura un « repas du libre » le samedi soir, si certains veulent s’y
inscrire en ligne.
Il y aura même une conférence sur l’agriculture libre.

Merci de consulter le programme régulièrement pour plus d’infos.

Médailles Fields 2010

mardi 24 août 2010

Les noms des quatre médaillés Fields 2010 ont été dévoilés lors de la cérémonie d’ouverture du Congrès international des mathématiciens à Hyderabad :

- Elon Lindenstrauss
- Ngô Bào Châu
- Stanislas Smirnov
- Cédric Villani.

Sur le Web : ICM 2010

L’univers de Labomath sur Netvibes

dimanche 23 mai 2010

Quand on aime les maths et qu’en plus on est prof de maths, on ne peut pas passer à côté de cet univers mathématique créé par Kostrzewa Bruno, auteur de l’excellent site personnel Labomath.
Il vous donnera peut-être envie de vous créer votre propre espace sur Netvibes et votre propre univers mathématique.
Allez-voir, c’est hallucinant !
Nathalie Carrié

MathRider : L’outil ultime ?

mardi 24 novembre 2009

MathRider ressemble un peu à Maple (serveur de maths avec calcul formel). Mais il est plus léger (moins de fonctionnalités, on s’y retrouve donc mieux). Et il est conçu pour faire de la programmation...

Cette suite logicielle (dedans il y a 3d-Xplor, GeoGebra, LaTeX etc.) est multiplateforme et les exemples correspondent assez bien au programme actuel du Lycée. Le seul reproche qu’on puisse lui faire est que l’aide est en Anglais (mais de toute façon si on veut programmer on écrit souvent des « for » et des « while »). Le chapitre sur les branchements conditionnels fait appel à un vocabulaire assez original.

Le moteur de calcul formel, MathPiper, est celui qui a été incorporé à GeoGebra.

Le blog du prof geek

lundi 16 novembre 2009

Voici un blog publié sous licence Creative Commons à consommer sans modération pour les enseignants qui utilisent l’outil informatique (et les TICE).

J’ai adoré notamment la vidéo sur le cahier de textes en ligne.

Blog découvert dans le Café pédagogique de ce matin.

Nathalie Carrié

Sur le Web : Le blog du prof geek

Cours vidéo en ligne pour le collège

dimanche 30 août 2009

Philippe Mercier, professeur à Morhange (Moselle), a mis en ligne un cours vidéo couvrant l’ensemble du programme de mathématiques du collège, de la 6e à la 3e. Cet outil pédagogique peut être utile aux collégiens, aux parents d’élèves, aux personnes en formation continue et aux formateurs. Le cours est complété par un forum d’aide en mathématiques.

Un merveilleux travail mathématique et artistique

jeudi 25 juin 2009

Maria Carla Palmeri est professeur de mathématiques dans un collège de Florence (Italie). Cette année, elle a fait utiliser Cabri à ses élèves de 11 ans, une heure par semaine pendant toute l’année. Il en est résulté une magnifique vidéo mettant en scène quelques-unes de leurs constructions et animations : Le Fabuleux Monde de Cabri.

Statistiques

Dernière mise à jour

dimanche 22 octobre 2017

Publication

772 Articles
Aucun album photo
133 Brèves
11 Sites Web
132 Auteurs

Visites

673 aujourd'hui
782 hier
2133142 depuis le début
35 visiteurs actuellement connectés