update reports, last fixes, cleanup

[talweg.git] / reports / Experiments.gj

-----
# Résultats numériques

Cette partie montre les résultats obtenus avec des variantes de l'algorithme décrit à la
section 4, en utilisant le package présenté au chapitre précédent. Cet algorithme est
systématiquement comparé à deux approches naïves :

 * la moyenne des lendemains des jours "similaires" dans tout le passé, c'est-à-dire
prédiction = moyenne de tous les mardis passés si le jour courant est un lundi.
 * la persistence, reproduisant le jour courant ou allant chercher le lendemain de la
dernière journée "similaire" (même principe que ci-dessus ; argument "same\_day").

Concernant l'algorithme principal à voisins, deux variantes sont comparées dans cette
partie :

 * avec simtype="mix" et raccordement "Neighbors" dans le cas "non local", i.e. on va
chercher des voisins n'importe où du moment qu'ils correspondent au premier élément d'un
couple de deux jours consécutifs sans valeurs manquantes.
 * avec simtype="none" (moyenne simple) et raccordement=NULL (aucun ajustement après
moyenne des courbes) dans le cas "local" : voisins de même niveau de pollution et même
saison.

Pour chaque période retenue $-$ chauffage, épandage, semaine non polluée $-$ les erreurs
de prédiction sont d'abord affichées, puis quelques graphes de courbes réalisées/prévues
(sur le jour "en moyenne le plus facile" à gauche, et "en moyenne le plus difficile" à
droite). Ensuite plusieurs types de graphes apportant des précisions sur la nature et la
difficulté du problème viennent compléter ces premières courbes. Concernant les graphes
de filaments, la moitié droite du graphe correspond aux jours similaires au jour courant,
tandis que la moitié gauche affiche les jours précédents : ce sont donc les voisinages
tels qu'utilisés dans l'algorithme.
<%
list_titles = ['Pollution par chauffage','Pollution par épandage','Semaine non polluée']
list_indices = ['indices_ch', 'indices_ep', 'indices_np']
%>
-----r
library(talweg)

P = ${P} #première heure de prévision
H = ${H} #dernière heure de prévision

ts_data = read.csv(system.file("extdata","pm10_mesures_H_loc_report.csv",
	package="talweg"))
exo_data = read.csv(system.file("extdata","meteo_extra_noNAs.csv",
	package="talweg"))
# NOTE: 'GMT' because DST gaps are filled and multiple values merged in
# above dataset. Prediction from P+1 to P+H included.
data = getData(ts_data, exo_data)

indices_ch = seq(as.Date("2015-01-19"),as.Date("2015-01-25"),"days")
indices_ep = seq(as.Date("2015-03-16"),as.Date("2015-03-22"),"days")
indices_np = seq(as.Date("2015-04-27"),as.Date("2015-05-03"),"days")
% for i in range(3):
-----
##<h2 style="color:blue;font-size:2em">${list_titles[i]}</h2>
${"##"} ${list_titles[i]}
-----r
p1 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Neighbors", predict_from=P,
	horizon=H, simtype="mix", local=FALSE)
p2 = computeForecast(data, ${list_indices[i]}, "Neighbors", NULL, predict_from=P,
	horizon=H, simtype="none", local=TRUE)
p3 = computeForecast(data, ${list_indices[i]}, "Average", "Zero", predict_from=P,
	horizon=H)
p4 = computeForecast(data, ${list_indices[i]}, "Persistence", "Zero", predict_from=P,
	horizon=H, same_day=${'TRUE' if loop.index < 2 else 'FALSE'})
-----r
e1 = computeError(data, p1, P, H)
e2 = computeError(data, p2, P, H)
e3 = computeError(data, p3, P, H)
e4 = computeError(data, p4, P, H)
options(repr.plot.width=9, repr.plot.height=7)
plotError(list(e1, e4, e3, e2), cols=c(1,2,colors()[258],4))

# noir: Neighbors non-local (p1), bleu: Neighbors local (p2),
# vert: moyenne (p3), rouge: persistence (p4)

