----- # Package R "talweg" Le package $-$ Time-series sAmpLes forecasted With ExoGenous variables $-$ contient le code permettant de (re)lancer les expériences numériques décrites dans cette partie et la suivante. Les fonctions principales sont respectivement * **getData()** pour construire un objet R contenant les données à partir de fichiers CSV (extraits de bases de données). Le format choisi en R est une classe R6 (du package du même nom) exposant en particulier les méthodes *getSerie(i)* et *getExo(i)* qui renvoient respectivement la $i^{eme}$ série de 24h et les variables exogènes (mesurées) correspondantes. Voir ?Data pour plus d'information, une fois le package chargé. * **computeForecast()** pour calculer des prédictions sur une certaine plage temporelle contenue dans *data <- getData(...)* * **computeError()** pour évaluer les erreurs commises par différentes méthodes. Le package contient en outre diverses fonctions graphiques *plotXXX()*, utilisées dans la partie suivante. -----r # Chargement de la librairie (après compilation, "R CMD INSTALL .") library(talweg) # Acquisition des données (depuis les fichiers CSV) ts_data <- read.csv(system.file("extdata","pm10_mesures_H_loc.csv", package="talweg")) exo_data <- read.csv(system.file("extdata","meteo_extra_noNAs.csv", package="talweg")) data <- getData(ts_data, exo_data, input_tz="GMT", date_format="%d/%m/%Y %H:%M", working_tz="GMT", predict_at=7, limit=120) # Plus de détails à la section 1 ci-après. # Prédiction de 10 courbes (jours 102 à 111) pred <- computeForecast(data, 101:110, "Persistence", "Zero", memory=50, horizon=12, ncores=1) # Plus de détails à la section 2 ci-après. # Calcul des erreurs (sur un horizon arbitraire <= horizon de prédiction) err <- computeError(data, pred, horizon=6) # Plus de détails à la section 3 ci-après. # Puis voir ?plotError et les autres plot dans le paragraphe 'seealso' ----- ${"##"} getData() Les arguments de cette fonction sont, dans l'ordre : 1. **ts_data** : séries temporelles (fichier CSV avec entête ou data.frame) ; la première colonne contient les heures, la seconde les valeurs. 2. **exo_data** : variables exogènes (fichier CSV avec entête ou data.frame) ; la première colonne contient les jours, les $m$ suivantes les variables mesurées pour ce jour, et les $m$ dernières les variables prédites pour ce même jour. Dans notre cas $m=4$ : pression, température, gradient de température, vitesse du vent. 3. **input_tz** : zone horaire pour ts_data (défaut : "GMT"). 4. **date_format** : format des heures dans ts_data (défaut : "%d/%m/%Y %H:%M", format du fichier transmis par Michel). 5. **working_tz** : zone horaire dans laquelle on souhaite travailler avec les données (défaut : "GMT"). 6. **predict_at** : heure à laquelle s'effectue la prévision $-$ et donc dernière heure d'un bloc de 24h, relativement à working_tz. data`$`getSerie(3) renvoit ainsi les 24 valeurs de 8h à 7h pour le $3^{eme}$ bloc de 24h présent dans le jeu de données. -----r print(data) #?Data ----- ${"##"} computeForecast() Les arguments de cette fonction sont, dans l'ordre : 1. **data** : le jeu de données renvoyé par getData() 2. **indices** : l'ensemble de jours dont on veut prévoir les "lendemains" (prochains blocs de 24h) ; peut être donnée sous forme d'un vecteur de dates ou d'entiers (correspondants aux numéros des jours). 3. **forecaster** : le nom du prédicteur principal à utiliser ; voir ?computeForecast 4. **pjump** : le nom du prédicteur de saut d'une série à l'autre ; voir ?computeForecast 5. **memory** : le nombre de jours à prendre en compte dans le passé pour chaque prévision (par défaut : Inf, c'est-à-dire tout l'historique pris en compte). 6. **horizon** : le nombre d'heures à prédire ; par défaut "data`$`getStdHorizon()", c'est-à-dire le nombre d'heures restantes à partir de l'instant de prévision + 1 jusqu'à minuit (17 pour predict_at=7 par exemple). 7. **ncores** : le nombre de processus parallèles (utiliser 1 pour une exécution séquentielle) -----r print(pred) #?computeForecast ----- ${"##"} computeError() Les arguments de cette fonction sont, dans l'ordre : 1. **data** : le jeu de données renvoyé par getData() 2. **pred** : les prédictions renvoyées par computeForecast() 3. **horizon** : le nombre d'heures à considérer pour le calcul de l'erreur ; doit être inférieur ou égal à l'horizon utilisé pour la prédiction (même valeur par défaut : "data`$`getStdHorizon()") -----r summary(err) summary(err$abs) summary(err$MAPE) ----- ${"##"} Graphiques Voir ?plotError : les autres fonctions graphiques sont dans la section 'seealso' : ‘plotCurves’, ‘plotPredReal’, ‘plotSimils’, ‘plotFbox’, ‘computeFilaments’, ‘plotFilamentsBox’, ‘plotRelVar’ ?plotXXX, etc. ## $\clearpage$ How to do that? ----- # Expérimentations Cette partie montre les résultats obtenus via des variantes de l'algorithme décrit à la section 2, en utilisant le package présenté à la section 3. Cet algorithme est systématiquement comparé à deux approches naïves : * la moyenne des lendemains des jours "similaires" dans tout le passé, c'est-à-dire prédiction = moyenne de tous les mardis passé si le jour courant est un lundi par exemple. * la persistence, reproduisant le jour courant ou allant chercher le lendemain de la dernière journée "similaire" (même principe que ci-dessus ; argument "same\_day"). Concernant l'algorithme principal à voisins, trois variantes sont étudiées dans cette partie : * avec simtype="mix" et raccordement "Neighbors" dans le cas "non local", i.e. on va chercher des voisins n'importe où du moment qu'ils correspondent au premier élément d'un couple de deux jours consécutifs sans valeurs manquantes. * avec simtype="endo" + raccordement "Neighbors" puis simtype="none" + raccordement "Zero" (sans ajustement) dans le cas "local" : voisins de même niveau de pollution et même saison. Pour chaque période retenue $-$ chauffage, épandage, semaine non polluée $-$ les erreurs de prédiction sont d'abord affichées, puis quelques graphes de courbes réalisées/prévues (sur le jour "en moyenne le plus facile" à gauche, et "en moyenne le plus difficile" à droite). Ensuite plusieurs types de graphes apportant des précisions sur la nature et la difficulté du problème viennent compléter ces premières courbes. Concernant les graphes de filaments, la moitié gauche du graphe correspond aux jours similaires au jour courant, tandis que la moitié droite affiche les lendemains : ce sont donc les voisinages tels qu'utilisés dans l'algorithme. <% list_titles = ['Pollution par chauffage','Pollution par épandage','Semaine non polluée'] list_indices = ['indices_ch', 'indices_ep', 'indices_np'] %> -----r library(talweg) P = ${P} #instant de prévision H = ${H} #horizon (en heures) ts_data = read.csv(system.file("extdata","pm10_mesures_H_loc_report.csv", package="talweg")) exo_data = read.csv(system.file("extdata","meteo_extra_noNAs.csv", package="talweg")) # NOTE: 'GMT' because DST gaps are filled and multiple values merged in # above dataset. Prediction from P+1 to P+H included. data = getData(ts_data, exo_data, input_tz = "GMT", working_tz="GMT", predict_at=P) indices_ch = seq(as.Date("2015-01-18"),as.Date("2015-01-24"),"days") indices_ep = seq(as.Date("2015-03-15"),as.Date("2015-03-21"),"days") indices_np = seq(as.Date("2015-04-26"),as.Date("2015-05-02"),"days") % for i in range(3): ----- ##

