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 i_np = which.min(e1$abs$indices)
  67 i_p = which.max(e1$abs$indices)
  68 -----r
  69 options(repr.plot.width=9, repr.plot.height=4)
  70 par(mfrow=c(1,2))
  71
  72 plotPredReal(data, p1, i_np); title(paste("PredReal p1 day",i_np))
  73 plotPredReal(data, p1, i_p); title(paste("PredReal p1 day",i_p))
  74
  75 plotPredReal(data, p2, i_np); title(paste("PredReal p2 day",i_np))
  76 plotPredReal(data, p2, i_p); title(paste("PredReal p2 day",i_p))
  77
  78 plotPredReal(data, p3, i_np); title(paste("PredReal p3 day",i_np))
  79 plotPredReal(data, p3, i_p); title(paste("PredReal p3 day",i_p))
  80
  81 # Bleu: prévue, noir: réalisée
  82 -----r
  83 par(mfrow=c(1,2))
  84 f_np1 = computeFilaments(data, p1, i_np, plot=TRUE); title(paste("Filaments p1 day",i_np))
  85 f_p1 = computeFilaments(data, p1, i_p, plot=TRUE); title(paste("Filaments p1 day",i_p))
  86
  87 f_np2 = computeFilaments(data, p2, i_np, plot=TRUE); title(paste("Filaments p2 day",i_np))
  88 f_p2 = computeFilaments(data, p2, i_p, plot=TRUE); title(paste("Filaments p2 day",i_p))
  89 -----r
  90 par(mfrow=c(1,2))
  91 plotFilamentsBox(data, f_np1); title(paste("FilBox p1 day",i_np))
  92 plotFilamentsBox(data, f_p1); title(paste("FilBox p1 day",i_p))
  93
  94 # Too few neighbors in the local case for this plot
  95 -----r
  96 par(mfrow=c(1,2))
  97 plotRelVar(data, f_np1); title(paste("StdDev p1 day",i_np))
  98 plotRelVar(data, f_p1); title(paste("StdDev p1 day",i_p))
  99
 100 plotRelVar(data, f_np2); title(paste("StdDev p2 day",i_np))
 101 plotRelVar(data, f_p2); title(paste("StdDev p2 day",i_p))
 102
 103 # Variabilité globale en rouge ; sur les 60 voisins (+ lendemains) en noir
 104 -----r
 105 par(mfrow=c(1,2))
 106 plotSimils(p1, i_np); title(paste("Weights p1 day",i_np))
 107 plotSimils(p1, i_p); title(paste("Weights p1 day",i_p))
 108
 109 plotSimils(p2, i_np); title(paste("Weights p2 day",i_np))
 110 plotSimils(p2, i_p); title(paste("Weights p2 day",i_p))
 111
 112 # - pollué à gauche, + pollué à droite
 113 -----r
 114 # Fenêtres sélectionnées dans ]0,7] / non-loc à gauche, loc à droite
 115 p1$getParams(i_np)$window
 116 p1$getParams(i_p)$window
 117
 118 p2$getParams(i_np)$window
 119 p2$getParams(i_p)$window
 120 % endfor