fixes and improvements
[talweg.git] / reports / report_2017-01-13.rnw
diff --git a/reports/report_2017-01-13.rnw b/reports/report_2017-01-13.rnw
new file mode 100644 (file)
index 0000000..c2425af
--- /dev/null
@@ -0,0 +1,186 @@
+\documentclass[a4paper,12pt]{article}
+\usepackage[utf8]{inputenc}
+\usepackage[T1]{fontenc}
+
+\renewcommand*\familydefault{\sfdefault}
+
+\marginparwidth 0pt
+\oddsidemargin 0pt
+\evensidemargin 0pt
+\marginparsep 0pt
+\topmargin 0pt
+\textwidth 16cm
+\textheight 23cm
+\parindent 5mm
+
+\begin{document}
+
+\section{Package R "ppmfun"}
+
+Le package $-$ Predict PM10 with FUNctional methods $-$ contient le code permettant de (re)lancer
+les expériences numériques décrites dans ce document. La fonction principale \emph{predictPM10}
+se divise en trois parties, décrites successivement au cours des trois paragraphes suivants.\\
+
+<<setup, out.width='7cm', out.height='7cm'>>=
+#Chargement de la librairie (après compilation, "R CMD INSTALL ppmfun/")
+library(ppmfun)
+
+#Exemple d'appel principal (détaillé ci-après)
+p_mix = predictPM10("GMT", "7h", 400, 30, "Neighbors", NULL, 0, "direct",
+  simtype="mix")
+
+#Allure des courbes prédites
+yrange = range(p_mix$forecasts[401:430], na.rm=TRUE)
+plot(0,xaxt='n',yaxt='n',bty='n',pch='',ylab='Prédictions PM10',
+  xlab='Temps',ylim=yrange,main="Courbes PM10 prédites")
+for (i in 401:430)
+{
+  if (!any(is.na(p_mix$forecasts[[i]])))
+  {
+    par(new=TRUE)
+    plot(p_mix$forecasts[[i]], type="l", ylim=yrange, xlab="", ylab="")
+  }
+}
+@
+
+L'appel à \emph{predictPM10()} ci-dessus se traduit par :
+\begin{enumerate}
+  \item charger les données découpées selon le temps universel, en segments de $24h$ de $7h15$ à
+    $7h$ le lendemain ;
+  \item commencer la prédiction au jour $400$, terminer au jour $400+30-1 = 429$ ;
+  \item utiliser la méthode "Neighbors" qui place plus de poids sur les voisins de la courbe de
+    PM10 du jour courant, en tenant compte de tout l'historique ;
+  \item raccorder continûment la prévision centrée aux mesures sur le dernier bloc de $24h$.
+\end{enumerate}
+
+\subsection{Acquisition des données}
+
+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 organisées sous forme d'une liste
+\emph{data}, la $i^{eme}$ cellule correspondant aux données disponibles au $i^{eme}$ jour à
+l'heure $H$ de prédiction choisie (0h15, 7h15 ou 13h15) : c'est-à-dire les valeurs des PM10 de
+$H-24h$ à $H-15m$, ainsi que les variables météo mesurées du dernier jour complet avant l'heure
+$H$, et les variables météo prédites pour la période de $0h15$ à $0h$ du jour courant.\\
+
+Exemple :\\
+<<data>>=
+#Le premier argument indique la zone horaire souhaitée ; "GMT" ou "local"
+#pour l'heure française, ou tout autre fuseau horaire.
+data = getData("GMT", "7h")
+@
+
+\subsection{Prédiction}
+
+Deux types de prévisions du prochain bloc de $24h$ sont à distinguer :
+\begin{itemize}
+  \item prévision de la forme (centrée) ;
+  \item prévision du niveau.
+\end{itemize}
+
+\noindent Si l'on choisit de raccorder la prévision de la forme au dernier PM10 mesuré, alors le niveau n'a
+pas à être prédit (d'où l'argument \texttt{NULL} dans l'appel principal). Dans le cas contraire il faut
+préciser une méthode ; seule la persistance est actuellement implémentée. la méthode de prévision
+de forme "Neighbors" est détaillée ci-après (voir aussi le fichier S\_Neighbors.R).\\
+
+\begin{enumerate}
+  \item \textbf{Préparation des données} : calcul des niveaux sur 24h, fenêtrage si demandé
+    (paramètre "memory").
+  \item \textbf{Optimisation des paramètres d'échelle} : via la fonction \emph{optim()}
+    minimisant la somme des 45 dernières erreurs jounalières L2.
+  \item \textbf{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.
+\end{enumerate}
+
+Exemple :\\
+<<forecasts>>=
+forecasts = as.list(rep(NA, length(data)))
+for (i in 400:429)
+{
+  #forecast with data up to index i
+  forecasts[[i+1]] = getForecasts(data[1:i], "Neighbors", NULL, 0,
+    "direct", simtype="mix")
+}
+@
+
+\subsection{Calcul des erreurs}
+
+Pour chacun des instants à prévoir jusqu'à minuit du jour courant, on calcule l'erreur moyenne
+sur tous les instants similaires du passé (sur la plage prédite, dans l'exemple 401 à 430). Deux
+types d'erreurs sont considérées :
+\begin{itemize}
+  \item l'erreur "L1" égale à la valeur absolue entre la mesure et la prédiction ;
+  \item l'erreur "MAPE" égale à l'erreur L1 normalisée par la mesure.
+\end{itemize}
+
+Code :\\
+<<errors, out.width='7cm', out.height='7cm', fig.show='hold'>>=
+e = getErrors(data, forecasts)
+#Erreurs absolues, point par point, moyennées sur les 30 jours
+plot(e$L1, type="l", xlab="Temps", ylab="Erreur absolue")
+#Erreurs relatives, point par point, moyennées sur les 30 jours
+plot(e$MAPE, type="l", xlab="Temps", ylab="Erreur relative")
+@
+
+\subsection{Autres expériences numériques}
+
+<<others1, out.width='7cm', out.height='7cm', fig.show='hold'>>=
+p_endo = predictPM10("GMT", "7h", 400, 30, "Neighbors", NULL, 0,
+  "direct", simtype="endo")
+p_exo = predictPM10("GMT", "7h", 400, 30, "Neighbors", NULL, 0,
+  "direct", simtype="exo")
+yrange_L1 = range(p_mix$errors$L1, p_endo$errors$L1, p_exo$errors$L1)
+plot(p_mix$errors$L1, type="l", main="Erreur L1", xlab="Temps",
+  ylab="Erreur absolue", ylim=yrange_L1) ; par(new=TRUE)
+plot(p_endo$errors$L1, type="l", col=2, xlab="", ylab="",
+  ylim=yrange_L1) ; par(new=TRUE)
+plot(p_exo$errors$L1, type="l", col=3, xlab="", ylab="",
+  ylim=yrange_L1)
+yrange_MAPE =
+  range(p_mix$errors$MAPE, p_endo$errors$MAPE, p_exo$errors$MAPE)
+plot(p_mix$errors$MAPE, type="l", main="Erreur MAPE", xlab="Temps",
+  ylab="Erreur relative", ylim=yrange_MAPE) ; par(new=TRUE)
+plot(p_endo$errors$MAPE, type="l", col=2, xlab="", ylab="",
+  ylim=yrange_MAPE) ; par(new=TRUE)
+plot(p_exo$errors$MAPE, type="l", col=3, xlab="", ylab="",
+  ylim=yrange_MAPE)
+
+#Ne tenir compte que des similarités sur les variables exogènes semble
+#conduire à l'erreur la plus faible.
+@
+
+<<others2, out.width='7cm', out.height='7cm', fig.show='hold'>>=
+p_exo_h = predictPM10("GMT", "7h", 400, 30, "Neighbors", NULL, 0,
+  "direct", simtype="exo", h_window=0.25)
+plot(p_exo_h$errors$L1, type="l", main="Erreur L1", xlab="Temps",
+  ylab="Erreur absolue")
+plot(p_exo_h$errors$MAPE, type="l", main="Erreur MAPE", xlab="Temps",
+  ylab="Erreur relative")
+
+#Diminuer la fenêtre n'améliore pas les performances moyennes
+#(car les données individuelles sont très variables).
+@
+
+<<others3, out.width='7cm', out.height='7cm', fig.show='hold'>>=
+p_exo_s = predictPM10("GMT", "7h", 400, 30, "Neighbors", "Persistence",
+  0, "separate", simtype="exo")
+plot(p_exo_s$errors$L1, type="l", main="Erreur L1", xlab="Temps",
+  ylab="Erreur absolue")
+plot(p_exo_s$errors$MAPE, type="l", main="Erreur MAPE", xlab="Temps",
+  ylab="Erreur relative")
+
+#Prédire séparément forme et niveau mène à une erreur plus grande ;
+#d'autres méthodes de prévision du niveau doivent tout de même être testées.
+@
+
+\subsection{Suite du travail}
+
+Le type de jour n'est pas pris en compte dans la recherche de voisins ; cela diminuerait
+nettement le nombre de similarités retenues, mais pourrait significativement améliorer les
+prévisions.\\
+
+\noindent Il serait intéressant également de disposer de plusieurs méthodes de prédiction, pour
+par exemple les agréger à l'aide de méthodes similaires à celles du précédent contrat.
+
+\end{document}