X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=pkg%2FR%2FF_Neighbors.R;h=8eb1ddc240be802f50b52afeb63ac95a373676d1;hb=7c4b2952874de1d40a742e72efe51999b99050f5;hp=f140b0bb5559d1f51852e55d7ea7d98aa02c1be6;hpb=cf3bb00128ac8cb930996455faf7c99a3fc102fb;p=talweg.git diff --git a/pkg/R/F_Neighbors.R b/pkg/R/F_Neighbors.R index f140b0b..8eb1ddc 100644 --- a/pkg/R/F_Neighbors.R +++ b/pkg/R/F_Neighbors.R @@ -1,27 +1,23 @@ #' Neighbors Forecaster #' -#' Predict next serie as a weighted combination of "futures of the past" days, -#' where days in the past are chosen and weighted according to some similarity measures. +#' Predict next serie as a weighted combination of curves observed on "similar" days in +#' the past (and future if 'opera'=FALSE); the nature of the similarity is controlled by +#' the options 'simtype' and 'local' (see below). #' -#' The main method is \code{predictShape()}, taking arguments data, today, memory, -#' predict_from, horizon respectively for the dataset (object output of -#' \code{getData()}), the current index, the data depth (in days), the first predicted -#' hour and the last predicted hour. -#' In addition, optional arguments can be passed: +#' Optional arguments: #' \itemize{ -#' \item local : TRUE (default) to constrain neighbors to be "same days within same -#' season" -#' \item simtype : 'endo' for a similarity based on the series only, +#' \item local: TRUE (default) to constrain neighbors to be "same days in same season" +#' \item simtype: 'endo' for a similarity based on the series only, #' 'exo' for a similarity based on exogenous variables only, #' 'mix' for the product of 'endo' and 'exo', #' 'none' (default) to apply a simple average: no computed weights -#' \item window : A window for similarities computations; override cross-validation +#' \item window: A window for similarities computations; override cross-validation #' window estimation. #' } #' The method is summarized as follows: #' \enumerate{ -#' \item Determine N (=20) recent days without missing values, and followed by a -#' tomorrow also without missing values. +#' \item Determine N (=20) recent days without missing values, and preceded by a +#' curve also without missing values. #' \item Optimize the window parameters (if relevant) on the N chosen days. #' \item Considering the optimized window, compute the neighbors (with locality #' constraint or not), compute their similarities -- using a gaussian kernel if @@ -51,21 +47,26 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", return (NA) } - # Determine indices of no-NAs days preceded by no-NAs yerstedays - tdays = .getNoNA2(data, max(today-memory,2), today-1) - # Get optional args local = ifelse(hasArg("local"), list(...)$local, TRUE) #same level + season? simtype = ifelse(hasArg("simtype"), list(...)$simtype, "none") #or "endo", or "exo" + opera = ifelse(hasArg("opera"), list(...)$opera, FALSE) #operational mode? + + # Determine indices of no-NAs days preceded by no-NAs yerstedays + tdays = .getNoNA2(data, max(today-memory,2), ifelse(opera,today-1,data$getSize())) + if (!opera) + tdays = setdiff(tdays, today) #always exclude current day + + # Shortcut if window is known #TODO: cross-validation for number of days, on similar (yerste)days if (hasArg("window")) { - return ( private$.predictShapeAux(data, - tdays, today, predict_from, horizon, local, list(...)$window, simtype, TRUE) ) + return ( private$.predictShapeAux(data, tdays, today, predict_from, horizon, + local, list(...)$window, simtype, opera, TRUE) ) } # Indices of similar days for cross-validation; TODO: 20 = magic number cv_days = getSimilarDaysIndices(today, data, limit=20, same_season=FALSE, - days_in=tdays) + days_in=tdays, operational=opera) # Optimize h : h |--> sum of prediction errors on last N "similar" days errorOnLastNdays = function(window, simtype) @@ -76,7 +77,7 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", { # mix_strategy is never used here (simtype != "mix"), therefore left blank prediction = private$.predictShapeAux(data, tdays, cv_days[i], predict_from, - horizon, local, window, simtype, FALSE) + horizon, local, window, simtype, opera, FALSE) if (!is.na(prediction[1])) { nb_jours = nb_jours + 1 @@ -110,27 +111,25 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", 1 return( private$.predictShapeAux(data, tdays, today, predict_from, horizon, local, - best_window, simtype, TRUE) ) + best_window, simtype, opera, TRUE) ) } ), private = list( - # Precondition: "today" is full (no NAs) + # Precondition: "yersteday until predict_from-1" is full (no NAs) .predictShapeAux = function(data, tdays, today, predict_from, horizon, local, window, - simtype, final_call) + simtype, opera, final_call) { - tdays_cut = tdays[ tdays <= today-1 ] - if (length(tdays_cut) <= 1) + tdays_cut = tdays[ tdays != today ] + if (length(tdays_cut) == 0) return (NA) if (local) { - # TODO: 60 == magic number - tdays = getSimilarDaysIndices(today, data, limit=60, same_season=TRUE, - days_in=tdays_cut) - if (length(tdays) <= 1) - return (NA) - # TODO: 10, 12 == magic numbers - tdays = .getConstrainedNeighbs(today,data,tdays,min_neighbs=10,max_neighbs=12) + # limit=Inf to not censor any day (TODO: finite limit? 60?) + tdays = getSimilarDaysIndices(today, data, limit=Inf, same_season=TRUE, + days_in=tdays_cut, operational=opera) + # TODO: 10 == magic number + tdays = .getConstrainedNeighbs(today, data, tdays, min_neighbs=10) if (length(tdays) == 1) { if (final_call) @@ -141,6 +140,13 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", } return ( data$getSerie(tdays[1])[predict_from:horizon] ) } + max_neighbs = 12 #TODO: 10 or 12 or... ? + if (length(tdays) > max_neighbs) + { + distances2 <- .computeDistsEndo(data, today, tdays, predict_from) + ordering <- order(distances2) + tdays <- tdays[ ordering[1:max_neighbs] ] + } } else tdays = tdays_cut #no conditioning @@ -150,14 +156,9 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", # Compute endogen similarities using given window window_endo = ifelse(simtype=="mix", window[1], window) - # Distances from last observed day to days in the past - lastSerie = c( data$getSerie(today-1), - data$getSerie(today)[if (predict_from>=2) 1:(predict_from-1) else c()] ) - distances2 = sapply(tdays, function(i) { - delta = lastSerie - c(data$getSerie(i-1), - data$getSerie(i)[if (predict_from>=2) 1:(predict_from-1) else c()]) - sqrt(mean(delta^2)) - }) + # Distances from last observed day to selected days in the past + # TODO: redundant computation if local==TRUE + distances2 <- .computeDistsEndo(data, today, tdays, predict_from) simils_endo <- .computeSimils(distances2, window_endo) } @@ -167,24 +168,7 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", # Compute exogen similarities using given window window_exo = ifelse(simtype=="mix", window[2], window) - M = matrix( ncol=1+length(tdays), nrow=1+length(data$getExo(1)) ) - M[,1] = c( data$getLevelHat(today), as.double(data$getExoHat(today)) ) - for (i in seq_along(tdays)) - M[,i+1] = c( data$getLevel(tdays[i]), as.double(data$getExo(tdays[i])) ) - - sigma = cov(t(M)) #NOTE: robust covariance is way too slow - # TODO: 10 == magic number; more robust way == det, or always ginv() - sigma_inv = - if (length(tdays) > 10) - solve(sigma) - else - MASS::ginv(sigma) - - # Distances from last observed day to days in the past - distances2 = sapply(seq_along(tdays), function(i) { - delta = M[,1] - M[,i+1] - delta %*% sigma_inv %*% delta - }) + distances2 <- .computeDistsExo(data, today, tdays) simils_exo <- .computeSimils(distances2, window_exo) } @@ -237,13 +221,10 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", # @param min_neighbs Minimum number of points in a neighborhood # @param max_neighbs Maximum number of points in a neighborhood # -.getConstrainedNeighbs = function(today, data, tdays, min_neighbs=10, max_neighbs=12) +.getConstrainedNeighbs = function(today, data, tdays, min_neighbs=10) { levelToday = data$getLevelHat(today) - levelYersteday = data$getLevel(today-1) - distances = sapply(tdays, function(i) { - sqrt((data$getLevel(i-1)-levelYersteday)^2 + (data$getLevel(i)-levelToday)^2) - }) + distances = sapply( tdays, function(i) abs(data$getLevel(i) - levelToday) ) #TODO: 1, +1, +3 : magic numbers dist_thresh = 1 min_neighbs = min(min_neighbs,length(tdays)) @@ -255,14 +236,7 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", break dist_thresh = dist_thresh + ifelse(dist_thresh>1,3,1) } - tdays = tdays[same_pollution] - max_neighbs = 12 - if (nb_neighbs > max_neighbs) - { - # Keep only max_neighbs closest neighbors - tdays = tdays[ order(distances[same_pollution])[1:max_neighbs] ] - } - tdays + tdays[same_pollution] } # compute similarities @@ -282,3 +256,36 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", } exp(-distances2/(sd_dist*window^2)) } + +.computeDistsEndo <- function(data, today, tdays, predict_from) +{ + lastSerie = c( data$getSerie(today-1), + data$getSerie(today)[if (predict_from>=2) 1:(predict_from-1) else c()] ) + sapply(tdays, function(i) { + delta = lastSerie - c(data$getSerie(i-1), + data$getSerie(i)[if (predict_from>=2) 1:(predict_from-1) else c()]) + sqrt(mean(delta^2)) + }) +} + +.computeDistsExo <- function(data, today, tdays) +{ + M = matrix( ncol=1+length(tdays), nrow=1+length(data$getExo(1)) ) + M[,1] = c( data$getLevelHat(today), as.double(data$getExoHat(today)) ) + for (i in seq_along(tdays)) + M[,i+1] = c( data$getLevel(tdays[i]), as.double(data$getExo(tdays[i])) ) + + sigma = cov(t(M)) #NOTE: robust covariance is way too slow + # TODO: 10 == magic number; more robust way == det, or always ginv() + sigma_inv = + if (length(tdays) > 10) + solve(sigma) + else + MASS::ginv(sigma) + + # Distances from last observed day to days in the past + sapply(seq_along(tdays), function(i) { + delta = M[,1] - M[,i+1] + delta %*% sigma_inv %*% delta + }) +}