X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=pkg%2FR%2FF_Neighbors.R;h=238274bd5f8b1202ab601c8ecea5aad83fc4b578;hb=98e958cab563866f8e00886b54336018a2e8bc97;hp=43a6a13a1d82d7d15d7ec9c6519aeb1084f9309c;hpb=f17665c7d3da672163779da686d9f4d1ebad31f9;p=talweg.git diff --git a/pkg/R/F_Neighbors.R b/pkg/R/F_Neighbors.R index 43a6a13..238274b 100644 --- a/pkg/R/F_Neighbors.R +++ b/pkg/R/F_Neighbors.R @@ -1,62 +1,30 @@ #' @include Forecaster.R #' -#' @title Neighbors Forecaster +#' Neighbors Forecaster #' -#' @description Predict tomorrow as a weighted combination of "futures of the past" days. -#' Inherits \code{\link{Forecaster}} -NeighborsForecaster = setRefClass( - Class = "NeighborsForecaster", - contains = "Forecaster", - - methods = list( - initialize = function(...) - { - callSuper(...) - }, - predictShape = function(today, memory, horizon, ...) +#' Predict tomorrow as a weighted combination of "futures of the past" days. +#' Inherits \code{\link{Forecaster}} +NeighborsForecaster = R6::R6Class("NeighborsForecaster", + inherit = Forecaster, + + public = list( + predictShape = function(data, today, memory, horizon, ...) { # (re)initialize computed parameters - params <<- list("weights"=NA, "indices"=NA, "window"=NA) + private$.params <- list("weights"=NA, "indices"=NA, "window"=NA) + + # Determine indices of no-NAs days followed by no-NAs tomorrows + 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)) - return (.predictShapeAux(fdays,today,horizon,list(...)$h_window,kernel,simtype,TRUE)) - - # HACK for test reports: complete some days with a few NAs, for nicer graphics - nas_in_serie = is.na(data$getSerie(today)) - if (any(nas_in_serie)) { - if (sum(nas_in_serie) >= length(nas_in_serie) / 2) - return (NA) - for (i in seq_along(nas_in_serie)) - { - if (nas_in_serie[i]) - { - #look left - left = i-1 - while (left>=1 && nas_in_serie[left]) - left = left-1 - #look right - right = i+1 - while (right<=length(nas_in_serie) && nas_in_serie[right]) - right = right+1 - #HACK: modify by-reference Data object... - data$data[[today]]$serie[i] <<- - if (left==0) data$data[[today]]$serie[right] - else if (right==0) data$data[[today]]$serie[left] - else (data$data[[today]]$serie[left] + data$data[[today]]$serie[right]) / 2. - } - } + return ( private$.predictShapeAux(data, + fdays, today, horizon, list(...)$h_window, kernel, simtype, TRUE) ) } - # Determine indices of no-NAs days followed by no-NAs tomorrows - first_day = max(today - memory, 1) - fdays = (first_day:(today-1))[ sapply(first_day:(today-1), function(i) { - !any(is.na(data$getSerie(i)) | is.na(data$getSerie(i+1))) - }) ] - # Indices of similar days for cross-validation; TODO: 45 = magic number sdays = getSimilarDaysIndices(today, limit=45, same_season=FALSE) @@ -68,36 +36,51 @@ NeighborsForecaster = setRefClass( for (i in intersect(fdays,sdays)) { # mix_strategy is never used here (simtype != "mix"), therefore left blank - prediction = .predictShapeAux(fdays, i, horizon, h, kernel, simtype, FALSE) + prediction = private$.predictShapeAux(data, + fdays, i, horizon, h, kernel, simtype, FALSE) if (!is.na(prediction[1])) { nb_jours = nb_jours + 1 - error = error + mean((data$getCenteredSerie(i+1)[1:horizon] - prediction)^2) + error = error + + mean((data$getCenteredSerie(i+1)[1:horizon] - prediction)^2) } } return (error / nb_jours) } if (simtype != "endo") - h_best_exo = optimize(errorOnLastNdays, c(0,10), kernel=kernel, simtype="exo")$minimum + { + h_best_exo = optimize( + errorOnLastNdays, c(0,10), kernel=kernel, simtype="exo")$minimum + } if (simtype != "exo") - h_best_endo = optimize(errorOnLastNdays, c(0,10), kernel=kernel, simtype="endo")$minimum + { + h_best_endo = optimize( + errorOnLastNdays, c(0,10), kernel=kernel, simtype="endo")$minimum + } if (simtype == "endo") - return (.predictShapeAux(fdays, today, horizon, h_best_endo, kernel, "endo", TRUE)) + { + return (private$.predictShapeAux(data, + fdays, today, horizon, h_best_endo, kernel, "endo", TRUE)) + } if (simtype == "exo") - return (.predictShapeAux(fdays, today, horizon, h_best_exo, kernel, "exo", TRUE)) + { + 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 (.predictShapeAux(fdays, today, horizon, h_best_mix, kernel, "mix", TRUE)) + return(private$.predictShapeAux(data, + fdays, today, horizon, h_best_mix, kernel, "mix", TRUE)) } - }, + } + ), + private = list( # Precondition: "today" is full (no NAs) - .predictShapeAux = function(fdays, today, horizon, h, kernel, simtype, final_call) + .predictShapeAux = function(data, fdays, today, horizon, h, kernel, simtype, final_call) { - dat = data$data #HACK: faster this way... - fdays = fdays[ fdays < today ] # TODO: 3 = magic number if (length(fdays) < 3) @@ -111,7 +94,7 @@ NeighborsForecaster = setRefClass( distances2 = rep(NA, length(fdays)) for (i in seq_along(fdays)) { - delta = dat[[today]]$serie - dat[[ fdays[i] ]]$serie + delta = data$getCenteredSerie(today) - data$getCenteredSerie(fdays[i]) # Require at least half of non-NA common values to compute the distance if (sum(is.na(delta)) <= 0) #length(delta)/2) distances2[i] = mean(delta^2) #, na.rm=TRUE) @@ -134,10 +117,10 @@ NeighborsForecaster = setRefClass( { h_exo = ifelse(simtype=="mix", h[2], h) - M = matrix( nrow=1+length(fdays), ncol=1+length(dat[[today]]$exo) ) - M[1,] = c( dat[[today]]$level, as.double(dat[[today]]$exo) ) + M = matrix( nrow=1+length(fdays), ncol=1+length(data$getExo(today)) ) + M[1,] = c( data$getLevel(today), as.double(data$getExo(today)) ) for (i in seq_along(fdays)) - M[i+1,] = c( dat[[ fdays[i] ]]$level, as.double(dat[[ fdays[i] ]]$exo) ) + M[i+1,] = c( data$getLevel(fdays[i]), as.double(data$getExo(fdays[i])) ) sigma = cov(M) #NOTE: robust covariance is way too slow sigma_inv = solve(sigma) #TODO: use pseudo-inverse if needed? @@ -170,15 +153,15 @@ NeighborsForecaster = setRefClass( simils_endo * simils_exo prediction = rep(0, horizon) - for (i in seq_along(fdays_indices)) - prediction = prediction + similarities[i] * dat[[ fdays_indices[i]+1 ]]$serie[1:horizon] + for (i in seq_along(fdays)) + prediction = prediction + similarities[i] * data$getSerie(fdays[i]+1)[1:horizon] prediction = prediction / sum(similarities, na.rm=TRUE) if (final_call) { - params$weights <<- similarities - params$indices <<- fdays_indices - params$window <<- + private$.params$weights <- similarities + private$.params$indices <- fdays + private$.params$window <- if (simtype=="endo") { h_endo } else if (simtype=="exo") {