add part 3 of final report
[talweg.git] / reports / OLD / report_2017-02-02.ipynb
CommitLineData
3d69ff21
BA
1{
2 "cells": [
3 {
4 "cell_type": "markdown",
5 "metadata": {},
6 "source": [
7 "## Package R \"vortex\"\n",
8 "\n",
9 "using Vectorial exOgenous variables to foRecast Time-sErieX.\n",
10 "\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",
12 "\n",
13 "Fonctions principales :\n",
14 "\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",
17 "\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>."
19 ]
20 },
21 {
22 "cell_type": "code",
23 "execution_count": null,
24 "metadata": {
25 "collapsed": false
26 },
27 "outputs": [],
28 "source": [
29 "#Chargement de la librairie (après compilation, \"R CMD INSTALL ppmfun/\")\n",
30 "library(vortex)"
31 ]
32 },
33 {
34 "cell_type": "markdown",
35 "metadata": {},
36 "source": [
37 "Note : sur la base de nos dernières expériences, on considère que \n",
38 "\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",
41 "\n",
42 "### Acquisition des données\n",
43 "\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",
45 "\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",
47 "\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",
54 "\n",
55 "Exemple :"
56 ]
57 },
58 {
59 "cell_type": "code",
09cf9c19 60 "execution_count": null,
3d69ff21
BA
61 "metadata": {
62 "collapsed": false
63 },
09cf9c19 64 "outputs": [],
3d69ff21
BA
65 "source": [
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"
71 ]
72 },
73 {
74 "cell_type": "markdown",
75 "metadata": {},
76 "source": [
77 "### Prédiction\n",
78 "\n",
79 "Deux types de prévisions du prochain bloc de $24h$ sont à distinguer :\n",
80 "\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",
83 "\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",
85 "\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",
89 "\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",
91 "\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",
95 "\n",
96 "### Calcul des erreurs\n",
97 "\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",
99 "\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",
102 "\n",
103 "### Expériences numériques"
104 ]
105 },
106 {
107 "cell_type": "code",
108 "execution_count": null,
109 "metadata": {
110 "collapsed": false
111 },
112 "outputs": [],
113 "source": [
114 "options(repr.plot.width=9, repr.plot.height=3)"
115 ]
116 },
117 {
118 "cell_type": "code",
119 "execution_count": null,
120 "metadata": {
121 "collapsed": false
122 },
123 "outputs": [],
124 "source": [
125 "p_endo = predictPM10(data, 2200, 2230, 0,0, \"Neighbors\", \"Neighbors\", simtype=\"endo\")"
126 ]
127 },
128 {
129 "cell_type": "code",
130 "execution_count": null,
131 "metadata": {
132 "collapsed": false
133 },
134 "outputs": [],
135 "source": [
136 "p_exo = predictPM10(data, 2200, 2230, 0,0, \"Neighbors\", \"Neighbors\", simtype=\"exo\")"
137 ]
138 },
139 {
140 "cell_type": "code",
141 "execution_count": null,
142 "metadata": {
143 "collapsed": false
144 },
145 "outputs": [],
146 "source": [
147 "p_mix = predictPM10(data, 2200, 2230, 0,0, \"Neighbors\", \"Neighbors\", simtype=\"mix\")"
148 ]
149 },
150 {
151 "cell_type": "code",
152 "execution_count": null,
153 "metadata": {
154 "collapsed": false
155 },
156 "outputs": [],
157 "source": [
158 "p = list(p_endo, p_exo, p_mix)"
159 ]
160 },
161 {
162 "cell_type": "code",
163 "execution_count": null,
164 "metadata": {
165 "collapsed": false
166 },
167 "outputs": [],
168 "source": [
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",
173 "\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",
177 "{\n",
178 " for (j in 1:3) #model (mix,endo,exo)\n",
179 " {\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",
184 " if (j<3)\n",
185 " par(new=TRUE)\n",
186 " }\n",
187 "}"
188 ]
189 },
190 {
191 "cell_type": "markdown",
192 "metadata": {},
193 "source": [
194 "Ne tenir compte que des similarités sur les variables exogènes semble conduire à l'erreur la plus faible."
195 ]
196 },
197 {
198 "cell_type": "code",
199 "execution_count": null,
200 "metadata": {
201 "collapsed": false
202 },
203 "outputs": [],
204 "source": [
205 "p_nn = predictPM10(data, 2200, 2230, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE)"
206 ]
207 },
208 {
209 "cell_type": "code",
210 "execution_count": null,
211 "metadata": {
212 "collapsed": false
213 },
214 "outputs": [],
215 "source": [
216 "p_np = predictPM10(data, 2200, 2230, 0, 0, \"Neighbors\", \"Persistence\", sameSeaon=TRUE)"
217 ]
218 },
219 {
220 "cell_type": "code",
221 "execution_count": null,
222 "metadata": {
223 "collapsed": false
224 },
225 "outputs": [],
226 "source": [
227 "p_nz = predictPM10(data, 2200, 2230, 0, 0, \"Neighbors\", \"Zero\", sameSeaon=TRUE)"
228 ]
229 },
230 {
231 "cell_type": "code",
232 "execution_count": null,
233 "metadata": {
234 "collapsed": false
235 },
236 "outputs": [],
237 "source": [
238 "p_pp = predictPM10(data, 2200, 2230, 0, 0, \"Persistence\", \"Persistence\")"
239 ]
240 },
241 {
242 "cell_type": "code",
243 "execution_count": null,
244 "metadata": {
245 "collapsed": false
246 },
247 "outputs": [],
248 "source": [
249 "p_pz = predictPM10(data, 2200, 2230, 0, 0, \"Persistence\", \"Zero\")"
250 ]
251 },
252 {
253 "cell_type": "code",
254 "execution_count": null,
255 "metadata": {
256 "collapsed": false
257 },
258 "outputs": [],
259 "source": [
260 "p = list(p_nn, p_np, p_nz, p_pp, p_pz)"
261 ]
262 },
263 {
264 "cell_type": "code",
265 "execution_count": null,
266 "metadata": {
267 "collapsed": false
268 },
269 "outputs": [],
270 "source": [
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",
275 "\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",
279 "{\n",
280 " for (j in 1:5) #model (nn,np,nz,pp,pz)\n",
281 " {\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",
286 " if (j<5)\n",
287 " par(new=TRUE)\n",
288 " }\n",
289 "}\n",
290 " \n",
291 "\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",
295 " \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",
301 "\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",
305 "{\n",
306 " for (j in 1:3) #model (nn,np,nz,pp,pz)\n",
307 " {\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",
312 " if (j<3)\n",
313 " par(new=TRUE)\n",
314 " }\n",
315 "}\n",
316 "\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",
326 "{\n",
327 " plot(futurs[[i]], col=1, ylim=c(r_min,r_max), type=\"l\")\n",
328 " if (i<length(futurs))\n",
329 " par(new=TRUE)\n",
330 "}\n",
331 "\n",
332 "plot(data[[2290]]$pm10, ylim=c(r_min, r_max), col=1, type=\"l\")\n",
333 " par(new=TRUE)\n",
334 "plot(data[[2291]]$pm10, ylim=c(r_min, r_max), col=2, type=\"l\")\n"
335 ]
336 },
337 {
338 "cell_type": "markdown",
339 "metadata": {},
340 "source": [
341 "Meilleurs results: nn et nz (np moins bon)"
342 ]
343 },
344 {
345 "cell_type": "code",
346 "execution_count": null,
347 "metadata": {
348 "collapsed": false
349 },
350 "outputs": [],
351 "source": [
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\")"
355 ]
356 },
357 {
358 "cell_type": "code",
359 "execution_count": null,
360 "metadata": {
361 "collapsed": false
362 },
363 "outputs": [],
364 "source": [
365 "p_nn_epandage = predictPM10(data, 2287, 2293, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"mix\")"
366 ]
367 },
368 {
369 "cell_type": "code",
370 "execution_count": null,
371 "metadata": {
372 "collapsed": false
373 },
374 "outputs": [],
375 "source": [
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))"
379 ]
380 },
381 {
382 "cell_type": "code",
383 "execution_count": null,
384 "metadata": {
385 "collapsed": false
386 },
387 "outputs": [],
388 "source": [
389 "#19/01/2015 2232\n",
390 "p_nn_chauffage = predictPM10(data, 2231, 2237, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE)"
391 ]
392 },
393 {
394 "cell_type": "code",
395 "execution_count": null,
396 "metadata": {
397 "collapsed": false
398 },
399 "outputs": [],
400 "source": [
401 "#23/02/2015 2267\n",
402 "p_nn_nonpollue = predictPM10(data, 2266, 2272, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"mix\")"
403 ]
404 },
405 {
406 "cell_type": "code",
407 "execution_count": null,
408 "metadata": {
409 "collapsed": false
410 },
411 "outputs": [],
412 "source": [
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))"
415 ]
416 },
417 {
418 "cell_type": "code",
419 "execution_count": null,
420 "metadata": {
421 "collapsed": false
422 },
423 "outputs": [],
424 "source": [
425 "library(ppmfun)\n",
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\")"
430 ]
431 },
432 {
433 "cell_type": "code",
434 "execution_count": null,
435 "metadata": {
436 "collapsed": false
437 },
438 "outputs": [],
439 "source": [
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]])"
446 ]
447 },
448 {
449 "cell_type": "code",
450 "execution_count": null,
451 "metadata": {
452 "collapsed": false
453 },
454 "outputs": [],
455 "source": [
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",
459 "{\n",
460 " for (j in 1:5) #model (nn,np,nz,pp,pz)\n",
461 " {\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",
466 " if (j<5)\n",
467 " par(new=TRUE)\n",
468 " }\n",
469 "}"
470 ]
471 },
472 {
473 "cell_type": "markdown",
474 "metadata": {},
475 "source": [
476 "## Bilan\n",
477 "\n",
478 "TODO"
479 ]
480 }
481 ],
482 "metadata": {
483 "kernelspec": {
484 "display_name": "R",
485 "language": "R",
486 "name": "ir"
487 },
488 "language_info": {
489 "codemirror_mode": "r",
490 "file_extension": ".r",
491 "mimetype": "text/x-r-source",
492 "name": "R",
493 "pygments_lexer": "r",
494 "version": "3.3.2"
495 }
496 },
497 "nbformat": 4,
498 "nbformat_minor": 2
499}