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