Commit | Line | Data |
---|---|---|
25b75559 | 1 | #' Neighbors Forecaster |
3d69ff21 | 2 | #' |
4f3fdbb8 BA |
3 | #' Predict next serie as a weighted combination of curves observed on "similar" days in |
4 | #' the past (and future if 'opera'=FALSE); the nature of the similarity is controlled by | |
5 | #' the options 'simtype' and 'local' (see below). | |
c4c329f6 | 6 | #' |
4f3fdbb8 | 7 | #' Optional arguments: |
102bcfda | 8 | #' \itemize{ |
4f3fdbb8 BA |
9 | #' \item local: TRUE (default) to constrain neighbors to be "same days in same season" |
10 | #' \item simtype: 'endo' for a similarity based on the series only,<cr> | |
d2ab47a7 | 11 | #' 'exo' for a similarity based on exogenous variables only,<cr> |
102bcfda BA |
12 | #' 'mix' for the product of 'endo' and 'exo',<cr> |
13 | #' 'none' (default) to apply a simple average: no computed weights | |
4f3fdbb8 | 14 | #' \item window: A window for similarities computations; override cross-validation |
102bcfda BA |
15 | #' window estimation. |
16 | #' } | |
17 | #' The method is summarized as follows: | |
18 | #' \enumerate{ | |
4f3fdbb8 BA |
19 | #' \item Determine N (=20) recent days without missing values, and preceded by a |
20 | #' curve also without missing values. | |
102bcfda BA |
21 | #' \item Optimize the window parameters (if relevant) on the N chosen days. |
22 | #' \item Considering the optimized window, compute the neighbors (with locality | |
23 | #' constraint or not), compute their similarities -- using a gaussian kernel if | |
24 | #' simtype != "none" -- and average accordingly the "tomorrows of neigbors" to | |
25 | #' obtain the final prediction. | |
26 | #' } | |
c4c329f6 | 27 | #' |
4e821712 | 28 | #' @usage # NeighborsForecaster$new(pjump) |
689aa1d3 | 29 | #' |
102bcfda | 30 | #' @docType class |
c4c329f6 | 31 | #' @format R6 class, inherits Forecaster |
3ddf1c12 | 32 | #' @aliases F_Neighbors |
546b0cb6 | 33 | #' |
25b75559 | 34 | NeighborsForecaster = R6::R6Class("NeighborsForecaster", |
a66a84b5 | 35 | inherit = Forecaster, |
25b75559 BA |
36 | |
37 | public = list( | |
d2ab47a7 | 38 | predictShape = function(data, today, memory, predict_from, horizon, ...) |
3d69ff21 BA |
39 | { |
40 | # (re)initialize computed parameters | |
a66a84b5 | 41 | private$.params <- list("weights"=NA, "indices"=NA, "window"=NA) |
3d69ff21 | 42 | |
a5a3a294 | 43 | # Do not forecast on days with NAs (TODO: softer condition...) |
cf3bb001 BA |
44 | if (any(is.na(data$getSerie(today-1))) || |
45 | (predict_from>=2 && any(is.na(data$getSerie(today)[1:(predict_from-1)])))) | |
d2ab47a7 | 46 | { |
a5a3a294 | 47 | return (NA) |
d2ab47a7 | 48 | } |
a5a3a294 | 49 | |
f17665c7 | 50 | # Get optional args |
445e7bbc BA |
51 | local = ifelse(hasArg("local"), list(...)$local, TRUE) #same level + season? |
52 | simtype = ifelse(hasArg("simtype"), list(...)$simtype, "none") #or "endo", or "exo" | |
638f27f4 BA |
53 | opera = ifelse(hasArg("opera"), list(...)$opera, FALSE) #operational mode? |
54 | ||
55 | # Determine indices of no-NAs days preceded by no-NAs yerstedays | |
56 | tdays = .getNoNA2(data, max(today-memory,2), ifelse(opera,today-1,data$getSize())) | |
57 | if (!opera) | |
58 | tdays = setdiff(tdays, today) #always exclude current day | |
59 | ||
dca259e4 BA |
60 | # Shortcut if window is known or local==TRUE && simtype==none |
61 | if (hasArg("window") || (local && simtype=="none")) | |
a66a84b5 | 62 | { |
638f27f4 BA |
63 | return ( private$.predictShapeAux(data, tdays, today, predict_from, horizon, |
64 | local, list(...)$window, simtype, opera, TRUE) ) | |
a66a84b5 | 65 | } |
3d69ff21 | 66 | |
6774e53d | 67 | # Indices of similar days for cross-validation; TODO: 20 = magic number |
aa059de7 | 68 | cv_days = getSimilarDaysIndices(today, data, limit=20, same_season=FALSE, |
638f27f4 | 69 | days_in=tdays, operational=opera) |
5e838b3e | 70 | |
445e7bbc | 71 | # Optimize h : h |--> sum of prediction errors on last N "similar" days |
aa059de7 | 72 | errorOnLastNdays = function(window, simtype) |
3d69ff21 BA |
73 | { |
74 | error = 0 | |
75 | nb_jours = 0 | |
5e838b3e | 76 | for (i in seq_along(cv_days)) |
3d69ff21 | 77 | { |
f17665c7 | 78 | # mix_strategy is never used here (simtype != "mix"), therefore left blank |
cf3bb001 | 79 | prediction = private$.predictShapeAux(data, tdays, cv_days[i], predict_from, |
638f27f4 | 80 | horizon, local, window, simtype, opera, FALSE) |
f17665c7 | 81 | if (!is.na(prediction[1])) |
3d69ff21 BA |
82 | { |
83 | nb_jours = nb_jours + 1 | |
af3b84f4 | 84 | error = error + |
cf3bb001 | 85 | mean((data$getSerie(cv_days[i])[predict_from:horizon] - prediction)^2) |
3d69ff21 BA |
86 | } |
87 | } | |
88 | return (error / nb_jours) | |
89 | } | |
90 | ||
445e7bbc | 91 | # TODO: 7 == magic number |
eef54517 | 92 | if (simtype=="endo" || simtype=="mix") |
af3b84f4 | 93 | { |
aa059de7 BA |
94 | best_window_endo = optimize( |
95 | errorOnLastNdays, c(0,7), simtype="endo")$minimum | |
af3b84f4 | 96 | } |
eef54517 | 97 | if (simtype=="exo" || simtype=="mix") |
af3b84f4 | 98 | { |
eef54517 BA |
99 | best_window_exo = optimize( |
100 | errorOnLastNdays, c(0,7), simtype="exo")$minimum | |
3d69ff21 | 101 | } |
eef54517 BA |
102 | |
103 | best_window = | |
104 | if (simtype == "endo") | |
105 | best_window_endo | |
106 | else if (simtype == "exo") | |
107 | best_window_exo | |
108 | else if (simtype == "mix") | |
109 | c(best_window_endo,best_window_exo) | |
110 | else #none: value doesn't matter | |
111 | 1 | |
112 | ||
cf3bb001 | 113 | return( private$.predictShapeAux(data, tdays, today, predict_from, horizon, local, |
638f27f4 | 114 | best_window, simtype, opera, TRUE) ) |
25b75559 BA |
115 | } |
116 | ), | |
117 | private = list( | |
638f27f4 | 118 | # Precondition: "yersteday until predict_from-1" is full (no NAs) |
cf3bb001 | 119 | .predictShapeAux = function(data, tdays, today, predict_from, horizon, local, window, |
638f27f4 | 120 | simtype, opera, final_call) |
3d69ff21 | 121 | { |
638f27f4 BA |
122 | tdays_cut = tdays[ tdays != today ] |
123 | if (length(tdays_cut) == 0) | |
3d69ff21 BA |
124 | return (NA) |
125 | ||
aa059de7 BA |
126 | if (local) |
127 | { | |
dca259e4 BA |
128 | # limit=Inf to not censor any day (TODO: finite limit? 60?) |
129 | tdays = getSimilarDaysIndices(today, data, limit=Inf, same_season=TRUE, | |
638f27f4 | 130 | days_in=tdays_cut, operational=opera) |
638f27f4 BA |
131 | # TODO: 10 == magic number |
132 | tdays = .getConstrainedNeighbs(today, data, tdays, min_neighbs=10) | |
cf3bb001 | 133 | if (length(tdays) == 1) |
aa059de7 BA |
134 | { |
135 | if (final_call) | |
136 | { | |
137 | private$.params$weights <- 1 | |
cf3bb001 | 138 | private$.params$indices <- tdays |
aa059de7 BA |
139 | private$.params$window <- 1 |
140 | } | |
cf3bb001 | 141 | return ( data$getSerie(tdays[1])[predict_from:horizon] ) |
aa059de7 | 142 | } |
10886062 | 143 | max_neighbs = 12 #TODO: 10 or 12 or... ? |
dca259e4 BA |
144 | if (length(tdays) > max_neighbs) |
145 | { | |
146 | distances2 <- .computeDistsEndo(data, today, tdays, predict_from) | |
147 | ordering <- order(distances2) | |
148 | tdays <- tdays[ ordering[1:max_neighbs] ] | |
149 | } | |
aa059de7 BA |
150 | } |
151 | else | |
cf3bb001 | 152 | tdays = tdays_cut #no conditioning |
aa059de7 | 153 | |
445e7bbc | 154 | if (simtype == "endo" || simtype == "mix") |
3d69ff21 | 155 | { |
aa059de7 BA |
156 | # Compute endogen similarities using given window |
157 | window_endo = ifelse(simtype=="mix", window[1], window) | |
3d69ff21 | 158 | |
638f27f4 | 159 | # Distances from last observed day to selected days in the past |
dca259e4 | 160 | # TODO: redundant computation if local==TRUE |
638f27f4 BA |
161 | distances2 <- .computeDistsEndo(data, today, tdays, predict_from) |
162 | ||
3ddf1c12 | 163 | simils_endo <- .computeSimils(distances2, window_endo) |
3d69ff21 BA |
164 | } |
165 | ||
445e7bbc | 166 | if (simtype == "exo" || simtype == "mix") |
3d69ff21 | 167 | { |
aa059de7 | 168 | # Compute exogen similarities using given window |
445e7bbc | 169 | window_exo = ifelse(simtype=="mix", window[2], window) |
3d69ff21 | 170 | |
638f27f4 | 171 | distances2 <- .computeDistsExo(data, today, tdays) |
3d69ff21 | 172 | |
3ddf1c12 | 173 | simils_exo <- .computeSimils(distances2, window_exo) |
3d69ff21 BA |
174 | } |
175 | ||
3d69ff21 | 176 | similarities = |
f17665c7 | 177 | if (simtype == "exo") |
3d69ff21 | 178 | simils_exo |
f17665c7 BA |
179 | else if (simtype == "endo") |
180 | simils_endo | |
445e7bbc | 181 | else if (simtype == "mix") |
f17665c7 | 182 | simils_endo * simils_exo |
445e7bbc | 183 | else #none |
cf3bb001 | 184 | rep(1, length(tdays)) |
ea5c7e56 | 185 | similarities = similarities / sum(similarities) |
3d69ff21 | 186 | |
d2ab47a7 | 187 | prediction = rep(0, horizon-predict_from+1) |
cf3bb001 | 188 | for (i in seq_along(tdays)) |
d2ab47a7 BA |
189 | { |
190 | prediction = prediction + | |
cf3bb001 | 191 | similarities[i] * data$getSerie(tdays[i])[predict_from:horizon] |
d2ab47a7 | 192 | } |
99f83c9a | 193 | |
3d69ff21 BA |
194 | if (final_call) |
195 | { | |
a66a84b5 | 196 | private$.params$weights <- similarities |
cf3bb001 | 197 | private$.params$indices <- tdays |
a66a84b5 | 198 | private$.params$window <- |
546b0cb6 | 199 | if (simtype=="endo") |
aa059de7 | 200 | window_endo |
546b0cb6 | 201 | else if (simtype=="exo") |
aa059de7 | 202 | window_exo |
eef54517 | 203 | else if (simtype=="mix") |
aa059de7 | 204 | c(window_endo,window_exo) |
eef54517 BA |
205 | else #none |
206 | 1 | |
3d69ff21 | 207 | } |
99f83c9a | 208 | |
3d69ff21 BA |
209 | return (prediction) |
210 | } | |
211 | ) | |
212 | ) | |
3ddf1c12 | 213 | |
689aa1d3 BA |
214 | # getConstrainedNeighbs |
215 | # | |
216 | # Get indices of neighbors of similar pollution level (among same season + day type). | |
217 | # | |
218 | # @param today Index of current day | |
219 | # @param data Object of class Data | |
cf3bb001 | 220 | # @param tdays Current set of "second days" (no-NA pairs) |
689aa1d3 BA |
221 | # @param min_neighbs Minimum number of points in a neighborhood |
222 | # @param max_neighbs Maximum number of points in a neighborhood | |
223 | # | |
638f27f4 | 224 | .getConstrainedNeighbs = function(today, data, tdays, min_neighbs=10) |
3ddf1c12 | 225 | { |
d2ab47a7 | 226 | levelToday = data$getLevelHat(today) |
9b9bb2d4 | 227 | distances = sapply( tdays, function(i) abs(data$getLevel(i) - levelToday) ) |
d2ab47a7 BA |
228 | #TODO: 1, +1, +3 : magic numbers |
229 | dist_thresh = 1 | |
cf3bb001 | 230 | min_neighbs = min(min_neighbs,length(tdays)) |
3ddf1c12 BA |
231 | repeat |
232 | { | |
233 | same_pollution = (distances <= dist_thresh) | |
234 | nb_neighbs = sum(same_pollution) | |
235 | if (nb_neighbs >= min_neighbs) #will eventually happen | |
236 | break | |
d2ab47a7 | 237 | dist_thresh = dist_thresh + ifelse(dist_thresh>1,3,1) |
3ddf1c12 | 238 | } |
9b9bb2d4 | 239 | tdays[same_pollution] |
3ddf1c12 BA |
240 | } |
241 | ||
689aa1d3 BA |
242 | # compute similarities |
243 | # | |
244 | # Apply the gaussian kernel on computed squared distances. | |
245 | # | |
246 | # @param distances2 Squared distances | |
247 | # @param window Window parameter for the kernel | |
248 | # | |
3ddf1c12 BA |
249 | .computeSimils <- function(distances2, window) |
250 | { | |
251 | sd_dist = sd(distances2) | |
252 | if (sd_dist < .25 * sqrt(.Machine$double.eps)) | |
253 | { | |
254 | # warning("All computed distances are very close: stdev too small") | |
255 | sd_dist = 1 #mostly for tests... FIXME: | |
256 | } | |
257 | exp(-distances2/(sd_dist*window^2)) | |
258 | } | |
638f27f4 BA |
259 | |
260 | .computeDistsEndo <- function(data, today, tdays, predict_from) | |
261 | { | |
262 | lastSerie = c( data$getSerie(today-1), | |
263 | data$getSerie(today)[if (predict_from>=2) 1:(predict_from-1) else c()] ) | |
264 | sapply(tdays, function(i) { | |
265 | delta = lastSerie - c(data$getSerie(i-1), | |
266 | data$getSerie(i)[if (predict_from>=2) 1:(predict_from-1) else c()]) | |
9b9bb2d4 | 267 | sqrt(mean(delta^2)) |
638f27f4 BA |
268 | }) |
269 | } | |
270 | ||
271 | .computeDistsExo <- function(data, today, tdays) | |
272 | { | |
273 | M = matrix( ncol=1+length(tdays), nrow=1+length(data$getExo(1)) ) | |
274 | M[,1] = c( data$getLevelHat(today), as.double(data$getExoHat(today)) ) | |
275 | for (i in seq_along(tdays)) | |
276 | M[,i+1] = c( data$getLevel(tdays[i]), as.double(data$getExo(tdays[i])) ) | |
277 | ||
278 | sigma = cov(t(M)) #NOTE: robust covariance is way too slow | |
279 | # TODO: 10 == magic number; more robust way == det, or always ginv() | |
280 | sigma_inv = | |
281 | if (length(tdays) > 10) | |
282 | solve(sigma) | |
283 | else | |
284 | MASS::ginv(sigma) | |
285 | ||
286 | # Distances from last observed day to days in the past | |
287 | sapply(seq_along(tdays), function(i) { | |
288 | delta = M[,1] - M[,i+1] | |
289 | delta %*% sigma_inv %*% delta | |
290 | }) | |
291 | } |