reports/report_02-02-2017.ipynb

   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",
  60    "execution_count": 1,
  61    "metadata": {
  62     "collapsed": false
  63    },
  64    "outputs": [
  65     {
  66      "ename": "ERROR",
  67      "evalue": "Error in eval(expr, envir, enclos): could not find function \"getData\"\n",
  68      "output_type": "error",
  69      "traceback": [
  70       "Error in eval(expr, envir, enclos): could not find function \"getData\"\nTraceback:\n"
  71      ]
  72     }
  73    ],
  74    "source": [
  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"
  80    ]
  81   },
  82   {
  83    "cell_type": "markdown",
  84    "metadata": {},
  85    "source": [
  86     "### Prédiction\n",
  87     "\n",
  88     "Deux types de prévisions du prochain bloc de $24h$ sont à distinguer :\n",
  89     "\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",
  92     "\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",
  94     "\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",
  98     "\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",
 100     "\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",
 104     "\n",
 105     "### Calcul des erreurs\n",
 106     "\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",
 108     "\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",
 111     "\n",
 112     "### Expériences numériques"
 113    ]
 114   },
 115   {
 116    "cell_type": "code",
 117    "execution_count": null,
 118    "metadata": {
 119     "collapsed": false
 120    },
 121    "outputs": [],
 122    "source": [
 123     "options(repr.plot.width=9, repr.plot.height=3)"
 124    ]
 125   },
 126   {
 127    "cell_type": "code",
 128    "execution_count": null,
 129    "metadata": {
 130     "collapsed": false
 131    },
 132    "outputs": [],
 133    "source": [
 134     "p_endo = predictPM10(data, 2200, 2230, 0,0, \"Neighbors\", \"Neighbors\", simtype=\"endo\")"
 135    ]
 136   },
 137   {
 138    "cell_type": "code",
 139    "execution_count": null,
 140    "metadata": {
 141     "collapsed": false
 142    },
 143    "outputs": [],
 144    "source": [
 145     "p_exo = predictPM10(data, 2200, 2230, 0,0, \"Neighbors\", \"Neighbors\", simtype=\"exo\")"
 146    ]
 147   },
 148   {
 149    "cell_type": "code",
 150    "execution_count": null,
 151    "metadata": {
 152     "collapsed": false
 153    },
 154    "outputs": [],
 155    "source": [
 156     "p_mix = predictPM10(data, 2200, 2230, 0,0, \"Neighbors\", \"Neighbors\", simtype=\"mix\")"
 157    ]
 158   },
 159   {
 160    "cell_type": "code",
 161    "execution_count": null,
 162    "metadata": {
 163     "collapsed": false
 164    },
 165    "outputs": [],
 166    "source": [
 167     "p = list(p_endo, p_exo, p_mix)"
 168    ]
 169   },
 170   {
 171    "cell_type": "code",
 172    "execution_count": null,
 173    "metadata": {
 174     "collapsed": false
 175    },
 176    "outputs": [],
 177    "source": [
 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",
 182     "\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",
 186     "{\n",
 187     "  for (j in 1:3) #model (mix,endo,exo)\n",
 188     "  {\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",
 193     "    if (j<3)\n",
 194     "        par(new=TRUE)\n",
 195     "  }\n",
 196     "}"
 197    ]
 198   },
 199   {
 200    "cell_type": "markdown",
 201    "metadata": {},
 202    "source": [
 203     "Ne tenir compte que des similarités sur les variables exogènes semble conduire à l'erreur la plus faible."
 204    ]
 205   },
 206   {
 207    "cell_type": "code",
 208    "execution_count": null,
 209    "metadata": {
 210     "collapsed": false
 211    },
 212    "outputs": [],
 213    "source": [
 214     "p_nn = predictPM10(data, 2200, 2230, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE)"
 215    ]
 216   },
 217   {
 218    "cell_type": "code",
 219    "execution_count": null,
 220    "metadata": {
 221     "collapsed": false
 222    },
 223    "outputs": [],
 224    "source": [
 225     "p_np = predictPM10(data, 2200, 2230, 0, 0, \"Neighbors\", \"Persistence\", sameSeaon=TRUE)"
 226    ]
 227   },
 228   {
 229    "cell_type": "code",
 230    "execution_count": null,
 231    "metadata": {
 232     "collapsed": false
 233    },
 234    "outputs": [],
 235    "source": [
 236     "p_nz = predictPM10(data, 2200, 2230, 0, 0, \"Neighbors\", \"Zero\", sameSeaon=TRUE)"
 237    ]
 238   },
 239   {
 240    "cell_type": "code",
 241    "execution_count": null,
 242    "metadata": {
 243     "collapsed": false
 244    },
 245    "outputs": [],
 246    "source": [
 247     "p_pp = predictPM10(data, 2200, 2230, 0, 0, \"Persistence\", \"Persistence\")"
 248    ]
 249   },
 250   {
 251    "cell_type": "code",
 252    "execution_count": null,
 253    "metadata": {
 254     "collapsed": false
 255    },
 256    "outputs": [],
 257    "source": [
 258     "p_pz = predictPM10(data, 2200, 2230, 0, 0, \"Persistence\", \"Zero\")"
 259    ]
 260   },
 261   {
 262    "cell_type": "code",
 263    "execution_count": null,
 264    "metadata": {
 265     "collapsed": false
 266    },
 267    "outputs": [],
 268    "source": [
 269     "p = list(p_nn, p_np, p_nz, p_pp, p_pz)"
 270    ]
 271   },
 272   {
 273    "cell_type": "code",
 274    "execution_count": null,
 275    "metadata": {
 276     "collapsed": false
 277    },
 278    "outputs": [],
 279    "source": [
 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",
 284     "\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",
 288     "{\n",
 289     "  for (j in 1:5) #model (nn,np,nz,pp,pz)\n",
 290     "  {\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",
 295     "    if (j<5)\n",
 296     "      par(new=TRUE)\n",
 297     "  }\n",
 298     "}\n",
 299     "        \n",
 300     "\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",
 304     " \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",
 310     "\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",
 314     "{\n",
 315     "  for (j in 1:3) #model (nn,np,nz,pp,pz)\n",
 316     "  {\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",
 321     "    if (j<3)\n",
 322     "      par(new=TRUE)\n",
 323     "  }\n",
 324     "}\n",
 325     "\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",
 335     "{\n",
 336     "    plot(futurs[[i]], col=1, ylim=c(r_min,r_max), type=\"l\")\n",
 337     "    if (i<length(futurs))\n",
 338     "        par(new=TRUE)\n",
 339     "}\n",
 340     "\n",
 341     "plot(data[[2290]]$pm10, ylim=c(r_min, r_max), col=1, type=\"l\")\n",
 342     "    par(new=TRUE)\n",
 343     "plot(data[[2291]]$pm10, ylim=c(r_min, r_max), col=2, type=\"l\")\n"
 344    ]
 345   },
 346   {
 347    "cell_type": "markdown",
 348    "metadata": {},
 349    "source": [
 350     "Meilleurs results: nn et nz (np moins bon)"
 351    ]
 352   },
 353   {
 354    "cell_type": "code",
 355    "execution_count": null,
 356    "metadata": {
 357     "collapsed": false
 358    },
 359    "outputs": [],
 360    "source": [
 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\")"
 364    ]
 365   },
 366   {
 367    "cell_type": "code",
 368    "execution_count": null,
 369    "metadata": {
 370     "collapsed": false
 371    },
 372    "outputs": [],
 373    "source": [
 374     "p_nn_epandage = predictPM10(data, 2287, 2293, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"mix\")"
 375    ]
 376   },
 377   {
 378    "cell_type": "code",
 379    "execution_count": null,
 380    "metadata": {
 381     "collapsed": false
 382    },
 383    "outputs": [],
 384    "source": [
 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))"
 388    ]
 389   },
 390   {
 391    "cell_type": "code",
 392    "execution_count": null,
 393    "metadata": {
 394     "collapsed": false
 395    },
 396    "outputs": [],
 397    "source": [
 398     "#19/01/2015 2232\n",
 399     "p_nn_chauffage = predictPM10(data, 2231, 2237, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE)"
 400    ]
 401   },
 402   {
 403    "cell_type": "code",
 404    "execution_count": null,
 405    "metadata": {
 406     "collapsed": false
 407    },
 408    "outputs": [],
 409    "source": [
 410     "#23/02/2015 2267\n",
 411     "p_nn_nonpollue = predictPM10(data, 2266, 2272, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"mix\")"
 412    ]
 413   },
 414   {
 415    "cell_type": "code",
 416    "execution_count": null,
 417    "metadata": {
 418     "collapsed": false
 419    },
 420    "outputs": [],
 421    "source": [
 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))"
 424    ]
 425   },
 426   {
 427    "cell_type": "code",
 428    "execution_count": null,
 429    "metadata": {
 430     "collapsed": false
 431    },
 432    "outputs": [],
 433    "source": [
 434     "library(ppmfun)\n",
 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\")"
 439    ]
 440   },
 441   {
 442    "cell_type": "code",
 443    "execution_count": null,
 444    "metadata": {
 445     "collapsed": false
 446    },
 447    "outputs": [],
 448    "source": [
 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]])"
 455    ]
 456   },
 457   {
 458    "cell_type": "code",
 459    "execution_count": null,
 460    "metadata": {
 461     "collapsed": false
 462    },
 463    "outputs": [],
 464    "source": [
 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",
 468     "{\n",
 469     "  for (j in 1:5) #model (nn,np,nz,pp,pz)\n",
 470     "  {\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",
 475     "    if (j<5)\n",
 476     "      par(new=TRUE)\n",
 477     "  }\n",
 478     "}"
 479    ]
 480   },
 481   {
 482    "cell_type": "markdown",
 483    "metadata": {},
 484    "source": [
 485     "## Bilan\n",
 486     "\n",
 487     "TODO"
 488    ]
 489   }
 490  ],
 491  "metadata": {
 492   "kernelspec": {
 493    "display_name": "R",
 494    "language": "R",
 495    "name": "ir"
 496   },
 497   "language_info": {
 498    "codemirror_mode": "r",
 499    "file_extension": ".r",
 500    "mimetype": "text/x-r-source",
 501    "name": "R",
 502    "pygments_lexer": "r",
 503    "version": "3.3.2"
 504   }
 505  },
 506  "nbformat": 4,
 507  "nbformat_minor": 2
 508 }