pkg/R/F_Neighbors.R

   1 #' Neighbors Forecaster
   2 #'
   3 #' Predict next serie as a weighted combination of "futures of the past" days,
   4 #' where days in the past are chosen and weighted according to some similarity measures.
   5 #'
   6 #' The main method is \code{predictShape()}, taking arguments data, today, memory,
   7 #' horizon respectively for the dataset (object output of \code{getData()}), the current
   8 #' index, the data depth (in days) and the number of time steps to forecast.
   9 #' In addition, optional arguments can be passed:
  10 #' \itemize{
  11 #'   \item local : TRUE (default) to constrain neighbors to be "same days within same
  12 #'     season"
  13 #'   \item simtype : 'endo' for a similarity based on the series only,<cr>
  14 #'             'exo' for a similaruty based on exogenous variables only,<cr>
  15 #'             'mix' for the product of 'endo' and 'exo',<cr>
  16 #'             'none' (default) to apply a simple average: no computed weights
  17 #'   \item window : A window for similarities computations; override cross-validation
  18 #'     window estimation.
  19 #' }
  20 #' The method is summarized as follows:
  21 #' \enumerate{
  22 #'   \item Determine N (=20) recent days without missing values, and followed by a
  23 #'     tomorrow also without missing values.
  24 #'   \item Optimize the window parameters (if relevant) on the N chosen days.
  25 #'   \item Considering the optimized window, compute the neighbors (with locality
  26 #'     constraint or not), compute their similarities -- using a gaussian kernel if
  27 #'     simtype != "none" -- and average accordingly the "tomorrows of neigbors" to
  28 #'     obtain the final prediction.
  29 #' }
  30 #'
  31 #' @docType class
  32 #' @format R6 class, inherits Forecaster
  33 #' @alias F_Neighbors
  34 #'
  35 NeighborsForecaster = R6::R6Class("NeighborsForecaster",
  36     inherit = Forecaster,
  37
  38     public = list(
  39         predictShape = function(data, today, memory, horizon, ...)
  40         {
  41             # (re)initialize computed parameters
  42             private$.params <- list("weights"=NA, "indices"=NA, "window"=NA)
  43
  44             # Do not forecast on days with NAs (TODO: softer condition...)
  45             if (any(is.na(data$getCenteredSerie(today))))
  46                 return (NA)
  47
  48             # Determine indices of no-NAs days followed by no-NAs tomorrows
  49             fdays = getNoNA2(data, max(today-memory,1), today-1)
  50
  51             # Get optional args
  52             local = ifelse(hasArg("local"), list(...)$local, TRUE) #same level + season?
  53             simtype = ifelse(hasArg("simtype"), list(...)$simtype, "none") #or "endo", or "exo"
  54             if (hasArg("window"))
  55             {
  56                 return ( private$.predictShapeAux(data,
  57                     fdays, today, horizon, local, list(...)$window, simtype, TRUE) )
  58             }
  59
  60             # Indices of similar days for cross-validation; TODO: 20 = magic number
  61             cv_days = getSimilarDaysIndices(today, data, limit=20, same_season=FALSE,
  62                 days_in=fdays)
  63
  64             # Optimize h : h |--> sum of prediction errors on last N "similar" days
  65             errorOnLastNdays = function(window, simtype)
  66             {
  67                 error = 0
  68                 nb_jours = 0
  69                 for (i in seq_along(cv_days))
  70                 {
  71                     # mix_strategy is never used here (simtype != "mix"), therefore left blank
  72                     prediction = private$.predictShapeAux(data, fdays, cv_days[i], horizon, local,
  73                         window, simtype, FALSE)
  74                     if (!is.na(prediction[1]))
  75                     {
  76                         nb_jours = nb_jours + 1
  77                         error = error +
  78                             mean((data$getSerie(cv_days[i]+1)[1:horizon] - prediction)^2)
  79                     }
  80                 }
  81                 return (error / nb_jours)
  82             }
  83
  84             # TODO: 7 == magic number
  85             if (simtype=="endo" || simtype=="mix")
  86             {
  87                 best_window_endo = optimize(
  88                     errorOnLastNdays, c(0,7), simtype="endo")$minimum
  89             }
  90             if (simtype=="exo" || simtype=="mix")
  91             {
  92                 best_window_exo = optimize(
  93                     errorOnLastNdays, c(0,7), simtype="exo")$minimum
  94             }
  95
  96             best_window =
  97                 if (simtype == "endo")
  98                     best_window_endo
  99                 else if (simtype == "exo")
 100                     best_window_exo
 101                 else if (simtype == "mix")
 102                     c(best_window_endo,best_window_exo)
 103                 else #none: value doesn't matter
 104                     1
 105
 106             return(private$.predictShapeAux(data, fdays, today, horizon, local,
 107                 best_window, simtype, TRUE))
 108         }
 109     ),
 110     private = list(
 111         # Precondition: "today" is full (no NAs)
 112         .predictShapeAux = function(data, fdays, today, horizon, local, window, simtype,
 113             final_call)
 114         {
 115             fdays_cut = fdays[ fdays < today ]
 116             if (length(fdays_cut) <= 1)
 117                 return (NA)
 118
 119             if (local)
 120             {
 121                 # Neighbors: days in "same season"; TODO: 60 == magic number...
 122                 fdays = getSimilarDaysIndices(today, data, limit=60, same_season=TRUE,
 123                     days_in=fdays_cut)
 124                 if (length(fdays) <= 1)
 125                     return (NA)
 126                 levelToday = data$getLevel(today)
 127                 distances = sapply(fdays, function(i) abs(data$getLevel(i)-levelToday))
 128                 #TODO: 2, 10, 3, 12 magic numbers here...
 129                 dist_thresh = 2
 130                 min_neighbs = min(10,length(fdays))
 131                 repeat
 132                 {
 133                     same_pollution = (distances <= dist_thresh)
 134                     nb_neighbs = sum(same_pollution)
 135                     if (nb_neighbs >= min_neighbs) #will eventually happen
 136                         break
 137                     dist_thresh = dist_thresh + 3
 138                 }
 139                 fdays = fdays[same_pollution]
 140                 max_neighbs = 12
 141                 if (nb_neighbs > max_neighbs)
 142                 {
 143                     # Keep only max_neighbs closest neighbors
 144                     fdays = fdays[
 145                         sort(distances[same_pollution],index.return=TRUE)$ix[1:max_neighbs] ]
 146                 }
 147                 if (length(fdays) == 1) #the other extreme...
 148                 {
 149                     if (final_call)
 150                     {
 151                         private$.params$weights <- 1
 152                         private$.params$indices <- fdays
 153                         private$.params$window <- 1
 154                     }
 155                     return ( data$getSerie(fdays[1])[1:horizon] )
 156                 }
 157             }
 158             else
 159                 fdays = fdays_cut #no conditioning
 160
 161             if (simtype == "endo" || simtype == "mix")
 162             {
 163                 # Compute endogen similarities using given window
 164                 window_endo = ifelse(simtype=="mix", window[1], window)
 165
 166                 # Distances from last observed day to days in the past
 167                 serieToday = data$getSerie(today)
 168                 distances2 = sapply(fdays, function(i) {
 169                     delta = serieToday - data$getSerie(i)
 170                     mean(delta^2)
 171                 })
 172
 173                 sd_dist = sd(distances2)
 174                 if (sd_dist < .25 * sqrt(.Machine$double.eps))
 175                 {
 176 #                   warning("All computed distances are very close: stdev too small")
 177                     sd_dist = 1 #mostly for tests... FIXME:
 178                 }
 179                 simils_endo = exp(-distances2/(sd_dist*window_endo^2))
 180             }
 181
 182             if (simtype == "exo" || simtype == "mix")
 183             {
 184                 # Compute exogen similarities using given window
 185                 window_exo = ifelse(simtype=="mix", window[2], window)
 186
 187                 M = matrix( nrow=1+length(fdays), ncol=1+length(data$getExo(today)) )
 188                 M[1,] = c( data$getLevel(today), as.double(data$getExo(today)) )
 189                 for (i in seq_along(fdays))
 190                     M[i+1,] = c( data$getLevel(fdays[i]), as.double(data$getExo(fdays[i])) )
 191
 192                 sigma = cov(M) #NOTE: robust covariance is way too slow
 193                 # TODO: 10 == magic number; more robust way == det, or always ginv()
 194                 sigma_inv =
 195                     if (length(fdays) > 10)
 196                         solve(sigma)
 197                     else
 198                         MASS::ginv(sigma)
 199
 200                 # Distances from last observed day to days in the past
 201                 distances2 = sapply(seq_along(fdays), function(i) {
 202                     delta = M[1,] - M[i+1,]
 203                     delta %*% sigma_inv %*% delta
 204                 })
 205
 206                 sd_dist = sd(distances2)
 207                 if (sd_dist < .25 * sqrt(.Machine$double.eps))
 208                 {
 209 #                   warning("All computed distances are very close: stdev too small")
 210                     sd_dist = 1 #mostly for tests... FIXME:
 211                 }
 212                 simils_exo = exp(-distances2/(sd_dist*window_exo^2))
 213             }
 214
 215             similarities =
 216                 if (simtype == "exo")
 217                     simils_exo
 218                 else if (simtype == "endo")
 219                     simils_endo
 220                 else if (simtype == "mix")
 221                     simils_endo * simils_exo
 222                 else #none
 223                     rep(1, length(fdays))
 224             similarities = similarities / sum(similarities)
 225
 226             prediction = rep(0, horizon)
 227             for (i in seq_along(fdays))
 228                 prediction = prediction + similarities[i] * data$getSerie(fdays[i]+1)[1:horizon]
 229
 230             if (final_call)
 231             {
 232                 private$.params$weights <- similarities
 233                 private$.params$indices <- fdays
 234                 private$.params$window <-
 235                     if (simtype=="endo")
 236                         window_endo
 237                     else if (simtype=="exo")
 238                         window_exo
 239                     else if (simtype=="mix")
 240                         c(window_endo,window_exo)
 241                     else #none
 242                         1
 243             }
 244
 245             return (prediction)
 246         }
 247     )
 248 )