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