sum_p23 = e2$abs$indices + e3$abs$indices
i_np = which.min(sum_p23) #indice de jour "facile"
i_p = which.max(sum_p23) #indice de jour "difficile"
-----
% if i == 0:
L'erreur absolue -- en haut à gauche -- reste modérée pour les meilleurs modèles
(variantes à voisins), ne dépassant 10 que deux jours. Les deux modèles naïfs ont des
erreurs similaires sauf sur la période "difficile" (jours 4 à 6), sur laquelle on gagne
donc à chercher des jours similaires pour effectuer la prévision.
Le MAPE reste en général inférieur à 35% pour les meilleurs méthodes.
% elif i == 1:
Le modèle à voisins avec contrainte de localité obtient ici les meilleurs résultats, son
erreur étant clairement en dessous des autres à partir du jour 4 (graphe en haut à
droite). Le MAPE jour après jour est du même ordre que précédemment pour cette méthode
(35%, graphe en bas à droite) sauf un jour sur lequel le MAPE explose.
% else:
Dans ce cas plus favorable les intensité des erreurs absolues ont clairement diminué :
elles sont souvent en dessous de 5. En revanche le MAPE moyen reste en général au-delà de
20%. Comme dans le cas de l'épandage on constate une croissance globale de la courbe
journalière d'erreur absolue moyenne (en haut à gauche) -- sauf pour la méthode à voisins
"locale" ; ceci peut être dû au fait que l'on ajuste le niveau du jour à prédire en le
recollant sur la dernière valeur observée (sauf pour "Neighbors local").
% endif
-----r
options(repr.plot.width=9, repr.plot.height=4)
par(mfrow=c(1,2))

plotPredReal(data, p1, i_np); title(paste("PredReal p1 day",i_np))
plotPredReal(data, p1, i_p); title(paste("PredReal p1 day",i_p))

plotPredReal(data, p2, i_np); title(paste("PredReal p2 day",i_np))
plotPredReal(data, p2, i_p); title(paste("PredReal p2 day",i_p))

# Bleu : prévue ; noir : réalisée (confondues jusqu'à predict_from-1)
-----
% if i == 0:
Le jour "facile à prévoir", à gauche, se décompose en deux modes : un léger vers 10h
(7+3), puis un beaucoup plus marqué vers 19h (7+12). Ces deux modes sont retrouvés par
les deux variantes de l'algorithme à voisins, bien que l'amplitude soit mal prédite.
Concernant le jour "difficile à prévoir" (à droite) il y a deux pics en tout début et
toute fin de journée (à 9h et 23h), qui ne sont pas du tout anticipés par les méthodes ;
la grande amplitude de ces pics explique alors l'intensité de l'erreur observée.
% elif i == 1:
Dans le cas d'un jour "facile" à prédire $-$ à gauche $-$ la forme est plutôt bien
retrouvée, ainsi que le niveau moyen pour la méthode sans contrainte de localité
(dans l'autre, l'algorithme a probablement écarté trop de voisins potentiels).
Concernant le jour "difficile" à droite, non seulement la forme n'est pas anticipée mais
surtout le niveau prédit est largement supérieur au niveau de pollution observé -- dans
une moindre mesure toutefois pour la variante "locale".
% else:
L'impression visuelle est plutôt mauvaise dans ce cas, mais les écart étant minimes les
erreurs au final ne sont pas très importantes. De plus deux des quatres graphes sont
satisfaisants (en haut à droite et en bas à gauche : forme + niveau acceptables.
% endif
-----r
par(mfrow=c(1,2))

f_np1 = computeFilaments(data, p1, i_np, plot=TRUE)
title(paste("Filaments p1 day",i_np))

f_p1 = computeFilaments(data, p1, i_p, plot=TRUE)
title(paste("Filaments p1 day",i_p))

f_np2 = computeFilaments(data, p2, i_np, plot=TRUE)
title(paste("Filaments p2 day",i_np))

f_p2 = computeFilaments(data, p2, i_p, plot=TRUE)
title(paste("Filaments p2 day",i_p))
-----
% if i == 0:
Les voisins du jour courant (période de 24h allant de 8h à 7h le lendemain) sont affichés
avec un trait d'autant plus sombre qu'ils sont proches. On constate dans le cas non
contraint (en haut) une grande variabilité des lendemains, très nette sur le graphe en
haut à droite. Ceci indique une faible corrélation entre la forme d'une courbe sur une
période de 24h et la forme sur les 24h suivantes ; **cette observation est la source des
difficultés rencontrées par l'algorithme sur ce jeu de données.**
% elif i == 1:
Les observations sont les mêmes qu'au paragraphe précédent : trop de variabilité des
voisins (et ce même le jour précédent).
% else:
Les graphes de filaments ont encore la même allure, avec une assez grande variabilité
observée. Cette observation est cependant trompeuse, comme l'indique plus bas le graphe
de variabilité relative.
% endif
-----r
par(mfrow=c(1,2))

plotFilamentsBox(data, f_np1, predict_from=P)
title(paste("FilBox p1 day",i_np))

plotFilamentsBox(data, f_p1, predict_from=P)
title(paste("FilBox p1 day",i_p))

# En pointillés la courbe du jour courant (à prédire) + précédent
-----
% if i == 0:
Sur cette boxplot fonctionnelle (voir la fonction fboxplot() du package R "rainbow") on
constate essentiellement deux choses : le lendemain d'un voisin "normal" peut se révéler
être une courbe atypique, fort éloignée de ce que l'on souhaite prédire (courbes bleue et
rouge à gauche) ; et, dans le cas d'une courbe à prédire atypique (à droite) la plupart
des voisins sont trop éloignés de la forme à prédire et forcent ainsi un aplatissement de
la prédiction.
% elif i == 1:
Concernant le jour "difficile" on constate la présence de voisins au lendemains
complètement atypiques avec un pic en début de journée (courbes en vert et rouge à
droite). Ajouté au fait que le jour à prévoir est lui-même "hors norme", cela montre
l'impossibilité de bien prévoir une courbe en utilisant l'algorithme à voisins.
% else:
On peut réappliquer les mêmes remarques qu'auparavant sur les boxplots fonctionnels :
voisins atypiques, courbe à prévoir elle-même légèrement "hors norme".
% endif
-----r
par(mfrow=c(1,2))

