reports/report.gj

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