4 "cell_type": "markdown",
7 "## Package R \"vortex\"\n",
9 "using Vectorial exOgenous variables to foRecast Time-sErieX.\n",
11 "Ce package permet de prévoir des courbes de PM10 (par exemple), en se basant sur l'historique des valeurs mais aussi des variables exogènes (par exemple la météo).\n",
13 "Fonctions principales :\n",
15 " * <code>getData</code> : charge un jeu de données en mémoire\n",
16 " * <code>getForecast</code> : prédit les lendemains aux indices demandés\n",
18 "Diverses méthodes permettent ensuite d'analyser les performances : <code>getError</code>, <code>plotXYZ</code> : voir la section \"see also\" dans <code>?plotError</code>."
23 "execution_count": null,
29 "#Chargement de la librairie (après compilation, \"R CMD INSTALL ppmfun/\")\n",
34 "cell_type": "markdown",
37 "Note : sur la base de nos dernières expériences, on considère que \n",
39 " * on ne touche pas à la fenêtre obtenue par la fonction <code>optimize</code> ;\n",
40 " * on oublie la méthode consistant à prédire forme et niveau de manière complètement déconnectée : il faut relier les deux.\n",
42 "### Acquisition des données\n",
44 "Compte-tenu de la nature hétérogène des données utilisées $-$ fonctionnelles pour les PM10, vectorielles pour les variables exogènes $-$, celles-ci sont encapsulées (comme des listes) dans un objet de type *Data*. En interne, la $i^{eme}$ cellule correspondant aux données disponibles au $i^{eme}$ jour à l'heure $H$ de prédiction choisie (1h00, 8h00 ou 14h00) : c'est-à-dire les valeurs des PM10 de $H-24h$ à $H-1h$, ainsi que les variables météo prédites pour la période de 1h à 0h du jour courant (sauf si on prédit à 0h : on prend alors les valeurs mesurées de la veille).\n",
46 "Méthodes d'un objet de classe \"Data\" : elles prennent comme argument \"index\", qui est un index entier ; mais une fonction de conversion existe : <code>dateIndexToInteger</code>.\n",
48 " * <code>getTime</code> : suite des date+heure\n",
49 " * <code>getCenteredSerie</code> : série centrée\n",
50 " * <code>getLevel</code> : niveau\n",
51 " * <code>getSerie</code> : série *non* centrée\n",
52 " * <code>getExoHat</code> : variables exogènes prévues\n",
53 " * <code>getExoDm1</code> : variables exogènes mesurées la veille\n",
60 "execution_count": null,
66 "# Voir ?getData pour les arguments\n",
67 "data = getData(ts_data=\"data/pm10_mesures_H_loc.csv\", exo_data=\"data/meteo_extra_noNAs.csv\",\n",
68 " input_tz = \"Europe/Paris\", working_tz=\"Europe/Paris\", predict_at=\"07\")\n",
69 "data$getLevel(10) #niveau du jour 10\n",
70 "data$getExoHat(17) #météo prévue pour le jour 18"
74 "cell_type": "markdown",
79 "Deux types de prévisions du prochain bloc de $24h$ sont à distinguer :\n",
81 " * prévision de la forme (centrée) ;\n",
82 " * prévision du saut d'une fin de série au début de la suivante.\n",
84 "Il faut ainsi préciser à la fois une méthode de prévision de forme (\"Average\", \"Persistence\" et \"Neighbors\" sont implémentées), et une méthode de prédiction de saut (\"Zero\", \"Persistence\" ou \"Neighbors\"). On détaille surtout la méthode à voisins ci-après, les autres étant des approches naïves que l'on peut considérer comme des références à améliorer.\n",
86 " 1. **Préparation des données** : fenêtrage si demandé (paramètre \"memory\"), recherche des paires de jours consécutifs sans valeurs manquantes.\n",
87 " 2. **Optimisation des paramètres d'échelle** : via la fonction <code>optimize</code> minimisant la somme des 45 dernières erreurs jounalières RMSE, sur des jours similaires.\n",
88 " 3. **Prédiction finale** : une fois le (ou les, si \"simtype\" vaut \"mix\") paramètre d'échelle $h$ déterminé, les similarités sont évaluées sur les variables exogènes et/ou endogènes, sous la forme $s(i,j) = \\mbox{exp}\\left(-\\frac{\\mbox{dist}^2(i,j)}{h^2}\\right)$. La formule indiquée plus haut dans le rapport est alors appliquée.\n",
90 "Détail technique : la sortie de la méthode <code>getForecast</code> est un objet de type Forecast, encapsulant les séries prévues ainsi que tous les paramètres optimisés par la méthode \"Neighbors\". Fonctions disponibles (argument \"index\" comme pour les fonctions sur Data) :\n",
92 " * <code>getSerie</code> : série prévue (sans les information de temps)\n",
93 " * <code>getParams</code> : liste des paramètres (poids, fenêtre, ...)\n",
94 " * <code>getIndexInData</code> : indice du jour où s'effectue la prévision relativement au jeu de données\n",
96 "### Calcul des erreurs\n",
98 "Pour chacun des instants à prévoir jusqu'à minuit du jour courant (ou avant : argument *horizon*), on calcule l'erreur moyenne sur tous les instants similaires du passé. Deux types d'erreurs sont considérées :\n",
100 " * l'erreur \"abs\" égale à la valeur absolue moyenne entre la mesure et la prédiction ;\n",
101 " * l'erreur \"MAPE\" égale à l'erreur absolue normalisée par la mesure.\n",
103 "### Expériences numériques"
108 "execution_count": null,
114 "options(repr.plot.width=9, repr.plot.height=3)"
119 "execution_count": null,
125 "p_endo = predictPM10(data, 2200, 2230, 0,0, \"Neighbors\", \"Neighbors\", simtype=\"endo\")"
130 "execution_count": null,
136 "p_exo = predictPM10(data, 2200, 2230, 0,0, \"Neighbors\", \"Neighbors\", simtype=\"exo\")"
141 "execution_count": null,
147 "p_mix = predictPM10(data, 2200, 2230, 0,0, \"Neighbors\", \"Neighbors\", simtype=\"mix\")"
152 "execution_count": null,
158 "p = list(p_endo, p_exo, p_mix)"
163 "execution_count": null,
169 "yrange_MAPE = range(p_mix$errors$MAPE, p_endo$errors$MAPE, p_exo$errors$MAPE)\n",
170 "yrange_abs = range(p_mix$errors$abs, p_endo$errors$abs, p_exo$errors$abs)\n",
171 "yrange_RMSE = range(p_mix$errors$RMSE, p_endo$errors$RMSE, p_exo$errors$RMSE)\n",
172 "ranges = list(yrange_MAPE,yrange_abs,yrange_RMSE)\n",
174 "par(mfrow=c(1,3))\n",
175 "titles = paste(\"Erreur\",c(\"MAPE\",\"abs\",\"RMSE\"))\n",
176 "for (i in 1:3) #error type (MAPE,abs,RMSE)\n",
178 " for (j in 1:3) #model (mix,endo,exo)\n",
180 " xlab = if (j>1) \"\" else \"Temps\"\n",
181 " ylab = if (j>1) \"\" else \"Erreur\"\n",
182 " main = if (j>1) \"\" else titles[i]\n",
183 " plot(p[[j]]$errors[[i]], type=\"l\", col=j, main=main, xlab=xlab, ylab=ylab, ylim=ranges[[i]])\n",
191 "cell_type": "markdown",
194 "Ne tenir compte que des similarités sur les variables exogènes semble conduire à l'erreur la plus faible."
199 "execution_count": null,
205 "p_nn = predictPM10(data, 2200, 2230, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE)"
210 "execution_count": null,
216 "p_np = predictPM10(data, 2200, 2230, 0, 0, \"Neighbors\", \"Persistence\", sameSeaon=TRUE)"
221 "execution_count": null,
227 "p_nz = predictPM10(data, 2200, 2230, 0, 0, \"Neighbors\", \"Zero\", sameSeaon=TRUE)"
232 "execution_count": null,
238 "p_pp = predictPM10(data, 2200, 2230, 0, 0, \"Persistence\", \"Persistence\")"
243 "execution_count": null,
249 "p_pz = predictPM10(data, 2200, 2230, 0, 0, \"Persistence\", \"Zero\")"
254 "execution_count": null,
260 "p = list(p_nn, p_np, p_nz, p_pp, p_pz)"
265 "execution_count": null,
271 "yrange_MAPE = range(p_nn$errors$MAPE, p_nz$errors$MAPE, p_np$errors$MAPE, p_pp$errors$MAPE, p_pz$errors$MAPE)\n",
272 "yrange_abs = range(p_nn$errors$abs, p_nz$errors$abs, p_np$errors$abs, p_pp$errors$abs, p_pz$errors$abs)\n",
273 "yrange_RMSE = range(p_nn$errors$RMSE, p_nz$errors$RMSE, p_np$errors$RMSE, p_pp$errors$RMSE, p_pz$errors$RMSE)\n",
274 "ranges = list(yrange_MAPE,yrange_abs,yrange_RMSE)\n",
276 "par(mfrow=c(1,3))\n",
277 "titles = paste(\"Erreur\",c(\"MAPE\",\"abs\",\"RMSE\"))\n",
278 "for (i in 1:3) #error type (MAPE,abs,RMSE)\n",
280 " for (j in 1:5) #model (nn,np,nz,pp,pz)\n",
282 " xlab = if (j>1) \"\" else \"Temps\"\n",
283 " ylab = if (j>1) \"\" else \"Erreur\"\n",
284 " main = if (j>1) \"\" else titles[i]\n",
285 " plot(p[[j]]$errors[[i]], type=\"l\", col=j, main=titles[i], xlab=xlab, ylab=ylab, ylim=ranges[[i]])\n",
292 "p = list(p_nn_epandage, p_nn_nonpollue, p_nn_chauffage)\n",
293 "forecasts_2 = lapply(1:length(data), function(index) ( if (is.na(p[[2]]$forecasts[[index]][1])) rep(NA,24) else p[[2]]$forecasts[[index]]$pred ) )\n",
294 "e1 = getErrors(data, forecasts_1, 17)\n",
296 "e = list(e1,e2,e3)\n",
297 "yrange_MAPE = range(e1$MAPE, e2$MAPE, e3$MAPE)\n",
298 "yrange_abs = range(e1$abs, e2$abs, e3$abs)\n",
299 "yrange_RMSE = range(e1$RMSE, e2$RMSE, e3$RMSE)\n",
300 "ranges = list(yrange_MAPE,yrange_abs,yrange_RMSE)\n",
302 "par(mfrow=c(1,3))\n",
303 "titles = paste(\"Erreur\",c(\"MAPE\",\"abs\",\"RMSE\"))\n",
304 "for (i in 1:3) #error type (MAPE,abs,RMSE)\n",
306 " for (j in 1:3) #model (nn,np,nz,pp,pz)\n",
308 " xlab = if (j>1) \"\" else \"Temps\"\n",
309 " ylab = if (j>1) \"\" else \"Erreur\"\n",
310 " main = if (j>1) \"\" else titles[i]\n",
311 " plot(e[[j]][[i]], type=\"l\", col=j, main=titles[i], xlab=xlab, ylab=ylab, ylim=ranges[[i]])\n",
317 "par(mfrow=c(1,2))\n",
318 "#p[[i]]$forecasts[[index]]\n",
319 "#futurs des blocs du passé pour le jour 2290 ::\n",
320 "futurs = lapply(1:length(p[[1]]$forecasts[[2290]]$indices),\n",
321 " function(index) ( data[[ p[[1]]$forecasts[[2290]]$indices[index]+1 ]]$pm10 ) )\n",
322 "#vrai futur (en rouge), vrai jour (en noir)\n",
323 "r_min = min( sapply( 1:length(futurs), function(index) ( min(futurs[[index]] ) ) ) )\n",
324 "r_max = max( sapply( 1:length(futurs), function(index) ( max(futurs[[index]] ) ) ) )\n",
325 "for (i in 1:length(futurs))\n",
327 " plot(futurs[[i]], col=1, ylim=c(r_min,r_max), type=\"l\")\n",
328 " if (i<length(futurs))\n",
332 "plot(data[[2290]]$pm10, ylim=c(r_min, r_max), col=1, type=\"l\")\n",
334 "plot(data[[2291]]$pm10, ylim=c(r_min, r_max), col=2, type=\"l\")\n"
338 "cell_type": "markdown",
341 "Meilleurs results: nn et nz (np moins bon)"
346 "execution_count": null,
352 "#%%TODO: analyse sur les trois périodes indiquées par Michel ; simtype==\"exo\" par defaut\n",
353 "#16/03/2015 2288\n",
354 "p_nn_epandage = predictPM10(data, 2287, 2293, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"mix\")"
359 "execution_count": null,
365 "p_nn_epandage = predictPM10(data, 2287, 2293, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"mix\")"
370 "execution_count": null,
376 "options(repr.plot.width=9, repr.plot.height=6)\n",
377 "plot(p_nn_epandage$errors$abs, type=\"l\", col=1, main=\"Erreur absolue\", xlab=\"Temps\",\n",
378 " ylab=\"Erreur\", ylim=range(p_nn_epandage$errors$abs))"
383 "execution_count": null,
389 "#19/01/2015 2232\n",
390 "p_nn_chauffage = predictPM10(data, 2231, 2237, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE)"
395 "execution_count": null,
401 "#23/02/2015 2267\n",
402 "p_nn_nonpollue = predictPM10(data, 2266, 2272, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"mix\")"
407 "execution_count": null,
413 "plot(p_nn_nonpollue$errors$abs, type=\"l\", col=2, main=\"Erreur absolue\", xlab=\"Temps\",\n",
414 " ylab=\"Erreur\", ylim=range(p_nn_nonpollue$errors$abs))"
419 "execution_count": null,
426 "data = getData(\"local\", \"7h\")\n",
427 "p_nn_epandage = predictPM10(data, 2287, 2293, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"endo\")\n",
428 "p_nn_nonpollue = predictPM10(data, 2266, 2272, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"endo\")\n",
429 "p_nn_chauffage = predictPM10(data, 2231, 2237, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"endo\")"
434 "execution_count": null,
440 "p = list(p_nn_epandage, p_nn_nonpollue, p_nn_chauffage)\n",
441 "#yrange_MAPE = range(p[[1]]$errors$MAPE, p[[2]]$errors$MAPE, p[[3]]$errors$MAPE)\n",
442 "#yrange_abs = range(p[[1]]$errors$abs, p[[2]]$errors$abs, p[[3]]$errors$abs)\n",
443 "#yrange_RMSE = range(p[[1]]$errors$RMSE, p[[2]]$errors$RMSE, p[[3]]$errors$RMSE)\n",
444 "#ranges = list(yrange_MAPE,yrange_abs,yrange_RMSE)\n",
445 "print(p[[1]]$forecasts[[2290]])"
450 "execution_count": null,
456 "par(mfrow=c(1,3))\n",
457 "titles = paste(\"Erreur\",c(\"MAPE\",\"abs\",\"RMSE\"))\n",
458 "for (i in 1:3) #error type (MAPE,abs,RMSE)\n",
460 " for (j in 1:5) #model (nn,np,nz,pp,pz)\n",
462 " xlab = if (j>1) \"\" else \"Temps\"\n",
463 " ylab = if (j>1) \"\" else \"Erreur\"\n",
464 " main = if (j>1) \"\" else titles[i]\n",
465 " plot(p[[j]]$errors[[i]], type=\"l\", col=j, main=titles[i], xlab=xlab, ylab=ylab, ylim=ranges[[i]])\n",
473 "cell_type": "markdown",
489 "codemirror_mode": "r",
490 "file_extension": ".r",
491 "mimetype": "text/x-r-source",
493 "pygments_lexer": "r",