plotRelVar(data, f_np1, predict_from=P)
title(paste("StdDev p1 day",i_np))

plotRelVar(data, f_p1, predict_from=P)
title(paste("StdDev p1 day",i_p))

plotRelVar(data, f_np2, predict_from=P)
title(paste("StdDev p2 day",i_np))

plotRelVar(data, f_p2, predict_from=P)
title(paste("StdDev p2 day",i_p))

# Variabilité globale en rouge ; sur les voisins en noir
-----
% if i == 0:
Ces graphes viennent confirmer l'impression visuelle après observation des filaments. En
effet, la variabilité globale en rouge (écart-type heure par heure sur l'ensemble des
couples "hier/aujourd'hui" du passé) devrait rester nettement au-dessus de la
variabilité locale, calculée respectivement sur un voisinage d'une soixantaine de jours
(pour p1) et d'une dizaine de jours (pour p2). Or ce n'est pas du tout le cas sur la
moitié droite, sauf pour le jour "facile" avec l'algorithme "local".
% elif i == 1:
Comme précédemment les variabilités locales et globales sont trop proches dans les
parties droites des graphes pour le jour "difficile". L'allure des graphes est
raisonnable ppour l'autre jour, qui est d'ailleurs bien prédit.
% else:
Cette fois la situation idéale est observée : la variabilité globale est nettement
au-dessus de la variabilité locale. Bien que cela ne suffise pas à obtenir de bonnes
prédictions de forme, on constate au moins l'amélioration dans la prédiction du niveau.
% endif
-----r
plotSimils(p1, i_np)
title(paste("Weights p1 day",i_np))

plotSimils(p1, i_p)
title(paste("Weights p1 day",i_p))

# Poids < 1/N à gauche, >= 1/N à droite ; jour facile en haut, difficile en bas
-----
% if i == 0:
Les poids se concentrent près de 0 : c'est ce que l'on souhaite observer pour éviter
d'effectuer une simple moyenne.
% elif i == 1:
En comparaison avec le paragraphe précédent on retrouve le même (bon) comportement des
poids pour la version "non locale".
% else:
Concernant les poids en revanche, deux cas a priori mauvais se cumulent : ...
% endif
-----r
options(digits=2)

p1$getParams(i_np)$window
p1$getParams(i_p)$window

# Fenêtres sélectionnées dans ]0,7]
% endfor
-----
${"##"} Bilan

Nos algorithmes à voisins donnent de meilleurs résultats que les approches naïves
(persistence, moyenne sur tout le jeu de données). Les erreurs restent cependant assez
élevées, notamment en terme de MAPE. Une possible poste d'amélioration consisterait à
aggréger les courbes spatialement (sur plusieurs stations situées dans la même
agglomération ou dans une même zone).