--- /dev/null
+*.pdf filter=fat
+*.zip filter=fat
+*.tar.xz filter=fat
+*.docx filter=fat
+*.ipynb filter=nbstripout
--- /dev/null
+[rsync]
+remote = gitfat@auder.net:~/files/rnormand
--- /dev/null
+#NOTE: knitr generated files in ignored folder under reports (?!)
+.RData
+.Rhistory
+.ipynb_checkpoints/
+*.swp
+*~
+/data/*.csv
+/man/*
+!/man/*-package.Rd
+*.tar.gz
+NAMESPACE
--- /dev/null
+Package: talweg
+Title: talweg : Time-series sAmpLes forecasted With ExoGenous variables
+Version: 0.1-0
+Description: Forecast a curve sampled within the day (seconds, minutes, hours...),
+ using past measured curves, history of some exogenous variables measurements
+ and the exogenous prediction for tomorrow. Main method is getForecast()
+Authors: Benjamin Auder <Benjamin.Auder@math.u-psud.fr> [aut,cre],
+ Jean-Michel Poggi <Jean-Michel.Poggi@parisdescartes.fr> [ctb],
+ Bruno Portier <Bruno.Portier@insa-rouen.fr>, [ctb]
+Maintainer: Benjamin Auder <Benjamin.Auder@math.u-psud.fr>
+Depends:
+ R (>= 3.0)
+Suggests:
+ roxygen2,
+ testthat,
+ rmarkdown,
+ rainbow
+LazyData: yes
+URL: http://git.auder.net/?p=talweg.git
+License: MIT + file LICENSE
+RoxygenNote: 5.0.1
+Collate:
+ 'D_Neighbors.R'
+ 'D_Persistence.R'
+ 'D_Zero.R'
+ 'Data.R'
+ 'Forecast.R'
+ 'ShapeForecaster.R'
+ 'S_Average.R'
+ 'S_Neighbors.R'
+ 'S_Persistence.R'
+ 'getData.R'
+ 'getError.R'
+ 'getForecast.R'
+ 'plot.R'
+ 'utils.R'
--- /dev/null
+Copyright (c) 2016-2017, Benjamin AUDER
+ 2016-2017, Jean-Michel Poggi
+ 2016-2017, Bruno Portier
+
+Permission is hereby granted, free of charge, to any person obtaining
+a copy of this software and associated documentation files (the
+"Software"), to deal in the Software without restriction, including
+without limitation the rights to use, copy, modify, merge, publish,
+distribute, sublicense, and/or sell copies of the Software, and to
+permit persons to whom the Software is furnished to do so, subject to
+the following conditions:
+
+The above copyright notice and this permission notice shall be
+included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
--- /dev/null
+NOTE: (Michel)
+La prévision ARPEGE pour le jour J est disponible dès 4h10 TU, soit 6h10 en été et 5h10 en hiver.
+ --> donc prévis utilisables à 7 et 13h mais pas à 0h
+Donc si 7h ou 13h on utilise les prévis du jour courant (à compléter),
+sinon (0h) les mesures du dernier jour (la veille)
+
+predict heure != moyenner 4 predict 1/4 heure
+plusieurs capteurs, moyennés pour indice ATMO
+moyenne spatiale sur les capteurs d'une ville ?
+règle pr calculer moyenne horaire : >= 75% de données
+
+Dates typiques (Michel):
+19/01 : pollution chauffage
+23/02 : pas pollué
+16/03 : pollution épandage
+
+#TODO:
+#indexer axe des temps avec vraies dates
+#==> tricher en regardant le lendemain, puis chercher paires consécutives similaires dans le
+#passé, ,,,, si aucune semblable : on va se planter !
+#--> retourner un indice de confiance ?!
+#fichier texte qui affiche le top 20 des + proches (genre sam 20 janvier 2010........)
+#fonction indices <---> dates
+systématiquement horizon == 12h comme plumelabs ?
+Courbes considérées comme objets ts (package R) ? (Sans doute inutile)
+Sur exo on doit prédire moyenne du jour au lieu de courbe --------> analyser (?!)
+
+--> javascript visualisation.................. + tard
+--> soft predict: ne pas virer les séries à NA
--- /dev/null
+#' Obtain delta forecast by the Neighbors method
+#'
+#' @inheritParams getForecast
+#' @inheritParams getZeroDeltaForecast
+getNeighborsDeltaForecast = function(data, today, memory, horizon, shape_params, ...)
+{
+ if (any(is.na(shape_params$weights) | is.na(shape_params$indices)))
+ return (NA)
+
+ gaps = sapply(shape_params$indices, function(i) {
+ data$getSerie(i+1)[1] - tail(data$getSerie(i), 1)
+ })
+
+ scal_product = shape_params$weights * gaps
+ norm_fact = sum( shape_params$weights[!is.na(scal_product)] )
+ sum(scal_product, na.rm=TRUE) / norm_fact
+}
--- /dev/null
+#' Obtain delta forecast by the Persistence method
+#'
+#' @inheritParams getForecast
+#' @inheritParams getZeroDeltaForecast
+getPersistenceDeltaForecast = function(data, today, memory, horizon, shape_params, ...)
+{
+ # Return level of last similar "tomorrow" (thus -6) computed on "horizon"
+ data$getSerie(today-6)[1] - tail(data$getSerie(today-7), 1)
+}
--- /dev/null
+#' Just predict zero deltas (for reference)
+#'
+#' @inheritParams getForecast
+#' @param today Index of the current day (predict tomorrow)
+#' @param shape_params Optional parameters returned by the shape forecaster
+getZeroDeltaForecast = function(data, today, memory, horizon, shape_params, ...)
+{
+ 0
+}
--- /dev/null
+#' @title Data
+#'
+#' @description Data encapsulation
+#'
+#' @field data List of
+#' \itemize{
+#' \item time: vector of times
+#' \item serie: centered series
+#' \item level: corresponding levels
+#' \item exo_hat: predicted exogenous variables
+#' \item exo_Dm1: List of measured exogenous variables at day minus 1
+#' }
+#'
+#' @exportClass Data
+#' @export Data
+Data = setRefClass(
+ Class = "Data",
+
+ fields = list(
+ data = "list"
+ ),
+
+ methods = list(
+ initialize = function(...)
+ {
+ "Initialize empty Data object"
+
+ callSuper(...)
+ },
+ getSize = function()
+ {
+ length(data)
+ },
+ append = function(new_time, new_serie, new_level, new_exo_hat, new_exo_Dm1)
+ {
+ "Acquire a new vector of lists (time, serie, level, exo_hat, exo_Dm1)"
+
+ data[[length(data)+1]] <<- list("time"=new_time,"serie"=new_serie,"level"=new_level,
+ "exo_hat"=new_exo_hat,"exo_Dm1"=new_exo_Dm1)
+ },
+ getTime = function(index)
+ {
+ "Get time values at specified index"
+
+ data[[index]]$time
+ },
+ getCenteredSerie = function(index)
+ {
+ "Get serie values (centered) at specified index"
+
+ data[[index]]$serie
+ },
+ getLevel = function(index)
+ {
+ "Get level for the serie at specified index"
+
+ data[[index]]$level
+ },
+ getSerie = function(index)
+ {
+ "Get serie values (centered+level) at specified index"
+
+ data[[index]]$serie + data[[index]]$level
+ },
+ getExoHat = function(index)
+ {
+ "Get exogeous predictions at specified index"
+
+ data[[index]]$exo_hat
+ },
+ getExoDm1 = function(index)
+ {
+ "Get exogenous measures the day before specified index"
+
+ data[[index]]$exo_Dm1
+ }
+ )
+)
--- /dev/null
+#' @title Forecast
+#'
+#' @description Forecast encapsulation
+#'
+#' @field pred List with
+#' \itemize{
+#' \item serie: forecasted serie
+#' \item params: corresponding list of parameters (weights, neighbors...)
+#' \item index: corresponding index in data object
+#' }
+#'
+#' @exportClass Forecast
+#' @export Forecast
+Forecast = setRefClass(
+ Class = "Forecast",
+
+ fields = list(
+ pred = "list"
+ ),
+
+ methods = list(
+ initialize = function(...)
+ {
+ "Initialize empty Forecast object"
+
+ callSuper(...)
+ },
+ append = function(new_serie, new_params, new_index)
+ {
+ "Obtain a new pair (serie, params)"
+
+ pred[[length(pred)+1]] <<- list("serie"=new_serie, "params"=new_params, "index"=new_index)
+ },
+ getSize = function()
+ {
+ length(pred)
+ },
+ getSerie = function(index)
+ {
+ "Get serie values at specified index"
+
+ pred[[index]]$serie
+ },
+ getParams = function(index)
+ {
+ "Get params at specified index"
+
+ pred[[index]]$params
+ },
+ getIndexInData = function(index)
+ {
+ "Get (day)index at specified index"
+
+ pred[[index]]$index
+ }
+ )
+)
--- /dev/null
+#' @include ShapeForecaster.R
+#'
+#' @title Average Shape Forecaster
+#'
+#' @description Return the (pointwise) average of the all the centered day curves in the past.
+#' Inherits \code{\link{ShapeForecaster}}
+AverageShapeForecaster = setRefClass(
+ Class = "AverageShapeForecaster",
+ contains = "ShapeForecaster",
+
+ methods = list(
+ initialize = function(...)
+ {
+ callSuper(...)
+ },
+ predict = function(today, memory, horizon, ...)
+ {
+ avg = rep(0., horizon)
+ for (i in (today-memory):today)
+ {
+ if (!any(is.na(data$getSerie(i))))
+ avg = avg + data$getCenteredSerie(i)
+ }
+ avg
+ }
+ )
+)
--- /dev/null
+#' @include ShapeForecaster.R
+#'
+#' @title Neighbors Shape Forecaster
+#'
+#' @description Predict tomorrow as a weighted combination of "futures of the past" days.
+#' Inherits \code{\link{ShapeForecaster}}
+NeighborsShapeForecaster = setRefClass(
+ Class = "NeighborsShapeForecaster",
+ contains = "ShapeForecaster",
+
+ methods = list(
+ initialize = function(...)
+ {
+ callSuper(...)
+ },
+ predict = function(today, memory, horizon, ...)
+ {
+ # (re)initialize computed parameters
+ params <<- list("weights"=NA, "indices"=NA, "window"=NA)
+
+ first_day = max(today - memory, 1)
+ # The first day is generally not complete:
+ if (length(data$getCenteredSerie(1)) < length(data$getCenteredSerie(2)))
+ first_day = 2
+
+ # Predict only on non-NAs days (TODO:)
+ if (any(is.na(data$getSerie(today))))
+ return (NA)
+
+ # Determine indices of no-NAs days followed by no-NAs tomorrows
+ fdays_indices = c()
+ for (i in first_day:(today-1))
+ {
+ if ( !any(is.na(data$getSerie(i)) | is.na(data$getSerie(i+1))) )
+ fdays_indices = c(fdays_indices, i)
+ }
+
+ #GET OPTIONAL PARAMS
+ # Similarity computed with exogenous variables ? endogenous ? both ? ("exo","endo","mix")
+ simtype = ifelse(hasArg("simtype"), list(...)$simtype, "exo")
+ simthresh = ifelse(hasArg("simthresh"), list(...)$simthresh, 0.)
+ kernel = ifelse(hasArg("kernel"), list(...)$kernel, "Gauss")
+ mix_strategy = ifelse(hasArg("mix_strategy"), list(...)$mix_strategy, "neighb") #or "mult"
+ same_season = ifelse(hasArg("same_season"), list(...)$same_season, TRUE)
+ if (hasArg(h_window))
+ return (.predictAux(fdays_indices, today, horizon, list(...)$h_window, kernel, simtype,
+ simthresh, mix_strategy, FALSE))
+ #END GET
+
+ # Indices for cross-validation; TODO: 45 = magic number
+ indices = getSimilarDaysIndices(today, limit=45, same_season=same_season)
+ #indices = (end_index-45):(end_index-1)
+
+ # Function to optimize h : h |--> sum of prediction errors on last 45 "similar" days
+ errorOnLastNdays = function(h, kernel, simtype)
+ {
+ error = 0
+ nb_jours = 0
+ for (i in indices)
+ {
+ # NOTE: predict only on non-NAs days followed by non-NAs (TODO:)
+ if (!any(is.na(data$getSerie(i)) | is.na(data$getSerie(i+1))))
+ {
+ nb_jours = nb_jours + 1
+ # mix_strategy is never used here (simtype != "mix"), therefore left blank
+ prediction = .predictAux(fdays_indices, i, horizon, h, kernel, simtype, simthresh,
+ "", FALSE)
+ if (!is.na(prediction[1]))
+ error = error + mean((data$getCenteredSerie(i+1)[1:horizon] - prediction)^2)
+ }
+ }
+ return (error / nb_jours)
+ }
+
+ h_best_exo = 1.
+ if (simtype != "endo" && !(simtype=="mix" && mix_strategy=="neighb"))
+ {
+ h_best_exo = optimize(errorOnLastNdays, interval=c(0,10), kernel=kernel,
+ simtype="exo")$minimum
+ }
+ if (simtype != "exo")
+ {
+ h_best_endo = optimize(errorOnLastNdays, interval=c(0,10), kernel=kernel,
+ simtype="endo")$minimum
+ }
+
+ if (simtype == "endo")
+ {
+ return (.predictAux(fdays_indices, today, horizon, h_best_endo, kernel, "endo",
+ simthresh, "", TRUE))
+ }
+ if (simtype == "exo")
+ {
+ return (.predictAux(fdays_indices, today, horizon, h_best_exo, kernel, "exo",
+ simthresh, "", TRUE))
+ }
+ if (simtype == "mix")
+ {
+ return (.predictAux(fdays_indices, today, horizon, c(h_best_endo,h_best_exo), kernel,
+ "mix", simthresh, mix_strategy, TRUE))
+ }
+ },
+ # Precondition: "today" is full (no NAs)
+ .predictAux = function(fdays_indices, today, horizon, h, kernel, simtype, simthresh,
+ mix_strategy, final_call)
+ {
+ dat = data$data #HACK: faster this way...
+
+ fdays_indices = fdays_indices[fdays_indices < today]
+ # TODO: 3 = magic number
+ if (length(fdays_indices) < 3)
+ return (NA)
+
+ if (simtype != "exo")
+ {
+ h_endo = ifelse(simtype=="mix", h[1], h)
+
+ # Distances from last observed day to days in the past
+ distances2 = rep(NA, length(fdays_indices))
+ for (i in seq_along(fdays_indices))
+ {
+ delta = dat[[today]]$serie - dat[[ fdays_indices[i] ]]$serie
+ # Require at least half of non-NA common values to compute the distance
+ if (sum(is.na(delta)) <= 0) #length(delta)/2)
+ distances2[i] = mean(delta^2) #, na.rm=TRUE)
+ }
+
+ sd_dist = sd(distances2)
+ simils_endo =
+ if (kernel=="Gauss") {
+ exp(-distances2/(sd_dist*h_endo^2))
+ } else { #Epanechnikov
+ u = 1 - distances2/(sd_dist*h_endo^2)
+ u[abs(u)>1] = 0.
+ u
+ }
+ }
+
+ if (simtype != "endo")
+ {
+ h_exo = ifelse(simtype=="mix", h[2], h)
+
+ # TODO: [rnormand] if predict_at == 0h then we should use measures from day minus 1
+ M = matrix( nrow=1+length(fdays_indices), ncol=1+length(dat[[today]]$exo_hat) )
+ M[1,] = c( dat[[today]]$level, as.double(dat[[today]]$exo_hat) )
+ for (i in seq_along(fdays_indices))
+ {
+ M[i+1,] = c( dat[[ fdays_indices[i] ]]$level,
+ as.double(dat[[ fdays_indices[i] ]]$exo_hat) )
+ }
+
+ sigma = cov(M) #NOTE: robust covariance is way too slow
+ sigma_inv = qr.solve(sigma)
+
+ # Distances from last observed day to days in the past
+ distances2 = rep(NA, nrow(M)-1)
+ for (i in 2:nrow(M))
+ {
+ delta = M[1,] - M[i,]
+ distances2[i-1] = delta %*% sigma_inv %*% delta
+ }
+
+ sd_dist = sd(distances2)
+ simils_exo =
+ if (kernel=="Gauss") {
+ exp(-distances2/(sd_dist*h_exo^2))
+ } else { #Epanechnikov
+ u = 1 - distances2/(sd_dist*h_exo^2)
+ u[abs(u)>1] = 0.
+ u
+ }
+ }
+
+ if (simtype=="mix")
+ {
+ if (mix_strategy == "neighb")
+ {
+ #Only (60) most similar days according to exogen variables are kept into consideration
+ #TODO: 60 = magic number
+ keep_indices = sort(simils_exo, index.return=TRUE)$ix[1:(min(60,length(simils_exo)))]
+ simils_endo[-keep_indices] = 0.
+ } else #mix_strategy == "mult"
+ {
+ simils_endo = simils_endo * simils_exo
+ }
+ }
+
+ similarities =
+ if (simtype != "exo") {
+ simils_endo
+ } else {
+ simils_exo
+ }
+
+ if (simthresh > 0.)
+ {
+ max_sim = max(similarities)
+ # Set to 0 all similarities s where s / max_sim < simthresh, but keep at least 60
+ ordering = sort(similarities / max_sim, index.return=TRUE)
+ if (ordering[60] < simthresh)
+ {
+ similarities[ ordering$ix[ - (1:60) ] ] = 0.
+ } else
+ {
+ limit = 61
+ while (limit < length(similarities) && ordering[limit] >= simthresh)
+ limit = limit + 1
+ similarities[ ordering$ix[ - 1:limit] ] = 0.
+ }
+ }
+
+ prediction = rep(0, horizon)
+ for (i in seq_along(fdays_indices))
+ prediction = prediction + similarities[i] * dat[[ fdays_indices[i]+1 ]]$serie[1:horizon]
+
+ prediction = prediction / sum(similarities, na.rm=TRUE)
+ if (final_call)
+ {
+ params$weights <<- similarities
+ params$indices <<- fdays_indices
+ params$window <<-
+ if (simtype=="endo") {
+ h_endo
+ } else if (simtype=="exo") {
+ h_exo
+ } else {
+ c(h_endo,h_exo)
+ }
+ }
+ return (prediction)
+ }
+ )
+)
--- /dev/null
+#' @include ShapeForecaster.R
+#'
+#' @title Persistence Shape Forecaster
+#'
+#' @description Return the last centered last (similar) day curve (may be a lot of NAs,
+#' but we cannot do better). Inherits \code{\link{ShapeForecaster}}
+PersistenceShapeForecaster = setRefClass(
+ Class = "PersistenceShapeForecaster",
+ contains = "ShapeForecaster",
+
+ methods = list(
+ initialize = function(...)
+ {
+ callSuper(...)
+ },
+ predict = function(today, memory, horizon, ...)
+ {
+ #return centered last (similar) day curve (may be a lot of NAs, but we cannot do better)
+ last_similar_serie = data$getCenteredSerie(today-6)[1:horizon]
+ if (any(is.na(last_similar_serie))) #TODO:
+ return (NA)
+ last_similar_serie
+ }
+ )
+)
--- /dev/null
+#' @title Shape Forecaster
+#'
+#' @description Generic class to represent a shape forecaster
+#'
+#' @field params List of computed parameters, potentially useful for some DeltaForecasters
+#' @field data Dataset, object of class Data
+ShapeForecaster = setRefClass(
+ Class = "ShapeForecaster",
+
+ fields = list(
+ params = "list",
+ data = "Data"
+ ),
+
+ methods = list(
+ initialize = function(...)
+ {
+ "Initialize (generic) ShapeForecaster object"
+
+ callSuper(...)
+ if (!hasArg(data))
+ stop("ShapeForecaster must be initialized with a Data object")
+ params <<- list()
+ },
+ predict = function(today, memory, horizon, ...)
+ {
+ "Obtain a new forecasted time-serie (+ side-effect: compute parameters)"
+
+ #empty default implementation: to implement in inherited classes
+ },
+ getParameters = function()
+ {
+ params
+ }
+ )
+)
+
--- /dev/null
+#' @title Acquire data in a clean format
+#'
+#' @description Take in input data frames and/or files containing raw data, and timezones, and
+#' output a Data object, roughly corresponding to a list where each cell contains all value
+#' for one day (see \code{?Data}). Current limitation: series (in working_tz) must start at
+#' right after midnight (to keep in sync with exogenous vars)
+#'
+#' @param ts_data Time-series, as a data frame (DB style: 2 columns, first is date/time,
+#' second is value) or a CSV file
+#' @param exo_data Exogenous variables, as a data frame or a CSV file; first comlumn is dates,
+#' next block are measurements for the day, and final block are exogenous forecasts
+#' @param input_tz Timezone in the input files ("GMT" or e.g. "Europe/Paris")
+#' @param date_format How date/time are stored (e.g. year/month/day hour:minutes;
+#' see \code{strptime})
+#' @param working_tz Timezone to work with ("GMT" or e.g. "Europe/Paris")
+#' @param predict_at When does the prediction take place ? ab[:cd][:ef] where a,b,c,d,e,f
+#' in (0,9) and define an hour[minute[second]]; time must be present in the file
+#'
+#' @return An object of class Data
+#'
+#' @export
+getData = function(ts_data, exo_data,
+ input_tz="GMT", date_format="%d/%m/%Y %H:%M", working_tz="GMT", predict_at="00")
+{
+ # Sanity checks (not full, but sufficient at this stage)
+ if (!is.character(input_tz) || !is.character(working_tz))
+ stop("Bad timezone (see ?timezone)")
+ input_tz = input_tz[1]
+ working_tz = working_tz[1]
+ if (!is.data.frame(ts_data) && !is.character(ts_data))
+ stop("Bad time-series input (data frame or CSV file)")
+ if (is.character(ts_data))
+ ts_data = ts_data[1]
+ pattern_index_in_predict_at = grep("^[0-9]{2}(:[0-9]{2}){0,2}$", predict_at)
+ if (!is.character(predict_at) || length(pattern_index_in_predict_at) == 0)
+ stop("Bad predict_at ( ^[0-9]{2}(:[0-9]{2}){0,2}$ )")
+ predict_at = predict_at[ pattern_index_in_predict_at[1] ]
+ if (!is.character(date_format))
+ stop("Bad date_format (character)")
+ date_format = date_format[1]
+
+ ts_df =
+ if (is.character(ts_data)) {
+ read.csv(ts_data)
+ } else {
+ ts_data
+ }
+ exo_df =
+ if (is.character(exo_data)) {
+ read.csv(exo_data)
+ } else {
+ exo_data
+ }
+ # Convert to the desired timezone (usually "GMT" or "Europe/Paris")
+ formatted_dates_POSIXlt = strptime(as.character(ts_df[,1]), date_format, tz=input_tz)
+ ts_df[,1] = format(as.POSIXct(formatted_dates_POSIXlt), tz=working_tz, usetz=TRUE)
+
+ if (nchar(predict_at) == 2)
+ predict_at = paste(predict_at,":00",sep="")
+ if (nchar(predict_at) == 5)
+ predict_at = paste(predict_at,":00",sep="")
+
+ line = 1 #index in PM10 file (24 lines for 1 cell)
+ nb_lines = nrow(ts_df)
+ nb_exos = ( ncol(exo_df) - 1 ) / 2
+ data = list() #new("Data")
+ i = 1 #index of a cell in data
+ while (line <= nb_lines)
+ {
+ time = c()
+ serie = c()
+ repeat {
+ {
+ time = c(time, ts_df[line,1])
+ serie = c(serie, ts_df[line,2])
+ line = line + 1
+ };
+ if (line >= nb_lines + 1
+ # NOTE: always second part of date/time, because it has been formatted
+ || strsplit(as.character(ts_df[line-1,1])," ")[[1]][2] == predict_at)
+ {
+ break
+ }}
+
+ # NOTE: if predict_at does not cut days at midnight,
+ # for the exogenous to be synchronized they need to be shifted
+ if (predict_at != "00:00:00")
+ {
+ exo_hat = as.data.frame(exo_df[max(1,i-1),(1+nb_exos+1):(1+2*nb_exos)])
+ exo_Dm1 = if (i>=3) as.data.frame(exo_df[i-1,2:(1+nb_exos)]) else NA
+ } else
+ {
+ exo_hat = as.data.frame(exo_df[i,(1+nb_exos+1):(1+2*nb_exos)])
+ exo_Dm1 = if (i>=2) as.data.frame(exo_df[i-1,2:(1+nb_exos)]) else NA
+ }
+ i = i + 1
+ #center data
+ level = mean(serie, na.rm=TRUE)
+ centered_serie = serie - level
+# data$append(time, centered_serie, level, exo_hat, exo_Jm1) #TODO: slow: why ?
+ data[[length(data)+1]] = list("time"=time, "serie"=centered_serie, "level"=level,
+ "exo_hat"=exo_hat, "exo_Dm1"=exo_Dm1)
+ }
+ new("Data",data=data)
+}
--- /dev/null
+#' @title Get error
+#'
+#' @description Obtain the errors between forecast and data
+#'
+#' @param data Dataset, object of class \code{Data} output of \code{getData}
+#' @param forecast Forecast object, class \code{Forecast} output of \code{getForecast}
+#' @param horizon Horizon where to compute the error (<= horizon used in \code{getForecast})
+#'
+#' @return A list (abs,MAPE) of lists (day,indices)
+#'
+#' @export
+getError = function(data, forecast, horizon)
+{
+ L = forecast$getSize()
+ mape_day = rep(0, horizon)
+ abs_day = rep(0, horizon)
+ mape_indices = rep(NA, L)
+ abs_indices = rep(NA, L)
+
+ nb_no_NA_data = 0
+ for (i in seq_len(L))
+ {
+ index = forecast$getIndexInData(i)
+ serie = data$getSerie(index+1)[1:horizon]
+ pred = forecast$getSerie(i)[1:horizon]
+ if (!any(is.na(serie)) && !any(is.na(pred)))
+ {
+ nb_no_NA_data = nb_no_NA_data + 1
+ mape_increment = abs(serie - pred) / serie
+ mape_increment[is.nan(mape_increment)] = 0. # 0 / 0
+ mape_increment[!is.finite(mape_increment)] = 1. # >0 / 0
+ mape_day = mape_day + mape_increment
+ abs_increment = abs(serie - pred)
+ abs_day = abs_day + abs_increment
+ mape_indices[i] = mean(mape_increment)
+ abs_indices[i] = mean(abs_increment)
+ }
+ }
+
+ list(
+ "abs" = list(
+ "day" = abs_day / nb_no_NA_data,
+ "indices" = abs_indices),
+ "MAPE" = list(
+ "day" = mape_day / nb_no_NA_data,
+ "indices" = mape_indices) )
+}
--- /dev/null
+#' @title get Forecast
+#'
+#' @description Predict time-series curves for the selected days indices (lines in data).
+#' Run the forecasting task described by \code{delta_forecaster_name} and
+#' \code{shape_forecaster_name} on data obtained with \code{getData}
+#'
+#' @param data Dataset, object of type \code{Data} output of \code{getData}
+#' @param indices Days indices where to forecast (the day after)
+#' @param memory Data depth (in days) to be used for prediction
+#' @param horizon Number of time steps to predict
+#' @param shape_forecaster_name Name of the shape forcaster
+#' \itemize{
+#' \item Persistence : use values of last (similar, next) day
+#' \item Neighbors : use PM10 from the k closest neighbors' tomorrows
+#' \item Average : global average of all the (similar) "tomorrow of past"
+#' }
+#' @param delta_forecaster_name Name of the delta forecaster
+#' \itemize{
+#' \item Persistence : use last (similar) day values
+#' \item Neighbors: re-use the weights optimized in corresponding shape forecaster
+#' \item Zero: just output 0 (no adjustment)
+#' }
+#' @param ... Additional parameters for the forecasting models
+#'
+#' @return An object of class Forecast
+#'
+#' @examples
+#' data = getData(ts_data="data/pm10_mesures_H_loc.csv", exo_data="data/meteo_extra_noNAs.csv",
+#' input_tz = "Europe/Paris", working_tz="Europe/Paris", predict_at="07")
+#' pred = getForecast(data, 2200:2230, Inf, 12, "Persistence", "Persistence")
+#' \dontrun{#Sketch for real-time mode:
+#' data = new("Data", ...)
+#' repeat {
+#' data$append(some_new_data)
+#' pred = getForecast(data, ...)
+#' #do_something_with_pred
+#' }}
+#' @export
+getForecast = function(data, indices, memory, horizon,
+ shape_forecaster_name, delta_forecaster_name, ...)
+{
+ horizon = as.integer(horizon)[1]
+ if (horizon<=0 || horizon>length(data$getCenteredSerie(2)))
+ stop("Horizon too short or too long")
+ indices = as.integer(indices)
+ if (any(indices<=0 | indices>data$getSize()))
+ stop("Indices out of range")
+ indices = sapply(indices, dateIndexToInteger, data)
+
+ #NOTE: some assymetry here...
+ shape_forecaster = new(paste(shape_forecaster_name,"ShapeForecaster",sep=""), data=data)
+ #A little bit strange, but match.fun() and get() fail
+ delta_forecaster = getFromNamespace(
+ paste("get",delta_forecaster_name,"DeltaForecast",sep=""), "talweg")
+
+ pred = list()
+ for (today in indices)
+ {
+ #NOTE: To estimate timing...
+# print(paste("Predict until index",today))
+
+ #shape always predicted first (on centered series, no scaling taken into account),
+ #with side-effect: optimize some parameters (h, weights, ...)
+ predicted_shape = shape_forecaster$predict(today, memory, horizon, ...)
+ #then, delta prediction can re-use some variables optimized previously (like neighbors infos)
+ predicted_delta = delta_forecaster(data, today, memory, horizon,
+ shape_forecaster$getParameters(), ...)
+
+ #TODO: this way is faster than a call to append(); why ?
+ pred[[length(pred)+1]] = list(
+ # Predict shape and align it on end of current day
+ serie = predicted_shape + tail( data$getSerie(today), 1 ) - predicted_shape[1],
+ predicted_delta, #add predicted jump
+ params = shape_forecaster$getParameters(),
+ index = today )
+ }
+ new("Forecast",pred=pred)
+}
--- /dev/null
+#' @title plot measured / predicted
+#'
+#' @description Plot measured curve (in black) and predicted curve (in red)
+#'
+#' @param data Object return by \code{getData}
+#' @param pred Object as returned by \code{getForecast}
+#' @param index Index in forecasts
+#'
+#' @export
+plotPredReal <- function(data, pred, index)
+{
+ horizon = length(pred$getSerie(1))
+ par(mar=c(4.7,5,1,1), cex.axis=2, cex.lab=2, lwd=2)
+ measure = data$getSerie(pred$getIndexInData(index)+1)[1:horizon]
+ yrange = range( pred$getSerie(index), measure )
+ plot(measure, type="l", ylim=yrange, lwd=3)
+ par(new=TRUE)
+ plot(pred$getSerie(index), type="l", col=2, ylim=yrange, lwd=3)
+}
+
+#' @title Plot filaments
+#'
+#' @description Plot similar days in the past + "past tomorrow", as black as distances are small
+#'
+#' @param data Object as returned by \code{getData}
+#' @param index Index in data
+#' @param limit Number of neighbors to consider
+#'
+#' @export
+plotFilaments <- function(data, index, limit=60)
+{
+ index = dateIndexToInteger(index, data)
+ ref_serie = data$getCenteredSerie(index)
+ if (any(is.na(ref_serie)))
+ stop("plotFilaments requires a serie without NAs")
+ L = length(ref_serie)
+ first_day = ifelse(length(data$getCenteredSerie(1)<L), 2, 1)
+ distances = sapply(first_day:(index-1), function(i) {
+ sqrt( sum( (ref_serie - data$getCenteredSerie(i))^2 ) / L )
+ })
+ # HACK to suppress NA effect while keeping indexation
+ distances[is.na(distances)] = max(distances,na.rm=TRUE) + 1
+ indices = sort(distances, index.return=TRUE)$ix[1:min(limit,index-first_day)]
+ yrange = range( ref_serie, sapply( indices, function(i) {
+ index = i - first_day + 1
+ serie = c(data$getCenteredSerie(index), data$getCenteredSerie(index+1))
+ if (!all(is.na(serie)))
+ return ( range(serie, na.rm=TRUE) )
+ return (0)
+ }) )
+ grays = gray.colors(20, 0.1, 0.9) #TODO: 20 == magic number
+ colors = c(
+ grays[ floor( 20.5 * distances[indices] / (1+max(distances[indices])) ) ], "#FF0000")
+ par(mar=c(4.7,5,1,1), cex.axis=2, cex.lab=2, lwd=2)
+ for (i in seq_len(length(indices)+1))
+ {
+ ind = ifelse(i<=length(indices), indices[i] - first_day + 1, index)
+ plot(c(data$getCenteredSerie(ind),data$getCenteredSerie(ind+1)),
+ ylim=yrange, type="l", col=colors[i],
+ xlab=ifelse(i==1,"Temps (en heures)",""), ylab=ifelse(i==1,"PM10 centré",""))
+ if (i <= length(indices))
+ par(new=TRUE)
+ }
+}
+
+#' @title Plot similarities
+#'
+#' @description Plot histogram of similarities (weights)
+#'
+#' @param pred Object as returned by \code{getForecast}
+#' @param index Index in forecasts (not in data)
+#'
+#' @export
+plotSimils <- function(pred, index)
+{
+ weights = pred$getParams(index)$weights
+ if (is.null(weights))
+ stop("plotSimils only works on 'Neighbors' forecasts")
+ par(mar=c(4.7,5,1,1))
+ hist(pred$getParams(index)$weights, nclass=20, xlab="Weight", ylab="Frequency")
+}
+
+#' @title Plot error
+#'
+#' @description Draw error graphs, potentially from several runs of \code{getForecast}
+#'
+#' @param err Error as returned by \code{getError}
+#'
+#' @seealso \code{\link{plotPredReal}}, \code{\link{plotFilaments}}, \code{\link{plotSimils}}
+#' \code{\link{plotFbox}}
+#'
+#' @export
+plotError <- function(err)
+{
+ par(mfrow=c(2,2), mar=c(4.7,5,1,1), cex.axis=2, cex.lab=2, lwd=2)
+ L = length(err)
+ yrange = range( sapply(1:L, function(index) ( err[[index]]$abs$day ) ), na.rm=TRUE )
+ for (i in seq_len(L))
+ {
+ plot(err[[i]]$abs$day, type="l", xlab=ifelse(i==1,"Temps (heures)",""),
+ ylab=ifelse(i==1,"Moyenne |y - y_hat|",""), ylim=yrange, col=i)
+ if (i < L)
+ par(new=TRUE)
+ }
+ yrange = range( sapply(1:L, function(index) ( err[[index]]$abs$indices ) ), na.rm=TRUE )
+ for (i in seq_len(L))
+ {
+ plot(err[[i]]$abs$indices, type="l", xlab=ifelse(i==1,"Temps (jours)",""),
+ ylab=ifelse(i==1,"Moyenne |y - y_hat|",""), ylim=yrange, col=i)
+ if (i < L)
+ par(new=TRUE)
+ }
+ yrange = range( sapply(1:L, function(index) ( err[[index]]$MAPE$day ) ), na.rm=TRUE )
+ for (i in seq_len(L))
+ {
+ plot(err[[i]]$MAPE$day, type="l", xlab=ifelse(i==1,"Temps (heures)",""),
+ ylab=ifelse(i==1,"MAPE moyen",""), ylim=yrange, col=i)
+ if (i < L)
+ par(new=TRUE)
+ }
+ yrange = range( sapply(1:L, function(index) ( err[[index]]$MAPE$indices ) ), na.rm=TRUE )
+ for (i in seq_len(L))
+ {
+ plot(err[[i]]$MAPE$indices, type="l", xlab=ifelse(i==1,"Temps (jours)",""),
+ ylab=ifelse(i==1,"MAPE moyen",""), ylim=yrange, col=i)
+ if (i < L)
+ par(new=TRUE)
+ }
+}
+
+#' @title Functional boxplot
+#'
+#' @description Draw the functional boxplot on the left, and bivariate plot on the right
+#'
+#' @param data Object return by \code{getData}
+#' @param fiter Optional filter: return TRUE on indices to process
+#'
+#' @export
+plotFbox <- function(data, filter=function(index) (TRUE))
+{
+ if (!requireNamespace("rainbow", quietly=TRUE))
+ stop("Functional boxplot requires the rainbow package")
+
+ start_index = 1
+ end_index = data$getSize()
+ if (length(data$getCenteredSerie(1)) < length(data$getCenteredSerie(2)))
+ {
+ # Shifted start (7am, or 1pm, or...)
+ start_index = 2
+ end_index = data$getSize() - 1
+ }
+
+ series_matrix = sapply(start_index:end_index, function(index) {
+ as.matrix(data$getSerie(index))
+ })
+ # Remove NAs. + filter TODO: merge with previous step: only one pass required...
+ nas_indices = seq_len(ncol(series_matrix))[ sapply( 1:ncol(series_matrix),
+ function(index) ( !filter(index) || any(is.na(series_matrix[,index])) ) ) ]
+ series_matrix = series_matrix[,-nas_indices]
+
+ series_fds = rainbow::fds(seq_len(nrow(series_matrix)), series_matrix)
+ par(mfrow=c(1,2), mar=c(4.7,5,1,1), cex.axis=2, cex.lab=2)
+ rainbow::fboxplot(series_fds, "functional", "hdr", xlab="Temps (heures)", ylab="PM10",
+ plotlegend=FALSE, lwd=2)
+ rainbow::fboxplot(series_fds, "bivariate", "hdr", plotlegend=FALSE)
+}
--- /dev/null
+#' @title dateIndexToInteger
+#'
+#' @description Transform a (potential) date index into an integer (relative to data)
+#'
+#' @param index Date (or integer) index
+#' @param data Object of class \code{Data}
+#'
+#' @export
+dateIndexToInteger = function(index, data)
+{
+ if (is.numeric(index))
+ index = as.integer(index)
+ #Undocumented "private" method to translate index --> date (is needed)
+ if (is.integer(index))
+ return (index[1])
+ if (is.character(index))
+ {
+ tryCatch(dt <- as.POSIXct(index[1]), finally=print("Unrecognized index format"))
+ #TODO: tz arg to difftime ?
+ integerIndex <- as.integer( difftime(dt, data$getTime(1)) ) + 1
+ if (integerIndex > 0 && integerIndex <= data$getSize())
+ return (integerIndex)
+ stop("Date outside data range")
+ }
+ stop("Unrecognized index format")
+}
+
+#' @title getSimilarDaysIndices
+#'
+#' @description Find similar days indices in the past
+#'
+#' @param index Day index (numeric or date)
+#' @param limit Maximum number of indices to return
+#' @param same_seaon Should the indices correspond to day in same season?
+#'
+#' @export
+getSimilarDaysIndices = function(index, limit, same_season)
+{
+ index = dateIndexToInteger(index)
+
+ #TODO: mardi similaire à lundi mercredi jeudi aussi ...etc
+ if (!same_season)
+ {
+ #take all similar days in recent past
+ nb_days = min( (index-1) %/% 7, limit)
+ return ( rep(index,nb_days) - 7*seq_len(nb_days) )
+ }
+
+ #Look for similar days in similar season (+/- 30 days)
+ days = c()
+ i = index
+ while (i >= 1 && length(days) < limit)
+ {
+ if (i < index)
+ {
+ days = c(days, i)
+ #look in the "future of the past"
+ for (j in 1:4)
+ days = c(days, i+7*j)
+ }
+ #...and in the "past of the past"
+ for (j in 1:4)
+ {
+ if (i - 7*j >= 1)
+ days = c(days, i-7*j)
+ }
+ # TODO: exact computation instead of -364
+ # 364 = closest multiple of 7 to 365 - drift along the years... but not so many years so OK
+ i = i - 364
+ }
+
+ return ( days[1:min(limit,length(days))] )
+}
--- /dev/null
+# Predict PM10 as functions of time
+
+Joint work with [Jean-Michel Poggi](http://www.math.u-psud.fr/~poggi/) and [Bruno Portier](http://lmi2.insa-rouen.fr/~bportier/)
+
+---
+
+TODO: ...
+
+---
+
+TODO: ...
+
+The final report can be found at [this location](http://www.airnormand.fr/Publications/Publications-telechargeables/Rapports-d-etudes)
--- /dev/null
+Données (compressées) dans data.tar.xz
+======================================
+
+pm10_mesures_[Q]H[_loc].csv : une mesure par ligne (quelques plages de NAs),
+du 10 décembre 2008 au 31 juillet 2016, en heure locale ou GMT.
+
+meteo[_extra[_noNA]].csv : les données météo initiales (mesures + prédictions).
+Mêmes dates ; noNA : avec NA remplacés par moyennes.
++ trois variables calculées,
+ * Season: 1 d'avril à août, 0 de septembre à mars
+ * Pollution: -1 si indéterminée (que des NAs), 0 si <30, 1 si entre 30 et 50, 2 si >50
+ * Week: 0 si week-end, 1 si jour de semaine
+(J'ai essayé d'appliquer en gros "0 = moins pollué, 1 = plus pollué")
--- /dev/null
+#$# git-fat fb2d21849524e743ce1b3e5589efb40d1182f468 564468
--- /dev/null
+meteo_df = read.csv("meteo.csv")
+
+meteo_df$Season = 0
+meteo_df$Week = 0
+meteo_df$Pollution = -1
+
+#Need to load and aggregate PM10 by days
+pm10_df = read.csv("pm10_mesures.csv")
+
+line_number = 1 #line number in pm10 file
+for (i in 1:(nrow(meteo_df)))
+{
+ pm10s = c()
+ repeat {
+ {
+ pm10s = c(pm10s, pm10_df[line_number,2])
+ line_number = line_number + 1
+ };
+ if (line_number >= nrow(pm10_df)+1
+ || strsplit(as.character(pm10_df[line_number,1])," ")[[1]][2] == '0:15')
+ {
+ break
+ }}
+ pm10_level = mean(pm10s, na.rm=TRUE)
+ #Fill Pollution column: -1 if no info, 0 to 2 for pollution level
+ if (!is.nan(pm10_level))
+ {
+ if (pm10_level < 30)
+ meteo_df$Pollution[i] = 0
+ else if (pm10_level <= 50)
+ meteo_df$Pollution[i] = 1
+ else #pm10 > 50
+ meteo_df$Pollution[i] = 2
+ }
+
+ #Also fill season + days of week variables
+ meteo_df$Season[i] = ifelse(
+ strsplit(as.character(meteo_df$Date[i]),'/')[[1]][1] %in% c("4","5","6","7","8"),
+ 1, 0)
+ current_datetime = strptime(as.character(meteo_df$Date[i]), "%m/%d/%Y", tz="GMT")
+ meteo_df$Week[i] = ifelse(current_datetime$wday %in% c(6,0), 0, 1)
+}
+
+#see also:
+#http://stackoverflow.com/questions/8214303/conditional-replacement-of-values-in-a-data-frame
+
+#Finally write new data
+write.csv(meteo_df, file="meteo_extra.csv", row.names=FALSE)
--- /dev/null
+fill_NAs = function(file)
+{
+ dat_with_date = read.csv(file)
+ dat = dat_with_date[,-1] #get rid of date column
+
+ line = 1
+ n = nrow(dat)
+ p = ncol(dat)
+ last_noNA_indices = rep(0, p)
+ while (line <= n)
+ {
+# print(line)
+ for (i in 1:p)
+ {
+# if (is.numeric(dat[line,i]))
+ if (!is.na(dat[line,i]))
+ {
+ if (last_noNA_indices[i] < line-1)
+ {
+ #do some completion
+ meanVal = ifelse(last_noNA_indices[i]>0,
+ 0.5*(dat[last_noNA_indices[i],i]+dat[line,i]),
+ dat[line,i])
+
+ for (j in (last_noNA_indices[i]+1):(line-1))
+ dat[j,i] = meanVal
+ }
+ last_noNA_indices[i] = line
+ }
+ }
+ line = line + 1
+ }
+
+ #complete until end if needed
+ for (i in 1:p)
+ {
+ if (last_noNA_indices[i] < n)
+ {
+ for (j in (last_noNA_indices[i]+1):n)
+ dat[j,i] = dat[last_noNA_indices[i],i]
+ }
+ }
+
+ dat_with_date[,2:(p+1)] = dat
+ write.csv(dat_with_date, "NONA.csv", row.names=FALSE)
+}
--- /dev/null
+\name{talweg-package}
+\alias{talweg-package}
+\alias{talweg}
+\docType{package}
+
+\title{
+ \packageTitle{talweg}
+}
+
+\description{
+ \packageDescription{talweg}
+}
+
+\details{
+ The package devtools should be useful in development stage, since we rely on testthat for
+ unit tests, and roxygen2 for documentation. knitr is used to generate the package vignette.
+ Concerning the other suggested packages:
+ \itemize{
+ \item{TODO...;}
+ \item{TODO...;}
+ }
+
+ The three main functions are located in R/main.R:
+ \itemize{
+ \item{TODO...;}
+ \item{TODO...;}
+ }
+}
+
+\author{
+ \packageAuthor{talweg}
+
+ Maintainer: \packageMaintainer{talweg}
+}
+
+%\references{
+% TODO: Literature or other references for background information
+%}
+
+%\examples{
+% TODO: simple examples of the most important functions
+%}
--- /dev/null
+#$# git-fat c45d30e178d516079466abb2b1d39ff10f1f0ade 96701
--- /dev/null
+\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
+
+%%%%%%% FROM KNITR
+
+\usepackage[pdftex]{graphicx}
+\usepackage[]{color}
+%% maxwidth is the original width if it is less than linewidth
+%% otherwise use linewidth (to make sure the graphics do not exceed the margin)
+\makeatletter
+\def\maxwidth{ %
+ \ifdim\Gin@nat@width>\linewidth
+ \linewidth
+ \else
+ \Gin@nat@width
+ \fi
+}
+\makeatother
+
+\definecolor{fgcolor}{rgb}{0.345, 0.345, 0.345}
+\newcommand{\hlnum}[1]{\textcolor[rgb]{0.686,0.059,0.569}{#1}}%
+\newcommand{\hlstr}[1]{\textcolor[rgb]{0.192,0.494,0.8}{#1}}%
+\newcommand{\hlcom}[1]{\textcolor[rgb]{0.678,0.584,0.686}{\textit{#1}}}%
+\newcommand{\hlopt}[1]{\textcolor[rgb]{0,0,0}{#1}}%
+\newcommand{\hlstd}[1]{\textcolor[rgb]{0.345,0.345,0.345}{#1}}%
+\newcommand{\hlkwa}[1]{\textcolor[rgb]{0.161,0.373,0.58}{\textbf{#1}}}%
+\newcommand{\hlkwb}[1]{\textcolor[rgb]{0.69,0.353,0.396}{#1}}%
+\newcommand{\hlkwc}[1]{\textcolor[rgb]{0.333,0.667,0.333}{#1}}%
+\newcommand{\hlkwd}[1]{\textcolor[rgb]{0.737,0.353,0.396}{\textbf{#1}}}%
+
+\usepackage{framed}
+\makeatletter
+\newenvironment{kframe}{%
+ \def\at@end@of@kframe{}%
+ \ifinner\ifhmode%
+ \def\at@end@of@kframe{\end{minipage}}%
+ \begin{minipage}{\columnwidth}%
+ \fi\fi%
+ \def\FrameCommand##1{\hskip\@totalleftmargin \hskip-\fboxsep
+ \colorbox{shadecolor}{##1}\hskip-\fboxsep
+ % There is no \\@totalrightmargin, so:
+ \hskip-\linewidth \hskip-\@totalleftmargin \hskip\columnwidth}%
+ \MakeFramed {\advance\hsize-\width
+ \@totalleftmargin\z@ \linewidth\hsize
+ \@setminipage}}%
+ {\par\unskip\endMakeFramed%
+ \at@end@of@kframe}
+\makeatother
+
+\definecolor{shadecolor}{rgb}{.97, .97, .97}
+\definecolor{messagecolor}{rgb}{0, 0, 0}
+\definecolor{warningcolor}{rgb}{1, 0, 1}
+\definecolor{errorcolor}{rgb}{1, 0, 0}
+\newenvironment{knitrout}{}{} % an empty environment to be redefined in TeX
+
+\usepackage{alltt}
+\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
+\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
+
+%%%%%%%% END FROM KNITR
+
+%%%%%%%%% page de garde
+
+%\usepackage[a4paper,left=2cm,right=2cm,top=2cm,bottom=2cm]{geometry}
+\usepackage[frenchb]{babel}
+%\usepackage[pdftex]{graphicx}
+
+\usepackage{Sweave}
+% \usepackage[dvips]{graphicx}
+\usepackage{float}
+%\usepackage{multido}
+
+\setlength{\parindent}{0cm}
+\setlength{\parskip}{1ex plus 0.5ex minus 0.2ex}
+%\newcommand{\hsp}{\hspace{20pt}}
+\newcommand{\HRule}{\rule{\linewidth}{0.5mm}}
+
+%%%%%%% end page de garde
+
+\include{PageDeGarde}
+\include{R_Package}
+
+\end{document}
--- /dev/null
+\begin{document}
+
+\begin{titlepage}
+ \begin{sffamily}
+ \begin{center}
+
+ % Upper part of the page. The '~' is needed because \\
+ % only works if a paragraph has started.
+ \includegraphics[width=3cm]{FIGURES/logo_upsud.png}\hspace{2cm}
+ \includegraphics[width=4cm]{FIGURES/logo_parisdescartes.png}~\\[0.5cm]
+ \includegraphics[width=3cm]{FIGURES/logo_insarouen.jpg}\hspace{2cm}
+ \includegraphics[width=4cm]{FIGURES/logo_airnormand.jpg}~\\[2cm]
+
+ \textsc{\Large Rapport intermédiaire 15/12/2016}\\[2cm]
+
+ % Title
+ \HRule \\[2cm]
+ { \LARGE \bfseries PM10 : prédiction de courbes journalières\\[2cm] }
+ {\normalsize
+ % Author and supervisor
+ \begin{minipage}{0.4\textwidth}
+ \begin{flushleft} \large
+ Benjamin AUDER\\
+ Jean-Michel POGGI\\
+ Bruno PORTIER\\
+ \end{flushleft}
+ \end{minipage}
+ \begin{minipage}{0.4\textwidth}
+ \begin{flushright} \large
+ Université Paris-Sud\\
+ Université Paris 5 \& Paris-Sud\\
+ INSA Rouen Normandie\\
+ \end{flushright}
+ \end{minipage}
+ }
+ \end{center}
+ \end{sffamily}
+\end{titlepage}
--- /dev/null
+#$# git-fat cdf287cbbb06e773d911a86fed964ba258f27814 5614008
--- /dev/null
+\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)
+@
+
+Note : sur la base de nos dernières expériences, on considère que
+\begin{itemize}
+ \item on ne touche pas à la fenêtre obtenue par optim() ;}
+ \item on oublie la méthode consistant à prédire forme et niveau de manière complètement
+ déconnectée : il faut relier les deux.
+\end{itemize}
+
+\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 (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 $1h00$ à $0h$ du
+jour courant (sauf si on prédit à 0h : on prend alors les valeurs mesurées de la veille).\\
+
+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("local", "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 saut d'une fin de série au début de la suivante.
+\end{itemize}
+
+\noindent Il faut ainsi préciser à la fois une méthode de prévision de forme ("Persistence" et
+"Neighbors" 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.\\
+
+\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}
+
+\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). Trois
+types d'erreurs sont considérées :
+\begin{itemize}
+ \item l'erreur "abs" égale à la valeur absolue moyenne entre la mesure et la prédiction ;
+ \item l'erreur "MAPE" égale à l'erreur absolue normalisée par la mesure.
+ \item l'erreur "RMSE" égale à la racine carrée de l'erreur quadratique moyenne.
+\end{itemize}
+
+\subsection{Expériences numériques}
+
+%, fig.show='hold'>>=
+<<xp1, out.width='18cm', out.height='6cm'>>=
+p_endo = predictPM10(data, 2200, 2230, 0,0, "Neighbors", "Neighbors", simtype="endo")
+p_exo = predictPM10(data, 2200, 2230, 0,0, "Neighbors", "Neighbors", simtype="exo")
+p_mix = predictPM10(data, 2200, 2230, 0,0, "Neighbors", "Neighbors", simtype="mix")
+p = c(p_endo, p_exo, p_mix)
+yrange_MAPE = range(p_mix$errors$MAPE, p_endo$errors$MAPE, p_exo$errors$MAPE)
+yrange_abs = range(p_mix$errors$abs, p_endo$errors$abs, p_exo$errors$abs)
+yrange_RMSE = range(p_mix$errors$RMSE, p_endo$errors$RMSE, p_exo$errors$RMSE)
+ranges = c(yrange_MAPE,yrange_abs,yrange_RMSE)
+par(mfrow=c(1,3))
+titles = paste("Erreur",c("MAPE","abs","RMSE"))
+for (i in 1:3) #error type (MAPE,abs,RMSE)
+{
+ for (j in 1:3) #model (mix,endo,exo)
+ {
+ plot(p[j]$errors[[i]], type="l", col=j, main=titles[i], xlab="Temps",
+ ylab="Erreur", ylim=ranges[i])
+ par(new=TRUE)
+ }
+}
+
+#Ne tenir compte que des similarités sur les variables exogènes semble
+#conduire à l'erreur la plus faible.
+@
+
+<<xp2, out.width='18cm', out.height='6cm'>>=
+p_nn = predictPM10(data, 2200, 2230, 0, 0, "Neighbors", "Neighbors", sameSeaon=TRUE)
+p_np = predictPM10(data, 2200, 2230, 0, 0, "Neighbors", "Persistence", sameSeaon=TRUE)
+p_nz = predictPM10(data, 2200, 2230, 0, 0, "Neighbors", "Zero", sameSeaon=TRUE)
+p_pp = predictPM10(data, 2200, 2230, 0, 0, "Persistence", "Persistence")
+p_pz = predictPM10(data, 2200, 2230, 0, 0, "Persistence", "Zero")
+p = c(p_nn, p_np, p_nz, p_pp, p_pz)
+yrange_MAPE = range(p_nn$errors$MAPE, p_nz$errors$MAPE, p_np$errors$MAPE, p_pp$errors$MAPE, p_pz$errors$MAPE)
+yrange_abs = range(p_nn$errors$abs, p_nz$errors$abs, p_np$errors$abs, p_pp$errors$abs, p_pz$errors$abs)
+yrange_RMSE = range(p_nn$errors$RMSE, p_nz$errors$RMSE, p_np$errors$RMSE, p_pp$errors$RMSE, p_pz$errors$RMSE)
+ranges = c(yrange_MAPE,yrange_abs,yrange_RMSE)
+par(mfrow=c(1,3))
+for (i in 1:3) #error type (MAPE,abs,RMSE)
+{
+ for (j in 1:5) #model (nn,np,nz,pp,pz)
+ {
+ plot(p[j]$errors[[i]], type="l", col=j, main=titles[i], xlab="Temps",
+ ylab="Erreur", ylim=ranges[i])
+ if (j<5)
+ par(new=TRUE)
+ }
+}
+
+#Meilleurs results: nn et nz (np moins bon)
+@
+
+%%TODO: analyse sur les trois périodes indiquées par Michel ; simtype=="exo" par defaut
+16/03/2015
+p_nn_epandage = predictPM10(data, 2200, 2200, 0, 0, "Neighbors", "Neighbors", sameSeaon=FALSE)
+19/01/2015
+p_nn_chauffage = predictPM10(data, 2200, 2200, 0, 0, "Neighbors", "Neighbors", sameSeaon=FALSE)
+23/02/2015
+p_nn_nonpollue = predictPM10(data, 2200, 2200, 0, 0, "Neighbors", "Neighbors", sameSeaon=FALSE)
+
+\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. \textcolor{blue}{OK : on le prend désormais en compte}\\
+
+\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.
+\textcolor{blue}{OK : on commence à en avoir quelques-unes}
+
+\end{document}
--- /dev/null
+{
+ "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
+}
--- /dev/null
+\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}
projects / talweg.git / commitdiff
