| 1 | # Package R "talweg" |
| 2 | |
| 3 | Le package $-$ Time-series sAmpLes forecasted With ExoGenous variables $-$ contient le |
| 4 | code permettant de (re)lancer les expériences numériques décrites dans cette partie et la |
| 5 | suivante. Les fonctions principales sont respectivement |
| 6 | |
| 7 | * **getData()** pour construire un objet R contenant les données à partir de fichiers |
| 8 | CSV (extraits de bases de données). Le format choisi en R est une classe R6 (du package |
| 9 | du même nom) exposant en particulier les méthodes *getSerie(i)* et *getExo(i)* qui |
| 10 | renvoient respectivement la $i^{eme}$ série de 24h et les variables exogènes (mesurées) |
| 11 | correspondantes. Voir ?Data pour plus d'information, une fois le package chargé. |
| 12 | * **computeForecast()** pour calculer des prédictions sur une certaine plage temporelle |
| 13 | contenue dans *data <- getData(...)* |
| 14 | * **computeError()** pour évaluer les erreurs commises par différentes méthodes. |
| 15 | |
| 16 | Le package contient en outre diverses fonctions graphiques *plotXXX()*, utilisées dans la |
| 17 | partie suivante. |
| 18 | -----r |
| 19 | # Chargement de la librairie (après compilation, "R CMD INSTALL .") |
| 20 | library(talweg) |
| 21 | |
| 22 | # Acquisition des données (depuis les fichiers CSV) |
| 23 | ts_data <- read.csv(system.file("extdata","pm10_mesures_H_loc.csv", |
| 24 | package="talweg")) |
| 25 | exo_data <- read.csv(system.file("extdata","meteo_extra_noNAs.csv", |
| 26 | package="talweg")) |
| 27 | data <- getData(ts_data, exo_data, input_tz="GMT", |
| 28 | date_format="%d/%m/%Y %H:%M", working_tz="GMT", predict_at=7, limit=120) |
| 29 | # Plus de détails à la section 1 ci-après. |
| 30 | |
| 31 | # Prédiction de 10 courbes (jours 102 à 111) |
| 32 | pred <- computeForecast(data, 101:110, "Persistence", "Zero", memory=50, |
| 33 | horizon=12, ncores=1) |
| 34 | # Plus de détails à la section 2 ci-après. |
| 35 | |
| 36 | # Calcul des erreurs (sur un horizon arbitraire <= horizon de prédiction) |
| 37 | err <- computeError(data, pred, horizon=6) |
| 38 | # Plus de détails à la section 3 ci-après. |
| 39 | |
| 40 | # Puis voir ?plotError et les autres plot dans le paragraphe 'seealso' |
| 41 | ----- |
| 42 | ## getData() |
| 43 | |
| 44 | Les arguments de cette fonction sont, dans l'ordre : |
| 45 | |
| 46 | 1. **ts_data** : séries temporelles (fichier CSV avec entête ou data.frame) ; la |
| 47 | première colonne contient les heures, la seconde les valeurs. |
| 48 | 2. **exo_data** : variables exogènes (fichier CSV avec entête ou data.frame) ; la |
| 49 | première colonne contient les jours, les $m$ suivantes les variables mesurées pour ce |
| 50 | jour, et les $m$ dernières les variables prédites pour ce même jour. Dans notre cas $m=4$ |
| 51 | : pression, température, gradient de température, vitesse du vent. |
| 52 | 3. **input_tz** : zone horaire pour ts_data (défaut : "GMT"). |
| 53 | 4. **date_format** : format des heures dans ts_data (défaut : "%d/%m/%Y %H:%M", format |
| 54 | du fichier transmis par Michel). |
| 55 | 5. **working_tz** : zone horaire dans laquelle on souhaite travailler avec les données |
| 56 | (défaut : "GMT"). |
| 57 | 6. **predict_at** : heure à laquelle s'effectue la prévision $-$ et donc dernière heure |
| 58 | d'un bloc de 24h, relativement à working_tz. data`$`getSerie(3) renvoit ainsi les 24 |
| 59 | valeurs de 8h à 7h pour le $3^{eme}$ bloc de 24h présent dans le jeu de données. |
| 60 | -----r |
| 61 | print(data) |
| 62 | #?Data |
| 63 | ----- |
| 64 | ## computeForecast() |
| 65 | |
| 66 | Les arguments de cette fonction sont, dans l'ordre : |
| 67 | |
| 68 | 1. **data** : le jeu de données renvoyé par getData() |
| 69 | 2. **indices** : l'ensemble de jours dont on veut prévoir les "lendemains" (prochains |
| 70 | blocs de 24h) ; peut être donnée sous forme d'un vecteur de dates ou d'entiers |
| 71 | (correspondants aux numéros des jours). |
| 72 | 3. **forecaster** : le nom du prédicteur principal à utiliser ; voir ?computeForecast |
| 73 | 4. **pjump** : le nom du prédicteur de saut d'une série à l'autre ; voir |
| 74 | ?computeForecast |
| 75 | 5. **memory** : le nombre de jours à prendre en compte dans le passé pour chaque |
| 76 | prévision (par défaut : Inf, c'est-à-dire tout l'historique pris en compte). |
| 77 | 6. **horizon** : le nombre d'heures à prédire ; par défaut "data`$`getStdHorizon()", |
| 78 | c'est-à-dire le nombre d'heures restantes à partir de l'instant de prévision + 1 jusqu'à |
| 79 | minuit (17 pour predict_at=7 par exemple). |
| 80 | 7. **ncores** : le nombre de processus parallèles (utiliser 1 pour une exécution |
| 81 | séquentielle) |
| 82 | -----r |
| 83 | print(pred) |
| 84 | #?computeForecast |
| 85 | ----- |
| 86 | ## computeError() |
| 87 | |
| 88 | Les arguments de cette fonction sont, dans l'ordre : |
| 89 | |
| 90 | 1. **data** : le jeu de données renvoyé par getData() |
| 91 | 2. **pred** : les prédictions renvoyées par computeForecast() |
| 92 | 3. **horizon** : le nombre d'heures à considérer pour le calcul de l'erreur ; doit être |
| 93 | inférieur ou égal à l'horizon utilisé pour la prédiction (même valeur par défaut : |
| 94 | "data`$`getStdHorizon()") |
| 95 | -----r |
| 96 | summary(err) |
| 97 | summary(err$abs) |
| 98 | summary(err$MAPE) |
| 99 | ----- |
| 100 | ## Graphiques |
| 101 | |
| 102 | Voir ?plotError : les autres fonctions graphiques sont dans la section 'seealso' : |
| 103 | |
| 104 | ‘plotCurves’, ‘plotPredReal’, ‘plotSimils’, ‘plotFbox’, |
| 105 | ‘computeFilaments’, ‘plotFilamentsBox’, ‘plotRelVar’ |
| 106 | |
| 107 | ?plotXXX, etc. |
| 108 | $\clearpage$ |
| 109 | ----- |
| 110 | # Expérimentations |
| 111 | |
| 112 | Cette partie montre les résultats obtenus via des variantes de l'algorithme décrit à la |
| 113 | section 2, en utilisant le package présenté à la section 3. Cet algorithme est |
| 114 | systématiquement comparé à deux approches naïves : |
| 115 | |
| 116 | * la moyenne des lendemains des jours "similaires" dans tout le passé, c'est-à-dire |
| 117 | prédiction = moyenne de tous les mardis passé si le jour courant est un lundi par |
| 118 | exemple. |
| 119 | * la persistence, reproduisant le jour courant ou allant chercher le lendemain de la |
| 120 | dernière journée "similaire" (même principe que ci-dessus ; argument "same\_day"). |
| 121 | |
| 122 | Concernant l'algorithme principal à voisins, trois variantes sont étudiées dans cette |
| 123 | partie : |
| 124 | |
| 125 | * avec simtype="mix" et raccordement "Neighbors" dans le cas "non local", i.e. on va |
| 126 | chercher des voisins n'importe où du moment qu'ils correspondent au premier élément d'un |
| 127 | couple de deux jours consécutifs sans valeurs manquantes. |
| 128 | * avec simtype="endo" + raccordement "Neighbors" puis simtype="none" + raccordement |
| 129 | "Zero" (sans ajustement) dans le cas "local" : voisins de même niveau de pollution et |
| 130 | même saison. |
| 131 | |
| 132 | Pour chaque période retenue $-$ chauffage, épandage, semaine non polluée $-$ les erreurs |
| 133 | de prédiction sont d'abord affichées, puis quelques graphes de courbes réalisées/prévues |
| 134 | (sur le jour "en moyenne le plus facile" à gauche, et "en moyenne le plus difficile" à |
| 135 | droite). Ensuite plusieurs types de graphes apportant des précisions sur la nature et la |
| 136 | difficulté du problème viennent compléter ces premières courbes. Concernant les graphes |
| 137 | de filaments, la moitié gauche du graphe correspond aux jours similaires au jour courant, |
| 138 | tandis que la moitié droite affiche les lendemains : ce sont donc les voisinages tels |
| 139 | qu'utilisés dans l'algorithme. |
| 140 | <% |
| 141 | list_titles = ['Pollution par chauffage','Pollution par épandage','Semaine non polluée'] |
| 142 | list_indices = ['indices_ch', 'indices_ep', 'indices_np'] |
| 143 | %> |
| 144 | -----r |
| 145 | library(talweg) |
| 146 | |
| 147 | P = ${P} #instant de prévision |
| 148 | H = ${H} #horizon (en heures) |
| 149 | |
| 150 | ts_data = read.csv(system.file("extdata","pm10_mesures_H_loc_report.csv", |
| 151 | package="talweg")) |
| 152 | exo_data = read.csv(system.file("extdata","meteo_extra_noNAs.csv", |
| 153 | package="talweg")) |
| 154 | # NOTE: 'GMT' because DST gaps are filled and multiple values merged in |
| 155 | # above dataset. Prediction from P+1 to P+H included. |
| 156 | data = getData(ts_data, exo_data, input_tz = "GMT", working_tz="GMT", |
| 157 | predict_at=P) |
| 158 | |
| 159 | indices_ch = seq(as.Date("2015-01-18"),as.Date("2015-01-24"),"days") |
| 160 | indices_ep = seq(as.Date("2015-03-15"),as.Date("2015-03-21"),"days") |
| 161 | indices_np = seq(as.Date("2015-04-26"),as.Date("2015-05-02"),"days") |
| 162 | |
| 163 | % for i in range(3): |
| 164 | ----- |
| 165 | % #<h2 style="color:blue;font-size:2em">${list_titles[i]}</h2> |
| 166 | ## ${list_titles[i]} |
| 167 | -----r |
| 168 | p1 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Neighbors", horizon=H, |
| 169 | simtype="mix", local=FALSE) |
| 170 | p2 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Neighbors", horizon=H, |
| 171 | simtype="endo", local=TRUE) |
| 172 | p3 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Zero", horizon=H, |
| 173 | simtype="none", local=TRUE) |
| 174 | p4 = computeForecast(data, ${list_indices[i]}, "Average", "Zero", horizon=H) |
| 175 | p5 = computeForecast(data, ${list_indices[i]}, "Persistence", "Zero", horizon=H, |
| 176 | same_day=${'TRUE' if loop.index < 2 else 'FALSE'}) |
| 177 | -----r |
| 178 | e1 = computeError(data, p1, H) |
| 179 | e2 = computeError(data, p2, H) |
| 180 | e3 = computeError(data, p3, H) |
| 181 | e4 = computeError(data, p4, H) |
| 182 | e5 = computeError(data, p5, H) |
| 183 | options(repr.plot.width=9, repr.plot.height=7) |
| 184 | plotError(list(e1, e5, e4, e2, e3), cols=c(1,2,colors()[258],4,6)) |
| 185 | |
| 186 | # noir: Neighbors non-local (p1), bleu: Neighbors local endo (p2), |
| 187 | # mauve: Neighbors local none (p3), vert: moyenne (p4), |
| 188 | # rouge: persistence (p5) |
| 189 | |
| 190 | sum_p123 = e1$abs$indices + e2$abs$indices + e3$abs$indices |
| 191 | i_np = which.min(sum_p123) #indice de (veille de) jour "facile" |
| 192 | i_p = which.max(sum_p123) #indice de (veille de) jour "difficile" |
| 193 | ----- |
| 194 | % if i == 1 |
| 195 | L'erreur absolue dépasse 20 sur 1 à 2 jours suivant les modèles (graphe en haut à |
| 196 | droite). C'est au-delà de ce que l'on aimerait voir (disons +/- 5 environ). Sur cet |
| 197 | exemple le modèle à voisins "contraint" (local=TRUE) utilisant des pondérations basées |
| 198 | sur les similarités de forme (simtype="endo") obtient en moyenne les meilleurs résultats, |
| 199 | avec un MAPE restant en général inférieur à 30% de 8h à 19h (7+1 à 7+12 : graphe en bas à |
| 200 | gauche). |
| 201 | % else if i == 2 |
| 202 | Il est difficile dans ce cas de déterminer une méthode meilleure que les autres : elles |
| 203 | donnent toutes de plutôt mauvais résultats, avec une erreur absolue moyennée sur la |
| 204 | journée dépassant presque toujours 15 (graphe en haut à droite). |
| 205 | % else |
| 206 | Dans ce cas plus favorable les intensité des erreurs absolues ont clairement diminué : |
| 207 | elles restent souvent en dessous de 5. En revanche le MAPE moyen reste au-delà de 20%, et |
| 208 | même souvent plus de 30%. Comme dans le cas de l'épandage on constate une croissance |
| 209 | globale de la courbe journalière d'erreur absolue moyenne (en haut à gauche) ; ceci peut |
| 210 | être dû au fait que l'on ajuste le niveau du jour à prédire en le recollant sur la |
| 211 | dernière valeur observée. |
| 212 | % endif |
| 213 | -----r |
| 214 | options(repr.plot.width=9, repr.plot.height=4) |
| 215 | par(mfrow=c(1,2)) |
| 216 | |
| 217 | plotPredReal(data, p1, i_np); title(paste("PredReal p1 day",i_np)) |
| 218 | plotPredReal(data, p1, i_p); title(paste("PredReal p1 day",i_p)) |
| 219 | |
| 220 | plotPredReal(data, p2, i_np); title(paste("PredReal p2 day",i_np)) |
| 221 | plotPredReal(data, p2, i_p); title(paste("PredReal p2 day",i_p)) |
| 222 | |
| 223 | plotPredReal(data, p3, i_np); title(paste("PredReal p3 day",i_np)) |
| 224 | plotPredReal(data, p3, i_p); title(paste("PredReal p3 day",i_p)) |
| 225 | |
| 226 | # Bleu : prévue ; noir : réalisée |
| 227 | ----- |
| 228 | % if i == 1 |
| 229 | Le jour "facile à prévoir", à gauche, se décompose en deux modes : un léger vers 10h |
| 230 | (7+3), puis un beaucoup plus marqué vers 19h (7+12). Ces deux modes sont retrouvés par |
| 231 | les trois variantes de l'algorithme à voisins, bien que l'amplitude soit mal prédite. |
| 232 | Concernant le jour "difficile à prévoir" il y a deux pics en tout début et toute fin de |
| 233 | journée (à 9h et 23h), qui ne sont pas du tout anticipés par le programme ; la grande |
| 234 | amplitude de ces pics explique alors l'intensité de l'erreur observée. |
| 235 | % else if i == 2 |
| 236 | Dans le cas d'un jour "facile" à prédire $-$ à gauche $-$ la forme est plus ou moins |
| 237 | retrouvée, mais le niveau moyen est trop bas (courbe en bleu). Concernant le jour |
| 238 | "difficile" à droite, non seulement la forme n'est pas anticipée mais surtout le niveau |
| 239 | prédit est très inférieur au niveau de pollution observé. Comme on le voit ci-dessous |
| 240 | cela découle d'un manque de voisins au comportement similaire. |
| 241 | % else |
| 242 | La forme est raisonnablement retrouvée pour les méthodes "locales", l'autre version |
| 243 | lissant trop les prédictions. Le biais reste cependant important, surtout en fin de |
| 244 | journée sur le jour "difficile". |
| 245 | % endif |
| 246 | -----r |
| 247 | par(mfrow=c(1,2)) |
| 248 | f_np1 = computeFilaments(data, p1, i_np, plot=TRUE) |
| 249 | title(paste("Filaments p1 day",i_np)) |
| 250 | f_p1 = computeFilaments(data, p1, i_p, plot=TRUE) |
| 251 | title(paste("Filaments p1 day",i_p)) |
| 252 | |
| 253 | f_np2 = computeFilaments(data, p2, i_np, plot=TRUE) |
| 254 | title(paste("Filaments p2 day",i_np)) |
| 255 | f_p2 = computeFilaments(data, p2, i_p, plot=TRUE) |
| 256 | title(paste("Filaments p2 day",i_p)) |
| 257 | ----- |
| 258 | % if i == 1 |
| 259 | Les voisins du jour courant (période de 24h allant de 8h à 7h le lendemain) sont affichés |
| 260 | avec un trait d'autant plus sombre qu'ils sont proches. On constate dans le cas non |
| 261 | contraint (en haut) une grande variabilité des lendemains, très nette sur le graphe en |
| 262 | haut à droite. Ceci indique une faible corrélation entre la forme d'une courbe sur une |
| 263 | période de 24h et la forme sur les 24h suivantes ; **cette observation est la source des |
| 264 | difficultés rencontrées par l'algorithme sur ce jeu de données.** |
| 265 | % else if i == 2 |
| 266 | Les observations sont les mêmes qu'au paragraphe précédent : trop de variabilité des |
| 267 | lendemains (et même des voisins du jour courant). |
| 268 | % else |
| 269 | Les graphes de filaments ont encore la même allure, avec une assez grande variabilité |
| 270 | observée. Cette observation est cependant trompeuse, comme l'indique plus bas le graphe |
| 271 | de variabilité relative. |
| 272 | % endif |
| 273 | -----r |
| 274 | par(mfrow=c(1,2)) |
| 275 | plotFilamentsBox(data, f_np1); title(paste("FilBox p1 day",i_np)) |
| 276 | plotFilamentsBox(data, f_p1); title(paste("FilBox p1 day",i_p)) |
| 277 | |
| 278 | # En pointillés la courbe du jour courant + lendemain (à prédire) |
| 279 | ----- |
| 280 | % if i == 1 |
| 281 | Sur cette boxplot fonctionnelle (voir la fonction fboxplot() du package R "rainbow") l'on |
| 282 | constate essentiellement deux choses : le lendemain d'un voisin "normal" peut se révéler |
| 283 | être une courbe atypique, fort éloignée de ce que l'on souhaite prédire (courbes bleue et |
| 284 | rouge à gauche) ; et, dans le cas d'une courbe à prédire atypique (à droite) la plupart |
| 285 | des voisins sont trop éloignés de la forme à prédire et forcent ainsi un aplatissement de |
| 286 | la prédiction. |
| 287 | % else if i == 2 |
| 288 | On constate la présence d'un voisin au lendemain complètement atypique avec un pic en |
| 289 | début de journée (courbe en vert à gauche), et d'un autre phénomène semblable avec la |
| 290 | courbe rouge sur le graphe de droite. Ajouté au fait que le lendemain à prévoir est |
| 291 | lui-même un jour "hors norme", cela montre l'impossibilité de bien prévoir une courbe en |
| 292 | utilisant l'algorithme à voisins. |
| 293 | % else |
| 294 | On peut réappliquer les mêmes remarques qu'auparavant sur les boxplots fonctionnels : |
| 295 | lendemains de voisins atypiques, courbe à prévoir elle-même légèrement "hors norme". |
| 296 | % endif |
| 297 | -----r |
| 298 | par(mfrow=c(1,2)) |
| 299 | plotRelVar(data, f_np1); title(paste("StdDev p1 day",i_np)) |
| 300 | plotRelVar(data, f_p1); title(paste("StdDev p1 day",i_p)) |
| 301 | |
| 302 | plotRelVar(data, f_np2); title(paste("StdDev p2 day",i_np)) |
| 303 | plotRelVar(data, f_p2); title(paste("StdDev p2 day",i_p)) |
| 304 | |
| 305 | # Variabilité globale en rouge ; sur les voisins (+ lendemains) en noir |
| 306 | ----- |
| 307 | % if i == 1 |
| 308 | Ces graphes viennent confirmer l'impression visuelle après observation des filaments. En |
| 309 | effet, la variabilité globale en rouge (écart-type heure par heure sur l'ensemble des |
| 310 | couples "aujourd'hui/lendemain"du passé) devrait rester nettement au-dessus de la |
| 311 | variabilité locale, calculée respectivement sur un voisinage d'une soixantaine de jours |
| 312 | (pour p1) et d'une dizaine de jours (pour p2). Or on constate que ce n'est pas du tout le |
| 313 | cas sur la période "lendemain", sauf en partie pour p2 le jour 4 $-$ mais ce n'est pas |
| 314 | suffisant. |
| 315 | % else if i == 2 |
| 316 | Comme précédemment les variabilités locales et globales sont confondues dans les parties |
| 317 | droites des graphes $-$ sauf pour la version "locale" sur le jour "facile"; mais cette |
| 318 | bonne propriété n'est pas suffisante si l'on ne trouve pas les bons poids à appliquer. |
| 319 | % else |
| 320 | Cette fois la situation idéale est observée : la variabilité globale est nettement |
| 321 | au-dessus de la variabilité locale. Bien que cela ne suffise pas à obtenir de bonnes |
| 322 | prédictions de forme, on constate au moins l'amélioration dans la prédiction du niveau. |
| 323 | % endif |
| 324 | -----r |
| 325 | par(mfrow=c(1,2)) |
| 326 | plotSimils(p1, i_np); title(paste("Weights p1 day",i_np)) |
| 327 | plotSimils(p1, i_p); title(paste("Weights p1 day",i_p)) |
| 328 | |
| 329 | plotSimils(p2, i_np); title(paste("Weights p2 day",i_np)) |
| 330 | plotSimils(p2, i_p); title(paste("Weights p2 day",i_p)) |
| 331 | ----- |
| 332 | % if i == 1 |
| 333 | Les poids se concentrent près de 0 dans le cas "non local" (p1), et se répartissent assez |
| 334 | uniformément dans [ 0, 0.2 ] dans le cas "local" (p2). C'est ce que l'on souhaite |
| 335 | observer pour éviter d'effectuer une simple moyenne. |
| 336 | % else if i == 2 |
| 337 | En comparaison avec le pragraphe précédent on retrouve le même (bon) comportement des |
| 338 | poids pour la version "non locale". En revanche la fenêtre optimisée est trop grande sur |
| 339 | le jour "facile" pour la méthode "locale" (voir affichage ci-dessous) : il en résulte des |
| 340 | poids tous semblables autour de 0.084, l'algorithme effectue donc une moyenne simple $-$ |
| 341 | expliquant pourquoi les courbes mauve et bleue sont très proches sur le graphe d'erreurs. |
| 342 | % else |
| 343 | Concernant les poids en revanche, deux cas a priori mauvais se cumulent : |
| 344 | |
| 345 | * les poids dans le cas "non local" ne sont pas assez concentrés autour de 0, menant à |
| 346 | un lissage trop fort $-$ comme observé sur les graphes des courbes réalisées/prévues ; |
| 347 | * les poids dans le cas "local" sont trop semblables (à cause de la trop grande fenêtre |
| 348 | optimisée par validation croisée, cf. ci-dessous), résultant encore en une moyenne simple |
| 349 | $-$ mais sur moins de jours, plus proches du jour courant. |
| 350 | % endif |
| 351 | -----r |
| 352 | # Fenêtres sélectionnées dans ]0,7] : |
| 353 | # "non-local" 2 premières lignes, "local" ensuite |
| 354 | p1$getParams(i_np)$window |
| 355 | p1$getParams(i_p)$window |
| 356 | |
| 357 | p2$getParams(i_np)$window |
| 358 | p2$getParams(i_p)$window |
| 359 | % endfor |
| 360 | ----- |
| 361 | ## Bilan |
| 362 | |
| 363 | Nos algorithmes à voisins ne sont pas adaptés à ce jeu de données où la forme varie |
| 364 | considérablement d'un jour à l'autre. Plus généralement cette décorrélation de forme rend |
| 365 | ardue la tâche de prévision pour toute autre méthode $-$ du moins, nous ne savons pas |
| 366 | comment procéder pour parvenir à une bonne précision. |
| 367 | |
| 368 | Toutefois, un espoir reste permis par exemple en aggréger les courbes spatialement (sur |
| 369 | plusieurs stations situées dans la même agglomération ou dans une même zone). |