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