'update'

[talweg.git] / reports / report.gj

-----
<h2>Introduction</h2>

J'ai fait quelques essais pour la méthode "Neighbors", qui incorpore désormais la "variante Bruno/Michel".

 * avec simtype="mix" et raccordement "Neighbors" (p1) 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="endo" + raccordement "Neighbors" (p2) puis "none" (p3, moyenne simple) + raccordement
   "Zero" (sans ajustement) dans le cas "local" : voisins de même niveau de pollution et même saison.

J'ai systématiquement comparé à une approche naïve : la moyenne des lendemains des jours
"similaires" dans tout le passé (p4), ainsi qu'à la persistence (p5) -- reproduisant le jour courant ou
allant chercher le futur similaire une semaine avant (argument "same_day").

Ensuite j'affiche les erreurs, quelques courbes prévues/mesurées, quelques filaments puis les
histogrammes de quelques poids. Concernant les graphes de filaments, la moitié gauche du graphe
correspond aux jours similaires au jour courant, tandis que la moitié droite affiche les
lendemains : 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} #instant de prévision
H = ${H} #horizon (en heures)

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, input_tz = "GMT", working_tz="GMT", predict_at=P)

indices_ch = seq(as.Date("2015-01-18"),as.Date("2015-01-24"),"days")
indices_ep = seq(as.Date("2015-03-15"),as.Date("2015-03-21"),"days")
indices_np = seq(as.Date("2015-04-26"),as.Date("2015-05-02"),"days")

% for i in range(3):
-----
<h2 style="color:blue;font-size:2em">${list_titles[i]}</h2>
-----r
p1 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Neighbors", horizon=H,
	simtype="mix", local=FALSE)
p2 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Neighbors", horizon=H,
	simtype="endo", local=TRUE)
p3 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Zero", horizon=H,
	simtype="none", local=TRUE)
p4 = computeForecast(data, ${list_indices[i]}, "Average", "Zero", horizon=H)
p5 = computeForecast(data, ${list_indices[i]}, "Persistence", "Zero", horizon=H,
	same_day=${'TRUE' if loop.index < 2 else 'FALSE'})
-----r
e1 = computeError(data, p1, H)
e2 = computeError(data, p2, H)
e3 = computeError(data, p3, H)
e4 = computeError(data, p4, H)
e5 = computeError(data, p5, H)
options(repr.plot.width=9, repr.plot.height=7)
plotError(list(e1, e5, e4, e2, e3), cols=c(1,2,colors()[258],4,6))

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

i_np = which.min(e1$abs$indices)
i_p = which.max(e1$abs$indices)
-----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))

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

# Bleu: prévue, noir: réalisée
-----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))
-----r
par(mfrow=c(1,2))
plotFilamentsBox(data, f_np1); title(paste("FilBox p1 day",i_np))
plotFilamentsBox(data, f_p1); title(paste("FilBox p1 day",i_p))

# Too few neighbors in the local case for this plot
-----r
par(mfrow=c(1,2))
plotRelVar(data, f_np1); title(paste("StdDev p1 day",i_np))
plotRelVar(data, f_p1); title(paste("StdDev p1 day",i_p))

plotRelVar(data, f_np2); title(paste("StdDev p2 day",i_np))
plotRelVar(data, f_p2); title(paste("StdDev p2 day",i_p))

# Variabilité globale en rouge ; sur les 60 voisins (+ lendemains) en noir
-----r
par(mfrow=c(1,2))
plotSimils(p1, i_np); title(paste("Weights p1 day",i_np))
plotSimils(p1, i_p); title(paste("Weights p1 day",i_p))

plotSimils(p2, i_np); title(paste("Weights p2 day",i_np))
plotSimils(p2, i_p); title(paste("Weights p2 day",i_p))

# - pollué à gauche, + pollué à droite
-----r
# Fenêtres sélectionnées dans ]0,7] / non-loc à gauche, loc à droite
p1$getParams(i_np)$window
p1$getParams(i_p)$window

p2$getParams(i_np)$window
p2$getParams(i_p)$window
% endfor