X-Git-Url: https://git.auder.net/?p=talweg.git;a=blobdiff_plain;f=pkg%2FR%2FF_Neighbors.R;h=43a6a13a1d82d7d15d7ec9c6519aeb1084f9309c;hp=ffb6d371df57a910451f0a79ffeed340886132ff;hb=f17665c7d3da672163779da686d9f4d1ebad31f9;hpb=42916467fa0d13aba2f345682d877adc11e1a77b diff --git a/pkg/R/F_Neighbors.R b/pkg/R/F_Neighbors.R index ffb6d37..43a6a13 100644 --- a/pkg/R/F_Neighbors.R +++ b/pkg/R/F_Neighbors.R @@ -18,16 +18,16 @@ NeighborsForecaster = setRefClass( # (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 + # 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)) - # Predict only on (almost) non-NAs days + # 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)) { - #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)) @@ -52,88 +52,55 @@ NeighborsForecaster = setRefClass( } # 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, "mix") - simthresh = ifelse(hasArg("simthresh"), list(...)$simthresh, 0.) - kernel = ifelse(hasArg("kernel"), list(...)$kernel, "Gauss") #or "Epan" - mix_strategy = ifelse(hasArg("mix_strategy"), list(...)$mix_strategy, "mult") #or "neighb" - same_season = ifelse(hasArg("same_season"), list(...)$same_season, FALSE) - if (hasArg(h_window)) - return (.predictShapeAux(fdays_indices, today, horizon, list(...)$h_window, kernel, - simtype, simthresh, mix_strategy, TRUE)) - #END GET + 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 for cross-validation; TODO: 45 = magic number - indices = getSimilarDaysIndices(today, limit=45, same_season=same_season) - if (tail(indices,1) == 1) - indices = head(indices,-1) + # Indices of similar days for cross-validation; TODO: 45 = magic number + sdays = getSimilarDaysIndices(today, limit=45, same_season=FALSE) # 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) + for (i in intersect(fdays,sdays)) { - # 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)))) + # mix_strategy is never used here (simtype != "mix"), therefore left blank + prediction = .predictShapeAux(fdays, i, horizon, h, kernel, simtype, FALSE) + if (!is.na(prediction[1])) { nb_jours = nb_jours + 1 - # mix_strategy is never used here (simtype != "mix"), therefore left blank - prediction = .predictShapeAux(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) + 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 != "endo") + h_best_exo = optimize(errorOnLastNdays, 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 - } + h_best_endo = optimize(errorOnLastNdays, c(0,10), kernel=kernel, simtype="endo")$minimum if (simtype == "endo") - { - return (.predictShapeAux(fdays_indices, today, horizon, h_best_endo, kernel, "endo", - simthresh, "", TRUE)) - } + return (.predictShapeAux(fdays, today, horizon, h_best_endo, kernel, "endo", TRUE)) if (simtype == "exo") - { - return (.predictShapeAux(fdays_indices, today, horizon, h_best_exo, kernel, "exo", - simthresh, "", TRUE)) - } + return (.predictShapeAux(fdays, today, horizon, h_best_exo, kernel, "exo", TRUE)) if (simtype == "mix") { - return (.predictShapeAux(fdays_indices, today, horizon, c(h_best_endo,h_best_exo), - kernel, "mix", simthresh, mix_strategy, TRUE)) + h_best_mix = c(h_best_endo,h_best_exo) + return (.predictShapeAux(fdays, today, horizon, h_best_mix, kernel, "mix", TRUE)) } }, # Precondition: "today" is full (no NAs) - .predictShapeAux = function(fdays_indices, today, horizon, h, kernel, simtype, simthresh, - mix_strategy, final_call) + .predictShapeAux = function(fdays, today, horizon, h, kernel, simtype, final_call) { dat = data$data #HACK: faster this way... - fdays_indices = fdays_indices[fdays_indices < today] + fdays = fdays[ fdays < today ] # TODO: 3 = magic number - if (length(fdays_indices) < 3) + if (length(fdays) < 3) return (NA) if (simtype != "exo") @@ -141,10 +108,10 @@ NeighborsForecaster = setRefClass( 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)) + distances2 = rep(NA, length(fdays)) + for (i in seq_along(fdays)) { - delta = dat[[today]]$serie - dat[[ fdays_indices[i] ]]$serie + delta = dat[[today]]$serie - dat[[ fdays[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) @@ -167,13 +134,10 @@ NeighborsForecaster = setRefClass( { h_exo = ifelse(simtype=="mix", h[2], h) - M = matrix( nrow=1+length(fdays_indices), ncol=1+length(dat[[today]]$exo) ) + M = matrix( nrow=1+length(fdays), ncol=1+length(dat[[today]]$exo) ) M[1,] = c( dat[[today]]$level, as.double(dat[[today]]$exo) ) - for (i in seq_along(fdays_indices)) - { - M[i+1,] = c( dat[[ fdays_indices[i] ]]$level, - as.double(dat[[ fdays_indices[i] ]]$exo) ) - } + for (i in seq_along(fdays)) + M[i+1,] = c( dat[[ fdays[i] ]]$level, as.double(dat[[ fdays[i] ]]$exo) ) sigma = cov(M) #NOTE: robust covariance is way too slow sigma_inv = solve(sigma) #TODO: use pseudo-inverse if needed? @@ -188,51 +152,22 @@ NeighborsForecaster = setRefClass( sd_dist = sd(distances2) simils_exo = - if (kernel=="Gauss") { + if (kernel=="Gauss") exp(-distances2/(sd_dist*h_exo^2)) - } else { #Epanechnikov + 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 { + if (simtype == "exo") 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. - } - } + else if (simtype == "endo") + simils_endo + else #mix + simils_endo * simils_exo prediction = rep(0, horizon) for (i in seq_along(fdays_indices)) @@ -248,7 +183,7 @@ NeighborsForecaster = setRefClass( h_endo } else if (simtype=="exo") { h_exo - } else { + } else { #mix c(h_endo,h_exo) } }