reports/report.gj

   1 -----
   2 <h2>Introduction</h2>
   3
   4 J'ai fait quelques essais dans différentes configurations pour la méthode "Neighbors"
   5 (la seule dont on a parlé) et sa variante récente appelée pour l'instant "Neighbors2",
   6 avec simtype="mix" : deux types de similarités prises en compte, puis multiplication des poids.
   7 Pour Neighbors on prédit le saut (par la moyenne pondérée des sauts passés), et pour Neighbors2
   8 on n'effectue aucun raccordement (prévision directe).
   9
  10 J'ai systématiquement comparé à une approche naïve : la moyenne des lendemains des jours
  11 "similaires" dans tout le passé, ainsi qu'à la persistence -- reproduisant le jour courant ou
  12 allant chercher le futur similaire une semaine avant.
  13
  14 Ensuite j'affiche les erreurs, quelques courbes prévues/mesurées, quelques filaments puis les
  15 histogrammes de quelques poids. Concernant les graphes de filaments, la moitié gauche du graphe
  16 correspond aux jours similaires au jour courant, tandis que la moitié droite affiche les
  17 lendemains : ce sont donc les voisinages tels qu'utilisés dans l'algorithme.
  18
  19 <%
  20 list_titles = ['Pollution par chauffage', 'Pollution par épandage', 'Semaine non polluée']
  21 list_indices = ['indices_ch', 'indices_ep', 'indices_np']
  22 %>
  23 -----r
  24 library(talweg)
  25
  26 P = ${P} #instant de prévision
  27 H = ${H} #horizon (en heures)
  28
  29 ts_data = read.csv(system.file("extdata","pm10_mesures_H_loc_report.csv",package="talweg"))
  30 exo_data = read.csv(system.file("extdata","meteo_extra_noNAs.csv",package="talweg"))
  31 # NOTE: 'GMT' because DST gaps are filled and multiple values merged in above dataset.
  32 # Prediction from P+1 to P+H included.
  33 data = getData(ts_data, exo_data, input_tz = "GMT", working_tz="GMT", predict_at=P)
  34
  35 indices_ch = seq(as.Date("2015-01-18"),as.Date("2015-01-24"),"days")
  36 indices_ep = seq(as.Date("2015-03-15"),as.Date("2015-03-21"),"days")
  37 indices_np = seq(as.Date("2015-04-26"),as.Date("2015-05-02"),"days")
  38
  39 % for i in range(3):
  40 -----
  41 <h2 style="color:blue;font-size:2em">${list_titles[i]}</h2>
  42 -----r
  43 p_nn = computeForecast(data, ${list_indices[i]}, "Neighbors", "Neighbors", horizon=H)
  44 p_nn2 = computeForecast(data, ${list_indices[i]}, "Neighbors2", "Zero", horizon=H)
  45 p_az = computeForecast(data, ${list_indices[i]}, "Average", "Zero", horizon=H)
  46 p_pz = computeForecast(data, ${list_indices[i]}, "Persistence", "Zero", horizon=H, same_day=${'TRUE' if loop.index < 2 else 'FALSE'})
  47 -----r
  48 e_nn = computeError(data, p_nn, H)
  49 e_nn2 = computeError(data, p_nn2, H)
  50 e_az = computeError(data, p_az, H)
  51 e_pz = computeError(data, p_pz, H)
  52 options(repr.plot.width=9, repr.plot.height=7)
  53 plotError(list(e_nn, e_pz, e_az, e_nn2), cols=c(1,2,colors()[258], 4))
  54
  55 # Noir: Neighbors, bleu: Neighbors2, vert: moyenne, rouge: persistence
  56
  57 i_np = which.min(e_nn$abs$indices)
  58 i_p = which.max(e_nn$abs$indices)
  59 -----r
  60 options(repr.plot.width=9, repr.plot.height=4)
  61 par(mfrow=c(1,2))
  62
  63 plotPredReal(data, p_nn, i_np); title(paste("PredReal nn day",i_np))
  64 plotPredReal(data, p_nn2, i_p); title(paste("PredReal nn day",i_p))
  65
  66 plotPredReal(data, p_nn2, i_np); title(paste("PredReal nn2 day",i_np))
  67 plotPredReal(data, p_nn2, i_p); title(paste("PredReal nn2 day",i_p))
  68
  69 plotPredReal(data, p_az, i_np); title(paste("PredReal az day",i_np))
  70 plotPredReal(data, p_az, i_p); title(paste("PredReal az day",i_p))
  71
  72 # Bleu: prévue, noir: réalisée
  73 -----r
  74 par(mfrow=c(1,2))
  75 f_np = computeFilaments(data, p_nn, i_np, plot=TRUE); title(paste("Filaments nn day",i_np))
  76 f_p = computeFilaments(data, p_nn, i_p, plot=TRUE); title(paste("Filaments nn day",i_p))
  77
  78 f_np2 = computeFilaments(data, p_nn2, i_np, plot=TRUE); title(paste("Filaments nn2 day",i_np))
  79 f_p2 = computeFilaments(data, p_nn2, i_p, plot=TRUE); title(paste("Filaments nn2 day",i_p))
  80 -----r
  81 par(mfrow=c(1,2))
  82 plotFilamentsBox(data, f_np); title(paste("FilBox nn day",i_np))
  83 plotFilamentsBox(data, f_p); title(paste("FilBox nn day",i_p))
  84
  85 # Generally too few neighbors:
  86 #plotFilamentsBox(data, f_np2); title(paste("FilBox nn2 day",i_np))
  87 #plotFilamentsBox(data, f_p2); title(paste("FilBox nn2 day",i_p))
  88 -----r
  89 par(mfrow=c(1,2))
  90 plotRelVar(data, f_np); title(paste("StdDev nn day",i_np))
  91 plotRelVar(data, f_p); title(paste("StdDev nn day",i_p))
  92
  93 plotRelVar(data, f_np2); title(paste("StdDev nn2 day",i_np))
  94 plotRelVar(data, f_p2); title(paste("StdDev nn2 day",i_p))
  95
  96 # Variabilité globale en rouge ; sur les 60 voisins (+ lendemains) en noir
  97 -----r
  98 par(mfrow=c(1,2))
  99 plotSimils(p_nn, i_np); title(paste("Weights nn day",i_np))
 100 plotSimils(p_nn, i_p); title(paste("Weights nn day",i_p))
 101
 102 plotSimils(p_nn2, i_np); title(paste("Weights nn2 day",i_np))
 103 plotSimils(p_nn2, i_p); title(paste("Weights nn2 day",i_p))
 104
 105 # - pollué à gauche, + pollué à droite
 106 -----r
 107 # Fenêtres sélectionnées dans ]0,7] / nn à gauche, nn2 à droite
 108 p_nn$getParams(i_np)$window
 109 p_nn$getParams(i_p)$window
 110
 111 p_nn2$getParams(i_np)$window
 112 p_nn2$getParams(i_p)$window
 113 % endfor