+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Package R \"vortex\"\n",
+ "\n",
+ "using Vectorial exOgenous variables to foRecast Time-sErieX.\n",
+ "\n",
+ "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",
+ "\n",
+ "Fonctions principales :\n",
+ "\n",
+ " * <code>getData</code> : charge un jeu de données en mémoire\n",
+ " * <code>getForecast</code> : prédit les lendemains aux indices demandés\n",
+ "\n",
+ "Diverses méthodes permettent ensuite d'analyser les performances : <code>getError</code>, <code>plotXYZ</code> : voir la section \"see also\" dans <code>?plotError</code>."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "#Chargement de la librairie (après compilation, \"R CMD INSTALL ppmfun/\")\n",
+ "library(vortex)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Note : sur la base de nos dernières expériences, on considère que \n",
+ "\n",
+ " * on ne touche pas à la fenêtre obtenue par la fonction <code>optimize</code> ;\n",
+ " * 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",
+ "\n",
+ "### Acquisition des données\n",
+ "\n",
+ "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",
+ "\n",
+ "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",
+ "\n",
+ " * <code>getTime</code> : suite des date+heure\n",
+ " * <code>getCenteredSerie</code> : série centrée\n",
+ " * <code>getLevel</code> : niveau\n",
+ " * <code>getSerie</code> : série *non* centrée\n",
+ " * <code>getExoHat</code> : variables exogènes prévues\n",
+ " * <code>getExoDm1</code> : variables exogènes mesurées la veille\n",
+ "\n",
+ "Exemple :"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "ename": "ERROR",
+ "evalue": "Error in eval(expr, envir, enclos): could not find function \"getData\"\n",
+ "output_type": "error",
+ "traceback": [
+ "Error in eval(expr, envir, enclos): could not find function \"getData\"\nTraceback:\n"
+ ]
+ }
+ ],
+ "source": [
+ "# Voir ?getData pour les arguments\n",
+ "data = getData(ts_data=\"data/pm10_mesures_H_loc.csv\", exo_data=\"data/meteo_extra_noNAs.csv\",\n",
+ " input_tz = \"Europe/Paris\", working_tz=\"Europe/Paris\", predict_at=\"07\")\n",
+ "data$getLevel(10) #niveau du jour 10\n",
+ "data$getExoHat(17) #météo prévue pour le jour 18"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Prédiction\n",
+ "\n",
+ "Deux types de prévisions du prochain bloc de $24h$ sont à distinguer :\n",
+ "\n",
+ " * prévision de la forme (centrée) ;\n",
+ " * prévision du saut d'une fin de série au début de la suivante.\n",
+ "\n",
+ "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",
+ "\n",
+ " 1. **Préparation des données** : fenêtrage si demandé (paramètre \"memory\"), recherche des paires de jours consécutifs sans valeurs manquantes.\n",
+ " 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",
+ " 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",
+ "\n",
+ "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",
+ "\n",
+ " * <code>getSerie</code> : série prévue (sans les information de temps)\n",
+ " * <code>getParams</code> : liste des paramètres (poids, fenêtre, ...)\n",
+ " * <code>getIndexInData</code> : indice du jour où s'effectue la prévision relativement au jeu de données\n",
+ "\n",
+ "### Calcul des erreurs\n",
+ "\n",
+ "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",
+ "\n",
+ " * l'erreur \"abs\" égale à la valeur absolue moyenne entre la mesure et la prédiction ;\n",
+ " * l'erreur \"MAPE\" égale à l'erreur absolue normalisée par la mesure.\n",
+ "\n",
+ "### Expériences numériques"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "options(repr.plot.width=9, repr.plot.height=3)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "p_endo = predictPM10(data, 2200, 2230, 0,0, \"Neighbors\", \"Neighbors\", simtype=\"endo\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "p_exo = predictPM10(data, 2200, 2230, 0,0, \"Neighbors\", \"Neighbors\", simtype=\"exo\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "p_mix = predictPM10(data, 2200, 2230, 0,0, \"Neighbors\", \"Neighbors\", simtype=\"mix\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "p = list(p_endo, p_exo, p_mix)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "yrange_MAPE = range(p_mix$errors$MAPE, p_endo$errors$MAPE, p_exo$errors$MAPE)\n",
+ "yrange_abs = range(p_mix$errors$abs, p_endo$errors$abs, p_exo$errors$abs)\n",
+ "yrange_RMSE = range(p_mix$errors$RMSE, p_endo$errors$RMSE, p_exo$errors$RMSE)\n",
+ "ranges = list(yrange_MAPE,yrange_abs,yrange_RMSE)\n",
+ "\n",
+ "par(mfrow=c(1,3))\n",
+ "titles = paste(\"Erreur\",c(\"MAPE\",\"abs\",\"RMSE\"))\n",
+ "for (i in 1:3) #error type (MAPE,abs,RMSE)\n",
+ "{\n",
+ " for (j in 1:3) #model (mix,endo,exo)\n",
+ " {\n",
+ " xlab = if (j>1) \"\" else \"Temps\"\n",
+ " ylab = if (j>1) \"\" else \"Erreur\"\n",
+ " main = if (j>1) \"\" else titles[i]\n",
+ " plot(p[[j]]$errors[[i]], type=\"l\", col=j, main=main, xlab=xlab, ylab=ylab, ylim=ranges[[i]])\n",
+ " if (j<3)\n",
+ " par(new=TRUE)\n",
+ " }\n",
+ "}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Ne tenir compte que des similarités sur les variables exogènes semble conduire à l'erreur la plus faible."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "p_nn = predictPM10(data, 2200, 2230, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "p_np = predictPM10(data, 2200, 2230, 0, 0, \"Neighbors\", \"Persistence\", sameSeaon=TRUE)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "p_nz = predictPM10(data, 2200, 2230, 0, 0, \"Neighbors\", \"Zero\", sameSeaon=TRUE)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "p_pp = predictPM10(data, 2200, 2230, 0, 0, \"Persistence\", \"Persistence\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "p_pz = predictPM10(data, 2200, 2230, 0, 0, \"Persistence\", \"Zero\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "p = list(p_nn, p_np, p_nz, p_pp, p_pz)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "yrange_MAPE = range(p_nn$errors$MAPE, p_nz$errors$MAPE, p_np$errors$MAPE, p_pp$errors$MAPE, p_pz$errors$MAPE)\n",
+ "yrange_abs = range(p_nn$errors$abs, p_nz$errors$abs, p_np$errors$abs, p_pp$errors$abs, p_pz$errors$abs)\n",
+ "yrange_RMSE = range(p_nn$errors$RMSE, p_nz$errors$RMSE, p_np$errors$RMSE, p_pp$errors$RMSE, p_pz$errors$RMSE)\n",
+ "ranges = list(yrange_MAPE,yrange_abs,yrange_RMSE)\n",
+ "\n",
+ "par(mfrow=c(1,3))\n",
+ "titles = paste(\"Erreur\",c(\"MAPE\",\"abs\",\"RMSE\"))\n",
+ "for (i in 1:3) #error type (MAPE,abs,RMSE)\n",
+ "{\n",
+ " for (j in 1:5) #model (nn,np,nz,pp,pz)\n",
+ " {\n",
+ " xlab = if (j>1) \"\" else \"Temps\"\n",
+ " ylab = if (j>1) \"\" else \"Erreur\"\n",
+ " main = if (j>1) \"\" else titles[i]\n",
+ " plot(p[[j]]$errors[[i]], type=\"l\", col=j, main=titles[i], xlab=xlab, ylab=ylab, ylim=ranges[[i]])\n",
+ " if (j<5)\n",
+ " par(new=TRUE)\n",
+ " }\n",
+ "}\n",
+ " \n",
+ "\n",
+ "p = list(p_nn_epandage, p_nn_nonpollue, p_nn_chauffage)\n",
+ "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",
+ "e1 = getErrors(data, forecasts_1, 17)\n",
+ " \n",
+ "e = list(e1,e2,e3)\n",
+ "yrange_MAPE = range(e1$MAPE, e2$MAPE, e3$MAPE)\n",
+ "yrange_abs = range(e1$abs, e2$abs, e3$abs)\n",
+ "yrange_RMSE = range(e1$RMSE, e2$RMSE, e3$RMSE)\n",
+ "ranges = list(yrange_MAPE,yrange_abs,yrange_RMSE)\n",
+ "\n",
+ "par(mfrow=c(1,3))\n",
+ "titles = paste(\"Erreur\",c(\"MAPE\",\"abs\",\"RMSE\"))\n",
+ "for (i in 1:3) #error type (MAPE,abs,RMSE)\n",
+ "{\n",
+ " for (j in 1:3) #model (nn,np,nz,pp,pz)\n",
+ " {\n",
+ " xlab = if (j>1) \"\" else \"Temps\"\n",
+ " ylab = if (j>1) \"\" else \"Erreur\"\n",
+ " main = if (j>1) \"\" else titles[i]\n",
+ " plot(e[[j]][[i]], type=\"l\", col=j, main=titles[i], xlab=xlab, ylab=ylab, ylim=ranges[[i]])\n",
+ " if (j<3)\n",
+ " par(new=TRUE)\n",
+ " }\n",
+ "}\n",
+ "\n",
+ "par(mfrow=c(1,2))\n",
+ "#p[[i]]$forecasts[[index]]\n",
+ "#futurs des blocs du passé pour le jour 2290 ::\n",
+ "futurs = lapply(1:length(p[[1]]$forecasts[[2290]]$indices),\n",
+ " function(index) ( data[[ p[[1]]$forecasts[[2290]]$indices[index]+1 ]]$pm10 ) )\n",
+ "#vrai futur (en rouge), vrai jour (en noir)\n",
+ "r_min = min( sapply( 1:length(futurs), function(index) ( min(futurs[[index]] ) ) ) )\n",
+ "r_max = max( sapply( 1:length(futurs), function(index) ( max(futurs[[index]] ) ) ) )\n",
+ "for (i in 1:length(futurs))\n",
+ "{\n",
+ " plot(futurs[[i]], col=1, ylim=c(r_min,r_max), type=\"l\")\n",
+ " if (i<length(futurs))\n",
+ " par(new=TRUE)\n",
+ "}\n",
+ "\n",
+ "plot(data[[2290]]$pm10, ylim=c(r_min, r_max), col=1, type=\"l\")\n",
+ " par(new=TRUE)\n",
+ "plot(data[[2291]]$pm10, ylim=c(r_min, r_max), col=2, type=\"l\")\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Meilleurs results: nn et nz (np moins bon)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "#%%TODO: analyse sur les trois périodes indiquées par Michel ; simtype==\"exo\" par defaut\n",
+ "#16/03/2015 2288\n",
+ "p_nn_epandage = predictPM10(data, 2287, 2293, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"mix\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "p_nn_epandage = predictPM10(data, 2287, 2293, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"mix\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "options(repr.plot.width=9, repr.plot.height=6)\n",
+ "plot(p_nn_epandage$errors$abs, type=\"l\", col=1, main=\"Erreur absolue\", xlab=\"Temps\",\n",
+ " ylab=\"Erreur\", ylim=range(p_nn_epandage$errors$abs))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "#19/01/2015 2232\n",
+ "p_nn_chauffage = predictPM10(data, 2231, 2237, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "#23/02/2015 2267\n",
+ "p_nn_nonpollue = predictPM10(data, 2266, 2272, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"mix\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "plot(p_nn_nonpollue$errors$abs, type=\"l\", col=2, main=\"Erreur absolue\", xlab=\"Temps\",\n",
+ " ylab=\"Erreur\", ylim=range(p_nn_nonpollue$errors$abs))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "library(ppmfun)\n",
+ "data = getData(\"local\", \"7h\")\n",
+ "p_nn_epandage = predictPM10(data, 2287, 2293, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"endo\")\n",
+ "p_nn_nonpollue = predictPM10(data, 2266, 2272, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"endo\")\n",
+ "p_nn_chauffage = predictPM10(data, 2231, 2237, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"endo\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "p = list(p_nn_epandage, p_nn_nonpollue, p_nn_chauffage)\n",
+ "#yrange_MAPE = range(p[[1]]$errors$MAPE, p[[2]]$errors$MAPE, p[[3]]$errors$MAPE)\n",
+ "#yrange_abs = range(p[[1]]$errors$abs, p[[2]]$errors$abs, p[[3]]$errors$abs)\n",
+ "#yrange_RMSE = range(p[[1]]$errors$RMSE, p[[2]]$errors$RMSE, p[[3]]$errors$RMSE)\n",
+ "#ranges = list(yrange_MAPE,yrange_abs,yrange_RMSE)\n",
+ "print(p[[1]]$forecasts[[2290]])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "par(mfrow=c(1,3))\n",
+ "titles = paste(\"Erreur\",c(\"MAPE\",\"abs\",\"RMSE\"))\n",
+ "for (i in 1:3) #error type (MAPE,abs,RMSE)\n",
+ "{\n",
+ " for (j in 1:5) #model (nn,np,nz,pp,pz)\n",
+ " {\n",
+ " xlab = if (j>1) \"\" else \"Temps\"\n",
+ " ylab = if (j>1) \"\" else \"Erreur\"\n",
+ " main = if (j>1) \"\" else titles[i]\n",
+ " plot(p[[j]]$errors[[i]], type=\"l\", col=j, main=titles[i], xlab=xlab, ylab=ylab, ylim=ranges[[i]])\n",
+ " if (j<5)\n",
+ " par(new=TRUE)\n",
+ " }\n",
+ "}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Bilan\n",
+ "\n",
+ "TODO"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "R",
+ "language": "R",
+ "name": "ir"
+ },
+ "language_info": {
+ "codemirror_mode": "r",
+ "file_extension": ".r",
+ "mimetype": "text/x-r-source",
+ "name": "R",
+ "pygments_lexer": "r",
+ "version": "3.3.2"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}