on the way to R6 class + remove truncated days (simplifications)

[talweg.git] / reports / report_2017-01-13.rnw

\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}