X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=pkg%2FR%2FF_Neighbors.R;h=c9eda053820615485e19079b493fefa7329b1172;hb=eef545170c5a76b710184db6b695c15b20759177;hp=d889a34ce64b82c51889a5a81322a1a7dec25b2c;hpb=ee8b1b4e3c13f8dcf13a2c8da6a3bef1520c8252;p=talweg.git diff --git a/pkg/R/F_Neighbors.R b/pkg/R/F_Neighbors.R index d889a34..c9eda05 100644 --- a/pkg/R/F_Neighbors.R +++ b/pkg/R/F_Neighbors.R @@ -22,82 +22,119 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", fdays = getNoNA2(data, max(today-memory,1), today-1) # Get optional args - simtype = ifelse(hasArg("simtype"), list(...)$simtype, "mix") #or "endo", or "exo" - kernel = ifelse(hasArg("kernel"), list(...)$kernel, "Gauss") #or "Epan" - if (hasArg(h_window)) + local = ifelse(hasArg("local"), list(...)$local, TRUE) #same level + season? + simtype = ifelse(hasArg("simtype"), list(...)$simtype, "none") #or "endo", or "exo" + if (hasArg("window")) { return ( private$.predictShapeAux(data, - fdays, today, horizon, list(...)$h_window, kernel, simtype, TRUE) ) + fdays, today, horizon, local, list(...)$window, simtype, TRUE) ) } - # Indices of similar days for cross-validation; TODO: 45 = magic number - sdays = getSimilarDaysIndices(today, data, limit=45, same_season=FALSE) + # Indices of similar days for cross-validation; TODO: 20 = magic number + cv_days = getSimilarDaysIndices(today, data, limit=20, same_season=FALSE, + days_in=fdays) - cv_days = intersect(fdays,sdays) - # Limit to 20 most recent matching days (TODO: 20 == magic number) - cv_days = sort(cv_days,decreasing=TRUE)[1:min(20,length(cv_days))] - - # Function to optimize h : h |--> sum of prediction errors on last 45 "similar" days - errorOnLastNdays = function(h, kernel, simtype) + # Optimize h : h |--> sum of prediction errors on last N "similar" days + errorOnLastNdays = function(window, simtype) { error = 0 nb_jours = 0 for (i in seq_along(cv_days)) { # mix_strategy is never used here (simtype != "mix"), therefore left blank - prediction = private$.predictShapeAux(data, - fdays, cv_days[i], horizon, h, kernel, simtype, FALSE) + prediction = private$.predictShapeAux(data, fdays, cv_days[i], horizon, local, + window, simtype, FALSE) if (!is.na(prediction[1])) { nb_jours = nb_jours + 1 error = error + - mean((data$getCenteredSerie(cv_days[i]+1)[1:horizon] - prediction)^2) + mean((data$getSerie(cv_days[i]+1)[1:horizon] - prediction)^2) } } return (error / nb_jours) } - if (simtype != "endo") + # TODO: 7 == magic number + if (simtype=="endo" || simtype=="mix") { - h_best_exo = optimize( - errorOnLastNdays, c(0,10), kernel=kernel, simtype="exo")$minimum + best_window_endo = optimize( + errorOnLastNdays, c(0,7), simtype="endo")$minimum } - if (simtype != "exo") + if (simtype=="exo" || simtype=="mix") { - h_best_endo = optimize( - errorOnLastNdays, c(0,10), kernel=kernel, simtype="endo")$minimum + best_window_exo = optimize( + errorOnLastNdays, c(0,7), simtype="exo")$minimum } - if (simtype == "endo") - { - return (private$.predictShapeAux(data, - fdays, today, horizon, h_best_endo, kernel, "endo", TRUE)) - } - if (simtype == "exo") - { - return (private$.predictShapeAux(data, - fdays, today, horizon, h_best_exo, kernel, "exo", TRUE)) - } - if (simtype == "mix") - { - h_best_mix = c(h_best_endo,h_best_exo) - return(private$.predictShapeAux(data, - fdays, today, horizon, h_best_mix, kernel, "mix", TRUE)) - } + best_window = + if (simtype == "endo") + best_window_endo + else if (simtype == "exo") + best_window_exo + else if (simtype == "mix") + c(best_window_endo,best_window_exo) + else #none: value doesn't matter + 1 + + return(private$.predictShapeAux(data, fdays, today, horizon, local, + best_window, simtype, TRUE)) } ), private = list( # Precondition: "today" is full (no NAs) - .predictShapeAux = function(data, fdays, today, horizon, h, kernel, simtype, final_call) + .predictShapeAux = function(data, fdays, today, horizon, local, window, simtype, + final_call) { - fdays = fdays[ fdays < today ] - # TODO: 3 = magic number - if (length(fdays) < 3) + fdays_cut = fdays[ fdays < today ] + if (length(fdays_cut) <= 1) return (NA) - if (simtype != "exo") + if (local) { - h_endo = ifelse(simtype=="mix", h[1], h) + # Neighbors: days in "same season"; TODO: 60 == magic number... + fdays = getSimilarDaysIndices(today, data, limit=60, same_season=TRUE, + days_in=fdays_cut) + if (length(fdays) <= 1) + return (NA) + levelToday = data$getLevel(today) + distances = sapply(fdays, function(i) abs(data$getLevel(i)-levelToday)) + #TODO: 2, 10, 3, 12 magic numbers here... + dist_thresh = 2 + min_neighbs = min(10,length(fdays)) + repeat + { + same_pollution = (distances <= dist_thresh) + nb_neighbs = sum(same_pollution) + if (nb_neighbs >= min_neighbs) #will eventually happen + break + dist_thresh = dist_thresh + 3 + } + fdays = fdays[same_pollution] + max_neighbs = 12 + if (nb_neighbs > max_neighbs) + { + # Keep only max_neighbs closest neighbors + fdays = fdays[ + sort(distances[same_pollution],index.return=TRUE)$ix[1:max_neighbs] ] + } + if (length(fdays) == 1) #the other extreme... + { + if (final_call) + { + private$.params$weights <- 1 + private$.params$indices <- fdays + private$.params$window <- 1 + } + return ( data$getSerie(fdays[1])[1:horizon] ) #what else?! + } + } + else + fdays = fdays_cut #no conditioning + + if (simtype == "endo" || simtype == "mix") + { + # Compute endogen similarities using given window + window_endo = ifelse(simtype=="mix", window[1], window) # Distances from last observed day to days in the past serieToday = data$getSerie(today) @@ -107,26 +144,18 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", }) sd_dist = sd(distances2) - if (sd_dist < .Machine$double.eps) + if (sd_dist < .25 * sqrt(.Machine$double.eps)) { # warning("All computed distances are very close: stdev too small") sd_dist = 1 #mostly for tests... FIXME: } - simils_endo = - if (kernel=="Gauss") - exp(-distances2/(sd_dist*h_endo^2)) - else - { - # Epanechnikov - u = 1 - distances2/(sd_dist*h_endo^2) - u[abs(u)>1] = 0. - u - } + simils_endo = exp(-distances2/(sd_dist*window_endo^2)) } - if (simtype != "endo") + if (simtype == "exo" || simtype == "mix") { - h_exo = ifelse(simtype=="mix", h[2], h) + # Compute exogen similarities using given window + window_exo = ifelse(simtype=="mix", window[2], window) M = matrix( nrow=1+length(fdays), ncol=1+length(data$getExo(today)) ) M[1,] = c( data$getLevel(today), as.double(data$getExo(today)) ) @@ -153,16 +182,7 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", # warning("All computed distances are very close: stdev too small") sd_dist = 1 #mostly for tests... FIXME: } - simils_exo = - if (kernel=="Gauss") - exp(-distances2/(sd_dist*h_exo^2)) - else - { - # Epanechnikov - u = 1 - distances2/(sd_dist*h_exo^2) - u[abs(u)>1] = 0. - u - } + simils_exo = exp(-distances2/(sd_dist*window_exo^2)) } similarities = @@ -170,13 +190,15 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", simils_exo else if (simtype == "endo") simils_endo - else #mix + else if (simtype == "mix") simils_endo * simils_exo + else #none + rep(1, length(fdays)) + similarities = similarities / sum(similarities) prediction = rep(0, horizon) for (i in seq_along(fdays)) - prediction = prediction + similarities[i] * data$getCenteredSerie(fdays[i]+1)[1:horizon] - prediction = prediction / sum(similarities, na.rm=TRUE) + prediction = prediction + similarities[i] * data$getSerie(fdays[i]+1)[1:horizon] if (final_call) { @@ -184,11 +206,13 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", private$.params$indices <- fdays private$.params$window <- if (simtype=="endo") - h_endo + window_endo else if (simtype=="exo") - h_exo - else #mix - c(h_endo,h_exo) + window_exo + else if (simtype=="mix") + c(window_endo,window_exo) + else #none + 1 } return (prediction)