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) !

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 !

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 Χ², 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 Χ², 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 Χ² 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 Χ² 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...

Faire un test de Χ² 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] :

• avec ces données, il calcule le critère du Χ² : (khi2Observe <- sum((histo$counts - EffTheo)^2 / EffTheo))
• Et là toute la puissance de R entre en jeu, avec l’application d’un test de Χ² à 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 Χ² é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.

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)