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