reports/report.gj

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