+++ /dev/null
-#$# git-fat cf4cf6c03e09a9aa2374a7250d21ce732c923e1e 49845703
+++ /dev/null
-#$# git-fat 9e89b67976812940a4329491351971ac550e629f 2781243
Package: talweg
Title: Time-Series Samples Forecasted With Exogenous Variables
-Version: 0.1-0
+Version: 0.2-0
Description: Forecast a curve sampled within the day (seconds, minutes,
hours...), using past measured curves + past exogenous informations, which
could be some aggregated measure on the past curves, the weather... Main
#' Data
#'
-#' Data encapsulation as a list (days) of lists (components).
+#' Data encapsulation, in the form of a few lists (time-series + exogenous variables).
#'
-#' The private field .data is a list where each cell contains the variables for a period
-#' of time of 24 hours, from time P+1 until P the next day where P is an integer between
-#' 0 and 23. .data[[1]] refers to the first measured data (serie and exogenous
-#' variables): the corresponding times can be accessed through the function
-#' \code{getTime()} below. Each cell .data[[i]] is itself a list containing five slots,
-#' as described in the 'field' section.
+#' The private field .tvp is a list where each cell contains the hourly variables for a
+#' period of time of 24 hours, from 1am to next midnight. The other lists contain
+#' informations on series' levels and exogenous variables (both measured and predicted).
#'
#' @usage # Data$new()
#'
-#' @field .data[[i]] List of
+#' @field .tvp List of "time-values"; in each cell:
#' \itemize{
#' \item time: vector of times
-#' \item centered_serie: centered serie
-#' \item level: corresponding level
-#' \item exo: exogenous variables
-#' \item exo_hat: predicted exogenous variables (for next day)
+#' \item serie: (measured) serie
+#' \item level_hat: predicted level for current day
#' }
+#' @field .level Vector of measured levels
+#' @field .exo List of measured exogenous variables, cell = numerical vector.
+#' @field .exo_hat List of predicted exogenous variables, cell = numerical vector.
#'
#' @section Methods:
#' \describe{
#' \item{\code{getSize()}}{
#' Number of series in dataset.}
-#' \item{\code{getStdHorizon()}}{
-#' Number of time steps from serie[1] until midnight}
-#' \item{\code{append(time, serie, exo, exo_hat)}}{
+#' \item{\code{append(time, value, level_hat, exo, exo_hat)}}{
#' Measured data for given vector of times + exogenous predictions from
#' last midgnight.}
#' \item{\code{getTime(index)}}{
#' Times (vector) at specified index.}
-#' \item{\code{getCenteredSerie(index)}}{
-#' Centered serie at specified index.}
-#' \item{\code{getCenteredSeries(indices)}}{
-#' Centered series at specified indices (in columns).}
-#' \item{\code{getLevel(index)}}{
-#' Level at specified index.}
#' \item{\code{getSerie(index)}}{
#' Serie (centered+level) at specified index.}
#' \item{\code{getSeries(indices)}}{
#' Series at specified indices (in columns).}
-#' \item{\code{getExoHat(index)}}{
-#' Predicted exogenous variables at specified index.}
+#' \item{\code{getLevel(index)}}{
+#' Measured level at specified index.}
+#' \item{\code{getLevelHat(index)}}{
+#' Predicted level vector at specified index (by hour).}
+#' \item{\code{getCenteredSerie(index)}}{
+#' Centered serie at specified index.}
+#' \item{\code{getCenteredSeries(indices)}}{
+#' Centered series at specified indices (in columns).}
#' \item{\code{getExo(index)}}{
#' Measured exogenous variables at specified index.}
-#' \item{\code{removeFirst()}}{
-#' Remove first list element (if truncated).}
-#' \item{\code{removeLast()}}{
-#' Remove last list element (if truncated).}
+#' \item{\code{getExoHat(index)}}{
+#' Predicted exogenous variables at specified index.}
#' }
#'
#' @docType class
#'
Data = R6::R6Class("Data",
private = list(
- .data = list()
+ .tvp = list(),
+ .level = vector("double",0),
+ .exo = list(),
+ .exo_hat = list()
),
public = list(
getSize = function()
- length(private$.data)
- ,
- getStdHorizon = function()
- 24 - as.POSIXlt( private$.data[[1]]$time[1] )$hour + 1
+ length(private$.tvp)
,
- append = function(time, serie, exo, exo_hat)
+ append = function(time=NULL, value=NULL, level_hat=NULL, exo=NULL, exo_hat=NULL)
{
- level = mean(serie, na.rm=TRUE)
- centered_serie = serie - level
- private$.data[[length(private$.data)+1]] <- list(
- "time"=time, #H-24 --> H-1
- "centered_serie"=centered_serie, #at 'time'
- "level"=level, #at 'time'
- "exo"=exo, #at 'time' (yersteday 0am to last midnight)
- "exo_hat"=exo_hat) #today 0am to next midnight
+ if (!is.null(time) && !is.null(value) && !is.null(level_hat))
+ {
+ L = length(private$.tvp)
+ if (L == 0 || strftime( tail(private$.tvp[[L]]$time,1),
+ format="%H:%M:%S", tz="GMT" ) == "00:00:00")
+ {
+ # Append a new cell
+ private$.tvp[[L+1]] <- list("time"=time, "serie"=value, "level_hat"=level_hat)
+ }
+ else
+ {
+ # Complete current cell
+ private$.tvp[[L]]$time <- c(private$.tvp[[L]]$time, time)
+ private$.tvp[[L]]$serie <- c(private$.tvp[[L]]$serie, value)
+ private$.tvp[[L]]$level_hat <- c(private$.tvp[[L]]$levem_hat, level_hat)
+ }
+ }
+ if (strftime( tail(private$.tvp[[length(private$.tvp)]]$time,1),
+ format="%H:%M:%S", tz="GMT" ) == "00:00:00")
+ {
+ private$.level = c(private$.level,
+ mean(private$.tvp[[length(private$.tvp)]]$serie, na.rm=TRUE))
+ }
+ if (!is.null(exo))
+ private$.exo[[length(private$.exo)+1]] = exo
+ if (!is.null(exo_hat))
+ private$.exo_hat[[length(private$.exo_hat)+1]] = exo_hat
},
getTime = function(index)
{
index = dateIndexToInteger(index, self)
- private$.data[[index]]$time
+ private$.tvp[[index]]$time
},
- getCenteredSerie = function(index)
+ getSerie = function(index)
{
index = dateIndexToInteger(index, self)
- private$.data[[index]]$centered_serie
+ private$.tvp[[index]]$serie
},
- getCenteredSeries = function(indices)
- sapply(indices, function(i) self$getCenteredSerie(i))
+ getSeries = function(indices)
+ sapply(indices, function(i) self$getSerie(i))
,
getLevel = function(index)
{
index = dateIndexToInteger(index, self)
- private$.data[[index]]$level
+ private$.level[index]
},
- getSerie = function(index)
+ getLevelHat = function(index)
{
index = dateIndexToInteger(index, self)
- private$.data[[index]]$centered_serie + private$.data[[index]]$level
+ private$.tvp[[index]]$level_hat
},
- getSeries = function(indices)
- sapply(indices, function(i) self$getSerie(i))
- ,
- getExoHat = function(index)
+ getCenteredSerie = function(index)
{
index = dateIndexToInteger(index, self)
- private$.data[[index]]$exo_hat
+ private$.tvp[[index]]$serie - private$.level[index]
},
+ getCenteredSeries = function(indices)
+ sapply(indices, function(i) self$getCenteredSerie(i))
+ ,
getExo = function(index)
{
index = dateIndexToInteger(index, self)
- private$.data[[index]]$exo
+ private$.exo[[index]]
},
- removeFirst = function()
- private$.data <- private$.data[2:length(private$.data)]
- ,
- removeLast = function()
- private$.data <- private$.data[1:(length(private$.data)-1)]
+ getExoHat = function(index)
+ {
+ index = dateIndexToInteger(index, self)
+ private$.exo_hat[[index]]
+ }
)
)
#' Pointwise average of all the series of the same day of week in the past.
#'
#' For example, if the current day (argument "today") is a tuesday, then all series
-#' corresponding to wednesdays in the past (until the beginning or memory limit) are
-#' averaged to provide a smooth prediction. This forecast will most of the time be wrong,
-#' but will also look plausible enough.
+#' corresponding to tuesday in the past (until the beginning or memory limit) -- and in
+#' the future if 'opera' is FALSE -- are averaged to provide a smooth prediction. This
+#' forecast will most of the time be wrong, but will also look plausible enough.
#'
#' @usage # AverageForecaster$new(pjump)
#'
inherit = Forecaster,
public = list(
- predictShape = function(data, today, memory, horizon, ...)
+ predictShape = function(data, today, memory, predict_from, horizon, ...)
{
- avg = rep(0., horizon)
+ avg = rep(0., (horizon-predict_from+1))
first_day = max(1, today-memory)
- index = today-7 + 1
+ index <- today
nb_no_na_series = 0
+ opera = ifelse(hasArg("opera"), list(...)$opera, FALSE)
repeat
{
- {
- serie_on_horizon = data$getCenteredSerie(index)[1:horizon]
- index = index - 7
- };
+ index = index - 7
+ if (index < first_day)
+ break
+ serie_on_horizon = data$getCenteredSerie(index)[predict_from:horizon]
if (!any(is.na(serie_on_horizon)))
{
avg = avg + serie_on_horizon
nb_no_na_series = nb_no_na_series + 1
- };
- if (index < first_day)
- break
+ }
+ }
+ if (!opera)
+ {
+ # The same, in the future
+ index <- today
+ repeat
+ {
+ index = index + 7
+ if (index > data$getSize())
+ break
+ serie_on_horizon = data$getCenteredSerie(index)[predict_from:horizon]
+ if (!any(is.na(serie_on_horizon)))
+ {
+ avg = avg + serie_on_horizon
+ nb_no_na_series = nb_no_na_series + 1
+ }
+ }
}
avg / nb_no_na_series
}
#' Neighbors Forecaster
#'
-#' Predict next serie as a weighted combination of "futures of the past" days,
-#' where days in the past are chosen and weighted according to some similarity measures.
+#' Predict next serie as a weighted combination of curves observed on "similar" days in
+#' the past (and future if 'opera'=FALSE); the nature of the similarity is controlled by
+#' the options 'simtype' and 'local' (see below).
#'
-#' The main method is \code{predictShape()}, taking arguments data, today, memory,
-#' horizon respectively for the dataset (object output of \code{getData()}), the current
-#' index, the data depth (in days) and the number of time steps to forecast.
-#' In addition, optional arguments can be passed:
+#' Optional arguments:
#' \itemize{
-#' \item local : TRUE (default) to constrain neighbors to be "same days within same
-#' season"
-#' \item simtype : 'endo' for a similarity based on the series only,<cr>
-#' 'exo' for a similaruty based on exogenous variables only,<cr>
+#' \item local: TRUE (default) to constrain neighbors to be "same days in same season"
+#' \item simtype: 'endo' for a similarity based on the series only,<cr>
+#' 'exo' for a similarity based on exogenous variables only,<cr>
#' 'mix' for the product of 'endo' and 'exo',<cr>
#' 'none' (default) to apply a simple average: no computed weights
-#' \item window : A window for similarities computations; override cross-validation
+#' \item window: A window for similarities computations; override cross-validation
#' window estimation.
#' }
#' The method is summarized as follows:
#' \enumerate{
-#' \item Determine N (=20) recent days without missing values, and followed by a
-#' tomorrow also without missing values.
+#' \item Determine N (=20) recent days without missing values, and preceded by a
+#' curve also without missing values.
#' \item Optimize the window parameters (if relevant) on the N chosen days.
#' \item Considering the optimized window, compute the neighbors (with locality
#' constraint or not), compute their similarities -- using a gaussian kernel if
inherit = Forecaster,
public = list(
- predictShape = function(data, today, memory, horizon, ...)
+ predictShape = function(data, today, memory, predict_from, horizon, ...)
{
# (re)initialize computed parameters
private$.params <- list("weights"=NA, "indices"=NA, "window"=NA)
# Do not forecast on days with NAs (TODO: softer condition...)
- if (any(is.na(data$getCenteredSerie(today))))
+ if (any(is.na(data$getSerie(today-1))) ||
+ (predict_from>=2 && any(is.na(data$getSerie(today)[1:(predict_from-1)]))))
+ {
return (NA)
-
- # Determine indices of no-NAs days followed by no-NAs tomorrows
- fdays = .getNoNA2(data, max(today-memory,1), today-1)
+ }
# Get optional args
local = ifelse(hasArg("local"), list(...)$local, TRUE) #same level + season?
simtype = ifelse(hasArg("simtype"), list(...)$simtype, "none") #or "endo", or "exo"
+ opera = ifelse(hasArg("opera"), list(...)$opera, FALSE) #operational mode?
+
+ # Determine indices of no-NAs days preceded by no-NAs yerstedays
+ tdays = .getNoNA2(data, max(today-memory,2), ifelse(opera,today-1,data$getSize()))
+ if (!opera)
+ tdays = setdiff(tdays, today) #always exclude current day
+
+ # Shortcut if window is known #TODO: cross-validation for number of days, on similar (yerste)days
if (hasArg("window"))
{
- return ( private$.predictShapeAux(data,
- fdays, today, horizon, local, list(...)$window, simtype, TRUE) )
+ return ( private$.predictShapeAux(data, tdays, today, predict_from, horizon,
+ local, list(...)$window, simtype, opera, TRUE) )
}
# Indices of similar days for cross-validation; TODO: 20 = magic number
cv_days = getSimilarDaysIndices(today, data, limit=20, same_season=FALSE,
- days_in=fdays)
+ days_in=tdays, operational=opera)
# Optimize h : h |--> sum of prediction errors on last N "similar" days
errorOnLastNdays = function(window, simtype)
for (i in seq_along(cv_days))
{
# mix_strategy is never used here (simtype != "mix"), therefore left blank
- prediction = private$.predictShapeAux(data, fdays, cv_days[i], horizon, local,
- window, simtype, FALSE)
+ prediction = private$.predictShapeAux(data, tdays, cv_days[i], predict_from,
+ horizon, local, window, simtype, opera, FALSE)
if (!is.na(prediction[1]))
{
nb_jours = nb_jours + 1
error = error +
- mean((data$getSerie(cv_days[i]+1)[1:horizon] - prediction)^2)
+ mean((data$getSerie(cv_days[i])[predict_from:horizon] - prediction)^2)
}
}
return (error / nb_jours)
else #none: value doesn't matter
1
- return(private$.predictShapeAux(data, fdays, today, horizon, local,
- best_window, simtype, TRUE))
+ return( private$.predictShapeAux(data, tdays, today, predict_from, horizon, local,
+ best_window, simtype, opera, TRUE) )
}
),
private = list(
- # Precondition: "today" is full (no NAs)
- .predictShapeAux = function(data, fdays, today, horizon, local, window, simtype,
- final_call)
+ # Precondition: "yersteday until predict_from-1" is full (no NAs)
+ .predictShapeAux = function(data, tdays, today, predict_from, horizon, local, window,
+ simtype, opera, final_call)
{
- fdays_cut = fdays[ fdays < today ]
- if (length(fdays_cut) <= 1)
+ tdays_cut = tdays[ tdays != today ]
+ if (length(tdays_cut) == 0)
return (NA)
if (local)
{
- # TODO: 60 == magic number
- fdays = getSimilarDaysIndices(today, data, limit=60, same_season=TRUE,
- days_in=fdays_cut)
- if (length(fdays) <= 1)
- return (NA)
- # TODO: 10, 12 == magic numbers
- fdays = .getConstrainedNeighbs(today,data,fdays,min_neighbs=10,max_neighbs=12)
- if (length(fdays) == 1)
+ # limit=Inf to not censor any day (TODO: finite limit? 60?)
+ tdays = getSimilarDaysIndices(today, data, limit=Inf, same_season=TRUE,
+ days_in=tdays_cut, operational=opera)
+ # TODO: 10 == magic number
+ tdays = .getConstrainedNeighbs(today, data, tdays, min_neighbs=10)
+ if (length(tdays) == 1)
{
if (final_call)
{
private$.params$weights <- 1
- private$.params$indices <- fdays
+ private$.params$indices <- tdays
private$.params$window <- 1
}
- return ( data$getSerie(fdays[1])[1:horizon] )
+ return ( data$getSerie(tdays[1])[predict_from:horizon] )
+ }
+ max_neighbs = 12 #TODO: 10 or 12 or... ?
+ if (length(tdays) > max_neighbs)
+ {
+ distances2 <- .computeDistsEndo(data, today, tdays, predict_from)
+ ordering <- order(distances2)
+ tdays <- tdays[ ordering[1:max_neighbs] ]
}
}
else
- fdays = fdays_cut #no conditioning
+ tdays = tdays_cut #no conditioning
if (simtype == "endo" || simtype == "mix")
{
# Compute endogen similarities using given window
window_endo = ifelse(simtype=="mix", window[1], window)
- # Distances from last observed day to days in the past
- serieToday = data$getSerie(today)
- distances2 = sapply(fdays, function(i) {
- delta = serieToday - data$getSerie(i)
- mean(delta^2)
- })
+ # Distances from last observed day to selected days in the past
+ # TODO: redundant computation if local==TRUE
+ distances2 <- .computeDistsEndo(data, today, tdays, predict_from)
simils_endo <- .computeSimils(distances2, window_endo)
}
# Compute exogen similarities using given window
window_exo = ifelse(simtype=="mix", window[2], window)
- M = matrix( nrow=1+length(fdays), ncol=1+length(data$getExo(today)) )
- M[1,] = c( data$getLevel(today), as.double(data$getExo(today)) )
- for (i in seq_along(fdays))
- M[i+1,] = c( data$getLevel(fdays[i]), as.double(data$getExo(fdays[i])) )
-
- sigma = cov(M) #NOTE: robust covariance is way too slow
- # TODO: 10 == magic number; more robust way == det, or always ginv()
- sigma_inv =
- if (length(fdays) > 10)
- solve(sigma)
- else
- MASS::ginv(sigma)
-
- # Distances from last observed day to days in the past
- distances2 = sapply(seq_along(fdays), function(i) {
- delta = M[1,] - M[i+1,]
- delta %*% sigma_inv %*% delta
- })
+ distances2 <- .computeDistsExo(data, today, tdays)
simils_exo <- .computeSimils(distances2, window_exo)
}
else if (simtype == "mix")
simils_endo * simils_exo
else #none
- rep(1, length(fdays))
+ rep(1, length(tdays))
similarities = similarities / sum(similarities)
- prediction = rep(0, horizon)
- for (i in seq_along(fdays))
- prediction = prediction + similarities[i] * data$getSerie(fdays[i]+1)[1:horizon]
+ prediction = rep(0, horizon-predict_from+1)
+ for (i in seq_along(tdays))
+ {
+ prediction = prediction +
+ similarities[i] * data$getSerie(tdays[i])[predict_from:horizon]
+ }
if (final_call)
{
private$.params$weights <- similarities
- private$.params$indices <- fdays
+ private$.params$indices <- tdays
private$.params$window <-
if (simtype=="endo")
window_endo
#
# @param today Index of current day
# @param data Object of class Data
-# @param fdays Current set of "first days" (no-NA pairs)
+# @param tdays Current set of "second days" (no-NA pairs)
# @param min_neighbs Minimum number of points in a neighborhood
# @param max_neighbs Maximum number of points in a neighborhood
#
-.getConstrainedNeighbs = function(today, data, fdays, min_neighbs=10, max_neighbs=12)
+.getConstrainedNeighbs = function(today, data, tdays, min_neighbs=10)
{
- levelToday = data$getLevel(today)
- distances = sapply(fdays, function(i) abs(data$getLevel(i)-levelToday))
- #TODO: 2, +3 : magic numbers
- dist_thresh = 2
- min_neighbs = min(min_neighbs,length(fdays))
+ levelToday = data$getLevelHat(today)
+ distances = sapply( tdays, function(i) abs(data$getLevel(i) - levelToday) )
+ #TODO: 1, +1, +3 : magic numbers
+ dist_thresh = 1
+ min_neighbs = min(min_neighbs,length(tdays))
repeat
{
same_pollution = (distances <= dist_thresh)
nb_neighbs = sum(same_pollution)
if (nb_neighbs >= min_neighbs) #will eventually happen
break
- dist_thresh = dist_thresh + 3
+ dist_thresh = dist_thresh + ifelse(dist_thresh>1,3,1)
}
- fdays = fdays[same_pollution]
- max_neighbs = 12
- if (nb_neighbs > max_neighbs)
- {
- # Keep only max_neighbs closest neighbors
- fdays = fdays[
- sort(distances[same_pollution],index.return=TRUE)$ix[1:max_neighbs] ]
- }
- fdays
+ tdays[same_pollution]
}
# compute similarities
}
exp(-distances2/(sd_dist*window^2))
}
+
+.computeDistsEndo <- function(data, today, tdays, predict_from)
+{
+ lastSerie = c( data$getSerie(today-1),
+ data$getSerie(today)[if (predict_from>=2) 1:(predict_from-1) else c()] )
+ sapply(tdays, function(i) {
+ delta = lastSerie - c(data$getSerie(i-1),
+ data$getSerie(i)[if (predict_from>=2) 1:(predict_from-1) else c()])
+ sqrt(mean(delta^2))
+ })
+}
+
+.computeDistsExo <- function(data, today, tdays)
+{
+ M = matrix( ncol=1+length(tdays), nrow=1+length(data$getExo(1)) )
+ M[,1] = c( data$getLevelHat(today), as.double(data$getExoHat(today)) )
+ for (i in seq_along(tdays))
+ M[,i+1] = c( data$getLevel(tdays[i]), as.double(data$getExo(tdays[i])) )
+
+ sigma = cov(t(M)) #NOTE: robust covariance is way too slow
+ # TODO: 10 == magic number; more robust way == det, or always ginv()
+ sigma_inv =
+ if (length(tdays) > 10)
+ solve(sigma)
+ else
+ MASS::ginv(sigma)
+
+ # Distances from last observed day to days in the past
+ sapply(seq_along(tdays), function(i) {
+ delta = M[,1] - M[,i+1]
+ delta %*% sigma_inv %*% delta
+ })
+}
inherit = Forecaster,
public = list(
- predictShape = function(data, today, memory, horizon, ...)
+ predictShape = function(data, today, memory, predict_from, horizon, ...)
{
# Return centered last (similar) day curve, avoiding NAs until memory is run
first_day = max(1, today-memory)
same_day = ifelse(hasArg("same_day"), list(...)$same_day, TRUE)
- # If 'same_day', get the last known future of similar day: -7 + 1 == -6
- index = today - ifelse(same_day,6,0)
+ index <- today
repeat
{
- {
- last_serie = data$getCenteredSerie(index)[1:horizon]
- index = index - ifelse(same_day,7,1)
- };
- if (!any(is.na(last_serie)))
- return (last_serie);
+ # If 'same_day', get the last known future of similar day
+ index = index - ifelse(same_day,7,1)
if (index < first_day)
return (NA)
+ last_serie = data$getCenteredSerie(index)[predict_from:horizon]
+ if (!any(is.na(last_serie)))
+ return (last_serie)
}
}
)
#' Zero Forecaster
#'
-#' Flat prediction: always forecast a serie of zeros.
-#' This serie is then adjusted using the ".pjump" function (see \code{Forecaster} class).
+#' Flat prediction: always forecast a serie of zeros (adjusted by 'pjump' function).
#'
#' @usage # ZeroForecaster$new(pjump)
#'
inherit = Forecaster,
public = list(
- predictShape = function(data, today, memory, horizon, ...)
- rep(0, horizon)
+ predictShape = function(data, today, memory, predict_from, horizon, ...)
+ rep(0, (horizon-predict_from+1))
)
)
#' Forecast encapsulation as a list (days where prediction occur) of lists (components).
#'
#' The private field .pred is a list where each cell contains the predicted variables for
-#' a period of time of H<=24 hours, from hour P+1 until P+H, where P+1 is taken right
-#' after the end of the period designated by \code{getIndexInData()}. In other terms,
-#' \code{forecast$getForecast(i)} return forecasts for \code{data$getSerie(i+1)}.
-#' Each cell .pred[[i]] is itself a list containing three slots, as described in the
-#' 'field' section.
+#' a period of time of H-P+1<=24 hours, from hour P until H, where P == predict_from.
+#' \code{forecast$getForecast(i)} output forecasts for
+#' \code{data$getSerie(forecast$getIndexInData(i))}.
#'
#' @usage # Forecast$new(dates)
#'
#' example "Neighbors" method stores informations about the considered neighborhood for
#' the current prediction task) and one main function: \code{predictSerie()}. This last
#' function (by default) calls \code{predictShape()} to get a forecast of a centered
-#' serie, and then calls the "jump prediction" function -- see "field" section -- to
-#' adjust it based on the last observed values.
+#' serie, and then calls the "jump prediction" function if it's provided -- see "field"
+#' section -- to adjust it based on the last observed values. The main method in derived
+#' forecasters is \code{predictShape()}; see 'Methods' section.
#'
#' @usage # Forecaster$new(pjump) #warning: predictShape() is unimplemented
#'
#' @field .pjump Function: how to predict the jump at day interface? The arguments of
#' this function are -- in this order:
#' \itemize{
-#' \item data : object output of \code{getData()},
-#' \item today : index (integer or date) of the last known day in data,
-#' \item memory : number of days to use in the past (including today),
-#' \item horizon : number of time steps to predict,
-#' \item params : optimized parameters in the main method \code{predictShape()},
-#' \item ... : additional arguments.
+#' \item data: object output of \code{getData()},
+#' \item today: index of the current day in data (known until predict_from-1),
+#' \item memory: number of days to use in the past (including today),
+#' \item predict_from: first time step to predict (in [1,24])
+#' \item horizon: last time step to predict (in [predict_from,24]),
+#' \item params: optimized parameters in the main method \code{predictShape()},
+#' \item ...: additional arguments.
#' }
#' .pjump returns an estimation of the jump after the last observed value.
#'
#' @section Methods:
#' \describe{
#' \item{\code{initialize(data, pjump)}}{
-#' Initialize a Forecaster object with a Data object and a jump prediction function.}
-#' \item{\code{predictSerie(today,memory,horizon,...)}}{
-#' Predict a new serie of \code{horizon} values at day index \code{today} using
-#' \code{memory} days in the past.}
-#' \item{\code{predictShape(today,memory,horizon,...)}}{
-#' Predict a new shape of \code{horizon} values at day index \code{today} using
+#' Initialize a Forecaster object with a Data object and a jump prediction function,
+#' or NULL if \code{predictShape()} returns an adjusted curve.}
+#' \item{\code{predictSerie(data,today,memory,predict_from,horizon,...)}}{
+#' Predict the next curve (at index today) from predict_from to horizon (hours), using
#' \code{memory} days in the past.}
+#' \item{\code{predictShape(data,today,memory,predict_from,horizon,...)}}{
+#' Predict the shape of the next curve (at index today) from predict_from to horizon
+#' (hours), using \code{memory} days in the past.}
#' \item{\code{getParameters()}}{
#' Return (internal) parameters.}
#' }
private$.pjump <- pjump
invisible(self)
},
- predictSerie = function(data, today, memory, horizon, ...)
+ predictSerie = function(data, today, memory, predict_from, horizon, ...)
{
# Parameters (potentially) computed during shape prediction stage
- predicted_shape = self$predictShape(data, today, memory, horizon, ...)
- predicted_delta = private$.pjump(data,today,memory,horizon,private$.params,...)
- # Predicted shape is aligned it on the end of current day + jump
- predicted_shape+tail(data$getSerie(today),1)-predicted_shape[1]+predicted_delta
+ predicted_shape <- self$predictShape(data,today,memory,predict_from,horizon,...)
+
+ if (is.na(predicted_shape))
+ return (NA)
+
+ predicted_delta <-
+ if (is.null(private$.pjump))
+ NULL
+ else
+ private$.pjump(data,today,memory,predict_from,horizon,private$.params,...)
+
+ # Predicted shape is aligned on the end of current day + jump (if jump!=NULL)
+ c( data$getSerie(today)[if (predict_from>=2) 1:(predict_from-1) else c()],
+ predicted_shape + ifelse( is.null(private$.pjump),
+ 0,
+ predicted_delta - predicted_shape[1] +
+ ifelse(predict_from>=2,
+ data$getSerie(today)[predict_from-1], tail(data$getSerie(today-1),1)) ) )
},
- predictShape = function(data, today, memory, horizon, ...)
+ predictShape = function(data, today, memory, predict_from, horizon, ...)
NULL #empty default implementation: to implement in inherited classes
,
getParameters = function()
#'
#' @aliases J_Neighbors
#'
-getNeighborsJumpPredict = function(data, today, memory, horizon, params, ...)
+getNeighborsJumpPredict = function(data, today, memory, predict_from, horizon,
+ params, ...)
{
first_day = max(1, today-memory)
filter = (params$indices >= first_day)
indices = params$indices[filter]
weights = params$weights[filter]
+ if (is.na(indices[1]))
+ return (NA)
+
gaps = sapply(indices, function(i) {
- head( data$getSerie(i+1),1 ) - tail( data$getSerie(i),1 )
+ if (predict_from >= 2)
+ data$getSerie(i)[predict_from] - data$getSerie(i)[predict_from-1]
+ else
+ head(data$getSerie(i),1) - tail(data$getSerie(i-1),1)
})
scal_product = weights * gaps
norm_fact = sum( weights[!is.na(scal_product)] )
#'
#' @aliases J_Persistence
#'
-getPersistenceJumpPredict = function(data, today, memory, horizon, params, ...)
+getPersistenceJumpPredict = function(data, today, memory, predict_from,
+ horizon, params, ...)
{
#return gap between end of similar day curve and first day of tomorrow (in the past)
first_day = max(1, today-memory)
same_day = ifelse(hasArg("same_day"), list(...)$same_day, TRUE)
- index = today - ifelse(same_day,7,1)
+ index <- today
repeat
{
- {
- last_serie_end = tail( data$getSerie(index), 1)
- last_tomorrow_begin = head( data$getSerie(index+1), 1)
- index = index - ifelse(same_day,7,1)
- };
- if (!is.na(last_serie_end) && !is.na(last_tomorrow_begin))
- return (last_tomorrow_begin - last_serie_end);
+ # If 'same_day', get the last known future of similar day
+ index = index - ifelse(same_day,7,1)
if (index < first_day)
return (NA)
+ gap <-
+ if (predict_from >= 2)
+ data$getSerie(index)[predict_from] - data$getSerie(index)[predict_from-1]
+ else
+ head(data$getSerie(index),1) - tail(data$getSerie(index-1),1)
+ if (!is.na(gap))
+ return (gap)
}
}
#'
#' @param data Object of class \code{Data} output of \code{getData}
#' @param pred Object of class \code{Forecast} output of \code{computeForecast}
-#' @param horizon Horizon where to compute the error
-#' (<= horizon used in \code{computeForecast})
+#' @param predict_from First time step to consider (>= predict_from used in
+#' \code{computeForecast()})
+#' @param horizon Horizon where to compute the error (<= horizon used in
+#' \code{computeForecast})
#'
#' @return A list (abs,MAPE) of lists (day,indices). The "indices" slots contain series
#' of size L where L is the number of predicted days; i-th value is the averaged error
#' step, averaged on the L forecasting days.
#'
#' @export
-computeError = function(data, pred, horizon=data$getStdHorizon())
+computeError = function(data, pred, predict_from, horizon=length(data$getSerie(1)))
{
L = pred$getSize()
- mape_day = rep(0, horizon)
- abs_day = rep(0, horizon)
+ mape_day = rep(0, horizon-predict_from+1)
+ abs_day = rep(0, horizon-predict_from+1)
mape_indices = rep(NA, L)
abs_indices = rep(NA, L)
for (i in seq_len(L))
{
index = pred$getIndexInData(i)
- serie = data$getSerie(index+1)[1:horizon]
- forecast = pred$getForecast(i)[1:horizon]
+ serie = data$getSerie(index)[predict_from:horizon]
+ forecast = pred$getForecast(i)[predict_from:horizon]
if (!any(is.na(serie)) && !any(is.na(forecast)))
{
nb_no_NA_data = nb_no_NA_data + 1
#' Compute forecast
#'
-#' Predict time-series curves ("tomorrows") at the selected days indices ("todays").
-#' This function just runs a loop over all requested indices, and stores the individual
-#' forecasts into a list which is then turned into a Forecast object.
+#' Predict time-series curves ("today" from predict_from to horizon) at the selected days
+#' indices ("today" from 1am to predict_from-1). This function just runs a loop over all
+#' requested indices, and stores the individual forecasts into a Forecast object.
#'
#' @param data Object of class Data, output of \code{getData()}.
#' @param indices Indices where to forecast (the day after); integers relative to the
#' beginning of data, or (convertible to) Date objects.
#' @param forecaster Name of the main forecaster; more details: ?F_<forecastername>
-#' \itemize{
-#' \item Persistence : use last (similar, next) day
-#' \item Neighbors : weighted tomorrows of similar days
-#' \item Average : average tomorrow of all same day-in-week
-#' \item Zero : just output 0 (benchmarking purpose)
-#' }
+#' \itemize{
+#' \item Persistence : use last (similar) day
+#' \item Neighbors : weighted similar days
+#' \item Average : average curve of all same day-in-week
+#' \item Zero : just output 0 (benchmarking purpose)
+#' }
#' @param pjump Function to predict the jump at the interface between two days;
#' more details: ?J_<functionname>
-#' \itemize{
-#' \item Persistence : use last (similar, next) day
-#' \item Neighbors: re-use the weights from F_Neighbors
-#' \item Zero: just output 0 (no adjustment)
-#' }
+#' \itemize{
+#' \item Persistence : use last (similar) day
+#' \item Neighbors: re-use the weights from F_Neighbors
+#' \item Zero: just output 0 (no adjustment)
+#' }
+#' If pjump=NULL, then no adjustment is performed (output of \code{predictShape()} is
+#' used directly).
+#' @param predict_from First time step to predict.
#' @param memory Data depth (in days) to be used for prediction.
-#' @param horizon Number of time steps to predict.
+#' @param horizon Last time step to predict.
#' @param ncores Number of cores for parallel execution (1 to disable).
+#' @param verbose TRUE to print basic traces (runs beginnings)
#' @param ... Additional parameters for the forecasting models.
#'
#' @return An object of class Forecast
#' @examples
#' ts_data <- system.file("extdata","pm10_mesures_H_loc.csv",package="talweg")
#' exo_data <- system.file("extdata","meteo_extra_noNAs.csv",package="talweg")
-#' data <- getData(ts_data, exo_data, input_tz="GMT", working_tz="GMT",
-#' predict_at=7, limit=200)
+#' data <- getData(ts_data, exo_data, limit=200)
#' pred <- computeForecast(data, 100:130, "Persistence", "Zero",
-#' memory=50, horizon=12, ncores=1)
-#' \dontrun{#Sketch for real-time mode:
+#' predict_from=8, memory=50, horizon=12, ncores=1)
+#' \dontrun{
+#' #Sketch for real-time mode:
#' data <- Data$new()
#' forecaster <- MyForecaster$new(myJumpPredictFunc)
#' repeat {
-#' # In the morning 7am+ or afternoon 1pm+:
+#' # As soon as daily predictions are available:
#' data$append(
-#' times_from_H+1_yersteday_to_Hnow,
-#' PM10_values_of_last_24h,
-#' exogenous_measures_of_last_24h,
-#' exogenous_predictions_for_next_24h)
+#' level_hat=predicted_level,
+#' exo_hat=predicted_exogenous)
+#' # When a day ends:
+#' data$append(
+#' level=observed_level,
+#' exo=observed_exogenous)
+#' # And, at every hour:
+#' data$append(
+#' time=current_hour,
+#' value=current_PM10)
+#' # Finally, a bit before predict_from hour:
#' pred <- forecaster$predictSerie(data, data$getSize(), ...)
#' #do_something_with_pred
-#' }}
+#' } }
#' @export
-computeForecast = function(data, indices, forecaster, pjump,
- memory=Inf, horizon=data$getStdHorizon(), ncores=3, ...)
+computeForecast = function(data, indices, forecaster, pjump, predict_from,
+ memory=Inf, horizon=length(data$getSerie(1)), ncores=3, verbose=FALSE, ...)
{
# (basic) Arguments sanity checks
+ predict_from = as.integer(predict_from)[1]
+ if (! predict_from %in% 1:length(data$getSerie(1)))
+ stop("predict_from in [1,24] (hours)")
+ if (hasArg("opera") && !list(...)$opera && memory < Inf)
+ memory <- Inf #finite memory in training mode makes no sense
horizon = as.integer(horizon)[1]
- if (horizon<=0 || horizon>length(data$getCenteredSerie(1)))
- stop("Horizon too short or too long")
+ if (horizon<=predict_from || horizon>length(data$getSerie(1)))
+ stop("Horizon in [predict_from+1,24] (hours)")
integer_indices = sapply(indices, function(i) dateIndexToInteger(i,data))
if (any(integer_indices<=0 | integer_indices>data$getSize()))
stop("Indices out of range")
- if (!is.character(forecaster) || !is.character(pjump))
- stop("forecaster (name) and pjump (function) should be of class character")
+ if (!is.character(forecaster))
+ stop("forecaster (name): character")
+ if (!is.null(pjump) && !is.character(pjump))
+ stop("pjump (function): character or NULL")
pred = Forecast$new( sapply(indices, function(i) integerIndexToDate(i,data)) )
forecaster_class_name = getFromNamespace(
paste(forecaster,"Forecaster",sep=""), "talweg")
- forecaster = forecaster_class_name$new( #.pjump =
- getFromNamespace(paste("get",pjump,"JumpPredict",sep=""), "talweg"))
- if (ncores > 1 && requireNamespace("parallel",quietly=TRUE))
- {
- p <- parallel::mclapply(seq_along(integer_indices), function(i) {
- list(
- "forecast" = forecaster$predictSerie(
- data, integer_indices[i], memory, horizon, ...),
- "params"= forecaster$getParameters(),
- "index" = integer_indices[i] )
- }, mc.cores=ncores)
- }
- else
+ if (!is.null(pjump))
+ pjump <- getFromNamespace(paste("get",pjump,"JumpPredict",sep=""), "talweg")
+ forecaster = forecaster_class_name$new(pjump)
+
+ computeOneForecast <- function(i)
{
- p <- lapply(seq_along(integer_indices), function(i) {
- list(
- "forecast" = forecaster$predictSerie(
- data, integer_indices[i], memory, horizon, ...),
- "params"= forecaster$getParameters(),
- "index" = integer_indices[i] )
- })
+ if (verbose)
+ print(paste("Index",i))
+ list(
+ "forecast" = forecaster$predictSerie(data,i,memory,predict_from,horizon,...),
+ "params" = forecaster$getParameters(),
+ "index" = i )
}
+ p <-
+ if (ncores > 1 && requireNamespace("parallel",quietly=TRUE))
+ parallel::mclapply(integer_indices, computeOneForecast, mc.cores=ncores)
+ else
+ lapply(integer_indices, computeOneForecast)
+
# TODO: find a way to fill pred in //...
for (i in seq_along(integer_indices))
{
#' @param exo_data Exogenous variables, as a data frame or a CSV file; first column is
#' dates, next block are measurements for the day, and final block are exogenous
#' forecasts (for the same day).
-#' @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 ?strptime)
-#' @param working_tz Timezone to work with ("GMT" or e.g. "Europe/Paris")
-#' @param predict_at When does the prediction take place? Integer, in hours. Default: 0
#' @param limit Number of days to extract (default: Inf, for "all")
#'
#' @return An object of class Data
#' exo_data = read.csv(system.file("extdata","meteo_extra_noNAs.csv",package="talweg"))
#' data = getData(ts_data, exo_data, limit=120)
#' @export
-getData = function(ts_data, exo_data, input_tz="GMT", date_format="%d/%m/%Y %H:%M",
- working_tz="GMT", predict_at=0, limit=Inf)
+getData = function(ts_data, exo_data, date_format="%d/%m/%Y %H:%M", limit=Inf)
{
# 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)) ||
(!is.data.frame(exo_data) && !is.character(exo_data)) )
stop("Bad time-series / exogenous input (data frame or CSV file)")
ts_data = ts_data[1]
if (is.character(exo_data))
exo_data = exo_data[1]
- predict_at = as.integer(predict_at)[1]
- if (predict_at<0 || predict_at>23)
- stop("Bad predict_at (0-23)")
if (!is.character(date_format))
stop("Bad date_format (character)")
date_format = date_format[1]
+ if (!is.numeric(limit) || limit < 0)
+ stop("limit: positive integer")
ts_df =
if (is.character(ts_data))
read.csv(ts_data)
else
ts_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=input_tz), tz=working_tz, usetz=TRUE)
+ # Convert to GMT (pretend it's GMT; no impact)
+ dates_POSIXlt = strptime(as.character(ts_df[,1]), date_format, tz="GMT")
+ ts_df[,1] = format(as.POSIXct(dates_POSIXlt, tz="GMT"), tz="GMT", usetz=TRUE)
exo_df =
if (is.character(exo_data))
line = line + 1
};
if (line >= nb_lines + 1
- || as.POSIXlt(ts_df[line-1,1],tz=working_tz)$hour == predict_at)
+ || as.POSIXlt(ts_df[line-1,1],tz="GMT")$hour == 0)
{
break
}
}
- exo = as.data.frame( exo_df[i,2:(1+nb_exos)] )
- exo_hat =
- if (i < nrow(exo_df))
- as.data.frame( exo_df[i+1,(1+nb_exos+1):(1+2*nb_exos)] )
- else
- NA #exogenous prediction for next day are useless on last day
- data$append(time, serie, exo, exo_hat)
+ # TODO: 2 modes, "operational" and "testing"; would need PM10 predictions
+ data$append(time=time, value=serie, level_hat=mean(serie,na.rm=TRUE),
+ exo=exo_df[i,2:(1+nb_exos)], exo_hat=exo_df[i,(1+nb_exos+1):(1+2*nb_exos)])
if (i >= limit)
break
i = i + 1
}
- if (length(data$getCenteredSerie(1)) < length(data$getCenteredSerie(2)))
- data$removeFirst()
- if (length(data$getCenteredSerie(data$getSize()))
- < length(data$getCenteredSerie(data$getSize()-1)))
- {
- data$removeLast()
- }
data
}
series = data$getSeries(indices)
yrange = quantile(series, probs=c(0.025,0.975), na.rm=TRUE)
par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5)
- for (i in seq_along(indices))
- {
- plot(series[,i], type="l", ylim=yrange,
- xlab=ifelse(i==1,"Time (hours)",""), ylab=ifelse(i==1,"PM10",""))
- if (i < length(indices))
- par(new=TRUE)
- }
+ matplot(series, type="l", ylim=yrange, xlab="Time (hours)", ylab="PM10")
}
#' Plot error
#' \code{\link{plotFilamentsBox}}, \code{\link{plotRelVar}}
#'
#' @export
-plotError <- function(err, cols=seq_along(err))
+plotError <- function(err, cols=seq_along(err), agg="day")
{
if (!is.null(err$abs))
err = list(err)
- par(mfrow=c(2,2), mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5, lwd=2)
+ par(mfrow=c(2,2), mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5)
L = length(err)
- yrange = range( sapply(1:L, function(i) ( err[[i]]$abs$day ) ), na.rm=TRUE )
- for (i in seq_len(L))
- {
- plot(err[[i]]$abs$day, type="l", xlab=ifelse(i==1,"Time (hours)",""),
- ylab=ifelse(i==1,"Mean |y - y_hat|",""), ylim=yrange, col=cols[i])
- if (i < L)
- par(new=TRUE)
- }
- yrange = range( sapply(1:L, function(i) ( err[[i]]$abs$indices ) ), na.rm=TRUE )
- for (i in seq_len(L))
- {
- plot(err[[i]]$abs$indices, type="l", xlab=ifelse(i==1,"Time (days)",""),
- ylab=ifelse(i==1,"Mean |y - y_hat|",""), ylim=yrange, col=cols[i])
- if (i < L)
- par(new=TRUE)
- }
- yrange = range( sapply(1:L, function(i) ( err[[i]]$MAPE$day ) ), na.rm=TRUE )
- for (i in seq_len(L))
- {
- plot(err[[i]]$MAPE$day, type="l", xlab=ifelse(i==1,"Time (hours)",""),
- ylab=ifelse(i==1,"Mean MAPE",""), ylim=yrange, col=cols[i])
- if (i < L)
- par(new=TRUE)
- }
- yrange = range( sapply(1:L, function(i) ( err[[i]]$MAPE$indices ) ), na.rm=TRUE )
- for (i in seq_len(L))
- {
- plot(err[[i]]$MAPE$indices, type="l", xlab=ifelse(i==1,"Time (days)",""),
- ylab=ifelse(i==1,"Mean MAPE",""), ylim=yrange, col=cols[i])
- if (i < L)
- par(new=TRUE)
- }
+
+ yrange = range( sapply(1:L, function(i) err[[i]]$abs$day), na.rm=TRUE )
+ matplot(sapply( seq_len(L), function(i) err[[i]]$abs$day ), type="l",
+ xlab="Time (hours)", ylab="Mean |y - y_hat|", ylim=yrange, col=cols, lwd=2, lty=1)
+
+ agg_curves <- sapply( seq_len(L), function(i) {
+ curve <- err[[i]]$abs$indices
+ delta <- if (agg=="day") 1 else if (agg=="week") 7 else if (agg=="month") 30
+ vapply( seq(1,length(curve),delta), function(i) {
+ mean(curve[i:(i+delta-1)], na.rm=TRUE)
+ }, vector("double",1), USE.NAMES=FALSE )
+ })
+ yrange = range(agg_curves, na.rm=TRUE)
+ matplot(agg_curves, type="l", xlab=paste("Time (",agg,"s)", sep=""),
+ ylab="Mean |y - y_hat|", ylim=yrange, col=cols, lwd=2, lty=1)
+
+ yrange = range( sapply(1:L, function(i) err[[i]]$MAPE$day), na.rm=TRUE )
+ matplot(sapply( seq_len(L), function(i) err[[i]]$MAPE$day ), type="l",
+ xlab="Time (hours)", ylab="Mean MAPE", ylim=yrange, col=cols, lwd=2, lty=1)
+
+ agg_curves <- sapply( seq_len(L), function(i) {
+ curve <- err[[i]]$MAPE$indices
+ delta <- if (agg=="day") 1 else if (agg=="week") 7 else if (agg=="month") 30
+ vapply( seq(1,length(curve),delta), function(i) {
+ mean(curve[i:(i+delta-1)], na.rm=TRUE)
+ }, vector("double",1), USE.NAMES=FALSE )
+ })
+ yrange = range(agg_curves, na.rm=TRUE)
+ matplot(agg_curves, type="l", xlab=paste("Time (",agg,"s)", sep=""),
+ ylab="Mean MAPE", ylim=yrange, col=cols, lwd=2, lty=1)
}
#' Plot measured / predicted
#' @export
plotPredReal <- function(data, pred, index)
{
- horizon = length(pred$getForecast(1))
- measure = data$getSerie( pred$getIndexInData(index)+1 )[1:horizon]
prediction = pred$getForecast(index)
+ measure = data$getSerie( pred$getIndexInData(index) )[1:length(pred$getForecast(1))]
+
+ # Remove the common part, where prediction == measure
+ dot_mark <- ifelse(prediction[1]==measure[1],
+ which.max(seq_along(prediction)[prediction==measure]), 0)
+ prediction = prediction[(dot_mark+1):length(prediction)]
+ measure = measure[(dot_mark+1):length(measure)]
+
yrange = range(measure, prediction)
par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5, lwd=3)
plot(measure, type="l", ylim=yrange, xlab="Time (hours)", ylab="PM10")
weights = pred$getParams(index)$weights
if (is.null(weights))
stop("plotSimils only works on 'Neighbors' forecasts")
- par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5)
- hist(pred$getParams(index)$weights, nclass=20, main="", xlab="Weight", ylab="Count")
+ par(mfrow=c(1,2), mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5)
+ small_weights = weights[ weights < 1/length(weights) ]
+ large_weights = weights[ weights >= 1/length(weights) ]
+ hist(small_weights, nclass=25, main="", xlab="Weight < 1/N", ylab="Count")
+ hist(large_weights, nclass=25, main="", xlab="Weight >= 1/N", ylab="Count")
}
#' Functional boxplot
#' @export
computeFilaments <- function(data, pred, index, limit=60, plot=TRUE)
{
- ref_serie = data$getCenteredSerie( pred$getIndexInData(index) )
- if (any(is.na(ref_serie)))
+ weights <- pred$getParams(index)$weights
+ if (is.null(weights) || is.na(pred$getParams(index)$weights[1]))
stop("computeFilaments requires a serie without NAs")
- # Compute colors for each neighbor (from darkest to lightest)
- sorted_dists = sort(-log(pred$getParams(index)$weights), index.return=TRUE)
- nn = min(limit, length(sorted_dists$x))
- min_dist = min(sorted_dists$x[1:nn])
- max_dist = max(sorted_dists$x[1:nn])
- color_values = floor(19.5*(sorted_dists$x[1:nn]-min_dist)/(max_dist-min_dist)) + 1
- colors = gray.colors(20,0.1,0.9)[color_values] #TODO: 20 == magic number
+ nn <- min(limit, length(weights))
+ sorted_dists = sort(-log(weights), index.return=TRUE)
+ # Compute colors for each neighbor (from darkest to lightest), if weights differ
+ if ( any( weights != weights[1] ) )
+ {
+ min_dist = min(sorted_dists$x[1:nn])
+ max_dist = max(sorted_dists$x[1:nn])
+ color_values = floor(19.5*(sorted_dists$x[1:nn]-min_dist)/(max_dist-min_dist)) + 1
+ colors = gray.colors(20,0.1,0.9)[color_values] #TODO: 20 == magic number
+ }
+ else
+ colors <- rep(colors()[17], length(weights))
if (plot)
{
# Complete series with (past and present) tomorrows
- ref_serie = c(ref_serie, data$getCenteredSerie( pred$getIndexInData(index)+1 ))
+ ref_serie = c( data$getCenteredSerie( pred$getIndexInData(index)-1 ),
+ data$getCenteredSerie( pred$getIndexInData(index) ) )
centered_series = rbind(
- data$getCenteredSeries( pred$getParams(index)$indices ),
- data$getCenteredSeries( pred$getParams(index)$indices+1 ) )
+ data$getCenteredSeries( pred$getParams(index)$indices-1 ),
+ data$getCenteredSeries( pred$getParams(index)$indices ) )
yrange = range( ref_serie,
quantile(centered_series, probs=c(0.025,0.975), na.rm=TRUE) )
par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5, lwd=2)
}
# Also plot ref curve, in red
plot(ref_serie, ylim=yrange, type="l", col="#FF0000", xlab="", ylab="")
- abline(v=24, lty=2, col=colors()[56], lwd=1)
+ dot_mark <- 0.5 + which.max( pred$getForecast(1) ==
+ data$getSerie( pred$getIndexInData(1) )[1:length(pred$getForecast(1))] )
+ abline(v=24+dot_mark, lty=2, col=colors()[56], lwd=1)
}
list(
#'
#' @inheritParams computeError
#' @param fil Output of \code{computeFilaments}
+#' @param predict_from First predicted time step
#'
#' @export
-plotFilamentsBox = function(data, fil)
+plotFilamentsBox = function(data, fil, predict_from)
{
if (!requireNamespace("rainbow", quietly=TRUE))
stop("Functional boxplot requires the rainbow package")
series_matrix = rbind(
- data$getSeries(fil$neighb_indices), data$getSeries(fil$neighb_indices+1) )
+ data$getSeries(fil$neighb_indices-1), data$getSeries(fil$neighb_indices) )
series_fds = rainbow::fds(seq_len(nrow(series_matrix)), series_matrix)
+
par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5)
rainbow::fboxplot(series_fds, "functional", "hdr", xlab="Time (hours)", ylab="PM10",
plotlegend=FALSE, lwd=2)
usr <- par("usr")
yr <- (usr[4] - usr[3]) / 27
par(new=TRUE)
- plot(c(data$getSerie(fil$index),data$getSerie(fil$index+1)), type="l", lwd=2, lty=2,
+ plot(c(data$getSerie(fil$index-1),data$getSerie(fil$index)), type="l", lwd=2, lty=2,
ylim=c(usr[3] + yr, usr[4] - yr), xlab="", ylab="")
- abline(v=24, lty=2, col=colors()[56])
+ abline(v=24+predict_from-0.5, lty=2, col=colors()[56])
}
#' Plot relative conditional variability / absolute variability
#' @inheritParams plotFilamentsBox
#'
#' @export
-plotRelVar = function(data, fil)
+plotRelVar = function(data, fil, predict_from)
{
- ref_var = c( apply(data$getSeries(fil$neighb_indices),1,sd),
- apply(data$getSeries(fil$neighb_indices+1),1,sd) )
- fdays = .getNoNA2(data, 1, fil$index-1)
+ ref_var = c( apply(data$getSeries(fil$neighb_indices-1),1,sd),
+ apply(data$getSeries(fil$neighb_indices),1,sd) )
+ tdays = .getNoNA2(data, 2, fil$index)
global_var = c(
- apply(data$getSeries(fdays),1,sd),
- apply(data$getSeries(fdays+1),1,sd) )
+ apply(data$getSeries(tdays-1),1,sd),
+ apply(data$getSeries(tdays),1,sd) )
yrange = range(ref_var, global_var)
par(mar=c(4.7,5,1,1), cex.axis=1.5, cex.lab=1.5)
xlab="Time (hours)", ylab="Standard deviation")
par(new=TRUE)
plot(global_var, type="l", col=2, lwd=3, ylim=yrange, xlab="", ylab="")
- abline(v=24, lty=2, col=colors()[56])
+ abline(v=24+predict_from-0.5, lty=2, col=colors()[56])
}
#' @param days_in Optional set to intersect with results (NULL to discard)
#'
#' @export
-getSimilarDaysIndices = function(index, data, limit, same_season, days_in=NULL)
+getSimilarDaysIndices = function(index, data, limit, same_season,
+ days_in=NULL, operational=TRUE)
{
index = dateIndexToInteger(index, data)
day_ref = dt_ref$wday #1=monday, ..., 6=saturday, 0=sunday
month_ref = as.POSIXlt(data$getTime(index)[1])$mon+1 #month in 1...12
i = index - 1
- while (i >= 1 && length(days) < limit)
+ if (!operational)
+ j = index + 1
+ while (length(days) < min( limit, ifelse(is.null(days_in),Inf,length(days_in)) ))
{
- dt = as.POSIXlt(data$getTime(i)[1])
- if ((is.null(days_in) || i %in% days_in) && .isSameDay(dt$wday, day_ref))
+ if (i < 1 && j > data$getSize())
+ break
+ if (i >= 1)
{
- if (!same_season || .isSameSeason(dt$mon+1, month_ref))
- days = c(days, i)
+ dt = as.POSIXlt(data$getTime(i)[1])
+ if ((is.null(days_in) || i %in% days_in) && .isSameDay(dt$wday, day_ref))
+ {
+ if (!same_season || .isSameSeason(dt$mon+1, month_ref))
+ days = c(days, i)
+ }
+ i = i - 1
+ }
+ if (!operational && j <= data$getSize())
+ {
+ dt = as.POSIXlt(data$getTime(j)[1])
+ if ((is.null(days_in) || j %in% days_in) && .isSameDay(dt$wday, day_ref))
+ {
+ if (!same_season || .isSameSeason(dt$mon+1, month_ref))
+ days = c(days, j)
+ }
+ j = j + 1
}
- i = i - 1
}
return ( days )
}
#
.isSameSeason = function(month, month_ref)
{
+# if (month_ref == 3) #TODO: same as Bruno (but weird)
+# return (month %in% c(2,3,4,9,10))
if (month_ref %in% c(11,12,1,2)) #~= mid-polluted
return (month %in% c(11,12,1,2))
if (month_ref %in% c(3,4,9,10)) #~= high-polluted
# getNoNA2
#
-# Get indices in data of no-NA series followed by no-NA, within [first,last] range.
+# Get indices in data of no-NA series preceded by no-NA, within [first,last] range.
#
# @inheritParams dateIndexToInteger
# @param first First index (included)
.getNoNA2 = function(data, first, last)
{
(first:last)[ sapply(first:last, function(i)
- !any( is.na(data$getCenteredSerie(i)) | is.na(data$getCenteredSerie(i+1)) )
+ !any( is.na(data$getSerie(i-1)) | is.na(data$getSerie(i)) )
) ]
}
time = as.POSIXct((i-1)*60*60*24+15*60*(1:96), origin="2007-01-01", tz="GMT")
exo = runif(4)
exo_hat = runif(4)
- data$append(time, serie, exo, exo_hat)
+ data$append(time=time, value=serie, exo=exo, exo_hat=exo_hat)
}
data
}
ts_data = system.file("testdata","ts_test.csv",package="talweg")
exo_data = system.file("testdata","exo_test.csv",package="talweg")
-data0 <<- getData(ts_data, exo_data, input_tz="GMT", date_format="%Y-%m-%d %H:%M",
- working_tz="GMT", predict_at=0, limit=Inf)
-data7 <<- getData(ts_data, exo_data, input_tz="GMT", date_format="%Y-%m-%d %H:%M",
- working_tz="GMT", predict_at=7, limit=Inf)
+data0 <<- getData(ts_data, exo_data, date_format="%Y-%m-%d %H:%M", limit=Inf)
+data7 <<- getData(ts_data, exo_data, date_format="%Y-%m-%d %H:%M", limit=Inf)
test_that("dateIndexToInteger works as expected",
{
ts_data = system.file("testdata","ts_test.csv",package="talweg")
exo_data = system.file("testdata","exo_test.csv",package="talweg")
-data00 <<- getData(ts_data, exo_data, input_tz="GMT", date_format="%Y-%m-%d %H:%M",
- working_tz="GMT", predict_at=0, limit=Inf)
-data13 <<- getData(ts_data, exo_data, input_tz="GMT", date_format="%Y-%m-%d %H:%M",
- working_tz="GMT", predict_at=13, limit=Inf)
-#Forecast at sunday to saturday (series 7 to 1), for monday to sunday (series 1 to 7)
-indices <<- seq(as.Date("2007-04-01"),as.Date("2007-04-07"),"days")
+data_p <<- getData(ts_data, exo_data, date_format="%Y-%m-%d %H:%M", limit=Inf)
+#Forecasts from monday to sunday (series 1 to 7)
+indices <<- seq(as.Date("2007-04-02"),as.Date("2007-04-08"),"days")
pred_order = c(7,1:6) #will facilitate tests
test_that("Average method behave as expected",
{
- pred00_z = computeForecast(data00, indices, "Average", "Zero", Inf, 24)
- pred00_p = computeForecast(data00, indices, "Average", "Persistence", Inf, 24)
+ pred00_z = computeForecast(data_p, indices, "Average", "Zero", 1, Inf, 24, ncores=1)
+ pred00_p = computeForecast(data_p, indices, "Average", "Persistence", 1, Inf, 24)
for (i in 1:7)
{
#zero jump: should predict true values minus 1
expect_equal( pred00_p$getForecast(i), rep(i,24) )
}
- #NOTE: days become
+ #NOTE: 24h-block become
#1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 (14h-->0h then 1h-->13h)
#No jump between days, thus zero and persistence are equivalent (and correct)
- pred13_z = computeForecast(data13, indices, "Average", "Zero", Inf, 24)
- pred13_p = computeForecast(data13, indices, "Average", "Persistence", Inf, 24)
+ pred13_z = computeForecast(data_p, indices, "Average", "Zero", 14, Inf, 24)
+ pred13_p = computeForecast(data_p, indices, "Average", "Persistence", 14, Inf, 24)
for (i in 1:7)
{
- expect_equal( pred13_z$getForecast(i), c( rep(i,11), rep(i%%7+1,13) ) )
- expect_equal( pred13_p$getForecast(i), c( rep(i,11), rep(i%%7+1,13) ) )
+ expect_equal( pred13_z$getForecast(i), rep(i,24) )
+ expect_equal( pred13_p$getForecast(i), rep(i,24) )
}
#A few extra checks
- expect_equal( pred00_p$getIndexInData(2), dateIndexToInteger("2007-04-02",data00) )
- expect_equal( pred00_z$getIndexInData(2), dateIndexToInteger("2007-04-02",data00) )
- expect_equal( pred13_p$getIndexInData(5), dateIndexToInteger("2007-04-05",data13) )
- expect_equal( pred13_z$getIndexInData(5), dateIndexToInteger("2007-04-05",data13) )
+ expect_equal( pred00_p$getIndexInData(2), dateIndexToInteger("2007-04-03",data_p) )
+ expect_equal( pred00_z$getIndexInData(2), dateIndexToInteger("2007-04-03",data_p) )
+ expect_equal( pred13_p$getIndexInData(5), dateIndexToInteger("2007-04-06",data_p) )
+ expect_equal( pred13_z$getIndexInData(5), dateIndexToInteger("2007-04-06",data_p) )
})
test_that("Persistence method behave as expected",
{
#Situation A: +Zero; (generally) correct if jump, wrong otherwise
- pred00_sd = computeForecast(data00, indices, "Persistence", "Zero", Inf, 24,
+ pred00_sd = computeForecast(data_p, indices, "Persistence", "Zero", 1, Inf, 24,
ncores=1, same_day=TRUE)
- pred00_dd = computeForecast(data00, indices, "Persistence", "Zero", Inf, 24,
+ pred00_dd = computeForecast(data_p, indices, "Persistence", "Zero", 1, Inf, 24,
ncores=1, same_day=FALSE)
for (i in 1:7)
{
expect_equal(pred00_dd$getForecast(i), rep(pred_order[i],24))
}
- pred13_sd = computeForecast(data13, indices, "Persistence", "Zero", Inf, 24,
+ pred13_sd = computeForecast(data_p, indices, "Persistence", "Zero", 14, Inf, 24,
ncores=1, same_day=TRUE)
- pred13_dd = computeForecast(data13, indices, "Persistence", "Zero", Inf, 24,
+ pred13_dd = computeForecast(data_p, indices, "Persistence", "Zero", 14, Inf, 24,
ncores=1, same_day=FALSE)
- for (i in 2:6)
+ for (i in 1:7)
{
- expect_equal(pred13_sd$getForecast(i), c( rep(i,11), rep(i%%7+1,13) ) )
- expect_equal(pred13_dd$getForecast(i), c( rep(i,11), rep(i%%7+1,13) ) )
+ expect_equal(pred13_sd$getForecast(i), rep(i,24) )
+ expect_equal(pred13_dd$getForecast(i), rep(i,24) )
}
- #boundaries are special cases: OK if same day, quite wrong otherwise
- expect_equal(pred13_sd$getForecast(1), c( rep(1,11), rep(2,13) ) )
- expect_equal(pred13_dd$getForecast(1), c( rep(1,11), rep(-5,13) ) )
- expect_equal(pred13_sd$getForecast(7), c( rep(7,11), rep(1,13) ) )
- expect_equal(pred13_dd$getForecast(7), c( rep(7,11), rep(8,13) ) )
#Situation B: +Persistence, (generally) correct
- pred00_sd = computeForecast(data00, indices, "Persistence", "Persistence", Inf, 24,
+ pred00_sd = computeForecast(data_p, indices, "Persistence", "Persistence", 1, Inf, 24,
ncores=1, same_day=TRUE)
- pred00_dd = computeForecast(data00, indices, "Persistence", "Persistence", Inf, 24,
+ pred00_dd = computeForecast(data_p, indices, "Persistence", "Persistence", 1, Inf, 24,
ncores=1, same_day=FALSE)
for (i in 3:7)
{
expect_equal(pred00_sd$getForecast(2), rep(2,24) )
expect_equal(pred00_dd$getForecast(2), rep(-5,24) )
- pred13_sd = computeForecast(data13, indices, "Persistence", "Persistence", Inf, 24,
+ pred13_sd = computeForecast(data_p, indices, "Persistence", "Persistence", 14, Inf, 24,
ncores=1, same_day=TRUE)
- pred13_dd = computeForecast(data13, indices, "Persistence", "Persistence", Inf, 24,
+ pred13_dd = computeForecast(data_p, indices, "Persistence", "Persistence", 14, Inf, 24,
ncores=1, same_day=FALSE)
- for (i in 2:6)
+ for (i in 1:7)
{
- expect_equal(pred13_sd$getForecast(i), c( rep(i,11), rep(i%%7+1,13) ) )
- expect_equal(pred13_dd$getForecast(i), c( rep(i,11), rep(i%%7+1,13) ) )
+ expect_equal(pred13_sd$getForecast(i), rep(i,24) )
+ expect_equal(pred13_dd$getForecast(i), rep(i,24) )
}
- #boundaries are special cases: OK if same day, quite wrong otherwise
- expect_equal(pred13_sd$getForecast(1), c( rep(1,11), rep(2,13) ) )
- expect_equal(pred13_dd$getForecast(1), c( rep(1,11), rep(-5,13) ) )
- expect_equal(pred13_sd$getForecast(7), c( rep(7,11), rep(1,13) ) )
- expect_equal(pred13_dd$getForecast(7), c( rep(7,11), rep(8,13) ) )
#A few extra checks
- expect_equal( pred00_sd$getIndexInData(3), dateIndexToInteger("2007-04-03",data00) )
- expect_equal( pred00_dd$getIndexInData(6), dateIndexToInteger("2007-04-06",data00) )
- expect_equal( pred13_sd$getIndexInData(3), dateIndexToInteger("2007-04-03",data13) )
- expect_equal( pred13_dd$getIndexInData(6), dateIndexToInteger("2007-04-06",data13) )
+ expect_equal( pred00_sd$getIndexInData(3), dateIndexToInteger("2007-04-04",data_p) )
+ expect_equal( pred00_dd$getIndexInData(6), dateIndexToInteger("2007-04-07",data_p) )
+ expect_equal( pred13_sd$getIndexInData(3), dateIndexToInteger("2007-04-04",data_p) )
+ expect_equal( pred13_dd$getIndexInData(6), dateIndexToInteger("2007-04-07",data_p) )
})
test_that("Neighbors method behave as expected",
{
#Situation A: +Zero; correct if jump, wrong otherwise
- pred00 = computeForecast(data00, indices, "Neighbors", "Zero", Inf, 24,
- simtype="mix", local=FALSE)
+ pred00 = computeForecast(data_p, indices, "Neighbors", "Zero", 1, Inf, 24,
+ simtype="mix", local=FALSE, window=c(1,1))
for (i in 1:7)
expect_equal(pred00$getForecast(i), rep(pred_order[i],24))
- pred13 = computeForecast(data13, indices, "Persistence", "Zero", Inf, 24,
- simtype="mix", local=FALSE)
+ pred13 = computeForecast(data_p, indices, "Persistence", "Zero", 14, Inf, 24,
+ simtype="mix", local=FALSE, window=c(1,1))
for (i in 1:7)
- expect_equal(pred13$getForecast(i), c( rep(i,11), rep(i%%7+1,13) ) )
+ expect_equal(pred13$getForecast(i), rep(i,24) )
#Situation B: +Neighbors == too difficult to eval in a unit test
-# pred00 = computeForecast(data00, indices, "Neighbors", "Neighbors", Inf, 24,
+# pred00 = computeForecast(data_p, indices, "Neighbors", "Neighbors", 1, Inf, 24,
# simtype="endo", local=FALSE)
# jumps = ...
# for (i in 1:7)
# expect_equal(pred00$getForecast(i), rep(pred_order[i]+jumps[i],24))
-# pred13 = computeForecast(data13, indices, "Neighbors", "Neighbors", Inf, 24,
+# pred13 = computeForecast(data_p, indices, "Neighbors", "Neighbors", 14, Inf, 24,
# simtype="endo", local=FALSE)
# for (i in 1:7)
# expect_equal(pred13$getForecast(i), c( rep(i,11), rep(i%%7+1,13) ) )
#A few extra checks
- expect_equal( pred00$getIndexInData(1), dateIndexToInteger("2007-04-01",data00) )
- expect_equal( pred00$getIndexInData(4), dateIndexToInteger("2007-04-04",data00) )
- expect_equal( pred13$getIndexInData(1), dateIndexToInteger("2007-04-01",data13) )
- expect_equal( pred13$getIndexInData(4), dateIndexToInteger("2007-04-04",data13) )
+ expect_equal( pred00$getIndexInData(1), dateIndexToInteger("2007-04-02",data_p) )
+ expect_equal( pred00$getIndexInData(4), dateIndexToInteger("2007-04-05",data_p) )
+ expect_equal( pred13$getIndexInData(1), dateIndexToInteger("2007-04-02",data_p) )
+ expect_equal( pred13$getIndexInData(4), dateIndexToInteger("2007-04-05",data_p) )
})
{
data = getDataTest(150)
- # index 143 : serie type 2
- pred = computeForecast(data, 143, "Neighbors", "Zero",
- horizon=length(data$getSerie(1)), simtype="endo", local=FALSE, h_window=1)
- f = computeFilaments(data, pred, 1, limit=60, plot=FALSE)
+ # index 144 : serie type 3, yersteday type 2
+ pred = computeForecast(data, 144, "Neighbors", "Zero", predict_from=1,
+ horizon=length(data$getSerie(1)), simtype="endo", local=FALSE, window=1, opera=TRUE)
+ f = computeFilaments(data, pred, 1, 8, limit=60, plot=FALSE)
- # Expected output: 50-3-10 series of type 2, then 23 series of type 3 (closest next)
+ # Expected output: 50-3-10 series of type 2+1 = 3,
+ # then 23 series of type 3+1 %% 3 = 1 (3 = closest next)
expect_identical(length(f$neighb_indices), as.integer(60))
expect_identical(length(f$colors), as.integer(60))
- expect_equal(f$index, 143)
- expect_true(all(I(f$neighb_indices) >= 2))
+ expect_equal(f$index, 144)
+ expect_true(all(I(f$neighb_indices) != 2))
for (i in 1:37)
{
- expect_equal(I(f$neighb_indices[i]), 2)
+ expect_equal(I(f$neighb_indices[i]), 3)
expect_match(f$colors[i], f$colors[1])
}
for (i in 38:60)
{
- expect_equal(I(f$neighb_indices[i]), 3)
+ expect_equal(I(f$neighb_indices[i]), 1)
expect_match(f$colors[i], f$colors[38])
}
expect_match(f$colors[1], "#1*")
expect_match(f$colors[38], "#E*")
- # index 142 : serie type 1
- pred = computeForecast(data, 142, "Neighbors", "Zero",
- horizon=length(data$getSerie(1)), simtype="endo", local=FALSE, h_window=1)
- f = computeFilaments(data, pred, 1, limit=50, plot=FALSE)
+ # index 143 : serie type 2
+ pred = computeForecast(data, 143, "Neighbors", "Zero", predict_from=1,
+ horizon=length(data$getSerie(1)), simtype="endo", local=FALSE, window=1, opera=TRUE)
+ f = computeFilaments(data, pred, 1, 8, limit=50, plot=FALSE)
- # Expected output: 50-10-3 series of type 1, then 13 series of type 3 (closest next)
- # NOTE: -10 because only past days with no-NAs tomorrow => exclude type 1 in [60,90[
+ # Expected output: 50-10-3 series of type 1+1=2,
+ # then 13 series of type 3+1 %% 3 = 1 (closest next)
+ # NOTE: -10 because only past tomorrows with no-NAs yerstedays
+ # => exclude type 2 in [60,90[
expect_identical(length(f$neighb_indices), as.integer(50))
expect_identical(length(f$colors), as.integer(50))
- expect_equal(f$index, 142)
- expect_true(all(I(f$neighb_indices) != 2))
+ expect_equal(f$index, 143)
+ expect_true(all(I(f$neighb_indices) <= 2))
for (i in 1:37)
{
- expect_equal(I(f$neighb_indices[i]), 1)
+ expect_equal(I(f$neighb_indices[i]), 2)
expect_match(f$colors[i], f$colors[1])
}
for (i in 38:50)
{
- expect_equal(I(f$neighb_indices[i]), 3)
+ expect_equal(I(f$neighb_indices[i]), 1)
expect_match(f$colors[i], f$colors[38])
}
expect_match(f$colors[1], "#1*")
-----
# Résultats numériques
-Cette partie montre les résultats obtenus avec des variantes de l'algorithme décrit au
-chapitre , en utilisant le package présenté à la section 3. Cet algorithme est
+% if P == 8:
+Cette partie montre les résultats obtenus avec des variantes de l'algorithme décrit à la
+section 4, en utilisant le package présenté au chapitre précédent. Cet algorithme est
systématiquement comparé à deux approches naïves :
* la moyenne des lendemains des jours "similaires" dans tout le passé, c'est-à-dire
* la persistence, reproduisant le jour courant ou allant chercher le lendemain de la
dernière journée "similaire" (même principe que ci-dessus ; argument "same\_day").
-Concernant l'algorithme principal à voisins, trois variantes sont étudiées dans cette
+Concernant l'algorithme principal à voisins, deux variantes sont comparées dans cette
partie :
* avec simtype="mix" et raccordement "Neighbors" dans le cas "non local", i.e. on va
chercher des voisins n'importe où du moment qu'ils correspondent au premier élément d'un
couple de deux jours consécutifs sans valeurs manquantes.
- * avec simtype="endo" + raccordement "Neighbors" puis simtype="none" + raccordement
-"Zero" (sans ajustement) dans le cas "local" : voisins de même niveau de pollution et
-même saison.
+ * avec simtype="none" (moyenne simple) et raccordement=NULL (aucun ajustement après
+moyenne des courbes) dans le cas "local" : voisins de même niveau de pollution et même
+saison.
Pour chaque période retenue $-$ chauffage, épandage, semaine non polluée $-$ les erreurs
de prédiction sont d'abord affichées, puis quelques graphes de courbes réalisées/prévues
(sur le jour "en moyenne le plus facile" à gauche, et "en moyenne le plus difficile" à
droite). Ensuite plusieurs types de graphes apportant des précisions sur la nature et la
difficulté du problème viennent compléter ces premières courbes. Concernant les graphes
-de filaments, la moitié gauche du graphe correspond aux jours similaires au jour courant,
-tandis que la moitié droite affiche les lendemains : ce sont donc les voisinages tels
-qu'utilisés dans l'algorithme.
+de filaments, la moitié droite du graphe correspond aux jours similaires au jour courant,
+tandis que la moitié gauche affiche les jours précédents : ce sont donc les voisinages
+tels qu'utilisés dans l'algorithme.
+% endif
<%
list_titles = ['Pollution par chauffage','Pollution par épandage','Semaine non polluée']
list_indices = ['indices_ch', 'indices_ep', 'indices_np']
-----r
library(talweg)
-P = ${P} #instant de prévision
-H = ${H} #horizon (en heures)
+P = ${P} #première heure de prévision
+H = ${H} #dernière heure de prévision
ts_data = read.csv(system.file("extdata","pm10_mesures_H_loc_report.csv",
package="talweg"))
exo_data = read.csv(system.file("extdata","meteo_extra_noNAs.csv",
package="talweg"))
-# NOTE: 'GMT' because DST gaps are filled and multiple values merged in
-# above dataset. Prediction from P+1 to P+H included.
-data = getData(ts_data, exo_data, input_tz = "GMT", working_tz="GMT",
- predict_at=P)
+data = getData(ts_data, exo_data)
-indices_ch = seq(as.Date("2015-01-18"),as.Date("2015-01-24"),"days")
-indices_ep = seq(as.Date("2015-03-15"),as.Date("2015-03-21"),"days")
-indices_np = seq(as.Date("2015-04-26"),as.Date("2015-05-02"),"days")
+indices_ch = seq(as.Date("2015-01-19"),as.Date("2015-01-25"),"days")
+indices_ep = seq(as.Date("2015-03-16"),as.Date("2015-03-22"),"days")
+indices_np = seq(as.Date("2015-04-27"),as.Date("2015-05-03"),"days")
% for i in range(3):
-----
##<h2 style="color:blue;font-size:2em">${list_titles[i]}</h2>
${"##"} ${list_titles[i]}
-----r
-p1 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Neighbors", horizon=H,
- simtype="mix", local=FALSE)
-p2 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Neighbors", horizon=H,
- simtype="endo", local=TRUE)
-p3 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Zero", horizon=H,
- simtype="none", local=TRUE)
-p4 = computeForecast(data, ${list_indices[i]}, "Average", "Zero", horizon=H)
-p5 = computeForecast(data, ${list_indices[i]}, "Persistence", "Zero", horizon=H,
- same_day=${'TRUE' if loop.index < 2 else 'FALSE'})
+p1 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Neighbors",
+ predict_from=P, horizon=H, simtype="mix", local=FALSE)
+p2 = computeForecast(data, ${list_indices[i]}, "Neighbors", NULL,
+ predict_from=P, horizon=H, simtype="none", local=TRUE)
+p3 = computeForecast(data, ${list_indices[i]}, "Average", "Zero",
+ predict_from=P, horizon=H)
+p4 = computeForecast(data, ${list_indices[i]}, "Persistence", "Zero",
+ predict_from=P, horizon=H, same_day=${'TRUE' if loop.index < 2 else 'FALSE'})
-----r
-e1 = computeError(data, p1, H)
-e2 = computeError(data, p2, H)
-e3 = computeError(data, p3, H)
-e4 = computeError(data, p4, H)
-e5 = computeError(data, p5, H)
+e1 = computeError(data, p1, P, H)
+e2 = computeError(data, p2, P, H)
+e3 = computeError(data, p3, P, H)
+e4 = computeError(data, p4, P, H)
options(repr.plot.width=9, repr.plot.height=7)
-plotError(list(e1, e5, e4, e2, e3), cols=c(1,2,colors()[258],4,6))
+plotError(list(e1, e4, e3, e2), cols=c(1,2,colors()[258],4))
-# noir: Neighbors non-local (p1), bleu: Neighbors local endo (p2),
-# mauve: Neighbors local none (p3), vert: moyenne (p4),
-# rouge: persistence (p5)
+# noir: Neighbors non-local (p1), bleu: Neighbors local (p2),
+# vert: moyenne (p3), rouge: persistence (p4)
-sum_p123 = e1$abs$indices + e2$abs$indices + e3$abs$indices
-i_np = which.min(sum_p123) #indice de (veille de) jour "facile"
-i_p = which.max(sum_p123) #indice de (veille de) jour "difficile"
+sum_p23 = e2$abs$indices + e3$abs$indices
+i_np = which.min(sum_p23) #indice de jour "facile"
+i_p = which.max(sum_p23) #indice de jour "difficile"
+% if P == 8:
-----
% if i == 0:
-L'erreur absolue dépasse 20 sur 1 à 2 jours suivant les modèles (graphe en haut à
-droite). Sur cet exemple le modèle à voisins "contraint" (local=TRUE) utilisant des
-pondérations basées sur les similarités de forme (simtype="endo") obtient en moyenne les
-meilleurs résultats, avec un MAPE restant en général inférieur à 30% de 8h à 19h (7+1 à
-7+12 : graphe en bas à gauche).
+L'erreur absolue $-$ en haut à droite $-$ reste modérée pour les meilleurs modèles
+(variantes à voisins), ne dépassant 10 que deux jours. Les deux modèles naïfs ont des
+erreurs similaires sauf sur la période "difficile" (jours 4 à 6), sur laquelle on gagne
+donc à chercher des jours semblables pour effectuer la prévision.
+Le MAPE reste en général inférieur à 35% pour les meilleurs méthodes.
% elif i == 1:
-Il est difficile dans ce cas de déterminer une méthode meilleure que les autres : elles
-donnent toutes de plutôt mauvais résultats, avec une erreur absolue moyennée sur la
-journée dépassant presque toujours 15 (graphe en haut à droite).
+Le modèle à voisins avec contrainte de localité obtient ici les meilleurs résultats, son
+erreur étant clairement en dessous des autres à partir du jour 4 (graphe en haut à
+droite). Le MAPE jour après jour est du même ordre que précédemment pour cette méthode
+(35%, graphe en bas à droite) sauf un jour sur lequel le MAPE explose.
% else:
Dans ce cas plus favorable les intensité des erreurs absolues ont clairement diminué :
-elles restent souvent en dessous de 5. En revanche le MAPE moyen reste au-delà de 20%, et
-même souvent plus de 30%. Comme dans le cas de l'épandage on constate une croissance
-globale de la courbe journalière d'erreur absolue moyenne (en haut à gauche) ; ceci peut
-être dû au fait que l'on ajuste le niveau du jour à prédire en le recollant sur la
-dernière valeur observée.
+elles sont souvent en dessous de 5. En revanche le MAPE moyen reste en général au-delà de
+20%. Comme dans le cas de l'épandage on constate une croissance globale de la courbe
+journalière d'erreur absolue moyenne (en haut à gauche) $-$ sauf pour la méthode à
+voisins "locale" ; ceci peut être dû au fait que l'on ajuste le niveau du jour à prédire
+en le recollant sur la dernière valeur observée (sauf pour "Neighbors local").
+% endif
% endif
-----r
options(repr.plot.width=9, repr.plot.height=4)
plotPredReal(data, p2, i_np); title(paste("PredReal p2 day",i_np))
plotPredReal(data, p2, i_p); title(paste("PredReal p2 day",i_p))
-plotPredReal(data, p3, i_np); title(paste("PredReal p3 day",i_np))
-plotPredReal(data, p3, i_p); title(paste("PredReal p3 day",i_p))
-
-# Bleu : prévue ; noir : réalisée
+# Bleu : prévue ; noir : réalisée (confondues jusqu'à predict_from-1)
+% if P == 8:
-----
% if i == 0:
+<<<<<<< HEAD
La courbe du jour "facile à prévoir", à gauche, se décompose en deux modes : un léger
vers 10h (7+3), puis un beaucoup plus marqué vers 19h (7+12). Ces deux modes sont
retrouvés par les trois variantes de l'algorithme à voisins, bien que l'amplitude soit
début et toute fin de journée (à 9h et 23h), qui ne sont pas du tout anticipés par les
méthodes ; la grande amplitude de ces pics explique alors l'intensité de l'erreur
observée.
+=======
+Le jour "facile à prévoir", à gauche, se décompose en deux modes : un léger vers 10h
+(7+3), puis un beaucoup plus marqué vers 19h (7+12). Ces deux modes sont retrouvés par
+les deux variantes de l'algorithme à voisins, bien que l'amplitude soit mal prédite.
+Concernant le jour "difficile à prévoir" (à droite) il y a deux pics en tout début et
+toute fin de journée (à 9h et 23h), qui ne sont pas du tout anticipés par les méthodes ;
+la grande amplitude de ces pics explique alors l'intensité de l'erreur observée.
+>>>>>>> 7c4b2952874de1d40a742e72efe51999b99050f5
% elif i == 1:
-Dans le cas d'un jour "facile" à prédire $-$ à gauche $-$ la forme est plus ou moins
-retrouvée, mais le niveau moyen est trop bas (courbe en bleu). Concernant le jour
-"difficile" à droite, non seulement la forme n'est pas anticipée mais surtout le niveau
-prédit est très inférieur au niveau de pollution observé. Comme on le voit ci-dessous
-cela découle d'un manque de voisins au comportement similaire.
+Dans le cas d'un jour "facile" à prédire $-$ à gauche $-$ la forme est plutôt bien
+retrouvée, ainsi que le niveau moyen pour la méthode sans contrainte de localité
+(dans l'autre, l'algorithme a probablement écarté trop de voisins potentiels).
+Concernant le jour "difficile" à droite, non seulement la forme n'est pas anticipée mais
+surtout le niveau prédit est largement supérieur au niveau de pollution observé $-$ dans
+une moindre mesure toutefois pour la variante "locale".
% else:
-La forme est raisonnablement retrouvée pour les méthodes "locales", l'autre version
-lissant trop les prédictions. Le biais reste cependant important, surtout en fin de
-journée sur la courbes "difficile à prévoir".
+L'impression visuelle est plutôt mauvaise dans ce cas, mais les écart étant minimes les
+erreurs au final ne sont pas très importantes. De plus deux des quatres graphes sont
+satisfaisants (en haut à droite et en bas à gauche : forme + niveau acceptables.
+% endif
% endif
-----r
par(mfrow=c(1,2))
+
f_np1 = computeFilaments(data, p1, i_np, plot=TRUE)
- title(paste("Filaments p1 day",i_np))
+title(paste("Filaments p1 day",i_np))
+
f_p1 = computeFilaments(data, p1, i_p, plot=TRUE)
- title(paste("Filaments p1 day",i_p))
+title(paste("Filaments p1 day",i_p))
f_np2 = computeFilaments(data, p2, i_np, plot=TRUE)
- title(paste("Filaments p2 day",i_np))
+title(paste("Filaments p2 day",i_np))
+
f_p2 = computeFilaments(data, p2, i_p, plot=TRUE)
- title(paste("Filaments p2 day",i_p))
+title(paste("Filaments p2 day",i_p))
+% if P == 8:
-----
% if i == 0:
Les voisins du jour courant (période de 24h allant de 8h à 7h le lendemain) sont affichés
difficultés rencontrées par l'algorithme sur ce jeu de données.**
% elif i == 1:
Les observations sont les mêmes qu'au paragraphe précédent : trop de variabilité des
-lendemains (et même des voisins du jour courant).
+voisins (et ce même le jour précédent).
% else:
Les graphes de filaments ont encore la même allure, avec une assez grande variabilité
observée. Cette observation est cependant trompeuse, comme l'indique plus bas le graphe
de variabilité relative.
% endif
+% endif
-----r
par(mfrow=c(1,2))
-plotFilamentsBox(data, f_np1); title(paste("FilBox p1 day",i_np))
-plotFilamentsBox(data, f_p1); title(paste("FilBox p1 day",i_p))
-# En pointillés la courbe du jour courant + lendemain (à prédire)
+plotFilamentsBox(data, f_np1, predict_from=P)
+title(paste("FilBox p1 day",i_np))
+
+plotFilamentsBox(data, f_p1, predict_from=P)
+title(paste("FilBox p1 day",i_p))
+
+# En pointillés la courbe du jour courant (à prédire) + précédent
+% if P == 8:
-----
% if i == 0:
Sur cette boxplot fonctionnelle (voir la fonction fboxplot() du package R "rainbow") on
des voisins sont trop éloignés de la forme à prédire et forcent ainsi un aplatissement de
la prédiction.
% elif i == 1:
-On constate la présence d'un voisin au lendemain complètement atypique avec un pic en
-début de journée (courbe en vert à gauche), et d'un autre phénomène semblable avec la
-courbe rouge sur le graphe de droite. Ajouté au fait que le lendemain à prévoir est
-lui-même un jour "hors norme", cela montre l'impossibilité de bien prévoir une courbe en
-utilisant l'algorithme à voisins.
+Concernant le jour "difficile" on constate la présence de voisins au lendemains
+complètement atypiques avec un pic en début de journée (courbes en vert et rouge à
+droite). Ajouté au fait que le jour à prévoir est lui-même "hors norme", cela montre
+l'impossibilité de bien prévoir une courbe en utilisant l'algorithme à voisins.
% else:
On peut réappliquer les mêmes remarques qu'auparavant sur les boxplots fonctionnels :
-lendemains de voisins atypiques, courbe à prévoir elle-même légèrement "hors norme".
+voisins atypiques, courbe à prévoir elle-même légèrement "hors norme".
+% endif
% endif
-----r
par(mfrow=c(1,2))
-plotRelVar(data, f_np1); title(paste("StdDev p1 day",i_np))
-plotRelVar(data, f_p1); title(paste("StdDev p1 day",i_p))
-plotRelVar(data, f_np2); title(paste("StdDev p2 day",i_np))
-plotRelVar(data, f_p2); title(paste("StdDev p2 day",i_p))
+plotRelVar(data, f_np1, predict_from=P)
+title(paste("StdDev p1 day",i_np))
-# Variabilité globale en rouge ; sur les voisins (+ lendemains) en noir
+plotRelVar(data, f_p1, predict_from=P)
+title(paste("StdDev p1 day",i_p))
+
+plotRelVar(data, f_np2, predict_from=P)
+title(paste("StdDev p2 day",i_np))
+
+plotRelVar(data, f_p2, predict_from=P)
+title(paste("StdDev p2 day",i_p))
+
+# Variabilité globale en rouge ; sur les voisins en noir
+% if P == 8:
-----
% if i == 0:
Ces graphes viennent confirmer l'impression visuelle après observation des filaments. En
effet, la variabilité globale en rouge (écart-type heure par heure sur l'ensemble des
-couples "aujourd'hui/lendemain"du passé) devrait rester nettement au-dessus de la
+couples "hier/aujourd'hui" du passé) devrait rester nettement au-dessus de la
variabilité locale, calculée respectivement sur un voisinage d'une soixantaine de jours
-(pour p1) et d'une dizaine de jours (pour p2). Or on constate que ce n'est pas du tout le
-cas sur la période "lendemain", sauf en partie pour p2 le jour 4 $-$ mais ce n'est pas
-suffisant.
+(pour p1) et d'une dizaine de jours (pour p2). Or ce n'est pas du tout le cas sur la
+moitié droite, sauf pour le jour "facile" avec l'algorithme "local".
% elif i == 1:
-Comme précédemment les variabilités locales et globales sont confondues dans les parties
-droites des graphes $-$ sauf pour la version "locale" sur le jour "facile"; mais cette
-bonne propriété n'est pas suffisante si l'on ne trouve pas les bons poids à appliquer.
+Comme précédemment les variabilités locales et globales sont trop proches dans les
+parties droites des graphes pour le jour "difficile". L'allure des graphes est
+raisonnable ppour l'autre jour, qui est d'ailleurs bien prédit.
% else:
Cette fois la situation idéale est observée : la variabilité globale est nettement
au-dessus de la variabilité locale. Bien que cela ne suffise pas à obtenir de bonnes
prédictions de forme, on constate au moins l'amélioration dans la prédiction du niveau.
% endif
+% endif
-----r
-par(mfrow=c(1,2))
-plotSimils(p1, i_np); title(paste("Weights p1 day",i_np))
-plotSimils(p1, i_p); title(paste("Weights p1 day",i_p))
+plotSimils(p1, i_np)
+title(paste("Weights p1 day",i_np))
-plotSimils(p2, i_np); title(paste("Weights p2 day",i_np))
-plotSimils(p2, i_p); title(paste("Weights p2 day",i_p))
+plotSimils(p1, i_p)
+title(paste("Weights p1 day",i_p))
+
+# Poids < 1/N à gauche, >= 1/N à droite ; jour facile en haut, difficile en bas
+% if P == 8:
-----
% if i == 0:
-Les poids se concentrent près de 0 dans le cas "non local" (p1), et se répartissent assez
-uniformément dans [ 0, 0.2 ] dans le cas "local" (p2). C'est ce que l'on souhaite
-observer pour éviter d'effectuer une simple moyenne.
+Les poids se concentrent près de 0 : c'est ce que l'on souhaite observer pour éviter
+d'effectuer une simple moyenne.
% elif i == 1:
-En comparaison avec le pragraphe précédent on retrouve le même (bon) comportement des
-poids pour la version "non locale". En revanche la fenêtre optimisée est trop grande sur
-le jour "facile" pour la méthode "locale" (voir affichage ci-dessous) : il en résulte des
-poids tous semblables autour de 0.084, l'algorithme effectue donc une moyenne simple $-$
-expliquant pourquoi les courbes mauve et bleue sont très proches sur le graphe d'erreurs.
+On retrouve le même (bon) comportement des poids : concentration vers 0, quelques poids
+non négligeables (presque trop peu pour le jour "difficile").
% else:
-Concernant les poids en revanche, deux cas a priori mauvais se cumulent :
-
- * les poids dans le cas "non local" ne sont pas assez concentrés autour de 0, menant à
-un lissage trop fort $-$ comme observé sur les graphes des courbes réalisées/prévues ;
- * les poids dans le cas "local" sont trop semblables (à cause de la trop grande fenêtre
-optimisée par validation croisée, cf. ci-dessous), résultant encore en une moyenne simple
-$-$ mais sur moins de jours, plus proches du jour courant.
+Les poids sont répartis comme souhaité : concentrés vers 0 avec quelques valeurs non
+négligeables.
+% endif
% endif
-----r
-# Fenêtres sélectionnées dans ]0,7] :
-# "non-local" 2 premières lignes, "local" ensuite
-p1$getParams(i_np)$window
-p1$getParams(i_p)$window
+options(digits=2)
-p2$getParams(i_np)$window
-p2$getParams(i_p)$window
+print(p1$getParams(i_np)$window)
+print(p1$getParams(i_p)$window)
+
+# Fenêtres sélectionnées dans ]0,7]
% endfor
+% if P == 8:
-----
${"##"} Bilan
-Nos algorithmes à voisins ne sont pas adaptés à ce jeu de données où la forme varie
-considérablement d'un jour à l'autre. Toutefois, un espoir reste permis par exemple en
-aggrégeant les courbes spatialement (sur plusieurs stations situées dans la même
+Nos algorithmes à voisins donnent de meilleurs résultats que les approches naïves
+(persistence, moyenne sur tout le jeu de données). Les erreurs restent cependant assez
+élevées, notamment en terme de MAPE. Une possible poste d'amélioration consisterait à
+aggréger les courbes spatialement (sur plusieurs stations situées dans la même
agglomération ou dans une même zone).
+% endif
+++ /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
-
-\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}
+++ /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": null,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "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
------
-# Résultats numériques
-
-Cette partie montre les résultats obtenus avec des variantes de l'algorithme décrit au
-chapitre 5, en utilisant le package présenté au chapitre 6.
-Les ........... options ...........
-Cet algorithme est
-systématiquement comparé à deux approches naïves :
-
- * la moyenne des lendemains des jours de même type dans tout le passé, c'est-à-dire
-prédiction = moyenne de tous les mardis passés si le jour courant est un lundi.
- * la persistence, reproduisant le jour courant ou allant chercher le lendemain de la
-dernière journée de même type (même principe que ci-dessus ; argument "same\_day").
-
-Concernant l'algorithme principal à voisins, trois variantes sont étudiées dans cette
-partie :
-
- * avec simtype="mix" et raccordement "Neighbors" dans le cas "non local", i.e. on va
-chercher des voisins n'importe où du moment qu'ils correspondent au premier élément d'un
-couple de deux jours consécutifs sans valeurs manquantes.
- * avec simtype="endo" + raccordement "Neighbors" puis simtype="none" + raccordement
-"Zero" (sans ajustement) dans le cas "local" : voisins de même niveau de pollution et
-même saison.
-
-Pour chaque période retenue $-$ chauffage, épandage, semaine non polluée $-$ les erreurs
-de prédiction sont d'abord affichées, puis quelques graphes de courbes réalisées/prévues
-(sur le jour "en moyenne le plus facile" à gauche, et "en moyenne le plus difficile" à
-droite). Ensuite plusieurs types de graphes apportant des précisions sur la nature et la
-difficulté du problème viennent compléter ces premières courbes. Concernant les graphes
-de filaments, la moitié gauche du graphe correspond aux jours similaires au jour courant,
-tandis que la moitié droite affiche les lendemains : ce sont donc les voisinages tels
-qu'utilisés dans l'algorithme.
-<%
-list_titles = ['Pollution par chauffage','Pollution par épandage','Semaine non polluée']
-list_indices = ['indices_ch', 'indices_ep', 'indices_np']
-%>
------r
-library(talweg)
-
-P = ${P} #instant de prévision
-H = ${H} #horizon (en heures)
-
-ts_data = read.csv(system.file("extdata","pm10_mesures_H_loc_report.csv",
- package="talweg"))
-exo_data = read.csv(system.file("extdata","meteo_extra_noNAs.csv",
- package="talweg"))
-# NOTE: 'GMT' because DST gaps are filled and multiple values merged in
-# above dataset. Prediction from P+1 to P+H included.
-data = getData(ts_data, exo_data, input_tz = "GMT", working_tz="GMT",
- predict_at=P)
-
-indices_ch = seq(as.Date("2015-01-18"),as.Date("2015-01-24"),"days")
-indices_ep = seq(as.Date("2015-03-15"),as.Date("2015-03-21"),"days")
-indices_np = seq(as.Date("2015-04-26"),as.Date("2015-05-02"),"days")
-% for i in range(3):
------
-##<h2 style="color:blue;font-size:2em">${list_titles[i]}</h2>
-${"##"} ${list_titles[i]}
------r
-p1 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Neighbors", horizon=H,
- simtype="mix", local=FALSE)
-p2 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Neighbors", horizon=H,
- simtype="endo", local=TRUE)
-p3 = computeForecast(data, ${list_indices[i]}, "Neighbors", "Zero", horizon=H,
- simtype="none", local=TRUE)
-p4 = computeForecast(data, ${list_indices[i]}, "Average", "Zero", horizon=H)
-p5 = computeForecast(data, ${list_indices[i]}, "Persistence", "Zero", horizon=H,
- same_day=${'TRUE' if loop.index < 2 else 'FALSE'})
------r
-e1 = computeError(data, p1, H)
-e2 = computeError(data, p2, H)
-e3 = computeError(data, p3, H)
-e4 = computeError(data, p4, H)
-e5 = computeError(data, p5, H)
-options(repr.plot.width=9, repr.plot.height=7)
-plotError(list(e1, e5, e4, e2, e3), cols=c(1,2,colors()[258],4,6))
-
-# noir: Neighbors non-local (p1), bleu: Neighbors local endo (p2),
-# mauve: Neighbors local none (p3), vert: moyenne (p4),
-# rouge: persistence (p5)
-
-##############TODO: expliquer "endo" "none"......etc
-## ajouter fenêtres essais dans rapport. --> dans chapitre actuel.
-## re-ajouter annexe sur ancienne méthode exo/endo/mix
-## ---------> fenetres comment elles sont optimisées
-#--------> ajouter à la fin quelques graphes montrant/comparant autres méthodes
-#chapitre résumé avec différents essais conclusions. ---> synthèse des essais réalisés,
-#avec sous-paragraphes avec conclusions H3/H17 sans surprises on améliore les choses,
-#mais il y a des situations où c'est pas mieux.
-#---------> fichier tex réinsérer synthèse de l'ensemble des essais réalisés.
-#++++++++ ajouter à 13h
-
-sum_p123 = e1$abs$indices + e2$abs$indices + e3$abs$indices
-i_np = which.min(sum_p123) #indice de (veille de) jour "facile"
-i_p = which.max(sum_p123) #indice de (veille de) jour "difficile"
------
-% if i == 0:
-L'erreur absolue dépasse 20 sur 1 à 2 jours suivant les modèles (graphe en haut à
-droite). ##C'est au-delà de ce que l'on aimerait voir (disons +/- 5 environ).
-Sur cet
-exemple le modèle à voisins "contraint" (local=TRUE) utilisant des pondérations basées
-sur les similarités de forme (simtype="endo") obtient en moyenne les meilleurs résultats,
-avec un MAPE restant en général inférieur à 30% de 8h à 19h (7+1 à 7+12 : graphe en bas à
-gauche).
-% elif i == 1:
-Il est difficile dans ce cas de déterminer une méthode meilleure que les autres : elles
-donnent toutes de plutôt mauvais résultats, avec une erreur absolue moyennée sur la
-journée dépassant presque toujours 15 (graphe en haut à droite).
-% else:
-Dans ce cas plus favorable les intensité des erreurs absolues ont clairement diminué :
-elles restent souvent en dessous de 5. En revanche le MAPE moyen reste au-delà de 20%, et
-même souvent plus de 30%. Comme dans le cas de l'épandage on constate une croissance
-globale de la courbe journalière d'erreur absolue moyenne (en haut à gauche) ; ceci peut
-être dû au fait que l'on ajuste le niveau du jour à prédire en le recollant sur la
-dernière valeur observée.
-% endif
------r
-options(repr.plot.width=9, repr.plot.height=4)
-par(mfrow=c(1,2))
-
-plotPredReal(data, p1, i_np); title(paste("PredReal p1 day",i_np))
-plotPredReal(data, p1, i_p); title(paste("PredReal p1 day",i_p))
-
-plotPredReal(data, p2, i_np); title(paste("PredReal p2 day",i_np))
-plotPredReal(data, p2, i_p); title(paste("PredReal p2 day",i_p))
-
-plotPredReal(data, p3, i_np); title(paste("PredReal p3 day",i_np))
-plotPredReal(data, p3, i_p); title(paste("PredReal p3 day",i_p))
-
-# Bleu : prévue ; noir : réalisée
------
-% if i == 0:
-La courbe non centrée du jour facile à prévoir (en noir),
-##Le jour "facile à prévoir",
-à gauche, se décompose en deux modes : un léger vers 10h
-(7+3), puis un beaucoup plus marqué vers 19h (7+12). Ces deux modes sont retrouvés par
-les trois variantes de l'algorithme à voisins, bien que l'amplitude soit mal prédite.
-Concernant le jour "difficile à prévoir" (à droite) il y a deux pics en tout début et toute fin de
-journée (à 9h et 23h), qui ne sont pas du tout anticipés par les méthodes ; la grande
-amplitude de ces pics explique alors l'intensité de l'erreur observée.
-% elif i == 1:
-Dans le cas d'un jour "facile" à prédire $-$ à gauche $-$ la forme est plus ou moins
-retrouvée, mais le niveau moyen est trop bas (courbe en bleu). Concernant le jour
-"difficile" à droite, non seulement la forme n'est pas anticipée mais surtout le niveau
-prédit est très inférieur au niveau de pollution observé. Comme on le voit ci-dessous
-cela découle d'un manque de voisins au comportement similaire.
-% else:
-La forme est raisonnablement retrouvée pour les méthodes "locales", l'autre version
-lissant trop les prédictions. Le biais reste cependant important, surtout en fin de
-journée sur le jour "difficile".
-% endif
------r
-par(mfrow=c(1,2))
-f_np1 = computeFilaments(data, p1, i_np, plot=TRUE)
- title(paste("Filaments p1 day",i_np))
-f_p1 = computeFilaments(data, p1, i_p, plot=TRUE)
- title(paste("Filaments p1 day",i_p))
-
-f_np2 = computeFilaments(data, p2, i_np, plot=TRUE)
- title(paste("Filaments p2 day",i_np))
-f_p2 = computeFilaments(data, p2, i_p, plot=TRUE)
- title(paste("Filaments p2 day",i_p))
------
-% if i == 0:
-Les voisins du jour courant (période de 24h allant de 8h à 7h le lendemain) sont affichés
-avec un trait d'autant plus sombre qu'ils sont proches. On constate dans le cas non
-contraint (en haut) une grande variabilité des lendemains, très nette sur le graphe en
-haut à droite. Ceci indique une faible corrélation entre la forme d'une courbe sur une
-période de 24h et la forme sur les 24h suivantes ; **cette observation est la source des
-difficultés rencontrées par l'algorithme sur ce jeu de données.**
-% elif i == 1:
-Les observations sont les mêmes qu'au paragraphe précédent : trop de variabilité des
-lendemains (et même des voisins du jour courant).
-% else:
-Les graphes de filaments ont encore la même allure, avec une assez grande variabilité
-observée. Cette observation est cependant trompeuse, comme l'indique plus bas le graphe
-de variabilité relative.
-% endif
------r
-par(mfrow=c(1,2))
-plotFilamentsBox(data, f_np1); title(paste("FilBox p1 day",i_np))
-plotFilamentsBox(data, f_p1); title(paste("FilBox p1 day",i_p))
-
-## Questions :
-#7h VS 13h
-#est-ce que prévoir 24h ou 13 ou 3 facilite.
-#amplitude erreur raisonnable ? probleme facile difficile ?
-#place des exogènes ?
-#H = ?
-#épandage > chauffage > np
-
-# En pointillés la courbe du jour courant + lendemain (à prédire)
------
-% if i == 0:
-Sur cette boxplot fonctionnelle (voir la fonction fboxplot() du package R "rainbow") on
-constate essentiellement deux choses : le lendemain d'un voisin "normal" peut se révéler
-être une courbe atypique, fort éloignée de ce que l'on souhaite prédire (courbes bleue et
-rouge à gauche) ; et, dans le cas d'une courbe à prédire atypique (à droite) la plupart
-des voisins sont trop éloignés de la forme à prédire et forcent ainsi un aplatissement de
-la prédiction.
-% elif i == 1:
-On constate la présence d'un voisin au lendemain complètement atypique avec un pic en
-début de journée (courbe en vert à gauche), et d'un autre phénomène semblable avec la
-courbe rouge sur le graphe de droite. Ajouté au fait que le lendemain à prévoir est
-lui-même un jour "hors norme", cela montre l'impossibilité de bien prévoir une courbe en
-utilisant l'algorithme à voisins.
-% else:
-On peut réappliquer les mêmes remarques qu'auparavant sur les boxplots fonctionnels :
-lendemains de voisins atypiques, courbe à prévoir elle-même légèrement "hors norme".
-% endif
------r
-par(mfrow=c(1,2))
-plotRelVar(data, f_np1); title(paste("StdDev p1 day",i_np))
-plotRelVar(data, f_p1); title(paste("StdDev p1 day",i_p))
-
-plotRelVar(data, f_np2); title(paste("StdDev p2 day",i_np))
-plotRelVar(data, f_p2); title(paste("StdDev p2 day",i_p))
-
-# Variabilité globale en rouge ; sur les voisins (+ lendemains) en noir
------
-% if i == 0:
-Ces graphes viennent confirmer l'impression visuelle après observation des filaments. En
-effet, la variabilité globale en rouge (écart-type heure par heure sur l'ensemble des
-couples "aujourd'hui/lendemain"du passé) devrait rester nettement au-dessus de la
-variabilité locale, calculée respectivement sur un voisinage d'une soixantaine de jours
-(pour p1) et d'une dizaine de jours (pour p2). Or on constate que ce n'est pas du tout le
-cas sur la période "lendemain", sauf en partie pour p2 le jour 4 $-$ mais ce n'est pas
-suffisant.
-% elif i == 1:
-Comme précédemment les variabilités locales et globales sont confondues dans les parties
-droites des graphes $-$ sauf pour la version "locale" sur le jour "facile"; mais cette
-bonne propriété n'est pas suffisante si l'on ne trouve pas les bons poids à appliquer.
-% else:
-Cette fois la situation idéale est observée : la variabilité globale est nettement
-au-dessus de la variabilité locale. Bien que cela ne suffise pas à obtenir de bonnes
-prédictions de forme, on constate au moins l'amélioration dans la prédiction du niveau.
-% endif
------r
-par(mfrow=c(1,2))
-plotSimils(p1, i_np); title(paste("Weights p1 day",i_np))
-plotSimils(p1, i_p); title(paste("Weights p1 day",i_p))
-
-plotSimils(p2, i_np); title(paste("Weights p2 day",i_np))
-plotSimils(p2, i_p); title(paste("Weights p2 day",i_p))
------
-% if i == 0:
-Les poids se concentrent près de 0 dans le cas "non local" (p1), et se répartissent assez
-uniformément dans [ 0, 0.2 ] dans le cas "local" (p2). C'est ce que l'on souhaite
-observer pour éviter d'effectuer une simple moyenne.
-% elif i == 1:
-En comparaison avec le pragraphe précédent on retrouve le même (bon) comportement des
-poids pour la version "non locale". En revanche la fenêtre optimisée est trop grande sur
-le jour "facile" pour la méthode "locale" (voir affichage ci-dessous) : il en résulte des
-poids tous semblables autour de 0.084, l'algorithme effectue donc une moyenne simple $-$
-expliquant pourquoi les courbes mauve et bleue sont très proches sur le graphe d'erreurs.
-% else:
-Concernant les poids en revanche, deux cas a priori mauvais se cumulent :
-
- * les poids dans le cas "non local" ne sont pas assez concentrés autour de 0, menant à
-un lissage trop fort $-$ comme observé sur les graphes des courbes réalisées/prévues ;
- * les poids dans le cas "local" sont trop semblables (à cause de la trop grande fenêtre
-optimisée par validation croisée, cf. ci-dessous), résultant encore en une moyenne simple
-$-$ mais sur moins de jours, plus proches du jour courant.
-% endif
------r
-# Fenêtres sélectionnées dans ]0,7] :
-# "non-local" 2 premières lignes, "local" ensuite
-p1$getParams(i_np)$window
-p1$getParams(i_p)$window
-
-p2$getParams(i_np)$window
-p2$getParams(i_p)$window
-% endfor
------
-${"##"} Bilan
-
-Nos algorithmes à voisins ne sont pas adaptés à ce jeu de données où la forme varie
-considérablement d'un jour à l'autre.
-Toutefois, un espoir reste permis par exemple en aggrégeant les courbes spatialement (sur
-plusieurs stations situées dans la même agglomération ou dans une même zone).
-##Plus généralement cette décorrélation de forme rend
-##ardue la tâche de prévision pour toute autre méthode $-$ du moins, nous ne savons pas
-##comment procéder pour parvenir à une bonne précision.
library(talweg)
# Acquisition des données (depuis les fichiers CSV)
-ts_data <- read.csv(system.file("extdata","pm10_mesures_H_loc.csv",
- package="talweg"))
-exo_data <- read.csv(system.file("extdata","meteo_extra_noNAs.csv",
- package="talweg"))
-data <- getData(ts_data, exo_data, input_tz="GMT",
- date_format="%d/%m/%Y %H:%M", working_tz="GMT",
- predict_at=7, limit=120)
+ts_data <- read.csv(system.file("extdata","pm10_mesures_H_loc.csv", package="talweg"))
+exo_data <- read.csv(system.file("extdata","meteo_extra_noNAs.csv", package="talweg"))
+data <- getData(ts_data, exo_data, date_format="%d/%m/%Y %H:%M", limit=120)
# Plus de détails à la section 1 ci-après.
# Prédiction de 10 courbes (jours 102 à 111)
-pred <- computeForecast(data, 101:110, "Persistence", "Zero", memory=50,
- horizon=12, ncores=1)
+pred <- computeForecast(data, 101:110, "Persistence", "Zero", predict_from=8, memory=50,
+ horizon=24, ncores=1)
# Plus de détails à la section 2 ci-après.
# Calcul des erreurs (sur un horizon arbitraire <= horizon de prédiction)
-err <- computeError(data, pred, horizon=6)
+err <- computeError(data, pred, predict_from=8, horizon=20)
# Plus de détails à la section 3 ci-après.
# Puis voir ?plotError et les autres plot dans le paragraphe 'seealso'
première colonne contient les jours, les $m$ suivantes les variables mesurées pour ce
jour, et les $m$ dernières les variables prédites pour ce même jour. Dans notre cas $m=4$
: pression, température, gradient de température, vitesse du vent.
- 3. **input_tz** : zone horaire pour ts_data (défaut : "GMT").
- 4. **date_format** : format des heures dans ts_data (défaut : "%d/%m/%Y %H:%M", format
-du fichier transmis par Michel).
- 5. **working_tz** : zone horaire dans laquelle on souhaite travailler avec les données
-(défaut : "GMT").
- 6. **predict_at** : heure à laquelle s'effectue la prévision $-$ et donc dernière heure
-d'un bloc de 24h, relativement à working_tz. data`$`getSerie(3) renvoit ainsi les 24
-valeurs de 8h à 7h pour le $3^{eme}$ bloc de 24h présent dans le jeu de données.
+ 3. **date_format** : format des heures dans ts_data (défaut : "%d/%m/%Y %H:%M", format
+du fichier transmis par Michel Bobbia).
+ 4. **limit** : nombre de séries à récupérer.
-----r
print(data)
#?Data
?computeForecast
5. **memory** : le nombre de jours à prendre en compte dans le passé pour chaque
prévision (par défaut : Inf, c'est-à-dire tout l'historique pris en compte).
- 6. **horizon** : le nombre d'heures à prédire ; par défaut "data`$`getStdHorizon()",
-c'est-à-dire le nombre d'heures restantes à partir de l'instant de prévision + 1 jusqu'à
-minuit (17 pour predict_at=7 par exemple).
- 7. **ncores** : le nombre de processus parallèles (utiliser 1 pour une exécution
+ 6. **predict_from** : première heure de prévision. Les séries prévues démarrent
+cependant toutes à 1h du matin (en reprenant les premières valeurs connues).
+ 7. **horizon** : dernière heure de prévision ; maximum 24 == minuit (valeur par défaut).
+pred`$`getForecast(i) retourne une journée complète de 01:00 à 00:00 si horizon=24.
+ 8. **ncores** : le nombre de processus parallèles (utiliser 1 pour une exécution
séquentielle)
-----r
print(pred)
1. **data** : le jeu de données renvoyé par getData()
2. **pred** : les prédictions renvoyées par computeForecast()
- 3. **horizon** : le nombre d'heures à considérer pour le calcul de l'erreur ; doit être
-inférieur ou égal à l'horizon utilisé pour la prédiction (même valeur par défaut :
-"data`$`getStdHorizon()")
+ 3. **predict_from** : première heure de prévision ; peut être différente de l'analogue
+dans l'appel à *computeForecast()*.
+ 4. **horizon** : dernière heure de prévision à considérer pour le calcul de l'erreur ;
+inférieure ou égale à la valeur de l'argument analogue dans *computeForecast()*
-----r
summary(err)
summary(err$abs)
+++ /dev/null
-
-% Default to the notebook output style
-
-
-
-
-% Inherit from the specified cell style.
-
-
-
-
-
-\documentclass[11pt]{article}
-
-
-
- \usepackage[T1]{fontenc}
- % Nicer default font (+ math font) than Computer Modern for most use cases
- \usepackage{mathpazo}
-
- % Basic figure setup, for now with no caption control since it's done
- % automatically by Pandoc (which extracts  syntax from Markdown).
- \usepackage{graphicx}
- % We will generate all images so they have a width \maxwidth. This means
- % that they will get their normal width if they fit onto the page, but
- % are scaled down if they would overflow the margins.
- \makeatletter
- \def\maxwidth{\ifdim\Gin@nat@width>\linewidth\linewidth
- \else\Gin@nat@width\fi}
- \makeatother
- \let\Oldincludegraphics\includegraphics
- % Set max figure width to be 80% of text width, for now hardcoded.
- \renewcommand{\includegraphics}[1]{\Oldincludegraphics[width=.8\maxwidth]{#1}}
- % Ensure that by default, figures have no caption (until we provide a
- % proper Figure object with a Caption API and a way to capture that
- % in the conversion process - todo).
- \usepackage{caption}
- \DeclareCaptionLabelFormat{nolabel}{}
- \captionsetup{labelformat=nolabel}
-
- \usepackage{adjustbox} % Used to constrain images to a maximum size
- \usepackage{xcolor} % Allow colors to be defined
- \usepackage{enumerate} % Needed for markdown enumerations to work
- \usepackage{geometry} % Used to adjust the document margins
- \usepackage{amsmath} % Equations
- \usepackage{amssymb} % Equations
- \usepackage{textcomp} % defines textquotesingle
- % Hack from http://tex.stackexchange.com/a/47451/13684:
- \AtBeginDocument{%
- \def\PYZsq{\textquotesingle}% Upright quotes in Pygmentized code
- }
- \usepackage{upquote} % Upright quotes for verbatim code
- \usepackage{eurosym} % defines \euro
- \usepackage[mathletters]{ucs} % Extended unicode (utf-8) support
- \usepackage[utf8x]{inputenc} % Allow utf-8 characters in the tex document
- \usepackage{fancyvrb} % verbatim replacement that allows latex
- \usepackage{grffile} % extends the file name processing of package graphics
- % to support a larger range
- % The hyperref package gives us a pdf with properly built
- % internal navigation ('pdf bookmarks' for the table of contents,
- % internal cross-reference links, web links for URLs, etc.)
- \usepackage{hyperref}
- \usepackage{longtable} % longtable support required by pandoc >1.10
- \usepackage{booktabs} % table support for pandoc > 1.12.2
- \usepackage[inline]{enumitem} % IRkernel/repr support (it uses the enumerate* environment)
- \usepackage[normalem]{ulem} % ulem is needed to support strikethroughs (\sout)
- % normalem makes italics be italics, not underlines
-
-
-
-
- % Colors for the hyperref package
- \definecolor{urlcolor}{rgb}{0,.145,.698}
- \definecolor{linkcolor}{rgb}{.71,0.21,0.01}
- \definecolor{citecolor}{rgb}{.12,.54,.11}
-
- % ANSI colors
- \definecolor{ansi-black}{HTML}{3E424D}
- \definecolor{ansi-black-intense}{HTML}{282C36}
- \definecolor{ansi-red}{HTML}{E75C58}
- \definecolor{ansi-red-intense}{HTML}{B22B31}
- \definecolor{ansi-green}{HTML}{00A250}
- \definecolor{ansi-green-intense}{HTML}{007427}
- \definecolor{ansi-yellow}{HTML}{DDB62B}
- \definecolor{ansi-yellow-intense}{HTML}{B27D12}
- \definecolor{ansi-blue}{HTML}{208FFB}
- \definecolor{ansi-blue-intense}{HTML}{0065CA}
- \definecolor{ansi-magenta}{HTML}{D160C4}
- \definecolor{ansi-magenta-intense}{HTML}{A03196}
- \definecolor{ansi-cyan}{HTML}{60C6C8}
- \definecolor{ansi-cyan-intense}{HTML}{258F8F}
- \definecolor{ansi-white}{HTML}{C5C1B4}
- \definecolor{ansi-white-intense}{HTML}{A1A6B2}
-
- % commands and environments needed by pandoc snippets
- % extracted from the output of `pandoc -s`
- \providecommand{\tightlist}{%
- \setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}
- \DefineVerbatimEnvironment{Highlighting}{Verbatim}{commandchars=\\\{\}}
- % Add ',fontsize=\small' for more characters per line
- \newenvironment{Shaded}{}{}
- \newcommand{\KeywordTok}[1]{\textcolor[rgb]{0.00,0.44,0.13}{\textbf{{#1}}}}
- \newcommand{\DataTypeTok}[1]{\textcolor[rgb]{0.56,0.13,0.00}{{#1}}}
- \newcommand{\DecValTok}[1]{\textcolor[rgb]{0.25,0.63,0.44}{{#1}}}
- \newcommand{\BaseNTok}[1]{\textcolor[rgb]{0.25,0.63,0.44}{{#1}}}
- \newcommand{\FloatTok}[1]{\textcolor[rgb]{0.25,0.63,0.44}{{#1}}}
- \newcommand{\CharTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{{#1}}}
- \newcommand{\StringTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{{#1}}}
- \newcommand{\CommentTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textit{{#1}}}}
- \newcommand{\OtherTok}[1]{\textcolor[rgb]{0.00,0.44,0.13}{{#1}}}
- \newcommand{\AlertTok}[1]{\textcolor[rgb]{1.00,0.00,0.00}{\textbf{{#1}}}}
- \newcommand{\FunctionTok}[1]{\textcolor[rgb]{0.02,0.16,0.49}{{#1}}}
- \newcommand{\RegionMarkerTok}[1]{{#1}}
- \newcommand{\ErrorTok}[1]{\textcolor[rgb]{1.00,0.00,0.00}{\textbf{{#1}}}}
- \newcommand{\NormalTok}[1]{{#1}}
-
- % Additional commands for more recent versions of Pandoc
- \newcommand{\ConstantTok}[1]{\textcolor[rgb]{0.53,0.00,0.00}{{#1}}}
- \newcommand{\SpecialCharTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{{#1}}}
- \newcommand{\VerbatimStringTok}[1]{\textcolor[rgb]{0.25,0.44,0.63}{{#1}}}
- \newcommand{\SpecialStringTok}[1]{\textcolor[rgb]{0.73,0.40,0.53}{{#1}}}
- \newcommand{\ImportTok}[1]{{#1}}
- \newcommand{\DocumentationTok}[1]{\textcolor[rgb]{0.73,0.13,0.13}{\textit{{#1}}}}
- \newcommand{\AnnotationTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{{#1}}}}}
- \newcommand{\CommentVarTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{{#1}}}}}
- \newcommand{\VariableTok}[1]{\textcolor[rgb]{0.10,0.09,0.49}{{#1}}}
- \newcommand{\ControlFlowTok}[1]{\textcolor[rgb]{0.00,0.44,0.13}{\textbf{{#1}}}}
- \newcommand{\OperatorTok}[1]{\textcolor[rgb]{0.40,0.40,0.40}{{#1}}}
- \newcommand{\BuiltInTok}[1]{{#1}}
- \newcommand{\ExtensionTok}[1]{{#1}}
- \newcommand{\PreprocessorTok}[1]{\textcolor[rgb]{0.74,0.48,0.00}{{#1}}}
- \newcommand{\AttributeTok}[1]{\textcolor[rgb]{0.49,0.56,0.16}{{#1}}}
- \newcommand{\InformationTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{{#1}}}}}
- \newcommand{\WarningTok}[1]{\textcolor[rgb]{0.38,0.63,0.69}{\textbf{\textit{{#1}}}}}
-
-
- % Define a nice break command that doesn't care if a line doesn't already
- % exist.
- \def\br{\hspace*{\fill} \\* }
- % Math Jax compatability definitions
- \def\gt{>}
- \def\lt{<}
- % Document parameters
- \title{report\_P7\_H17}
-
-
-
-
- % Pygments definitions
-
-\makeatletter
-\def\PY@reset{\let\PY@it=\relax \let\PY@bf=\relax%
- \let\PY@ul=\relax \let\PY@tc=\relax%
- \let\PY@bc=\relax \let\PY@ff=\relax}
-\def\PY@tok#1{\csname PY@tok@#1\endcsname}
-\def\PY@toks#1+{\ifx\relax#1\empty\else%
- \PY@tok{#1}\expandafter\PY@toks\fi}
-\def\PY@do#1{\PY@bc{\PY@tc{\PY@ul{%
- \PY@it{\PY@bf{\PY@ff{#1}}}}}}}
-\def\PY#1#2{\PY@reset\PY@toks#1+\relax+\PY@do{#2}}
-
-\expandafter\def\csname PY@tok@w\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.73,0.73}{##1}}}
-\expandafter\def\csname PY@tok@c\endcsname{\let\PY@it=\textit\def\PY@tc##1{\textcolor[rgb]{0.25,0.50,0.50}{##1}}}
-\expandafter\def\csname PY@tok@cp\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.74,0.48,0.00}{##1}}}
-\expandafter\def\csname PY@tok@k\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
-\expandafter\def\csname PY@tok@kp\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
-\expandafter\def\csname PY@tok@kt\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.69,0.00,0.25}{##1}}}
-\expandafter\def\csname PY@tok@o\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}}
-\expandafter\def\csname PY@tok@ow\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.67,0.13,1.00}{##1}}}
-\expandafter\def\csname PY@tok@nb\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
-\expandafter\def\csname PY@tok@nf\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,1.00}{##1}}}
-\expandafter\def\csname PY@tok@nc\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,1.00}{##1}}}
-\expandafter\def\csname PY@tok@nn\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,1.00}{##1}}}
-\expandafter\def\csname PY@tok@ne\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.82,0.25,0.23}{##1}}}
-\expandafter\def\csname PY@tok@nv\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.10,0.09,0.49}{##1}}}
-\expandafter\def\csname PY@tok@no\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.53,0.00,0.00}{##1}}}
-\expandafter\def\csname PY@tok@nl\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.63,0.63,0.00}{##1}}}
-\expandafter\def\csname PY@tok@ni\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.60,0.60,0.60}{##1}}}
-\expandafter\def\csname PY@tok@na\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.49,0.56,0.16}{##1}}}
-\expandafter\def\csname PY@tok@nt\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
-\expandafter\def\csname PY@tok@nd\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.67,0.13,1.00}{##1}}}
-\expandafter\def\csname PY@tok@s\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}}
-\expandafter\def\csname PY@tok@sd\endcsname{\let\PY@it=\textit\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}}
-\expandafter\def\csname PY@tok@si\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.73,0.40,0.53}{##1}}}
-\expandafter\def\csname PY@tok@se\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.73,0.40,0.13}{##1}}}
-\expandafter\def\csname PY@tok@sr\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.40,0.53}{##1}}}
-\expandafter\def\csname PY@tok@ss\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.10,0.09,0.49}{##1}}}
-\expandafter\def\csname PY@tok@sx\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
-\expandafter\def\csname PY@tok@m\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}}
-\expandafter\def\csname PY@tok@gh\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,0.50}{##1}}}
-\expandafter\def\csname PY@tok@gu\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.50,0.00,0.50}{##1}}}
-\expandafter\def\csname PY@tok@gd\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.63,0.00,0.00}{##1}}}
-\expandafter\def\csname PY@tok@gi\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.63,0.00}{##1}}}
-\expandafter\def\csname PY@tok@gr\endcsname{\def\PY@tc##1{\textcolor[rgb]{1.00,0.00,0.00}{##1}}}
-\expandafter\def\csname PY@tok@ge\endcsname{\let\PY@it=\textit}
-\expandafter\def\csname PY@tok@gs\endcsname{\let\PY@bf=\textbf}
-\expandafter\def\csname PY@tok@gp\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,0.50}{##1}}}
-\expandafter\def\csname PY@tok@go\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.53,0.53,0.53}{##1}}}
-\expandafter\def\csname PY@tok@gt\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.27,0.87}{##1}}}
-\expandafter\def\csname PY@tok@err\endcsname{\def\PY@bc##1{\setlength{\fboxsep}{0pt}\fcolorbox[rgb]{1.00,0.00,0.00}{1,1,1}{\strut ##1}}}
-\expandafter\def\csname PY@tok@kc\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
-\expandafter\def\csname PY@tok@kd\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
-\expandafter\def\csname PY@tok@kn\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
-\expandafter\def\csname PY@tok@kr\endcsname{\let\PY@bf=\textbf\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
-\expandafter\def\csname PY@tok@bp\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.50,0.00}{##1}}}
-\expandafter\def\csname PY@tok@fm\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.00,0.00,1.00}{##1}}}
-\expandafter\def\csname PY@tok@vc\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.10,0.09,0.49}{##1}}}
-\expandafter\def\csname PY@tok@vg\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.10,0.09,0.49}{##1}}}
-\expandafter\def\csname PY@tok@vi\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.10,0.09,0.49}{##1}}}
-\expandafter\def\csname PY@tok@vm\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.10,0.09,0.49}{##1}}}
-\expandafter\def\csname PY@tok@sa\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}}
-\expandafter\def\csname PY@tok@sb\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}}
-\expandafter\def\csname PY@tok@sc\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}}
-\expandafter\def\csname PY@tok@dl\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}}
-\expandafter\def\csname PY@tok@s2\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}}
-\expandafter\def\csname PY@tok@sh\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}}
-\expandafter\def\csname PY@tok@s1\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.73,0.13,0.13}{##1}}}
-\expandafter\def\csname PY@tok@mb\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}}
-\expandafter\def\csname PY@tok@mf\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}}
-\expandafter\def\csname PY@tok@mh\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}}
-\expandafter\def\csname PY@tok@mi\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}}
-\expandafter\def\csname PY@tok@il\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}}
-\expandafter\def\csname PY@tok@mo\endcsname{\def\PY@tc##1{\textcolor[rgb]{0.40,0.40,0.40}{##1}}}
-\expandafter\def\csname PY@tok@ch\endcsname{\let\PY@it=\textit\def\PY@tc##1{\textcolor[rgb]{0.25,0.50,0.50}{##1}}}
-\expandafter\def\csname PY@tok@cm\endcsname{\let\PY@it=\textit\def\PY@tc##1{\textcolor[rgb]{0.25,0.50,0.50}{##1}}}
-\expandafter\def\csname PY@tok@cpf\endcsname{\let\PY@it=\textit\def\PY@tc##1{\textcolor[rgb]{0.25,0.50,0.50}{##1}}}
-\expandafter\def\csname PY@tok@c1\endcsname{\let\PY@it=\textit\def\PY@tc##1{\textcolor[rgb]{0.25,0.50,0.50}{##1}}}
-\expandafter\def\csname PY@tok@cs\endcsname{\let\PY@it=\textit\def\PY@tc##1{\textcolor[rgb]{0.25,0.50,0.50}{##1}}}
-
-\def\PYZbs{\char`\\}
-\def\PYZus{\char`\_}
-\def\PYZob{\char`\{}
-\def\PYZcb{\char`\}}
-\def\PYZca{\char`\^}
-\def\PYZam{\char`\&}
-\def\PYZlt{\char`\<}
-\def\PYZgt{\char`\>}
-\def\PYZsh{\char`\#}
-\def\PYZpc{\char`\%}
-\def\PYZdl{\char`\$}
-\def\PYZhy{\char`\-}
-\def\PYZsq{\char`\'}
-\def\PYZdq{\char`\"}
-\def\PYZti{\char`\~}
-% for compatibility with earlier versions
-\def\PYZat{@}
-\def\PYZlb{[}
-\def\PYZrb{]}
-\makeatother
-
-
- % Exact colors from NB
- \definecolor{incolor}{rgb}{0.0, 0.0, 0.5}
- \definecolor{outcolor}{rgb}{0.545, 0.0, 0.0}
-
-
-
-
- % Prevent overflowing lines due to hard-to-break entities
- \sloppy
- % Setup hyperref package
- \hypersetup{
- breaklinks=true, % so long urls are correctly broken across lines
- colorlinks=true,
- urlcolor=urlcolor,
- linkcolor=linkcolor,
- citecolor=citecolor,
- }
- % Slightly bigger margins than the latex defaults
-
- \geometry{verbose,tmargin=1in,bmargin=1in,lmargin=1in,rmargin=1in}
-
-
-\graphicspath{{./figs/}}
-
- \begin{document}
-
-
- \maketitle
-
-
-
-
- \subsection*{Introduction}
-
-Cette partie montre les résultats obtenus via des variantes de
-l'algorithme décrit à la section 2, en utilisant le package présenté à
-la section 3. Cet algorithme est systématiquement comparé à deux
-approches naïves : * la moyenne des lendemains des jours "similaires"
-dans tout le passé, c'est-à-dire prédiction = moyenne de tous les mardis
-passé si le jour courant est un lndi par exemple. * la persistence,
-reproduisant le jour courant ou allant chercher le lendemain de la
-dernière journée "similaire" (même principe que ci-dessus ; argument
-"same\_day").
-
-Concernant l'algorithme principal à voisins, trois variantes sont
-étudiées dans cette partie : * avec simtype="mix" et raccordement
-"Neighbors" dans le cas "non local", i.e. on va chercher des voisins
-n'importe où du moment qu'ils correspondent au premier élément d'un
-couple de deux jours consécutifs sans valeurs manquantes. * avec
-simtype="endo" + raccordement "Neighbors" puis simtype="none" +
-raccordement "Zero" (sans ajustement) dans le cas "local" : voisins de
-même niveau de pollution et même saison.
-
-Pour chaque période retenue -\/- chauffage, épandage, semaine non
-polluée -\/- les erreurs de prédiction sont d'abord affichées, puis
-quelques graphes de courbes réalisées/prévues (sur le jour "en moyenne
-le plus facile" à gauche, et "en moyenne le plus difficile" à droite).
-Ensuite plusieurs types de graphes apportant des précisions sur la
-nature et la difficulté du problème viennent compléter ces premières
-courbes. Concernant les graphes de filaments, la moitié gauche du graphe
-correspond aux jours similaires au jour courant, tandis que la moitié
-droite affiche les lendemains : ce sont donc les voisinages tels
-qu'utilisés dans l'algorithme.
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}1}]:} \PY{n}{library}\PY{p}{(}\PY{n}{talweg}\PY{p}{)}
-
- \PY{n}{P} \PY{o}{=} \PY{l+m+mi}{7} \PY{c+c1}{\PYZsh{}instant de prévision}
- \PY{n}{H} \PY{o}{=} \PY{l+m+mi}{17} \PY{c+c1}{\PYZsh{}horizon (en heures)}
-
- \PY{n}{ts\PYZus{}data} \PY{o}{=} \PY{n}{read}\PY{o}{.}\PY{n}{csv}\PY{p}{(}\PY{n}{system}\PY{o}{.}\PY{n}{file}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{extdata}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{pm10\PYZus{}mesures\PYZus{}H\PYZus{}loc\PYZus{}report.csv}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}
- \PY{n}{package}\PY{o}{=}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{talweg}\PY{l+s+s2}{\PYZdq{}}\PY{p}{)}\PY{p}{)}
- \PY{n}{exo\PYZus{}data} \PY{o}{=} \PY{n}{read}\PY{o}{.}\PY{n}{csv}\PY{p}{(}\PY{n}{system}\PY{o}{.}\PY{n}{file}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{extdata}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{meteo\PYZus{}extra\PYZus{}noNAs.csv}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{package}\PY{o}{=}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{talweg}\PY{l+s+s2}{\PYZdq{}}\PY{p}{)}\PY{p}{)}
- \PY{c+c1}{\PYZsh{} NOTE: \PYZsq{}GMT\PYZsq{} because DST gaps are filled and multiple values merged in above dataset.}
- \PY{c+c1}{\PYZsh{} Prediction from P+1 to P+H included.}
- \PY{n}{data} \PY{o}{=} \PY{n}{getData}\PY{p}{(}\PY{n}{ts\PYZus{}data}\PY{p}{,} \PY{n}{exo\PYZus{}data}\PY{p}{,} \PY{n}{input\PYZus{}tz} \PY{o}{=} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{GMT}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{working\PYZus{}tz}\PY{o}{=}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{GMT}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{predict\PYZus{}at}\PY{o}{=}\PY{n}{P}\PY{p}{)}
-
- \PY{n}{indices\PYZus{}ch} \PY{o}{=} \PY{n}{seq}\PY{p}{(}\PY{k}{as}\PY{o}{.}\PY{n}{Date}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{2015\PYZhy{}01\PYZhy{}18}\PY{l+s+s2}{\PYZdq{}}\PY{p}{)}\PY{p}{,}\PY{k}{as}\PY{o}{.}\PY{n}{Date}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{2015\PYZhy{}01\PYZhy{}24}\PY{l+s+s2}{\PYZdq{}}\PY{p}{)}\PY{p}{,}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{days}\PY{l+s+s2}{\PYZdq{}}\PY{p}{)}
- \PY{n}{indices\PYZus{}ep} \PY{o}{=} \PY{n}{seq}\PY{p}{(}\PY{k}{as}\PY{o}{.}\PY{n}{Date}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{2015\PYZhy{}03\PYZhy{}15}\PY{l+s+s2}{\PYZdq{}}\PY{p}{)}\PY{p}{,}\PY{k}{as}\PY{o}{.}\PY{n}{Date}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{2015\PYZhy{}03\PYZhy{}21}\PY{l+s+s2}{\PYZdq{}}\PY{p}{)}\PY{p}{,}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{days}\PY{l+s+s2}{\PYZdq{}}\PY{p}{)}
- \PY{n}{indices\PYZus{}np} \PY{o}{=} \PY{n}{seq}\PY{p}{(}\PY{k}{as}\PY{o}{.}\PY{n}{Date}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{2015\PYZhy{}04\PYZhy{}26}\PY{l+s+s2}{\PYZdq{}}\PY{p}{)}\PY{p}{,}\PY{k}{as}\PY{o}{.}\PY{n}{Date}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{2015\PYZhy{}05\PYZhy{}02}\PY{l+s+s2}{\PYZdq{}}\PY{p}{)}\PY{p}{,}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{days}\PY{l+s+s2}{\PYZdq{}}\PY{p}{)}
-\end{Verbatim}
-
- \subsection*{Pollution par chauffage}
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}2}]:} \PY{n}{p1} \PY{o}{=} \PY{n}{computeForecast}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{indices\PYZus{}ch}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Neighbors}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Neighbors}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{horizon}\PY{o}{=}\PY{n}{H}\PY{p}{,}
- \PY{n}{simtype}\PY{o}{=}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{mix}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{local}\PY{o}{=}\PY{n}{FALSE}\PY{p}{)}
- \PY{n}{p2} \PY{o}{=} \PY{n}{computeForecast}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{indices\PYZus{}ch}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Neighbors}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Neighbors}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{horizon}\PY{o}{=}\PY{n}{H}\PY{p}{,}
- \PY{n}{simtype}\PY{o}{=}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{endo}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{local}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
- \PY{n}{p3} \PY{o}{=} \PY{n}{computeForecast}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{indices\PYZus{}ch}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Neighbors}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Zero}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{horizon}\PY{o}{=}\PY{n}{H}\PY{p}{,}
- \PY{n}{simtype}\PY{o}{=}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{none}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{local}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
- \PY{n}{p4} \PY{o}{=} \PY{n}{computeForecast}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{indices\PYZus{}ch}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Average}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Zero}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{horizon}\PY{o}{=}\PY{n}{H}\PY{p}{)}
- \PY{n}{p5} \PY{o}{=} \PY{n}{computeForecast}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{indices\PYZus{}ch}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Persistence}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Zero}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{horizon}\PY{o}{=}\PY{n}{H}\PY{p}{,}
- \PY{n}{same\PYZus{}day}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
-\end{Verbatim}
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}3}]:} \PY{n}{e1} \PY{o}{=} \PY{n}{computeError}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p1}\PY{p}{,} \PY{n}{H}\PY{p}{)}
- \PY{n}{e2} \PY{o}{=} \PY{n}{computeError}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p2}\PY{p}{,} \PY{n}{H}\PY{p}{)}
- \PY{n}{e3} \PY{o}{=} \PY{n}{computeError}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p3}\PY{p}{,} \PY{n}{H}\PY{p}{)}
- \PY{n}{e4} \PY{o}{=} \PY{n}{computeError}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p4}\PY{p}{,} \PY{n}{H}\PY{p}{)}
- \PY{n}{e5} \PY{o}{=} \PY{n}{computeError}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p5}\PY{p}{,} \PY{n}{H}\PY{p}{)}
- \PY{n}{options}\PY{p}{(}\PY{n+nb}{repr}\PY{o}{.}\PY{n}{plot}\PY{o}{.}\PY{n}{width}\PY{o}{=}\PY{l+m+mi}{9}\PY{p}{,} \PY{n+nb}{repr}\PY{o}{.}\PY{n}{plot}\PY{o}{.}\PY{n}{height}\PY{o}{=}\PY{l+m+mi}{7}\PY{p}{)}
- \PY{n}{plotError}\PY{p}{(}\PY{n+nb}{list}\PY{p}{(}\PY{n}{e1}\PY{p}{,} \PY{n}{e5}\PY{p}{,} \PY{n}{e4}\PY{p}{,} \PY{n}{e2}\PY{p}{,} \PY{n}{e3}\PY{p}{)}\PY{p}{,} \PY{n}{cols}\PY{o}{=}\PY{n}{c}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{,}\PY{n}{colors}\PY{p}{(}\PY{p}{)}\PY{p}{[}\PY{l+m+mi}{258}\PY{p}{]}\PY{p}{,}\PY{l+m+mi}{4}\PY{p}{,}\PY{l+m+mi}{6}\PY{p}{)}\PY{p}{)}
-
- \PY{c+c1}{\PYZsh{} noir: Neighbors non\PYZhy{}local (p1), bleu: Neighbors local endo (p2),}
- \PY{c+c1}{\PYZsh{} mauve: Neighbors local none (p3), vert: moyenne (p4),}
- \PY{c+c1}{\PYZsh{} rouge: persistence (p5)}
-
- \PY{n}{sum\PYZus{}p123} \PY{o}{=} \PY{n}{e1}\PY{err}{\PYZdl{}}\PY{n+nb}{abs}\PY{err}{\PYZdl{}}\PY{n}{indices} \PY{o}{+} \PY{n}{e2}\PY{err}{\PYZdl{}}\PY{n+nb}{abs}\PY{err}{\PYZdl{}}\PY{n}{indices} \PY{o}{+} \PY{n}{e3}\PY{err}{\PYZdl{}}\PY{n+nb}{abs}\PY{err}{\PYZdl{}}\PY{n}{indices}
- \PY{n}{i\PYZus{}np} \PY{o}{=} \PY{n}{which}\PY{o}{.}\PY{n}{min}\PY{p}{(}\PY{n}{sum\PYZus{}p123}\PY{p}{)} \PY{c+c1}{\PYZsh{}indice de (veille de) jour \PYZdq{}facile\PYZdq{}}
- \PY{n}{i\PYZus{}p} \PY{o}{=} \PY{n}{which}\PY{o}{.}\PY{n}{max}\PY{p}{(}\PY{n}{sum\PYZus{}p123}\PY{p}{)} \PY{c+c1}{\PYZsh{}indice de (veille de) jour \PYZdq{}difficile\PYZdq{}}
-\end{Verbatim}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_4_0.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- L'erreur absolue dépasse 20 sur 1 à 2 jours suivant les modèles (graphe
-en haut à droite). C'est au-delà de ce que l'on aimerait voir (disons
-+/- 5 environ). Sur cet exemple le modèle à voisins "contraint"
-(local=TRUE) utilisant des pondérations basées sur les similarités de
-forme (simtype="endo") obtient en moyenne les meilleurs résultats, avec
-un MAPE restant en général inférieur à 30\% de 8h à 19h (7+1 à 7+12 :
-graphe en bas à gauche).
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}4}]:} \PY{n}{options}\PY{p}{(}\PY{n+nb}{repr}\PY{o}{.}\PY{n}{plot}\PY{o}{.}\PY{n}{width}\PY{o}{=}\PY{l+m+mi}{9}\PY{p}{,} \PY{n+nb}{repr}\PY{o}{.}\PY{n}{plot}\PY{o}{.}\PY{n}{height}\PY{o}{=}\PY{l+m+mi}{4}\PY{p}{)}
- \PY{n}{par}\PY{p}{(}\PY{n}{mfrow}\PY{o}{=}\PY{n}{c}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{)}\PY{p}{)}
-
- \PY{n}{plotPredReal}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p1}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{PredReal p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotPredReal}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p1}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{PredReal p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{n}{plotPredReal}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p2}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{PredReal p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotPredReal}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p2}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{PredReal p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{n}{plotPredReal}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p3}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{PredReal p3 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotPredReal}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p3}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{PredReal p3 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{c+c1}{\PYZsh{} Bleu: prévue, noir: réalisée}
-\end{Verbatim}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_6_0.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_6_1.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_6_2.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- Le jour "facile à prévoir", à gauche, se décompose en deux modes : un
-léger vers 10h (7+3), puis un beaucoup plus marqué vers 19h (7+12). Ces
-deux modes sont retrouvés par les trois variantes de l'algorithme à
-voisins, bien que l'amplitude soit mal prédite. Concernant le jour
-"difficile à prévoir" il y a deux pics en tout début et toute fin de
-journée (à 9h et 23h), qui ne sont pas du tout anticipés par le
-programme ; la grande amplitude de ces pics explique alors l'intensité
-de l'erreur observée.
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}5}]:} \PY{n}{par}\PY{p}{(}\PY{n}{mfrow}\PY{o}{=}\PY{n}{c}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{)}\PY{p}{)}
- \PY{n}{f\PYZus{}np1} \PY{o}{=} \PY{n}{computeFilaments}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p1}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{,} \PY{n}{plot}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
- \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Filaments p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{f\PYZus{}p1} \PY{o}{=} \PY{n}{computeFilaments}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p1}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{,} \PY{n}{plot}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
- \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Filaments p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{n}{f\PYZus{}np2} \PY{o}{=} \PY{n}{computeFilaments}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p2}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{,} \PY{n}{plot}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
- \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Filaments p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{f\PYZus{}p2} \PY{o}{=} \PY{n}{computeFilaments}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p2}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{,} \PY{n}{plot}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
- \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Filaments p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-\end{Verbatim}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_8_0.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_8_1.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- Les voisins du jour courant (période de 24h allant de 8h à 7h le
-lendemain) sont affichés avec un trait d'autant plus sombre qu'ils sont
-proches. On constate dans le cas non contraint (en haut) une grande
-variabilité des lendemains, très nette sur le graphe en haut à droite.
-Ceci indique une faible corrélation entre la forme d'une courbe sur une
-période de 24h et la forme sur les 24h suivantes ; \textbf{cette
-observation est la source des difficultés rencontrées par l'algorithme
-sur ce jeu de données.}
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}6}]:} \PY{n}{par}\PY{p}{(}\PY{n}{mfrow}\PY{o}{=}\PY{n}{c}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotFilamentsBox}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{f\PYZus{}np1}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{FilBox p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotFilamentsBox}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{f\PYZus{}p1}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{FilBox p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{c+c1}{\PYZsh{} En pointillés la courbe du jour courant + lendemain (à prédire)}
-\end{Verbatim}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_10_0.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- Sur cette boxplot fonctionnelle (voir la fonction fboxplot() du package
-R "rainbow") l'on constate essentiellement deux choses : le lendemain
-d'un voisin "normal" peut se révéler être une courbe atypique, fort
-éloignée de ce que l'on souhaite prédire (courbes bleue et rouge à
-gauche) ; et, dans le cas d'une courbe à prédire atypique (à droite) la
-plupart des voisins sont trop éloignés de la forme à prédire et forcent
-ainsi un aplatissement de la prédiction.
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}7}]:} \PY{n}{par}\PY{p}{(}\PY{n}{mfrow}\PY{o}{=}\PY{n}{c}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotRelVar}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{f\PYZus{}np1}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{StdDev p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotRelVar}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{f\PYZus{}p1}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{StdDev p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{n}{plotRelVar}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{f\PYZus{}np2}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{StdDev p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotRelVar}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{f\PYZus{}p2}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{StdDev p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{c+c1}{\PYZsh{} Variabilité globale en rouge ; sur les 60 voisins (+ lendemains) en noir}
-\end{Verbatim}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_12_0.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_12_1.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- Ces graphes viennent confirmer l'impression visuelle après observation
-des filaments. En effet, la variabilité globale en rouge (écart-type
-heure par heure sur l'ensemble des couples "aujourd'hui/lendemain" du
-passé) devrait rester nettement au-dessus de la variabilité locale,
-calculée respectivement sur un voisinage d'une soixantaine de jours
-(pour p1) et d'une dizaine de jours (pour p2). Or on constate que ce
-n'est pas du tout le cas sur la période "lendemain", sauf en partie pour
-p2 le jour 4 \(-\) mais ce n'est pas suffisant.
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}8}]:} \PY{n}{par}\PY{p}{(}\PY{n}{mfrow}\PY{o}{=}\PY{n}{c}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotSimils}\PY{p}{(}\PY{n}{p1}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Weights p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotSimils}\PY{p}{(}\PY{n}{p1}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Weights p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{n}{plotSimils}\PY{p}{(}\PY{n}{p2}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Weights p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotSimils}\PY{p}{(}\PY{n}{p2}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Weights p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{c+c1}{\PYZsh{} \PYZhy{} pollué à gauche, + pollué à droite}
-\end{Verbatim}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_14_0.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_14_1.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- Les poids se concentrent près de 0 dans le cas "non local" (p1), et se
-répartissent assez uniformément dans \([0,0.2]\) dans le cas "local"
-(p2). C'est ce que l'on souhaite observer pour éviter d'effectuer une
-simple moyenne.
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}9}]:} \PY{c+c1}{\PYZsh{} Fenêtres sélectionnées dans ]0,7] / non\PYZhy{}loc 2 premières lignes, loc ensuite}
- \PY{n}{p1}\PY{err}{\PYZdl{}}\PY{n}{getParams}\PY{p}{(}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{err}{\PYZdl{}}\PY{n}{window}
- \PY{n}{p1}\PY{err}{\PYZdl{}}\PY{n}{getParams}\PY{p}{(}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{err}{\PYZdl{}}\PY{n}{window}
-
- \PY{n}{p2}\PY{err}{\PYZdl{}}\PY{n}{getParams}\PY{p}{(}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{err}{\PYZdl{}}\PY{n}{window}
- \PY{n}{p2}\PY{err}{\PYZdl{}}\PY{n}{getParams}\PY{p}{(}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{err}{\PYZdl{}}\PY{n}{window}
-\end{Verbatim}
-
- \begin{enumerate*}
-\item 0.168824188864717
-\item 0.336969608767438
-\end{enumerate*}
-
-
- \begin{enumerate*}
-\item 0.18004595760919
-\item 0.353963007643311
-\end{enumerate*}
-
-
- 1.16620655388085
-
-
- 1.18148881881259
-
-
- \subsection*{Pollution par épandage}
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}10}]:} \PY{n}{p1} \PY{o}{=} \PY{n}{computeForecast}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{indices\PYZus{}ep}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Neighbors}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Neighbors}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{horizon}\PY{o}{=}\PY{n}{H}\PY{p}{,}
- \PY{n}{simtype}\PY{o}{=}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{mix}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{local}\PY{o}{=}\PY{n}{FALSE}\PY{p}{)}
- \PY{n}{p2} \PY{o}{=} \PY{n}{computeForecast}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{indices\PYZus{}ep}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Neighbors}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Neighbors}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{horizon}\PY{o}{=}\PY{n}{H}\PY{p}{,}
- \PY{n}{simtype}\PY{o}{=}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{endo}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{local}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
- \PY{n}{p3} \PY{o}{=} \PY{n}{computeForecast}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{indices\PYZus{}ep}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Neighbors}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Zero}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{horizon}\PY{o}{=}\PY{n}{H}\PY{p}{,}
- \PY{n}{simtype}\PY{o}{=}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{none}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{local}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
- \PY{n}{p4} \PY{o}{=} \PY{n}{computeForecast}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{indices\PYZus{}ep}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Average}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Zero}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{horizon}\PY{o}{=}\PY{n}{H}\PY{p}{)}
- \PY{n}{p5} \PY{o}{=} \PY{n}{computeForecast}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{indices\PYZus{}ep}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Persistence}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Zero}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{horizon}\PY{o}{=}\PY{n}{H}\PY{p}{,}
- \PY{n}{same\PYZus{}day}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
-\end{Verbatim}
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}11}]:} \PY{n}{e1} \PY{o}{=} \PY{n}{computeError}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p1}\PY{p}{,} \PY{n}{H}\PY{p}{)}
- \PY{n}{e2} \PY{o}{=} \PY{n}{computeError}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p2}\PY{p}{,} \PY{n}{H}\PY{p}{)}
- \PY{n}{e3} \PY{o}{=} \PY{n}{computeError}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p3}\PY{p}{,} \PY{n}{H}\PY{p}{)}
- \PY{n}{e4} \PY{o}{=} \PY{n}{computeError}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p4}\PY{p}{,} \PY{n}{H}\PY{p}{)}
- \PY{n}{e5} \PY{o}{=} \PY{n}{computeError}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p5}\PY{p}{,} \PY{n}{H}\PY{p}{)}
- \PY{n}{options}\PY{p}{(}\PY{n+nb}{repr}\PY{o}{.}\PY{n}{plot}\PY{o}{.}\PY{n}{width}\PY{o}{=}\PY{l+m+mi}{9}\PY{p}{,} \PY{n+nb}{repr}\PY{o}{.}\PY{n}{plot}\PY{o}{.}\PY{n}{height}\PY{o}{=}\PY{l+m+mi}{7}\PY{p}{)}
- \PY{n}{plotError}\PY{p}{(}\PY{n+nb}{list}\PY{p}{(}\PY{n}{e1}\PY{p}{,} \PY{n}{e5}\PY{p}{,} \PY{n}{e4}\PY{p}{,} \PY{n}{e2}\PY{p}{,} \PY{n}{e3}\PY{p}{)}\PY{p}{,} \PY{n}{cols}\PY{o}{=}\PY{n}{c}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{,}\PY{n}{colors}\PY{p}{(}\PY{p}{)}\PY{p}{[}\PY{l+m+mi}{258}\PY{p}{]}\PY{p}{,}\PY{l+m+mi}{4}\PY{p}{,}\PY{l+m+mi}{6}\PY{p}{)}\PY{p}{)}
-
- \PY{c+c1}{\PYZsh{} noir: Neighbors non\PYZhy{}local (p1), bleu: Neighbors local endo (p2),}
- \PY{c+c1}{\PYZsh{} mauve: Neighbors local none (p3), vert: moyenne (p4),}
- \PY{c+c1}{\PYZsh{} rouge: persistence (p5)}
-
- \PY{n}{sum\PYZus{}p123} \PY{o}{=} \PY{n}{e1}\PY{err}{\PYZdl{}}\PY{n+nb}{abs}\PY{err}{\PYZdl{}}\PY{n}{indices} \PY{o}{+} \PY{n}{e2}\PY{err}{\PYZdl{}}\PY{n+nb}{abs}\PY{err}{\PYZdl{}}\PY{n}{indices} \PY{o}{+} \PY{n}{e3}\PY{err}{\PYZdl{}}\PY{n+nb}{abs}\PY{err}{\PYZdl{}}\PY{n}{indices}
- \PY{n}{i\PYZus{}np} \PY{o}{=} \PY{n}{which}\PY{o}{.}\PY{n}{min}\PY{p}{(}\PY{n}{sum\PYZus{}p123}\PY{p}{)} \PY{c+c1}{\PYZsh{}indice de (veille de) jour \PYZdq{}facile\PYZdq{}}
- \PY{n}{i\PYZus{}p} \PY{o}{=} \PY{n}{which}\PY{o}{.}\PY{n}{max}\PY{p}{(}\PY{n}{sum\PYZus{}p123}\PY{p}{)} \PY{c+c1}{\PYZsh{}indice de (veille de) jour \PYZdq{}difficile\PYZdq{}}
-\end{Verbatim}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_19_0.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- Il est difficile dans ce cas de déterminer une méthode meilleure que les
-autres : elles donnent toutes de plutôt mauvais résultats, avec une
-erreur absolue moyennée sur la journée dépassant presque toujours 15
-(graphe en haut à droite).
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}12}]:} \PY{n}{options}\PY{p}{(}\PY{n+nb}{repr}\PY{o}{.}\PY{n}{plot}\PY{o}{.}\PY{n}{width}\PY{o}{=}\PY{l+m+mi}{9}\PY{p}{,} \PY{n+nb}{repr}\PY{o}{.}\PY{n}{plot}\PY{o}{.}\PY{n}{height}\PY{o}{=}\PY{l+m+mi}{4}\PY{p}{)}
- \PY{n}{par}\PY{p}{(}\PY{n}{mfrow}\PY{o}{=}\PY{n}{c}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{)}\PY{p}{)}
-
- \PY{n}{plotPredReal}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p1}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{PredReal p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotPredReal}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p1}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{PredReal p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{n}{plotPredReal}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p2}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{PredReal p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotPredReal}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p2}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{PredReal p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{n}{plotPredReal}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p3}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{PredReal p3 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotPredReal}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p3}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{PredReal p3 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{c+c1}{\PYZsh{} Bleu: prévue, noir: réalisée}
-\end{Verbatim}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_21_0.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_21_1.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_21_2.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- Dans le cas d'un jour "facile" à prédire \(-\) à gauche \(-\) la forme
-est plus ou moins retrouvée, mais le niveau moyen est trop bas (courbe
-en bleu). Concernant le jour "difficile" à droite, non seulement la
-forme n'est pas anticipée mais surtout le niveau prédit est très
-inférieur au niveau de pollution observé. Comme on le voit ci-dessous
-cela découle d'un manque de voisins au comportement similaire.
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}13}]:} \PY{n}{par}\PY{p}{(}\PY{n}{mfrow}\PY{o}{=}\PY{n}{c}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{)}\PY{p}{)}
- \PY{n}{f\PYZus{}np1} \PY{o}{=} \PY{n}{computeFilaments}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p1}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{,} \PY{n}{plot}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
- \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Filaments p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{f\PYZus{}p1} \PY{o}{=} \PY{n}{computeFilaments}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p1}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{,} \PY{n}{plot}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
- \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Filaments p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{n}{f\PYZus{}np2} \PY{o}{=} \PY{n}{computeFilaments}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p2}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{,} \PY{n}{plot}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
- \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Filaments p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{f\PYZus{}p2} \PY{o}{=} \PY{n}{computeFilaments}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p2}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{,} \PY{n}{plot}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
- \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Filaments p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-\end{Verbatim}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_23_0.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_23_1.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- Les observations sont les mêmes qu'au paragraphe précédent : trop de
-variabilité des lendemains (et même des voisins du jour courant).
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}14}]:} \PY{n}{par}\PY{p}{(}\PY{n}{mfrow}\PY{o}{=}\PY{n}{c}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotFilamentsBox}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{f\PYZus{}np1}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{FilBox p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotFilamentsBox}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{f\PYZus{}p1}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{FilBox p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{c+c1}{\PYZsh{} En pointillés la courbe du jour courant + lendemain (à prédire)}
-\end{Verbatim}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_25_0.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- On constate la présence d'un voisin au lendemain complètement atypique
-avec un pic en début de journée (courbe en vert à gauche), et d'un autre
-phénomène semblable avec la courbe rouge sur le graphe de droite. Ajouté
-au fait que le lendemain à prévoir est lui-même un jour "hors norme",
-cela montre l'impossibilité de bien prévoir une courbe en utilisant
-l'algorithme à voisins.
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}15}]:} \PY{n}{par}\PY{p}{(}\PY{n}{mfrow}\PY{o}{=}\PY{n}{c}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotRelVar}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{f\PYZus{}np1}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{StdDev p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotRelVar}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{f\PYZus{}p1}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{StdDev p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{n}{plotRelVar}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{f\PYZus{}np2}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{StdDev p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotRelVar}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{f\PYZus{}p2}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{StdDev p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{c+c1}{\PYZsh{} Variabilité globale en rouge ; sur les 60 voisins (+ lendemains) en noir}
-\end{Verbatim}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_27_0.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_27_1.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- Comme précédemment les variabilités locales et globales sont confondues
-dans les parties droites des graphes \(-\) sauf pour la version "locale"
-sur le jour "facile" ; mais cette bonne propriété n'est pas suffisante
-si l'on ne trouve pas les bons poids à appliquer.
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}16}]:} \PY{n}{par}\PY{p}{(}\PY{n}{mfrow}\PY{o}{=}\PY{n}{c}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotSimils}\PY{p}{(}\PY{n}{p1}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Weights p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotSimils}\PY{p}{(}\PY{n}{p1}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Weights p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{n}{plotSimils}\PY{p}{(}\PY{n}{p2}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Weights p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotSimils}\PY{p}{(}\PY{n}{p2}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Weights p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{c+c1}{\PYZsh{} \PYZhy{} pollué à gauche, + pollué à droite}
-\end{Verbatim}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_29_0.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_29_1.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- En comparaison avec le pragraphe précédent on retrouve le même (bon)
-comportement des poids pour la version "non locale". En revanche la
-fenêtre optimisée est trop grande sur le jour "facile" pour la méthode
-"locale" (voir affichage ci-dessous) : il en résulte des poids tous
-semblables autour de 0.084, l'algorithme effectue donc une moyenne
-simple \(-\) expliquant pourquoi les courbes mauve et bleue sont très
-proches sur le graphe d'erreurs.
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}17}]:} \PY{c+c1}{\PYZsh{} Fenêtres sélectionnées dans ]0,7] / non\PYZhy{}loc 2 premières lignes, loc ensuite}
- \PY{n}{p1}\PY{err}{\PYZdl{}}\PY{n}{getParams}\PY{p}{(}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{err}{\PYZdl{}}\PY{n}{window}
- \PY{n}{p1}\PY{err}{\PYZdl{}}\PY{n}{getParams}\PY{p}{(}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{err}{\PYZdl{}}\PY{n}{window}
-
- \PY{n}{p2}\PY{err}{\PYZdl{}}\PY{n}{getParams}\PY{p}{(}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{err}{\PYZdl{}}\PY{n}{window}
- \PY{n}{p2}\PY{err}{\PYZdl{}}\PY{n}{getParams}\PY{p}{(}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{err}{\PYZdl{}}\PY{n}{window}
-\end{Verbatim}
-
- \begin{enumerate*}
-\item 0.206604806619633
-\item 0.661854053860987
-\end{enumerate*}
-
-
- \begin{enumerate*}
-\item 0.367945958636072
-\item 0.244429852740092
-\end{enumerate*}
-
-
- 6.99993248587025
-
-
- 1.24825506305085
-
-
- \subsection*{Semaine non polluée}
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}18}]:} \PY{n}{p1} \PY{o}{=} \PY{n}{computeForecast}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{indices\PYZus{}np}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Neighbors}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Neighbors}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{horizon}\PY{o}{=}\PY{n}{H}\PY{p}{,}
- \PY{n}{simtype}\PY{o}{=}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{mix}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{local}\PY{o}{=}\PY{n}{FALSE}\PY{p}{)}
- \PY{n}{p2} \PY{o}{=} \PY{n}{computeForecast}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{indices\PYZus{}np}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Neighbors}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Neighbors}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{horizon}\PY{o}{=}\PY{n}{H}\PY{p}{,}
- \PY{n}{simtype}\PY{o}{=}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{endo}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{local}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
- \PY{n}{p3} \PY{o}{=} \PY{n}{computeForecast}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{indices\PYZus{}np}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Neighbors}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Zero}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{horizon}\PY{o}{=}\PY{n}{H}\PY{p}{,}
- \PY{n}{simtype}\PY{o}{=}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{none}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{local}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
- \PY{n}{p4} \PY{o}{=} \PY{n}{computeForecast}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{indices\PYZus{}np}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Average}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Zero}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{horizon}\PY{o}{=}\PY{n}{H}\PY{p}{)}
- \PY{n}{p5} \PY{o}{=} \PY{n}{computeForecast}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{indices\PYZus{}np}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Persistence}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Zero}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,} \PY{n}{horizon}\PY{o}{=}\PY{n}{H}\PY{p}{,}
- \PY{n}{same\PYZus{}day}\PY{o}{=}\PY{n}{FALSE}\PY{p}{)}
-\end{Verbatim}
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}19}]:} \PY{n}{e1} \PY{o}{=} \PY{n}{computeError}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p1}\PY{p}{,} \PY{n}{H}\PY{p}{)}
- \PY{n}{e2} \PY{o}{=} \PY{n}{computeError}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p2}\PY{p}{,} \PY{n}{H}\PY{p}{)}
- \PY{n}{e3} \PY{o}{=} \PY{n}{computeError}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p3}\PY{p}{,} \PY{n}{H}\PY{p}{)}
- \PY{n}{e4} \PY{o}{=} \PY{n}{computeError}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p4}\PY{p}{,} \PY{n}{H}\PY{p}{)}
- \PY{n}{e5} \PY{o}{=} \PY{n}{computeError}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p5}\PY{p}{,} \PY{n}{H}\PY{p}{)}
- \PY{n}{options}\PY{p}{(}\PY{n+nb}{repr}\PY{o}{.}\PY{n}{plot}\PY{o}{.}\PY{n}{width}\PY{o}{=}\PY{l+m+mi}{9}\PY{p}{,} \PY{n+nb}{repr}\PY{o}{.}\PY{n}{plot}\PY{o}{.}\PY{n}{height}\PY{o}{=}\PY{l+m+mi}{7}\PY{p}{)}
- \PY{n}{plotError}\PY{p}{(}\PY{n+nb}{list}\PY{p}{(}\PY{n}{e1}\PY{p}{,} \PY{n}{e5}\PY{p}{,} \PY{n}{e4}\PY{p}{,} \PY{n}{e2}\PY{p}{,} \PY{n}{e3}\PY{p}{)}\PY{p}{,} \PY{n}{cols}\PY{o}{=}\PY{n}{c}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{,}\PY{n}{colors}\PY{p}{(}\PY{p}{)}\PY{p}{[}\PY{l+m+mi}{258}\PY{p}{]}\PY{p}{,}\PY{l+m+mi}{4}\PY{p}{,}\PY{l+m+mi}{6}\PY{p}{)}\PY{p}{)}
-
- \PY{c+c1}{\PYZsh{} noir: Neighbors non\PYZhy{}local (p1), bleu: Neighbors local endo (p2),}
- \PY{c+c1}{\PYZsh{} mauve: Neighbors local none (p3), vert: moyenne (p4),}
- \PY{c+c1}{\PYZsh{} rouge: persistence (p5)}
-
- \PY{n}{sum\PYZus{}p123} \PY{o}{=} \PY{n}{e1}\PY{err}{\PYZdl{}}\PY{n+nb}{abs}\PY{err}{\PYZdl{}}\PY{n}{indices} \PY{o}{+} \PY{n}{e2}\PY{err}{\PYZdl{}}\PY{n+nb}{abs}\PY{err}{\PYZdl{}}\PY{n}{indices} \PY{o}{+} \PY{n}{e3}\PY{err}{\PYZdl{}}\PY{n+nb}{abs}\PY{err}{\PYZdl{}}\PY{n}{indices}
- \PY{n}{i\PYZus{}np} \PY{o}{=} \PY{n}{which}\PY{o}{.}\PY{n}{min}\PY{p}{(}\PY{n}{sum\PYZus{}p123}\PY{p}{)} \PY{c+c1}{\PYZsh{}indice de (veille de) jour \PYZdq{}facile\PYZdq{}}
- \PY{n}{i\PYZus{}p} \PY{o}{=} \PY{n}{which}\PY{o}{.}\PY{n}{max}\PY{p}{(}\PY{n}{sum\PYZus{}p123}\PY{p}{)} \PY{c+c1}{\PYZsh{}indice de (veille de) jour \PYZdq{}difficile\PYZdq{}}
-\end{Verbatim}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_34_0.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- Dans ce cas plus favorable les intensité des erreurs absolues ont
-clairement diminué : elles restent souvent en dessous de 5. En revanche
-le MAPE moyen reste au-delà de 20\%, et même souvent plus de 30\%. Comme
-dans le cas de l'épandage on constate une croissance globale de la
-courbe journalière d'erreur absolue moyenne (en haut à gauche) ; ceci
-peut être dû au fait que l'on ajuste le niveau du jour à prédire en le
-recollant sur la dernière valeur observée.
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}20}]:} \PY{n}{options}\PY{p}{(}\PY{n+nb}{repr}\PY{o}{.}\PY{n}{plot}\PY{o}{.}\PY{n}{width}\PY{o}{=}\PY{l+m+mi}{9}\PY{p}{,} \PY{n+nb}{repr}\PY{o}{.}\PY{n}{plot}\PY{o}{.}\PY{n}{height}\PY{o}{=}\PY{l+m+mi}{4}\PY{p}{)}
- \PY{n}{par}\PY{p}{(}\PY{n}{mfrow}\PY{o}{=}\PY{n}{c}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{)}\PY{p}{)}
-
- \PY{n}{plotPredReal}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p1}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{PredReal p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotPredReal}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p1}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{PredReal p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{n}{plotPredReal}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p2}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{PredReal p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotPredReal}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p2}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{PredReal p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{n}{plotPredReal}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p3}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{PredReal p3 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotPredReal}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p3}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{PredReal p3 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{c+c1}{\PYZsh{} Bleu: prévue, noir: réalisée}
-\end{Verbatim}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_36_0.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_36_1.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_36_2.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- La forme est raisonnablement retrouvée pour les méthodes "locales",
-l'autre version lissant trop les prédictions. Le biais reste cependant
-important, surtout en fin de journée sur le jour "difficile".
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}21}]:} \PY{n}{par}\PY{p}{(}\PY{n}{mfrow}\PY{o}{=}\PY{n}{c}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{)}\PY{p}{)}
- \PY{n}{f\PYZus{}np1} \PY{o}{=} \PY{n}{computeFilaments}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p1}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{,} \PY{n}{plot}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
- \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Filaments p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{f\PYZus{}p1} \PY{o}{=} \PY{n}{computeFilaments}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p1}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{,} \PY{n}{plot}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
- \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Filaments p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{n}{f\PYZus{}np2} \PY{o}{=} \PY{n}{computeFilaments}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p2}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{,} \PY{n}{plot}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
- \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Filaments p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{f\PYZus{}p2} \PY{o}{=} \PY{n}{computeFilaments}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{p2}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{,} \PY{n}{plot}\PY{o}{=}\PY{n}{TRUE}\PY{p}{)}
- \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Filaments p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-\end{Verbatim}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_38_0.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_38_1.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- Les graphes de filaments ont encore la même allure, avec une assez
-grande variabilité observée. Cette observation est cependant trompeuse,
-comme l'indique plus bas le graphe de variabilité relative.
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}22}]:} \PY{n}{par}\PY{p}{(}\PY{n}{mfrow}\PY{o}{=}\PY{n}{c}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotFilamentsBox}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{f\PYZus{}np1}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{FilBox p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotFilamentsBox}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{f\PYZus{}p1}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{FilBox p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{c+c1}{\PYZsh{} En pointillés la courbe du jour courant + lendemain (à prédire)}
-\end{Verbatim}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_40_0.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- On peut réappliquer les mêmes remarques qu'auparavant sur les boxplots
-fonctionnels : lendemains de voisins atypiques, courbe à prévoir
-elle-même légèrement "hors norme".
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}23}]:} \PY{n}{par}\PY{p}{(}\PY{n}{mfrow}\PY{o}{=}\PY{n}{c}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotRelVar}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{f\PYZus{}np1}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{StdDev p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotRelVar}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{f\PYZus{}p1}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{StdDev p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{n}{plotRelVar}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{f\PYZus{}np2}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{StdDev p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotRelVar}\PY{p}{(}\PY{n}{data}\PY{p}{,} \PY{n}{f\PYZus{}p2}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{StdDev p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{c+c1}{\PYZsh{} Variabilité globale en rouge ; sur les 60 voisins (+ lendemains) en noir}
-\end{Verbatim}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_42_0.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_42_1.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- Cette fois la situation idéale est observée : la variabilité globale est
-nettement au-dessus de la variabilité locale. Bien que cela ne suffise
-pas à obtenir de bonnes prédictions de forme, on constate au moins
-l'amélioration dans la prédiction du niveau.
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}24}]:} \PY{n}{par}\PY{p}{(}\PY{n}{mfrow}\PY{o}{=}\PY{n}{c}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,}\PY{l+m+mi}{2}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotSimils}\PY{p}{(}\PY{n}{p1}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Weights p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotSimils}\PY{p}{(}\PY{n}{p1}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Weights p1 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{n}{plotSimils}\PY{p}{(}\PY{n}{p2}\PY{p}{,} \PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Weights p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{p}{)}
- \PY{n}{plotSimils}\PY{p}{(}\PY{n}{p2}\PY{p}{,} \PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{;} \PY{n}{title}\PY{p}{(}\PY{n}{paste}\PY{p}{(}\PY{l+s+s2}{\PYZdq{}}\PY{l+s+s2}{Weights p2 day}\PY{l+s+s2}{\PYZdq{}}\PY{p}{,}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{p}{)}
-
- \PY{c+c1}{\PYZsh{} \PYZhy{} pollué à gauche, + pollué à droite}
-\end{Verbatim}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_44_0.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- \begin{center}
- \adjustimage{max size={0.9\linewidth}{0.9\paperheight}}{output_44_1.png}
- \end{center}
- { \hspace*{\fill} \\}
-
- Concernant les poids en revanche, deux cas a priori mauvais se cumulent
-: * les poids dans le cas "non local" ne sont pas assez concentrés
-autour de 0, menant à un lissage trop fort \(-\) comme observé sur les
-graphes des courbes réalisées/prévues ; * les poids dans le cas "local"
-sont trop semblables (à cause de la trop grande fenêtre optimisée par
-validation croisée, cf. ci-dessous), résultant encore en une moyenne
-simple \(-\) mais sur moins de jours, plus proches du jour courant.
-
- \begin{Verbatim}[commandchars=\\\{\}]
-{\color{incolor}In [{\color{incolor}25}]:} \PY{c+c1}{\PYZsh{} Fenêtres sélectionnées dans ]0,7] / non\PYZhy{}loc 2 premières lignes, loc ensuite}
- \PY{n}{p1}\PY{err}{\PYZdl{}}\PY{n}{getParams}\PY{p}{(}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{err}{\PYZdl{}}\PY{n}{window}
- \PY{n}{p1}\PY{err}{\PYZdl{}}\PY{n}{getParams}\PY{p}{(}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{err}{\PYZdl{}}\PY{n}{window}
-
- \PY{n}{p2}\PY{err}{\PYZdl{}}\PY{n}{getParams}\PY{p}{(}\PY{n}{i\PYZus{}np}\PY{p}{)}\PY{err}{\PYZdl{}}\PY{n}{window}
- \PY{n}{p2}\PY{err}{\PYZdl{}}\PY{n}{getParams}\PY{p}{(}\PY{n}{i\PYZus{}p}\PY{p}{)}\PY{err}{\PYZdl{}}\PY{n}{window}
-\end{Verbatim}
-
- \begin{enumerate*}
-\item 0.205055690975915
-\item 0.703482647754766
-\end{enumerate*}
-
-
- \begin{enumerate*}
-\item 1.1038650998802
-\item 0.885155748316133
-\end{enumerate*}
-
-
- 3.64336124381868
-
-
- 6.99994501761361
-
-
- \subsection*{Bilan}
-
-Nos algorithmes à voisins ne sont pas adaptés à ce jeu de données où la
-forme varie considérablement d'un jour à l'autre. Plus généralement
-cette décorrélation de forme rend ardue la tâche de prévision pour toute
-autre méthode \(-\) du moins, nous ne savons pas comment procéder pour
-parvenir à une bonne précision.
-
-Toutefois, un espoir reste permis par exemple en aggréger les courbes
-spatialement (sur plusieurs stations situées dans la même agglomération
-ou dans une même zone).
-
-
- % Add a bibliography block to the postdoc
-
-
-
- \end{document}
+++ /dev/null
-#$# git-fat 9c8c710cadee05fb9695e943614093cb75e5c5cd 2744676
+++ /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
./ipynb_generator.py $1.gj - P=$2 H=$3
+# Timeout for 1 cell : 24*3600s = 24h
jupyter-nbconvert \
--ExecutePreprocessor.kernel_name='ir' \
- --ExecutePreprocessor.timeout=1800 \
+ --ExecutePreprocessor.timeout=86400 \
--execute $1.ipynb \
--to notebook --output $1.out.ipynb
# --to html --output=$1.html
--- /dev/null
+-----r
+library(talweg)
+
+P = ${P} #première heure de prévision
+H = ${H} #dernière heure de prévision
+
+ts_data = read.csv(system.file("extdata","pm10_mesures_H_loc.csv",
+ package="talweg"))
+exo_data = read.csv(system.file("extdata","meteo_extra_noNAs.csv",
+ package="talweg"))
+data = getData(ts_data, exo_data)
+
+indices = seq(as.Date("2015-01-01"),as.Date("2015-12-31"),"days")
+-----r
+p1 = computeForecast(data, indices, "Neighbors", "Neighbors",
+ predict_from=P, horizon=H, simtype="mix", local=FALSE)
+p2 = computeForecast(data, indices, "Neighbors", NULL,
+ predict_from=P, horizon=H, simtype="none", local=TRUE)
+p3 = computeForecast(data, indices, "Average", "Zero",
+ predict_from=P, horizon=H)
+p4 = computeForecast(data, indices, "Persistence", "Zero",
+ predict_from=P, horizon=H, same_day=TRUE)
+-----r
+e1 = computeError(data, p1, P, H)
+e2 = computeError(data, p2, P, H)
+e3 = computeError(data, p3, P, H)
+e4 = computeError(data, p4, P, H)
+options(repr.plot.width=9, repr.plot.height=7)
+plotError(list(e1, e4, e3, e2), cols=c(1,2,colors()[258],4), agg="month")
+
+# noir: Neighbors non-local (p1), bleu: Neighbors local (p2),
+# vert: moyenne (p3), rouge: persistence (p4)
projects / talweg.git / commitdiff