'update'
[talweg.git] / reports / report.gj
CommitLineData
63ff1ecb 1-----
ff5df8e3 2<h2>Introduction</h2>
63ff1ecb 3
aa059de7
BA
4J'ai fait quelques essais dans deux configurations pour la méthode "Neighbors"
5(la seule dont on a parlé, incorporant désormais la "variante Bruno/Michel").
6
445e7bbc
BA
7 * avec simtype="mix" et raccordement "Neighbors" (p1) dans le cas "non local", i.e. on va
8 chercher des voisins n'importe où du moment qu'ils correspondent au premier élément d'un
9 couple de deux jours consécutifs sans valeurs manquantes.
10 * avec simtype="endo" + raccordement "Neighbors" (p2) puis "none" (p3, moyenne simple) + raccordement
11 "Zero" (sans ajustement) dans le cas "local" : voisins de même niveau de pollution et même saison.
63ff1ecb 12
ee8b1b4e 13J'ai systématiquement comparé à une approche naïve : la moyenne des lendemains des jours
445e7bbc
BA
14"similaires" dans tout le passé (p4), ainsi qu'à la persistence (p5) -- reproduisant le jour courant ou
15allant chercher le futur similaire une semaine avant (argument "same_day").
63ff1ecb
BA
16
17Ensuite j'affiche les erreurs, quelques courbes prévues/mesurées, quelques filaments puis les
18histogrammes de quelques poids. Concernant les graphes de filaments, la moitié gauche du graphe
19correspond aux jours similaires au jour courant, tandis que la moitié droite affiche les
20lendemains : ce sont donc les voisinages tels qu'utilisés dans l'algorithme.
21
22<%
23list_titles = ['Pollution par chauffage', 'Pollution par épandage', 'Semaine non polluée']
24list_indices = ['indices_ch', 'indices_ep', 'indices_np']
25%>
63ff1ecb 26-----r
63ff1ecb
BA
27library(talweg)
28
d09b09b0
BA
29P = ${P} #instant de prévision
30H = ${H} #horizon (en heures)
31
63ff1ecb
BA
32ts_data = read.csv(system.file("extdata","pm10_mesures_H_loc_report.csv",package="talweg"))
33exo_data = read.csv(system.file("extdata","meteo_extra_noNAs.csv",package="talweg"))
ee8b1b4e
BA
34# NOTE: 'GMT' because DST gaps are filled and multiple values merged in above dataset.
35# Prediction from P+1 to P+H included.
36data = getData(ts_data, exo_data, input_tz = "GMT", working_tz="GMT", predict_at=P)
63ff1ecb
BA
37
38indices_ch = seq(as.Date("2015-01-18"),as.Date("2015-01-24"),"days")
39indices_ep = seq(as.Date("2015-03-15"),as.Date("2015-03-21"),"days")
40indices_np = seq(as.Date("2015-04-26"),as.Date("2015-05-02"),"days")
d4841a3f 41
ff5df8e3 42% for i in range(3):
63ff1ecb 43-----
ff5df8e3 44<h2 style="color:blue;font-size:2em">${list_titles[i]}</h2>
63ff1ecb 45-----r
445e7bbc 46p1 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Neighbors", horizon=H,
aa059de7 47 simtype="mix", local=FALSE)
445e7bbc 48p2 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Neighbors", horizon=H,
aa059de7 49 simtype="endo", local=TRUE)
445e7bbc
BA
50p3 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Zero", horizon=H,
51 simtype="none", local=TRUE)
52p4 = computeForecast(data, ${list_indices[i]}, "Average", "Zero", horizon=H)
53p5 = computeForecast(data, ${list_indices[i]}, "Persistence", "Zero", horizon=H,
aa059de7 54 same_day=${'TRUE' if loop.index < 2 else 'FALSE'})
63ff1ecb 55-----r
445e7bbc
BA
56e1 = computeError(data, p_nz_mf, H)
57e2 = computeError(data, p_nz_mfl, H)
58e3 = computeError(data, p_a, H)
59e4 = computeError(data, p_p, H)
63ff1ecb 60options(repr.plot.width=9, repr.plot.height=7)
aa059de7 61plotError(list(e_n, e_p, e_a, e_l), cols=c(1,2,colors()[258], 4))
63ff1ecb 62
445e7bbc
BA
63# noir: Neighbors non-local (p1), bleu: Neighbors local endo (p2), mauve: Neighbors local none (p3),
64# vert: moyenne (p4), rouge: persistence (p5)
63ff1ecb 65
aa059de7
BA
66i_np = which.min(e_n$abs$indices)
67i_p = which.max(e_n$abs$indices)
63ff1ecb
BA
68-----r
69options(repr.plot.width=9, repr.plot.height=4)
70par(mfrow=c(1,2))
71
445e7bbc
BA
72plotPredReal(data, p1, i_np); title(paste("PredReal p1 day",i_np))
73plotPredReal(data, p1, i_p); title(paste("PredReal p1 day",i_p))
63ff1ecb 74
445e7bbc
BA
75plotPredReal(data, p2, i_np); title(paste("PredReal p2 day",i_np))
76plotPredReal(data, p2, i_p); title(paste("PredReal p2 day",i_p))
63ff1ecb 77
445e7bbc
BA
78plotPredReal(data, p3, i_np); title(paste("PredReal p3 day",i_np))
79plotPredReal(data, p3, i_p); title(paste("PredReal p3 day",i_p))
63ff1ecb 80
ff5df8e3 81# Bleu: prévue, noir: réalisée
63ff1ecb
BA
82-----r
83par(mfrow=c(1,2))
445e7bbc
BA
84f_np1 = computeFilaments(data, p1, i_np, plot=TRUE); title(paste("Filaments p1 day",i_np))
85f_p1 = computeFilaments(data, p1, i_p, plot=TRUE); title(paste("Filaments p1 day",i_p))
63ff1ecb 86
445e7bbc
BA
87f_np2 = computeFilaments(data, p2, i_np, plot=TRUE); title(paste("Filaments p2 day",i_np))
88f_p2 = computeFilaments(data, p2, i_p, plot=TRUE); title(paste("Filaments p2 day",i_p))
63ff1ecb
BA
89-----r
90par(mfrow=c(1,2))
445e7bbc
BA
91plotFilamentsBox(data, f_np1); title(paste("FilBox p1 day",i_np))
92plotFilamentsBox(data, f_p1); title(paste("FilBox p1 day",i_p))
63ff1ecb 93
445e7bbc 94# Too few neighbors in the local case for this plot
63ff1ecb
BA
95-----r
96par(mfrow=c(1,2))
445e7bbc
BA
97plotRelVar(data, f_np1); title(paste("StdDev p1 day",i_np))
98plotRelVar(data, f_p1); title(paste("StdDev p1 day",i_p))
63ff1ecb 99
445e7bbc
BA
100plotRelVar(data, f_np2); title(paste("StdDev p2 day",i_np))
101plotRelVar(data, f_p2); title(paste("StdDev p2 day",i_p))
63ff1ecb 102
ff5df8e3 103# Variabilité globale en rouge ; sur les 60 voisins (+ lendemains) en noir
63ff1ecb
BA
104-----r
105par(mfrow=c(1,2))
445e7bbc
BA
106plotSimils(p1, i_np); title(paste("Weights p1 day",i_np))
107plotSimils(p1, i_p); title(paste("Weights p1 day",i_p))
63ff1ecb 108
445e7bbc
BA
109plotSimils(p2, i_np); title(paste("Weights p2 day",i_np))
110plotSimils(p2, i_p); title(paste("Weights p2 day",i_p))
63ff1ecb 111
ff5df8e3 112# - pollué à gauche, + pollué à droite
63ff1ecb 113-----r
aa059de7 114# Fenêtres sélectionnées dans ]0,7] / non-loc à gauche, loc à droite
445e7bbc
BA
115p1$getParams(i_np)$window
116p1$getParams(i_p)$window
63ff1ecb 117
445e7bbc
BA
118p2$getParams(i_np)$window
119p2$getParams(i_p)$window
63ff1ecb 120% endfor