'update'
[talweg.git] / reports / Experiments.gj
CommitLineData
8eafefbc 1-----
4d376294 2# Résultats numériques
b6233fa6 3
882ae735
BA
4Cette partie montre les résultats obtenus avec des variantes de l'algorithme décrit à la
5section 4, en utilisant le package présenté au chapitre précédent. Cet algorithme est
6systématiquement comparé à deux approches naïves :
63ff1ecb 7
b6233fa6 8 * la moyenne des lendemains des jours "similaires" dans tout le passé, c'est-à-dire
4d376294 9prédiction = moyenne de tous les mardis passés si le jour courant est un lundi.
b6233fa6
BA
10 * la persistence, reproduisant le jour courant ou allant chercher le lendemain de la
11dernière journée "similaire" (même principe que ci-dessus ; argument "same\_day").
12
9b9bb2d4 13Concernant l'algorithme principal à voisins, deux variantes sont comparées dans cette
b6233fa6
BA
14partie :
15
16 * avec simtype="mix" et raccordement "Neighbors" dans le cas "non local", i.e. on va
17chercher des voisins n'importe où du moment qu'ils correspondent au premier élément d'un
18couple de deux jours consécutifs sans valeurs manquantes.
9b9bb2d4
BA
19 * avec simtype="none" (moyenne simple) et raccordement=NULL (aucun ajustement après
20moyenne des courbes) dans le cas "local" : voisins de même niveau de pollution et même
21saison.
b6233fa6
BA
22
23Pour chaque période retenue $-$ chauffage, épandage, semaine non polluée $-$ les erreurs
24de prédiction sont d'abord affichées, puis quelques graphes de courbes réalisées/prévues
25(sur le jour "en moyenne le plus facile" à gauche, et "en moyenne le plus difficile" à
26droite). Ensuite plusieurs types de graphes apportant des précisions sur la nature et la
27difficulté du problème viennent compléter ces premières courbes. Concernant les graphes
882ae735
BA
28de filaments, la moitié droite du graphe correspond aux jours similaires au jour courant,
29tandis que la moitié gauche affiche les jours précédents : ce sont donc les voisinages
30tels qu'utilisés dans l'algorithme.
63ff1ecb 31<%
b6233fa6 32list_titles = ['Pollution par chauffage','Pollution par épandage','Semaine non polluée']
63ff1ecb
BA
33list_indices = ['indices_ch', 'indices_ep', 'indices_np']
34%>
63ff1ecb 35-----r
63ff1ecb
BA
36library(talweg)
37
1e8327df
BA
38P = ${P} #première heure de prévision
39H = ${H} #dernière heure de prévision
d09b09b0 40
b6233fa6
BA
41ts_data = read.csv(system.file("extdata","pm10_mesures_H_loc_report.csv",
42 package="talweg"))
43exo_data = read.csv(system.file("extdata","meteo_extra_noNAs.csv",
44 package="talweg"))
45# NOTE: 'GMT' because DST gaps are filled and multiple values merged in
46# above dataset. Prediction from P+1 to P+H included.
1e8327df 47data = getData(ts_data, exo_data)
63ff1ecb 48
882ae735
BA
49indices_ch = seq(as.Date("2015-01-19"),as.Date("2015-01-25"),"days")
50indices_ep = seq(as.Date("2015-03-16"),as.Date("2015-03-22"),"days")
51indices_np = seq(as.Date("2015-04-27"),as.Date("2015-05-03"),"days")
ff5df8e3 52% for i in range(3):
63ff1ecb 53-----
8eafefbc
BA
54##<h2 style="color:blue;font-size:2em">${list_titles[i]}</h2>
55${"##"} ${list_titles[i]}
63ff1ecb 56-----r
1e8327df
BA
57p1 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Neighbors", predict_from=P,
58 horizon=H, simtype="mix", local=FALSE)
882ae735 59p2 = computeForecast(data, ${list_indices[i]}, "Neighbors", NULL, predict_from=P,
1e8327df 60 horizon=H, simtype="none", local=TRUE)
9b9bb2d4 61p3 = computeForecast(data, ${list_indices[i]}, "Average", "Zero", predict_from=P,
1e8327df 62 horizon=H)
9b9bb2d4 63p4 = computeForecast(data, ${list_indices[i]}, "Persistence", "Zero", predict_from=P,
1e8327df 64 horizon=H, same_day=${'TRUE' if loop.index < 2 else 'FALSE'})
63ff1ecb 65-----r
1e8327df
BA
66e1 = computeError(data, p1, P, H)
67e2 = computeError(data, p2, P, H)
68e3 = computeError(data, p3, P, H)
69e4 = computeError(data, p4, P, H)
63ff1ecb 70options(repr.plot.width=9, repr.plot.height=7)
9b9bb2d4 71plotError(list(e1, e4, e3, e2), cols=c(1,2,colors()[258],4))
63ff1ecb 72
9b9bb2d4
BA
73# noir: Neighbors non-local (p1), bleu: Neighbors local (p2),
74# vert: moyenne (p3), rouge: persistence (p4)
63ff1ecb 75
882ae735 76sum_p23 = e2$abs$indices + e3$abs$indices
9b9bb2d4
BA
77i_np = which.min(sum_p23) #indice de jour "facile"
78i_p = which.max(sum_p23) #indice de jour "difficile"
b6233fa6 79-----
12119d21 80% if i == 0:
9b9bb2d4
BA
81L'erreur absolue -- en haut à gauche -- reste modérée pour les meilleurs modèles
82(variantes à voisins), ne dépassant 10 que deux jours. Les deux modèles naïfs ont des
83erreurs similaires sauf sur la période "difficile" (jours 4 à 6), sur laquelle on gagne
84donc à chercher des jours similaires pour effectuer la prévision.
85Le MAPE reste en général inférieur à 35% pour les meilleurs méthodes.
12119d21 86% elif i == 1:
9b9bb2d4
BA
87Le modèle à voisins avec contrainte de localité obtient ici les meilleurs résultats, son
88erreur étant clairement en dessous des autres à partir du jour 4 (graphe en haut à
89droite). Le MAPE jour après jour est du même ordre que précédemment pour cette méthode
90(35%, graphe en bas à droite) sauf un jour sur lequel le MAPE explose.
8eafefbc 91% else:
b6233fa6 92Dans ce cas plus favorable les intensité des erreurs absolues ont clairement diminué :
9b9bb2d4
BA
93elles sont souvent en dessous de 5. En revanche le MAPE moyen reste en général au-delà de
9420%. Comme dans le cas de l'épandage on constate une croissance globale de la courbe
95journalière d'erreur absolue moyenne (en haut à gauche) -- sauf pour la méthode à voisins
96"locale" ; ceci peut être dû au fait que l'on ajuste le niveau du jour à prédire en le
97recollant sur la dernière valeur observée (sauf pour "Neighbors local").
b6233fa6 98% endif
63ff1ecb
BA
99-----r
100options(repr.plot.width=9, repr.plot.height=4)
101par(mfrow=c(1,2))
102
445e7bbc
BA
103plotPredReal(data, p1, i_np); title(paste("PredReal p1 day",i_np))
104plotPredReal(data, p1, i_p); title(paste("PredReal p1 day",i_p))
63ff1ecb 105
445e7bbc
BA
106plotPredReal(data, p2, i_np); title(paste("PredReal p2 day",i_np))
107plotPredReal(data, p2, i_p); title(paste("PredReal p2 day",i_p))
63ff1ecb 108
9b9bb2d4 109# Bleu : prévue ; noir : réalisée (confondues jusqu'à predict_from-1)
b6233fa6 110-----
12119d21 111% if i == 0:
b6233fa6
BA
112Le jour "facile à prévoir", à gauche, se décompose en deux modes : un léger vers 10h
113(7+3), puis un beaucoup plus marqué vers 19h (7+12). Ces deux modes sont retrouvés par
9b9bb2d4 114les deux variantes de l'algorithme à voisins, bien que l'amplitude soit mal prédite.
4d376294
BA
115Concernant le jour "difficile à prévoir" (à droite) il y a deux pics en tout début et
116toute fin de journée (à 9h et 23h), qui ne sont pas du tout anticipés par les méthodes ;
117la grande amplitude de ces pics explique alors l'intensité de l'erreur observée.
12119d21 118% elif i == 1:
9b9bb2d4
BA
119Dans le cas d'un jour "facile" à prédire $-$ à gauche $-$ la forme est plutôt bien
120retrouvée, ainsi que le niveau moyen pour la méthode sans contrainte de localité
121(dans l'autre, l'algorithme a probablement écarté trop de voisins potentiels).
122Concernant le jour "difficile" à droite, non seulement la forme n'est pas anticipée mais
123surtout le niveau prédit est largement supérieur au niveau de pollution observé -- dans
124une moindre mesure toutefois pour la variante "locale".
8eafefbc 125% else:
9b9bb2d4
BA
126L'impression visuelle est plutôt mauvaise dans ce cas, mais les écart étant minimes les
127erreurs au final ne sont pas très importantes. De plus deux des quatres graphes sont
128satisfaisants (en haut à droite et en bas à gauche : forme + niveau acceptables.
b6233fa6 129% endif
63ff1ecb
BA
130-----r
131par(mfrow=c(1,2))
9b9bb2d4
BA
132
133f_np1 = computeFilaments(data, p1, i_np, plot=TRUE)
134title(paste("Filaments p1 day",i_np))
135
136f_p1 = computeFilaments(data, p1, i_p, plot=TRUE)
137title(paste("Filaments p1 day",i_p))
138
139f_np2 = computeFilaments(data, p2, i_np, plot=TRUE)
140title(paste("Filaments p2 day",i_np))
141
142f_p2 = computeFilaments(data, p2, i_p, plot=TRUE)
143title(paste("Filaments p2 day",i_p))
b6233fa6 144-----
12119d21 145% if i == 0:
b6233fa6
BA
146Les voisins du jour courant (période de 24h allant de 8h à 7h le lendemain) sont affichés
147avec un trait d'autant plus sombre qu'ils sont proches. On constate dans le cas non
148contraint (en haut) une grande variabilité des lendemains, très nette sur le graphe en
149haut à droite. Ceci indique une faible corrélation entre la forme d'une courbe sur une
150période de 24h et la forme sur les 24h suivantes ; **cette observation est la source des
151difficultés rencontrées par l'algorithme sur ce jeu de données.**
12119d21 152% elif i == 1:
b6233fa6 153Les observations sont les mêmes qu'au paragraphe précédent : trop de variabilité des
9b9bb2d4 154voisins (et ce même le jour précédent).
8eafefbc 155% else:
b6233fa6
BA
156Les graphes de filaments ont encore la même allure, avec une assez grande variabilité
157observée. Cette observation est cependant trompeuse, comme l'indique plus bas le graphe
158de variabilité relative.
159% endif
63ff1ecb
BA
160-----r
161par(mfrow=c(1,2))
63ff1ecb 162
9b9bb2d4
BA
163plotFilamentsBox(data, f_np1, predict_from=P)
164title(paste("FilBox p1 day",i_np))
165
166plotFilamentsBox(data, f_p1, predict_from=P)
167title(paste("FilBox p1 day",i_p))
168
169# En pointillés la courbe du jour courant (à prédire) + précédent
b6233fa6 170-----
12119d21 171% if i == 0:
4d376294 172Sur cette boxplot fonctionnelle (voir la fonction fboxplot() du package R "rainbow") on
b6233fa6
BA
173constate essentiellement deux choses : le lendemain d'un voisin "normal" peut se révéler
174être une courbe atypique, fort éloignée de ce que l'on souhaite prédire (courbes bleue et
175rouge à gauche) ; et, dans le cas d'une courbe à prédire atypique (à droite) la plupart
176des voisins sont trop éloignés de la forme à prédire et forcent ainsi un aplatissement de
177la prédiction.
12119d21 178% elif i == 1:
9b9bb2d4
BA
179Concernant le jour "difficile" on constate la présence de voisins au lendemains
180complètement atypiques avec un pic en début de journée (courbes en vert et rouge à
181droite). Ajouté au fait que le jour à prévoir est lui-même "hors norme", cela montre
182l'impossibilité de bien prévoir une courbe en utilisant l'algorithme à voisins.
8eafefbc 183% else:
b6233fa6 184On peut réappliquer les mêmes remarques qu'auparavant sur les boxplots fonctionnels :
9b9bb2d4 185voisins atypiques, courbe à prévoir elle-même légèrement "hors norme".
b6233fa6 186% endif
63ff1ecb
BA
187-----r
188par(mfrow=c(1,2))
63ff1ecb 189
9b9bb2d4
BA
190plotRelVar(data, f_np1, predict_from=P)
191title(paste("StdDev p1 day",i_np))
63ff1ecb 192
9b9bb2d4
BA
193plotRelVar(data, f_p1, predict_from=P)
194title(paste("StdDev p1 day",i_p))
195
196plotRelVar(data, f_np2, predict_from=P)
197title(paste("StdDev p2 day",i_np))
198
199plotRelVar(data, f_p2, predict_from=P)
200title(paste("StdDev p2 day",i_p))
201
202# Variabilité globale en rouge ; sur les voisins en noir
b6233fa6 203-----
12119d21 204% if i == 0:
b6233fa6
BA
205Ces graphes viennent confirmer l'impression visuelle après observation des filaments. En
206effet, la variabilité globale en rouge (écart-type heure par heure sur l'ensemble des
9b9bb2d4 207couples "hier/aujourd'hui" du passé) devrait rester nettement au-dessus de la
b6233fa6 208variabilité locale, calculée respectivement sur un voisinage d'une soixantaine de jours
9b9bb2d4
BA
209(pour p1) et d'une dizaine de jours (pour p2). Or ce n'est pas du tout le cas sur la
210moitié droite, sauf pour le jour "facile" avec l'algorithme "local".
12119d21 211% elif i == 1:
9b9bb2d4
BA
212Comme précédemment les variabilités locales et globales sont trop proches dans les
213parties droites des graphes pour le jour "difficile". L'allure des graphes est
214raisonnable ppour l'autre jour, qui est d'ailleurs bien prédit.
8eafefbc 215% else:
b6233fa6
BA
216Cette fois la situation idéale est observée : la variabilité globale est nettement
217au-dessus de la variabilité locale. Bien que cela ne suffise pas à obtenir de bonnes
218prédictions de forme, on constate au moins l'amélioration dans la prédiction du niveau.
219% endif
63ff1ecb 220-----r
9b9bb2d4
BA
221plotSimils(p1, i_np)
222title(paste("Weights p1 day",i_np))
63ff1ecb 223
9b9bb2d4
BA
224plotSimils(p1, i_p)
225title(paste("Weights p1 day",i_p))
226
227# Poids < 1/N à gauche, >= 1/N à droite ; jour facile en haut, difficile en bas
b6233fa6 228-----
12119d21 229% if i == 0:
9b9bb2d4
BA
230Les poids se concentrent près de 0 : c'est ce que l'on souhaite observer pour éviter
231d'effectuer une simple moyenne.
12119d21 232% elif i == 1:
9b9bb2d4
BA
233En comparaison avec le paragraphe précédent on retrouve le même (bon) comportement des
234poids pour la version "non locale".
8eafefbc 235% else:
9b9bb2d4 236Concernant les poids en revanche, deux cas a priori mauvais se cumulent : ...
b6233fa6 237% endif
63ff1ecb 238-----r
9b9bb2d4
BA
239options(digits=2)
240
445e7bbc
BA
241p1$getParams(i_np)$window
242p1$getParams(i_p)$window
63ff1ecb 243
9b9bb2d4 244# Fenêtres sélectionnées dans ]0,7]
63ff1ecb 245% endfor
b6233fa6 246-----
8eafefbc 247${"##"} Bilan
b6233fa6 248
9b9bb2d4
BA
249Nos algorithmes à voisins donnent de meilleurs résultats que les approches naïves
250(persistence, moyenne sur tout le jeu de données). Les erreurs restent cependant assez
251élevées, notamment en terme de MAPE. Une possible poste d'amélioration consisterait à
252aggréger les courbes spatialement (sur plusieurs stations situées dans la même
4d376294 253agglomération ou dans une même zone).