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