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