X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=pkg%2FR%2FF_Neighbors.R;h=c55291aa996838a2a02535725f72c4fb8709b922;hb=aa059de77cbcd28a3a66c7ff29ebe0346882867b;hp=d889a34ce64b82c51889a5a81322a1a7dec25b2c;hpb=ee8b1b4e3c13f8dcf13a2c8da6a3bef1520c8252;p=talweg.git diff --git a/pkg/R/F_Neighbors.R b/pkg/R/F_Neighbors.R index d889a34..c55291a 100644 --- a/pkg/R/F_Neighbors.R +++ b/pkg/R/F_Neighbors.R @@ -22,36 +22,33 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", fdays = getNoNA2(data, max(today-memory,1), today-1) # Get optional args + local = ifelse(hasArg("local"), list(...)$local, FALSE) #same level + season? simtype = ifelse(hasArg("simtype"), list(...)$simtype, "mix") #or "endo", or "exo" - kernel = ifelse(hasArg("kernel"), list(...)$kernel, "Gauss") #or "Epan" - if (hasArg(h_window)) + 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 45 "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) @@ -59,45 +56,87 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", if (simtype != "endo") { - h_best_exo = optimize( - errorOnLastNdays, c(0,10), kernel=kernel, simtype="exo")$minimum + best_window_exo = optimize( + errorOnLastNdays, c(0,7), simtype="exo")$minimum } if (simtype != "exo") { - h_best_endo = optimize( - errorOnLastNdays, c(0,10), kernel=kernel, simtype="endo")$minimum + best_window_endo = optimize( + errorOnLastNdays, c(0,7), simtype="endo")$minimum } if (simtype == "endo") { - return (private$.predictShapeAux(data, - fdays, today, horizon, h_best_endo, kernel, "endo", TRUE)) + return (private$.predictShapeAux(data, fdays, today, horizon, local, + best_window_endo, "endo", TRUE)) } if (simtype == "exo") { - return (private$.predictShapeAux(data, - fdays, today, horizon, h_best_exo, kernel, "exo", TRUE)) + return (private$.predictShapeAux(data, fdays, today, horizon, local, + best_window_exo, "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)) + return(private$.predictShapeAux(data, fdays, today, horizon, local, + c(best_window_endo,best_window_exo), "mix", 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 (local) + { + # 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, 3, 5, 10 magic numbers here... + dist_thresh = 2 + min_neighbs = min(3,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 = 10 + 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 != "exo") { - h_endo = ifelse(simtype=="mix", h[1], h) + # 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 +146,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") { - h_exo = ifelse(simtype=="mix", h[2], h) + # Compute exogen similarities using given window + h_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 +184,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 = @@ -172,11 +194,11 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", simils_endo else #mix simils_endo * simils_exo + 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,11 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", private$.params$indices <- fdays private$.params$window <- if (simtype=="endo") - h_endo + window_endo else if (simtype=="exo") - h_exo + window_exo else #mix - c(h_endo,h_exo) + c(window_endo,window_exo) } return (prediction)