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 BA |
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 | #' } | |
c4c329f6 | 30 | #' |
102bcfda | 31 | #' @docType class |
c4c329f6 BA |
32 | #' @format R6 class, inherits Forecaster |
33 | #' @alias F_Neighbors | |
546b0cb6 | 34 | #' |
25b75559 | 35 | NeighborsForecaster = R6::R6Class("NeighborsForecaster", |
a66a84b5 | 36 | inherit = Forecaster, |
25b75559 BA |
37 | |
38 | public = list( | |
98e958ca | 39 | predictShape = function(data, today, memory, horizon, ...) |
3d69ff21 BA |
40 | { |
41 | # (re)initialize computed parameters | |
a66a84b5 | 42 | private$.params <- list("weights"=NA, "indices"=NA, "window"=NA) |
3d69ff21 | 43 | |
a5a3a294 BA |
44 | # Do not forecast on days with NAs (TODO: softer condition...) |
45 | if (any(is.na(data$getCenteredSerie(today)))) | |
46 | return (NA) | |
47 | ||
af3b84f4 | 48 | # Determine indices of no-NAs days followed by no-NAs tomorrows |
98e958ca | 49 | fdays = getNoNA2(data, max(today-memory,1), today-1) |
af3b84f4 | 50 | |
f17665c7 | 51 | # Get optional args |
445e7bbc BA |
52 | local = ifelse(hasArg("local"), list(...)$local, TRUE) #same level + season? |
53 | simtype = ifelse(hasArg("simtype"), list(...)$simtype, "none") #or "endo", or "exo" | |
aa059de7 | 54 | if (hasArg("window")) |
a66a84b5 | 55 | { |
98e958ca | 56 | return ( private$.predictShapeAux(data, |
aa059de7 | 57 | fdays, today, horizon, local, list(...)$window, simtype, TRUE) ) |
a66a84b5 | 58 | } |
3d69ff21 | 59 | |
6774e53d | 60 | # Indices of similar days for cross-validation; TODO: 20 = magic number |
aa059de7 BA |
61 | cv_days = getSimilarDaysIndices(today, data, limit=20, same_season=FALSE, |
62 | days_in=fdays) | |
5e838b3e | 63 | |
445e7bbc | 64 | # Optimize h : h |--> sum of prediction errors on last N "similar" days |
aa059de7 | 65 | errorOnLastNdays = function(window, simtype) |
3d69ff21 BA |
66 | { |
67 | error = 0 | |
68 | nb_jours = 0 | |
5e838b3e | 69 | for (i in seq_along(cv_days)) |
3d69ff21 | 70 | { |
f17665c7 | 71 | # mix_strategy is never used here (simtype != "mix"), therefore left blank |
aa059de7 BA |
72 | prediction = private$.predictShapeAux(data, fdays, cv_days[i], horizon, local, |
73 | window, simtype, FALSE) | |
f17665c7 | 74 | if (!is.na(prediction[1])) |
3d69ff21 BA |
75 | { |
76 | nb_jours = nb_jours + 1 | |
af3b84f4 | 77 | error = error + |
aa059de7 | 78 | mean((data$getSerie(cv_days[i]+1)[1:horizon] - prediction)^2) |
3d69ff21 BA |
79 | } |
80 | } | |
81 | return (error / nb_jours) | |
82 | } | |
83 | ||
445e7bbc | 84 | # TODO: 7 == magic number |
eef54517 | 85 | if (simtype=="endo" || simtype=="mix") |
af3b84f4 | 86 | { |
aa059de7 BA |
87 | best_window_endo = optimize( |
88 | errorOnLastNdays, c(0,7), simtype="endo")$minimum | |
af3b84f4 | 89 | } |
eef54517 | 90 | if (simtype=="exo" || simtype=="mix") |
af3b84f4 | 91 | { |
eef54517 BA |
92 | best_window_exo = optimize( |
93 | errorOnLastNdays, c(0,7), simtype="exo")$minimum | |
3d69ff21 | 94 | } |
eef54517 BA |
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)) | |
25b75559 BA |
108 | } |
109 | ), | |
110 | private = list( | |
3d69ff21 | 111 | # Precondition: "today" is full (no NAs) |
aa059de7 BA |
112 | .predictShapeAux = function(data, fdays, today, horizon, local, window, simtype, |
113 | final_call) | |
3d69ff21 | 114 | { |
aa059de7 BA |
115 | fdays_cut = fdays[ fdays < today ] |
116 | if (length(fdays_cut) <= 1) | |
3d69ff21 BA |
117 | return (NA) |
118 | ||
aa059de7 BA |
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)) | |
445e7bbc | 128 | #TODO: 2, 10, 3, 12 magic numbers here... |
aa059de7 | 129 | dist_thresh = 2 |
445e7bbc | 130 | min_neighbs = min(10,length(fdays)) |
aa059de7 BA |
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] | |
445e7bbc | 140 | max_neighbs = 12 |
aa059de7 BA |
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 | } | |
2057c793 | 155 | return ( data$getSerie(fdays[1])[1:horizon] ) |
aa059de7 BA |
156 | } |
157 | } | |
158 | else | |
159 | fdays = fdays_cut #no conditioning | |
160 | ||
445e7bbc | 161 | if (simtype == "endo" || simtype == "mix") |
3d69ff21 | 162 | { |
aa059de7 BA |
163 | # Compute endogen similarities using given window |
164 | window_endo = ifelse(simtype=="mix", window[1], window) | |
3d69ff21 BA |
165 | |
166 | # Distances from last observed day to days in the past | |
5e838b3e BA |
167 | serieToday = data$getSerie(today) |
168 | distances2 = sapply(fdays, function(i) { | |
169 | delta = serieToday - data$getSerie(i) | |
170 | mean(delta^2) | |
171 | }) | |
3d69ff21 BA |
172 | |
173 | sd_dist = sd(distances2) | |
aa059de7 | 174 | if (sd_dist < .25 * sqrt(.Machine$double.eps)) |
546b0cb6 | 175 | { |
fa5b7bfc | 176 | # warning("All computed distances are very close: stdev too small") |
99f83c9a | 177 | sd_dist = 1 #mostly for tests... FIXME: |
546b0cb6 | 178 | } |
aa059de7 | 179 | simils_endo = exp(-distances2/(sd_dist*window_endo^2)) |
3d69ff21 BA |
180 | } |
181 | ||
445e7bbc | 182 | if (simtype == "exo" || simtype == "mix") |
3d69ff21 | 183 | { |
aa059de7 | 184 | # Compute exogen similarities using given window |
445e7bbc | 185 | window_exo = ifelse(simtype=="mix", window[2], window) |
3d69ff21 | 186 | |
25b75559 BA |
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)) ) | |
f17665c7 | 189 | for (i in seq_along(fdays)) |
25b75559 | 190 | M[i+1,] = c( data$getLevel(fdays[i]), as.double(data$getExo(fdays[i])) ) |
3d69ff21 BA |
191 | |
192 | sigma = cov(M) #NOTE: robust covariance is way too slow | |
ee8b1b4e BA |
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) | |
3d69ff21 BA |
199 | |
200 | # Distances from last observed day to days in the past | |
5e838b3e BA |
201 | distances2 = sapply(seq_along(fdays), function(i) { |
202 | delta = M[1,] - M[i+1,] | |
203 | delta %*% sigma_inv %*% delta | |
204 | }) | |
3d69ff21 BA |
205 | |
206 | sd_dist = sd(distances2) | |
ee8b1b4e | 207 | if (sd_dist < .25 * sqrt(.Machine$double.eps)) |
546b0cb6 | 208 | { |
fa5b7bfc | 209 | # warning("All computed distances are very close: stdev too small") |
546b0cb6 BA |
210 | sd_dist = 1 #mostly for tests... FIXME: |
211 | } | |
aa059de7 | 212 | simils_exo = exp(-distances2/(sd_dist*window_exo^2)) |
3d69ff21 BA |
213 | } |
214 | ||
3d69ff21 | 215 | similarities = |
f17665c7 | 216 | if (simtype == "exo") |
3d69ff21 | 217 | simils_exo |
f17665c7 BA |
218 | else if (simtype == "endo") |
219 | simils_endo | |
445e7bbc | 220 | else if (simtype == "mix") |
f17665c7 | 221 | simils_endo * simils_exo |
445e7bbc BA |
222 | else #none |
223 | rep(1, length(fdays)) | |
ea5c7e56 | 224 | similarities = similarities / sum(similarities) |
3d69ff21 BA |
225 | |
226 | prediction = rep(0, horizon) | |
a66a84b5 | 227 | for (i in seq_along(fdays)) |
aa059de7 | 228 | prediction = prediction + similarities[i] * data$getSerie(fdays[i]+1)[1:horizon] |
99f83c9a | 229 | |
3d69ff21 BA |
230 | if (final_call) |
231 | { | |
a66a84b5 BA |
232 | private$.params$weights <- similarities |
233 | private$.params$indices <- fdays | |
234 | private$.params$window <- | |
546b0cb6 | 235 | if (simtype=="endo") |
aa059de7 | 236 | window_endo |
546b0cb6 | 237 | else if (simtype=="exo") |
aa059de7 | 238 | window_exo |
eef54517 | 239 | else if (simtype=="mix") |
aa059de7 | 240 | c(window_endo,window_exo) |
eef54517 BA |
241 | else #none |
242 | 1 | |
3d69ff21 | 243 | } |
99f83c9a | 244 | |
3d69ff21 BA |
245 | return (prediction) |
246 | } | |
247 | ) | |
248 | ) |