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