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",
67 "evalue": "Error in eval(expr, envir, enclos): could not find function \"getData\"\n",
68 "output_type": "error",
70 "Error in eval(expr, envir, enclos): could not find function \"getData\"\nTraceback:\n"
75 "# Voir ?getData pour les arguments\n",
76 "data = getData(ts_data=\"data/pm10_mesures_H_loc.csv\", exo_data=\"data/meteo_extra_noNAs.csv\",\n",
77 " input_tz = \"Europe/Paris\", working_tz=\"Europe/Paris\", predict_at=\"07\")\n",
78 "data$getLevel(10) #niveau du jour 10\n",
79 "data$getExoHat(17) #météo prévue pour le jour 18"
83 "cell_type": "markdown",
88 "Deux types de prévisions du prochain bloc de $24h$ sont à distinguer :\n",
90 " * prévision de la forme (centrée) ;\n",
91 " * prévision du saut d'une fin de série au début de la suivante.\n",
93 "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",
95 " 1. **Préparation des données** : fenêtrage si demandé (paramètre \"memory\"), recherche des paires de jours consécutifs sans valeurs manquantes.\n",
96 " 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",
97 " 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",
99 "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",
101 " * <code>getSerie</code> : série prévue (sans les information de temps)\n",
102 " * <code>getParams</code> : liste des paramètres (poids, fenêtre, ...)\n",
103 " * <code>getIndexInData</code> : indice du jour où s'effectue la prévision relativement au jeu de données\n",
105 "### Calcul des erreurs\n",
107 "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",
109 " * l'erreur \"abs\" égale à la valeur absolue moyenne entre la mesure et la prédiction ;\n",
110 " * l'erreur \"MAPE\" égale à l'erreur absolue normalisée par la mesure.\n",
112 "### Expériences numériques"
117 "execution_count": null,
123 "options(repr.plot.width=9, repr.plot.height=3)"
128 "execution_count": null,
134 "p_endo = predictPM10(data, 2200, 2230, 0,0, \"Neighbors\", \"Neighbors\", simtype=\"endo\")"
139 "execution_count": null,
145 "p_exo = predictPM10(data, 2200, 2230, 0,0, \"Neighbors\", \"Neighbors\", simtype=\"exo\")"
150 "execution_count": null,
156 "p_mix = predictPM10(data, 2200, 2230, 0,0, \"Neighbors\", \"Neighbors\", simtype=\"mix\")"
161 "execution_count": null,
167 "p = list(p_endo, p_exo, p_mix)"
172 "execution_count": null,
178 "yrange_MAPE = range(p_mix$errors$MAPE, p_endo$errors$MAPE, p_exo$errors$MAPE)\n",
179 "yrange_abs = range(p_mix$errors$abs, p_endo$errors$abs, p_exo$errors$abs)\n",
180 "yrange_RMSE = range(p_mix$errors$RMSE, p_endo$errors$RMSE, p_exo$errors$RMSE)\n",
181 "ranges = list(yrange_MAPE,yrange_abs,yrange_RMSE)\n",
183 "par(mfrow=c(1,3))\n",
184 "titles = paste(\"Erreur\",c(\"MAPE\",\"abs\",\"RMSE\"))\n",
185 "for (i in 1:3) #error type (MAPE,abs,RMSE)\n",
187 " for (j in 1:3) #model (mix,endo,exo)\n",
189 " xlab = if (j>1) \"\" else \"Temps\"\n",
190 " ylab = if (j>1) \"\" else \"Erreur\"\n",
191 " main = if (j>1) \"\" else titles[i]\n",
192 " plot(p[[j]]$errors[[i]], type=\"l\", col=j, main=main, xlab=xlab, ylab=ylab, ylim=ranges[[i]])\n",
200 "cell_type": "markdown",
203 "Ne tenir compte que des similarités sur les variables exogènes semble conduire à l'erreur la plus faible."
208 "execution_count": null,
214 "p_nn = predictPM10(data, 2200, 2230, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE)"
219 "execution_count": null,
225 "p_np = predictPM10(data, 2200, 2230, 0, 0, \"Neighbors\", \"Persistence\", sameSeaon=TRUE)"
230 "execution_count": null,
236 "p_nz = predictPM10(data, 2200, 2230, 0, 0, \"Neighbors\", \"Zero\", sameSeaon=TRUE)"
241 "execution_count": null,
247 "p_pp = predictPM10(data, 2200, 2230, 0, 0, \"Persistence\", \"Persistence\")"
252 "execution_count": null,
258 "p_pz = predictPM10(data, 2200, 2230, 0, 0, \"Persistence\", \"Zero\")"
263 "execution_count": null,
269 "p = list(p_nn, p_np, p_nz, p_pp, p_pz)"
274 "execution_count": null,
280 "yrange_MAPE = range(p_nn$errors$MAPE, p_nz$errors$MAPE, p_np$errors$MAPE, p_pp$errors$MAPE, p_pz$errors$MAPE)\n",
281 "yrange_abs = range(p_nn$errors$abs, p_nz$errors$abs, p_np$errors$abs, p_pp$errors$abs, p_pz$errors$abs)\n",
282 "yrange_RMSE = range(p_nn$errors$RMSE, p_nz$errors$RMSE, p_np$errors$RMSE, p_pp$errors$RMSE, p_pz$errors$RMSE)\n",
283 "ranges = list(yrange_MAPE,yrange_abs,yrange_RMSE)\n",
285 "par(mfrow=c(1,3))\n",
286 "titles = paste(\"Erreur\",c(\"MAPE\",\"abs\",\"RMSE\"))\n",
287 "for (i in 1:3) #error type (MAPE,abs,RMSE)\n",
289 " for (j in 1:5) #model (nn,np,nz,pp,pz)\n",
291 " xlab = if (j>1) \"\" else \"Temps\"\n",
292 " ylab = if (j>1) \"\" else \"Erreur\"\n",
293 " main = if (j>1) \"\" else titles[i]\n",
294 " plot(p[[j]]$errors[[i]], type=\"l\", col=j, main=titles[i], xlab=xlab, ylab=ylab, ylim=ranges[[i]])\n",
301 "p = list(p_nn_epandage, p_nn_nonpollue, p_nn_chauffage)\n",
302 "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",
303 "e1 = getErrors(data, forecasts_1, 17)\n",
305 "e = list(e1,e2,e3)\n",
306 "yrange_MAPE = range(e1$MAPE, e2$MAPE, e3$MAPE)\n",
307 "yrange_abs = range(e1$abs, e2$abs, e3$abs)\n",
308 "yrange_RMSE = range(e1$RMSE, e2$RMSE, e3$RMSE)\n",
309 "ranges = list(yrange_MAPE,yrange_abs,yrange_RMSE)\n",
311 "par(mfrow=c(1,3))\n",
312 "titles = paste(\"Erreur\",c(\"MAPE\",\"abs\",\"RMSE\"))\n",
313 "for (i in 1:3) #error type (MAPE,abs,RMSE)\n",
315 " for (j in 1:3) #model (nn,np,nz,pp,pz)\n",
317 " xlab = if (j>1) \"\" else \"Temps\"\n",
318 " ylab = if (j>1) \"\" else \"Erreur\"\n",
319 " main = if (j>1) \"\" else titles[i]\n",
320 " plot(e[[j]][[i]], type=\"l\", col=j, main=titles[i], xlab=xlab, ylab=ylab, ylim=ranges[[i]])\n",
326 "par(mfrow=c(1,2))\n",
327 "#p[[i]]$forecasts[[index]]\n",
328 "#futurs des blocs du passé pour le jour 2290 ::\n",
329 "futurs = lapply(1:length(p[[1]]$forecasts[[2290]]$indices),\n",
330 " function(index) ( data[[ p[[1]]$forecasts[[2290]]$indices[index]+1 ]]$pm10 ) )\n",
331 "#vrai futur (en rouge), vrai jour (en noir)\n",
332 "r_min = min( sapply( 1:length(futurs), function(index) ( min(futurs[[index]] ) ) ) )\n",
333 "r_max = max( sapply( 1:length(futurs), function(index) ( max(futurs[[index]] ) ) ) )\n",
334 "for (i in 1:length(futurs))\n",
336 " plot(futurs[[i]], col=1, ylim=c(r_min,r_max), type=\"l\")\n",
337 " if (i<length(futurs))\n",
341 "plot(data[[2290]]$pm10, ylim=c(r_min, r_max), col=1, type=\"l\")\n",
343 "plot(data[[2291]]$pm10, ylim=c(r_min, r_max), col=2, type=\"l\")\n"
347 "cell_type": "markdown",
350 "Meilleurs results: nn et nz (np moins bon)"
355 "execution_count": null,
361 "#%%TODO: analyse sur les trois périodes indiquées par Michel ; simtype==\"exo\" par defaut\n",
362 "#16/03/2015 2288\n",
363 "p_nn_epandage = predictPM10(data, 2287, 2293, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"mix\")"
368 "execution_count": null,
374 "p_nn_epandage = predictPM10(data, 2287, 2293, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"mix\")"
379 "execution_count": null,
385 "options(repr.plot.width=9, repr.plot.height=6)\n",
386 "plot(p_nn_epandage$errors$abs, type=\"l\", col=1, main=\"Erreur absolue\", xlab=\"Temps\",\n",
387 " ylab=\"Erreur\", ylim=range(p_nn_epandage$errors$abs))"
392 "execution_count": null,
398 "#19/01/2015 2232\n",
399 "p_nn_chauffage = predictPM10(data, 2231, 2237, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE)"
404 "execution_count": null,
410 "#23/02/2015 2267\n",
411 "p_nn_nonpollue = predictPM10(data, 2266, 2272, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"mix\")"
416 "execution_count": null,
422 "plot(p_nn_nonpollue$errors$abs, type=\"l\", col=2, main=\"Erreur absolue\", xlab=\"Temps\",\n",
423 " ylab=\"Erreur\", ylim=range(p_nn_nonpollue$errors$abs))"
428 "execution_count": null,
435 "data = getData(\"local\", \"7h\")\n",
436 "p_nn_epandage = predictPM10(data, 2287, 2293, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"endo\")\n",
437 "p_nn_nonpollue = predictPM10(data, 2266, 2272, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"endo\")\n",
438 "p_nn_chauffage = predictPM10(data, 2231, 2237, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"endo\")"
443 "execution_count": null,
449 "p = list(p_nn_epandage, p_nn_nonpollue, p_nn_chauffage)\n",
450 "#yrange_MAPE = range(p[[1]]$errors$MAPE, p[[2]]$errors$MAPE, p[[3]]$errors$MAPE)\n",
451 "#yrange_abs = range(p[[1]]$errors$abs, p[[2]]$errors$abs, p[[3]]$errors$abs)\n",
452 "#yrange_RMSE = range(p[[1]]$errors$RMSE, p[[2]]$errors$RMSE, p[[3]]$errors$RMSE)\n",
453 "#ranges = list(yrange_MAPE,yrange_abs,yrange_RMSE)\n",
454 "print(p[[1]]$forecasts[[2290]])"
459 "execution_count": null,
465 "par(mfrow=c(1,3))\n",
466 "titles = paste(\"Erreur\",c(\"MAPE\",\"abs\",\"RMSE\"))\n",
467 "for (i in 1:3) #error type (MAPE,abs,RMSE)\n",
469 " for (j in 1:5) #model (nn,np,nz,pp,pz)\n",
471 " xlab = if (j>1) \"\" else \"Temps\"\n",
472 " ylab = if (j>1) \"\" else \"Erreur\"\n",
473 " main = if (j>1) \"\" else titles[i]\n",
474 " plot(p[[j]]$errors[[i]], type=\"l\", col=j, main=titles[i], xlab=xlab, ylab=ylab, ylim=ranges[[i]])\n",
482 "cell_type": "markdown",
498 "codemirror_mode": "r",
499 "file_extension": ".r",
500 "mimetype": "text/x-r-source",
502 "pygments_lexer": "r",