'update'
[talweg.git] / pkg / R / F_Neighbors.R
CommitLineData
e030a6e3 1#' @include Forecaster.R
3d69ff21 2#'
25b75559 3#' Neighbors Forecaster
3d69ff21 4#'
25b75559
BA
5#' Predict tomorrow as a weighted combination of "futures of the past" days.
6#' Inherits \code{\link{Forecaster}}
7NeighborsForecaster = R6::R6Class("NeighborsForecaster",
a66a84b5 8 inherit = Forecaster,
25b75559
BA
9
10 public = list(
98e958ca 11 predictShape = function(data, today, memory, horizon, ...)
3d69ff21
BA
12 {
13 # (re)initialize computed parameters
a66a84b5 14 private$.params <- list("weights"=NA, "indices"=NA, "window"=NA)
3d69ff21 15
a5a3a294
BA
16 # Do not forecast on days with NAs (TODO: softer condition...)
17 if (any(is.na(data$getCenteredSerie(today))))
18 return (NA)
19
af3b84f4 20 # Determine indices of no-NAs days followed by no-NAs tomorrows
98e958ca 21 fdays = getNoNA2(data, max(today-memory,1), today-1)
af3b84f4 22
f17665c7
BA
23 # Get optional args
24 simtype = ifelse(hasArg("simtype"), list(...)$simtype, "mix") #or "endo", or "exo"
25 kernel = ifelse(hasArg("kernel"), list(...)$kernel, "Gauss") #or "Epan"
26 if (hasArg(h_window))
a66a84b5 27 {
98e958ca 28 return ( private$.predictShapeAux(data,
a66a84b5
BA
29 fdays, today, horizon, list(...)$h_window, kernel, simtype, TRUE) )
30 }
3d69ff21 31
f17665c7
BA
32 # Indices of similar days for cross-validation; TODO: 45 = magic number
33 sdays = getSimilarDaysIndices(today, limit=45, same_season=FALSE)
3d69ff21
BA
34
35 # Function to optimize h : h |--> sum of prediction errors on last 45 "similar" days
36 errorOnLastNdays = function(h, kernel, simtype)
37 {
38 error = 0
39 nb_jours = 0
f17665c7 40 for (i in intersect(fdays,sdays))
3d69ff21 41 {
f17665c7 42 # mix_strategy is never used here (simtype != "mix"), therefore left blank
98e958ca
BA
43 prediction = private$.predictShapeAux(data,
44 fdays, i, horizon, h, kernel, simtype, FALSE)
f17665c7 45 if (!is.na(prediction[1]))
3d69ff21
BA
46 {
47 nb_jours = nb_jours + 1
af3b84f4 48 error = error +
98e958ca 49 mean((data$getCenteredSerie(i+1)[1:horizon] - prediction)^2)
3d69ff21
BA
50 }
51 }
52 return (error / nb_jours)
53 }
54
f17665c7 55 if (simtype != "endo")
af3b84f4
BA
56 {
57 h_best_exo = optimize(
58 errorOnLastNdays, c(0,10), kernel=kernel, simtype="exo")$minimum
59 }
3d69ff21 60 if (simtype != "exo")
af3b84f4
BA
61 {
62 h_best_endo = optimize(
63 errorOnLastNdays, c(0,10), kernel=kernel, simtype="endo")$minimum
64 }
3d69ff21
BA
65
66 if (simtype == "endo")
af3b84f4 67 {
98e958ca 68 return (private$.predictShapeAux(data,
af3b84f4
BA
69 fdays, today, horizon, h_best_endo, kernel, "endo", TRUE))
70 }
3d69ff21 71 if (simtype == "exo")
af3b84f4 72 {
98e958ca 73 return (private$.predictShapeAux(data,
af3b84f4
BA
74 fdays, today, horizon, h_best_exo, kernel, "exo", TRUE))
75 }
3d69ff21
BA
76 if (simtype == "mix")
77 {
f17665c7 78 h_best_mix = c(h_best_endo,h_best_exo)
98e958ca 79 return(private$.predictShapeAux(data,
af3b84f4 80 fdays, today, horizon, h_best_mix, kernel, "mix", TRUE))
3d69ff21 81 }
25b75559
BA
82 }
83 ),
84 private = list(
3d69ff21 85 # Precondition: "today" is full (no NAs)
98e958ca 86 .predictShapeAux = function(data, fdays, today, horizon, h, kernel, simtype, final_call)
3d69ff21 87 {
f17665c7 88 fdays = fdays[ fdays < today ]
3d69ff21 89 # TODO: 3 = magic number
f17665c7 90 if (length(fdays) < 3)
3d69ff21
BA
91 return (NA)
92
93 if (simtype != "exo")
94 {
95 h_endo = ifelse(simtype=="mix", h[1], h)
96
97 # Distances from last observed day to days in the past
f17665c7
BA
98 distances2 = rep(NA, length(fdays))
99 for (i in seq_along(fdays))
3d69ff21 100 {
25b75559 101 delta = data$getCenteredSerie(today) - data$getCenteredSerie(fdays[i])
3d69ff21
BA
102 # Require at least half of non-NA common values to compute the distance
103 if (sum(is.na(delta)) <= 0) #length(delta)/2)
104 distances2[i] = mean(delta^2) #, na.rm=TRUE)
105 }
106
107 sd_dist = sd(distances2)
99f83c9a
BA
108 if (sd_dist < .Machine$double.eps)
109 sd_dist = 1 #mostly for tests... FIXME:
3d69ff21 110 simils_endo =
99f83c9a 111 if (kernel=="Gauss")
3d69ff21 112 exp(-distances2/(sd_dist*h_endo^2))
99f83c9a 113 else { #Epanechnikov
3d69ff21
BA
114 u = 1 - distances2/(sd_dist*h_endo^2)
115 u[abs(u)>1] = 0.
116 u
117 }
118 }
119
120 if (simtype != "endo")
121 {
122 h_exo = ifelse(simtype=="mix", h[2], h)
123
25b75559
BA
124 M = matrix( nrow=1+length(fdays), ncol=1+length(data$getExo(today)) )
125 M[1,] = c( data$getLevel(today), as.double(data$getExo(today)) )
f17665c7 126 for (i in seq_along(fdays))
25b75559 127 M[i+1,] = c( data$getLevel(fdays[i]), as.double(data$getExo(fdays[i])) )
3d69ff21
BA
128
129 sigma = cov(M) #NOTE: robust covariance is way too slow
613a986f 130 sigma_inv = solve(sigma) #TODO: use pseudo-inverse if needed?
3d69ff21
BA
131
132 # Distances from last observed day to days in the past
133 distances2 = rep(NA, nrow(M)-1)
134 for (i in 2:nrow(M))
135 {
136 delta = M[1,] - M[i,]
137 distances2[i-1] = delta %*% sigma_inv %*% delta
138 }
139
140 sd_dist = sd(distances2)
141 simils_exo =
f17665c7 142 if (kernel=="Gauss")
3d69ff21 143 exp(-distances2/(sd_dist*h_exo^2))
f17665c7 144 else { #Epanechnikov
3d69ff21
BA
145 u = 1 - distances2/(sd_dist*h_exo^2)
146 u[abs(u)>1] = 0.
147 u
148 }
149 }
150
3d69ff21 151 similarities =
f17665c7 152 if (simtype == "exo")
3d69ff21 153 simils_exo
f17665c7
BA
154 else if (simtype == "endo")
155 simils_endo
156 else #mix
157 simils_endo * simils_exo
3d69ff21
BA
158
159 prediction = rep(0, horizon)
a66a84b5
BA
160 for (i in seq_along(fdays))
161 prediction = prediction + similarities[i] * data$getSerie(fdays[i]+1)[1:horizon]
3d69ff21 162 prediction = prediction / sum(similarities, na.rm=TRUE)
99f83c9a 163
3d69ff21
BA
164 if (final_call)
165 {
a66a84b5
BA
166 private$.params$weights <- similarities
167 private$.params$indices <- fdays
168 private$.params$window <-
3d69ff21
BA
169 if (simtype=="endo") {
170 h_endo
171 } else if (simtype=="exo") {
172 h_exo
f17665c7 173 } else { #mix
3d69ff21
BA
174 c(h_endo,h_exo)
175 }
176 }
99f83c9a 177
3d69ff21
BA
178 return (prediction)
179 }
180 )
181)