X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=pkg%2FR%2FF_Neighbors2.R;h=fb63e4087f1b2aa4b426712d140132a2f2e83193;hb=5e838b3e17465c376ca075b766cf2543c82e9864;hp=e6adddea7a03cd5d0f21cf6ee84dfc26b40a1ab8;hpb=5c49f6cecd547358b327e9363e62bcc8219e9e33;p=talweg.git diff --git a/pkg/R/F_Neighbors2.R b/pkg/R/F_Neighbors2.R index e6addde..fb63e40 100644 --- a/pkg/R/F_Neighbors2.R +++ b/pkg/R/F_Neighbors2.R @@ -9,6 +9,15 @@ Neighbors2Forecaster = R6::R6Class("Neighbors2Forecaster", inherit = Forecaster, public = list( + predictSerie = function(data, today, memory, horizon, ...) + { + # Parameters (potentially) computed during shape prediction stage + predicted_shape = self$predictShape(data, today, memory, horizon, ...) +# predicted_delta = private$.pjump(data,today,memory,horizon,private$.params,...) + # Predicted shape is aligned it on the end of current day + jump +# predicted_shape+tail(data$getSerie(today),1)-predicted_shape[1]+predicted_delta + predicted_shape + }, predictShape = function(data, today, memory, horizon, ...) { # (re)initialize computed parameters @@ -22,101 +31,159 @@ Neighbors2Forecaster = R6::R6Class("Neighbors2Forecaster", 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 ( private$.predictShapeAux(data, - fdays, today, horizon, list(...)$h_window, kernel, TRUE) ) + fdays, today, horizon, list(...)$h_window, kernel, simtype, TRUE) ) } - # Indices of similar days for cross-validation; TODO: 45 = magic number - # TODO: ici faut une sorte de "same_season==TRUE" --> mois similaires epandage sdays = getSimilarDaysIndices(today, limit=45, same_season=FALSE) + 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) + errorOnLastNdays = function(h, kernel, simtype) { error = 0 nb_jours = 0 - for (i in intersect(fdays,sdays)) + for (i in seq_along(cv_days)) { # mix_strategy is never used here (simtype != "mix"), therefore left blank - prediction = private$.predictShapeAux(data, fdays, i, horizon, h, kernel, FALSE) + prediction = private$.predictShapeAux(data, + fdays, cv_days[i], horizon, h, kernel, simtype, FALSE) if (!is.na(prediction[1])) { nb_jours = nb_jours + 1 error = error + - mean((data$getSerie(i+1)[1:horizon] - prediction)^2) + mean((data$getSerie(cv_days[i]+1)[1:horizon] - prediction)^2) } } return (error / nb_jours) } - # h :: only for endo in this variation - h_best_endo = optimize(errorOnLastNdays, c(0,10), kernel=kernel)$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, c(0,10), kernel=kernel, simtype="endo")$minimum + } - return (private$.predictShapeAux(data, fdays, today, horizon, h_best, kernel, TRUE)) + if (simtype == "endo") + { + return (private$.predictShapeAux(data, + fdays, today, horizon, h_best_endo, kernel, "endo", TRUE)) + } + if (simtype == "exo") + { + 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(private$.predictShapeAux(data, + fdays, today, horizon, h_best_mix, kernel, "mix", TRUE)) + } } ), private = list( # Precondition: "today" is full (no NAs) - .predictShapeAux = function(data, fdays, today, horizon, h, kernel, final_call) + .predictShapeAux = function(data, fdays, today, horizon, h, kernel, simtype, final_call) { fdays = fdays[ fdays < today ] # TODO: 3 = magic number if (length(fdays) < 3) return (NA) - # ENDO:: Distances from last observed day to days in the past - distances2 = rep(NA, length(fdays)) - for (i in seq_along(fdays)) + # Neighbors: days in "same season" + sdays = getSimilarDaysIndices(today, limit=45, same_season=TRUE, data) + indices = intersect(fdays,sdays) + levelToday = data$getLevel(today) + distances = sapply(indices, function(i) abs(data$getLevel(i)-levelToday)) + same_pollution = (distances <= 2) + if (sum(same_pollution) < 3) #TODO: 3 == magic number { - delta = data$getSerie(today) - data$getSerie(fdays[i]) - # Require at least half of non-NA common values to compute the distance - if ( !any( is.na(delta) ) ) - distances2[i] = mean(delta^2) + same_pollution = (distances <= 5) + if (sum(same_pollution) < 3) + return (NA) } + indices = indices[same_pollution] - sd_dist = sd(distances2) - if (sd_dist < .Machine$double.eps) + if (simtype != "exo") { -# 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 - } + h_endo = ifelse(simtype=="mix", h[1], h) - # EXOGENS: distances computations are enough - # TODO: search among similar concentrations....... at this stage ?! - 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( data$getLevel(fdays[i]), as.double(data$getExo(fdays[i])) ) + # Distances from last observed day to days in the past + serieToday = data$getSerie(today) + distances2 = sapply(indices, function(i) { + delta = serieToday - data$getSerie(i) + mean(delta^2) + }) - sigma = cov(M) #NOTE: robust covariance is way too slow - sigma_inv = solve(sigma) #TODO: use pseudo-inverse if needed? + sd_dist = sd(distances2) + if (sd_dist < .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 + } + } - # Distances from last observed day to days in the past - distances2 = rep(NA, nrow(M)-1) - for (i in 2:nrow(M)) + if (simtype != "endo") { - delta = M[1,] - M[i,] - distances2[i-1] = delta %*% sigma_inv %*% delta + h_exo = ifelse(simtype=="mix", h[2], h) + + M = matrix( nrow=1+length(indices), ncol=1+length(data$getExo(today)) ) + M[1,] = c( data$getLevel(today), as.double(data$getExo(today)) ) + for (i in seq_along(indices)) + M[i+1,] = c( data$getLevel(indices[i]), as.double(data$getExo(indices[i])) ) + + sigma = cov(M) #NOTE: robust covariance is way too slow +# sigma_inv = solve(sigma) #TODO: use pseudo-inverse if needed? + sigma_inv = MASS::ginv(sigma) +#if (final_call) browser() + # Distances from last observed day to days in the past + distances2 = sapply(seq_along(indices), function(i) { + delta = M[1,] - M[i+1,] + delta %*% sigma_inv %*% delta + }) + + sd_dist = sd(distances2) + 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_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 + } } - ppv <- sort(distances2, index.return=TRUE)$ix[1:10] #.............. -#PPV pour endo ? - similarities = if (simtype == "exo") simils_exo @@ -126,8 +193,8 @@ Neighbors2Forecaster = R6::R6Class("Neighbors2Forecaster", simils_endo * simils_exo prediction = rep(0, horizon) - for (i in seq_along(fdays)) - prediction = prediction + similarities[i] * data$getSerie(fdays[i]+1)[1:horizon] + for (i in seq_along(indices)) + prediction = prediction + similarities[i] * data$getSerie(indices[i]+1)[1:horizon] prediction = prediction / sum(similarities, na.rm=TRUE) if (final_call)