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