${list_titles[i]}

${"##"} ${list_titles[i]} -----r p1 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Neighbors", horizon=H, simtype="mix", local=FALSE) p2 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Neighbors", horizon=H, simtype="endo", local=TRUE) p3 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Zero", horizon=H, simtype="none", local=TRUE) p4 = computeForecast(data, ${list_indices[i]}, "Average", "Zero", horizon=H) p5 = computeForecast(data, ${list_indices[i]}, "Persistence", "Zero", horizon=H, same_day=${'TRUE' if loop.index < 2 else 'FALSE'}) -----r e1 = computeError(data, p1, H) e2 = computeError(data, p2, H) e3 = computeError(data, p3, H) e4 = computeError(data, p4, H) e5 = computeError(data, p5, H) options(repr.plot.width=9, repr.plot.height=7) plotError(list(e1, e5, e4, e2, e3), cols=c(1,2,colors()[258],4,6)) # noir: Neighbors non-local (p1), bleu: Neighbors local endo (p2), # mauve: Neighbors local none (p3), vert: moyenne (p4), # rouge: persistence (p5) sum_p123 = e1$abs$indices + e2$abs$indices + e3$abs$indices i_np = which.min(sum_p123) #indice de (veille de) jour "facile" i_p = which.max(sum_p123) #indice de (veille de) jour "difficile" ----- % if i == 0: L'erreur absolue dépasse 20 sur 1 à 2 jours suivant les modèles (graphe en haut à droite). C'est au-delà de ce que l'on aimerait voir (disons +/- 5 environ). Sur cet exemple le modèle à voisins "contraint" (local=TRUE) utilisant des pondérations basées sur les similarités de forme (simtype="endo") obtient en moyenne les meilleurs résultats, avec un MAPE restant en général inférieur à 30% de 8h à 19h (7+1 à 7+12 : graphe en bas à gauche). % elif i == 1: Il est difficile dans ce cas de déterminer une méthode meilleure que les autres : elles donnent toutes de plutôt mauvais résultats, avec une erreur absolue moyennée sur la journée dépassant presque toujours 15 (graphe en haut à droite). % else: Dans ce cas plus favorable les intensité des erreurs absolues ont clairement diminué : elles restent souvent en dessous de 5. En revanche le MAPE moyen reste au-delà de 20%, et même souvent plus de 30%. Comme dans le cas de l'épandage on constate une croissance globale de la courbe journalière d'erreur absolue moyenne (en haut à gauche) ; ceci peut être dû au fait que l'on ajuste le niveau du jour à prédire en le recollant sur la dernière valeur observée. % endif -----r options(repr.plot.width=9, repr.plot.height=4) par(mfrow=c(1,2)) plotPredReal(data, p1, i_np); title(paste("PredReal p1 day",i_np)) plotPredReal(data, p1, i_p); title(paste("PredReal p1 day",i_p)) plotPredReal(data, p2, i_np); title(paste("PredReal p2 day",i_np)) plotPredReal(data, p2, i_p); title(paste("PredReal p2 day",i_p)) plotPredReal(data, p3, i_np); title(paste("PredReal p3 day",i_np)) plotPredReal(data, p3, i_p); title(paste("PredReal p3 day",i_p)) # Bleu : prévue ; noir : réalisée ----- % if i == 0: Le jour "facile à prévoir", à gauche, se décompose en deux modes : un léger vers 10h (7+3), puis un beaucoup plus marqué vers 19h (7+12). Ces deux modes sont retrouvés par les trois variantes de l'algorithme à voisins, bien que l'amplitude soit mal prédite. Concernant le jour "difficile à prévoir" il y a deux pics en tout début et toute fin de journée (à 9h et 23h), qui ne sont pas du tout anticipés par le programme ; la grande amplitude de ces pics explique alors l'intensité de l'erreur observée. % elif i == 1: Dans le cas d'un jour "facile" à prédire $-$ à gauche $-$ la forme est plus ou moins retrouvée, mais le niveau moyen est trop bas (courbe en bleu). Concernant le jour "difficile" à droite, non seulement la forme n'est pas anticipée mais surtout le niveau prédit est très inférieur au niveau de pollution observé. Comme on le voit ci-dessous cela découle d'un manque de voisins au comportement similaire. % else: La forme est raisonnablement retrouvée pour les méthodes "locales", l'autre version lissant trop les prédictions. Le biais reste cependant important, surtout en fin de journée sur le jour "difficile". % endif -----r par(mfrow=c(1,2)) f_np1 = computeFilaments(data, p1, i_np, plot=TRUE) title(paste("Filaments p1 day",i_np)) f_p1 = computeFilaments(data, p1, i_p, plot=TRUE) title(paste("Filaments p1 day",i_p)) f_np2 = computeFilaments(data, p2, i_np, plot=TRUE) title(paste("Filaments p2 day",i_np)) f_p2 = computeFilaments(data, p2, i_p, plot=TRUE) title(paste("Filaments p2 day",i_p)) ----- % if i == 0: Les voisins du jour courant (période de 24h allant de 8h à 7h le lendemain) sont affichés avec un trait d'autant plus sombre qu'ils sont proches. On constate dans le cas non contraint (en haut) une grande variabilité des lendemains, très nette sur le graphe en haut à droite. Ceci indique une faible corrélation entre la forme d'une courbe sur une période de 24h et la forme sur les 24h suivantes ; **cette observation est la source des difficultés rencontrées par l'algorithme sur ce jeu de données.** % elif i == 1: Les observations sont les mêmes qu'au paragraphe précédent : trop de variabilité des lendemains (et même des voisins du jour courant). % else: Les graphes de filaments ont encore la même allure, avec une assez grande variabilité observée. Cette observation est cependant trompeuse, comme l'indique plus bas le graphe de variabilité relative. % endif -----r par(mfrow=c(1,2)) plotFilamentsBox(data, f_np1); title(paste("FilBox p1 day",i_np)) plotFilamentsBox(data, f_p1); title(paste("FilBox p1 day",i_p)) # En pointillés la courbe du jour courant + lendemain (à prédire) ----- % if i == 0: Sur cette boxplot fonctionnelle (voir la fonction fboxplot() du package R "rainbow") l'on constate essentiellement deux choses : le lendemain d'un voisin "normal" peut se révéler être une courbe atypique, fort éloignée de ce que l'on souhaite prédire (courbes bleue et rouge à gauche) ; et, dans le cas d'une courbe à prédire atypique (à droite) la plupart des voisins sont trop éloignés de la forme à prédire et forcent ainsi un aplatissement de la prédiction. % elif i == 1: On constate la présence d'un voisin au lendemain complètement atypique avec un pic en début de journée (courbe en vert à gauche), et d'un autre phénomène semblable avec la courbe rouge sur le graphe de droite. Ajouté au fait que le lendemain à prévoir est lui-même un jour "hors norme", cela montre l'impossibilité de bien prévoir une courbe en utilisant l'algorithme à voisins. % else: On peut réappliquer les mêmes remarques qu'auparavant sur les boxplots fonctionnels : lendemains de voisins atypiques, courbe à prévoir elle-même légèrement "hors norme". % endif -----r par(mfrow=c(1,2)) plotRelVar(data, f_np1); title(paste("StdDev p1 day",i_np)) plotRelVar(data, f_p1); title(paste("StdDev p1 day",i_p)) plotRelVar(data, f_np2); title(paste("StdDev p2 day",i_np)) plotRelVar(data, f_p2); title(paste("StdDev p2 day",i_p)) # Variabilité globale en rouge ; sur les voisins (+ lendemains) en noir ----- % if i == 0: Ces graphes viennent confirmer l'impression visuelle après observation des filaments. En effet, la variabilité globale en rouge (écart-type heure par heure sur l'ensemble des couples "aujourd'hui/lendemain"du passé) devrait rester nettement au-dessus de la variabilité locale, calculée respectivement sur un voisinage d'une soixantaine de jours (pour p1) et d'une dizaine de jours (pour p2). Or on constate que ce n'est pas du tout le cas sur la période "lendemain", sauf en partie pour p2 le jour 4 $-$ mais ce n'est pas suffisant. % elif i == 1: Comme précédemment les variabilités locales et globales sont confondues dans les parties droites des graphes $-$ sauf pour la version "locale" sur le jour "facile"; mais cette bonne propriété n'est pas suffisante si l'on ne trouve pas les bons poids à appliquer. % else: Cette fois la situation idéale est observée : la variabilité globale est nettement au-dessus de la variabilité locale. Bien que cela ne suffise pas à obtenir de bonnes prédictions de forme, on constate au moins l'amélioration dans la prédiction du niveau. % endif -----r par(mfrow=c(1,2)) plotSimils(p1, i_np); title(paste("Weights p1 day",i_np)) plotSimils(p1, i_p); title(paste("Weights p1 day",i_p)) plotSimils(p2, i_np); title(paste("Weights p2 day",i_np)) plotSimils(p2, i_p); title(paste("Weights p2 day",i_p)) ----- % if i == 0: Les poids se concentrent près de 0 dans le cas "non local" (p1), et se répartissent assez uniformément dans [ 0, 0.2 ] dans le cas "local" (p2). C'est ce que l'on souhaite observer pour éviter d'effectuer une simple moyenne. % elif i == 1: En comparaison avec le pragraphe précédent on retrouve le même (bon) comportement des poids pour la version "non locale". En revanche la fenêtre optimisée est trop grande sur le jour "facile" pour la méthode "locale" (voir affichage ci-dessous) : il en résulte des poids tous semblables autour de 0.084, l'algorithme effectue donc une moyenne simple $-$ expliquant pourquoi les courbes mauve et bleue sont très proches sur le graphe d'erreurs. % else: Concernant les poids en revanche, deux cas a priori mauvais se cumulent : * les poids dans le cas "non local" ne sont pas assez concentrés autour de 0, menant à un lissage trop fort $-$ comme observé sur les graphes des courbes réalisées/prévues ; * les poids dans le cas "local" sont trop semblables (à cause de la trop grande fenêtre optimisée par validation croisée, cf. ci-dessous), résultant encore en une moyenne simple $-$ mais sur moins de jours, plus proches du jour courant. % endif -----r # Fenêtres sélectionnées dans ]0,7] : # "non-local" 2 premières lignes, "local" ensuite p1$getParams(i_np)$window p1$getParams(i_p)$window p2$getParams(i_np)$window p2$getParams(i_p)$window % endfor ----- ${"##"} Bilan Nos algorithmes à voisins ne sont pas adaptés à ce jeu de données où la forme varie considérablement d'un jour à l'autre. Plus généralement cette décorrélation de forme rend ardue la tâche de prévision pour toute autre méthode $-$ du moins, nous ne savons pas comment procéder pour parvenir à une bonne précision. Toutefois, un espoir reste permis par exemple en aggréger les courbes spatialement (sur plusieurs stations situées dans la même agglomération ou dans une même zone).