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