X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=R%2FS_Neighbors.R;fp=R%2FS_Neighbors.R;h=0000000000000000000000000000000000000000;hb=e030a6e31232332b73187eda25870e843152c174;hp=b8e32cc9df784852099e8596452e6cd60190981c;hpb=31f7d913d4a99d0a4db9bcfe40e31cebf90b22e6;p=talweg.git diff --git a/R/S_Neighbors.R b/R/S_Neighbors.R deleted file mode 100644 index b8e32cc..0000000 --- a/R/S_Neighbors.R +++ /dev/null @@ -1,257 +0,0 @@ -#' @include ShapeForecaster.R -#' -#' @title Neighbors Shape Forecaster -#' -#' @description Predict tomorrow as a weighted combination of "futures of the past" days. -#' Inherits \code{\link{ShapeForecaster}} -NeighborsShapeForecaster = setRefClass( - Class = "NeighborsShapeForecaster", - contains = "ShapeForecaster", - - methods = list( - initialize = function(...) - { - callSuper(...) - }, - predict = function(today, memory, horizon, ...) - { - # (re)initialize computed parameters - params <<- list("weights"=NA, "indices"=NA, "window"=NA) - - first_day = max(today - memory, 1) - # The first day is generally not complete: - if (length(data$getCenteredSerie(1)) < length(data$getCenteredSerie(2))) - first_day = 2 - - # Predict only on (almost) non-NAs days - nas_in_serie = is.na(data$getSerie(today)) - if (any(nas_in_serie)) - { - #TODO: better define "repairing" conditions (and method) - 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. - } - } - } - - # Determine indices of no-NAs days followed by no-NAs tomorrows - fdays_indices = c() - for (i in first_day:(today-1)) - { - if ( !any(is.na(data$getSerie(i)) | is.na(data$getSerie(i+1))) ) - fdays_indices = c(fdays_indices, i) - } - - #GET OPTIONAL PARAMS - # Similarity computed with exogenous variables ? endogenous ? both ? ("exo","endo","mix") - simtype = ifelse(hasArg("simtype"), list(...)$simtype, "exo") - simthresh = ifelse(hasArg("simthresh"), list(...)$simthresh, 0.) - kernel = ifelse(hasArg("kernel"), list(...)$kernel, "Gauss") - mix_strategy = ifelse(hasArg("mix_strategy"), list(...)$mix_strategy, "neighb") #or "mult" - same_season = ifelse(hasArg("same_season"), list(...)$same_season, TRUE) - if (hasArg(h_window)) - return (.predictAux(fdays_indices, today, horizon, list(...)$h_window, kernel, simtype, - simthresh, mix_strategy, FALSE)) - #END GET - - # Indices for cross-validation; TODO: 45 = magic number - indices = getSimilarDaysIndices(today, limit=45, same_season=same_season) - #indices = (end_index-45):(end_index-1) - - # Function to optimize h : h |--> sum of prediction errors on last 45 "similar" days - errorOnLastNdays = function(h, kernel, simtype) - { - error = 0 - nb_jours = 0 - for (i in indices) - { - # NOTE: predict only on non-NAs days followed by non-NAs (TODO:) - if (!any(is.na(data$getSerie(i)) | is.na(data$getSerie(i+1)))) - { - nb_jours = nb_jours + 1 - # mix_strategy is never used here (simtype != "mix"), therefore left blank - prediction = .predictAux(fdays_indices, i, horizon, h, kernel, simtype, simthresh, - "", FALSE) - if (!is.na(prediction[1])) - error = error + mean((data$getCenteredSerie(i+1)[1:horizon] - prediction)^2) - } - } - return (error / nb_jours) - } - - h_best_exo = 1. - if (simtype != "endo" && !(simtype=="mix" && mix_strategy=="neighb")) - { - h_best_exo = optimize(errorOnLastNdays, interval=c(0,10), kernel=kernel, - simtype="exo")$minimum - } - if (simtype != "exo") - { - h_best_endo = optimize(errorOnLastNdays, interval=c(0,10), kernel=kernel, - simtype="endo")$minimum - } - - if (simtype == "endo") - { - return (.predictAux(fdays_indices, today, horizon, h_best_endo, kernel, "endo", - simthresh, "", TRUE)) - } - if (simtype == "exo") - { - return (.predictAux(fdays_indices, today, horizon, h_best_exo, kernel, "exo", - simthresh, "", TRUE)) - } - if (simtype == "mix") - { - return (.predictAux(fdays_indices, today, horizon, c(h_best_endo,h_best_exo), kernel, - "mix", simthresh, mix_strategy, TRUE)) - } - }, - # Precondition: "today" is full (no NAs) - .predictAux = function(fdays_indices, today, horizon, h, kernel, simtype, simthresh, - mix_strategy, final_call) - { - dat = data$data #HACK: faster this way... - - fdays_indices = fdays_indices[fdays_indices < today] - # TODO: 3 = magic number - if (length(fdays_indices) < 3) - return (NA) - - if (simtype != "exo") - { - h_endo = ifelse(simtype=="mix", h[1], h) - - # Distances from last observed day to days in the past - distances2 = rep(NA, length(fdays_indices)) - for (i in seq_along(fdays_indices)) - { - delta = dat[[today]]$serie - dat[[ fdays_indices[i] ]]$serie - # 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) - } - - sd_dist = sd(distances2) - 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 - } - } - - if (simtype != "endo") - { - h_exo = ifelse(simtype=="mix", h[2], h) - - # TODO: [rnormand] if predict_at == 0h then we should use measures from day minus 1 - M = matrix( nrow=1+length(fdays_indices), ncol=1+length(dat[[today]]$exo_hat) ) - M[1,] = c( dat[[today]]$level, as.double(dat[[today]]$exo_hat) ) - for (i in seq_along(fdays_indices)) - { - M[i+1,] = c( dat[[ fdays_indices[i] ]]$level, - as.double(dat[[ fdays_indices[i] ]]$exo_hat) ) - } - - sigma = cov(M) #NOTE: robust covariance is way too slow - sigma_inv = qr.solve(sigma) - - # Distances from last observed day to days in the past - distances2 = rep(NA, nrow(M)-1) - for (i in 2:nrow(M)) - { - delta = M[1,] - M[i,] - distances2[i-1] = delta %*% sigma_inv %*% delta - } - - sd_dist = sd(distances2) - 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 - } - } - - if (simtype=="mix") - { - if (mix_strategy == "neighb") - { - #Only (60) most similar days according to exogen variables are kept into consideration - #TODO: 60 = magic number - keep_indices = sort(simils_exo, index.return=TRUE)$ix[1:(min(60,length(simils_exo)))] - simils_endo[-keep_indices] = 0. - } else #mix_strategy == "mult" - { - simils_endo = simils_endo * simils_exo - } - } - - similarities = - if (simtype != "exo") { - simils_endo - } else { - simils_exo - } - - if (simthresh > 0.) - { - max_sim = max(similarities) - # Set to 0 all similarities s where s / max_sim < simthresh, but keep at least 60 - ordering = sort(similarities / max_sim, index.return=TRUE) - if (ordering[60] < simthresh) - { - similarities[ ordering$ix[ - (1:60) ] ] = 0. - } else - { - limit = 61 - while (limit < length(similarities) && ordering[limit] >= simthresh) - limit = limit + 1 - similarities[ ordering$ix[ - 1:limit] ] = 0. - } - } - - prediction = rep(0, horizon) - for (i in seq_along(fdays_indices)) - prediction = prediction + similarities[i] * dat[[ fdays_indices[i]+1 ]]$serie[1:horizon] - - prediction = prediction / sum(similarities, na.rm=TRUE) - if (final_call) - { - params$weights <<- similarities - params$indices <<- fdays_indices - params$window <<- - if (simtype=="endo") { - h_endo - } else if (simtype=="exo") { - h_exo - } else { - c(h_endo,h_exo) - } - } - return (prediction) - } - ) -)