From: Benjamin Auder Date: Wed, 29 Mar 2017 17:22:21 +0000 (+0200) Subject: finished merging F_Neighbors.R; TODO: test X-Git-Url: https://git.auder.net/%7B%7B%20asset%28%27mixstore/css/static/current/%7B%7B%20targetUrl%20%7D%7D?a=commitdiff_plain;h=aa059de77cbcd28a3a66c7ff29ebe0346882867b;p=talweg.git finished merging F_Neighbors.R; TODO: test --- diff --git a/pkg/R/F_Neighbors.R b/pkg/R/F_Neighbors.R index 27cd23a..c55291a 100644 --- a/pkg/R/F_Neighbors.R +++ b/pkg/R/F_Neighbors.R @@ -22,32 +22,33 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", fdays = getNoNA2(data, max(today-memory,1), today-1) # Get optional args + local = ifelse(hasArg("local"), list(...)$local, FALSE) #same level + season? simtype = ifelse(hasArg("simtype"), list(...)$simtype, "mix") #or "endo", or "exo" - kernel = ifelse(hasArg("kernel"), list(...)$kernel, "Gauss") #or "Epan" - if (hasArg(h_window)) + if (hasArg("window")) { return ( private$.predictShapeAux(data, - fdays, today, horizon, list(...)$h_window, kernel, simtype, TRUE) ) + fdays, today, horizon, local, list(...)$window, simtype, TRUE) ) } # Indices of similar days for cross-validation; TODO: 20 = magic number - cv_days = getSimilarDaysIndices(today, data, limit=20, same_season=FALSE, days_in=fdays) + cv_days = getSimilarDaysIndices(today, data, limit=20, same_season=FALSE, + days_in=fdays) - # Function to optimize h : h |--> sum of prediction errors on last 45 "similar" days - errorOnLastNdays = function(h, kernel, simtype) + # Optimize h : h |--> sum of prediction errors on last 45 "similar" days + errorOnLastNdays = function(window, simtype) { error = 0 nb_jours = 0 for (i in seq_along(cv_days)) { # mix_strategy is never used here (simtype != "mix"), therefore left blank - prediction = private$.predictShapeAux(data, - fdays, cv_days[i], horizon, h, kernel, simtype, FALSE) + prediction = private$.predictShapeAux(data, fdays, cv_days[i], horizon, local, + window, simtype, FALSE) if (!is.na(prediction[1])) { nb_jours = nb_jours + 1 error = error + - mean((data$getCenteredSerie(cv_days[i]+1)[1:horizon] - prediction)^2) + mean((data$getSerie(cv_days[i]+1)[1:horizon] - prediction)^2) } } return (error / nb_jours) @@ -55,45 +56,87 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", if (simtype != "endo") { - h_best_exo = optimize( - errorOnLastNdays, c(0,7), kernel=kernel, simtype="exo")$minimum + best_window_exo = optimize( + errorOnLastNdays, c(0,7), simtype="exo")$minimum } if (simtype != "exo") { - h_best_endo = optimize( - errorOnLastNdays, c(0,7), kernel=kernel, simtype="endo")$minimum + best_window_endo = optimize( + errorOnLastNdays, c(0,7), simtype="endo")$minimum } if (simtype == "endo") { - return (private$.predictShapeAux(data, - fdays, today, horizon, h_best_endo, kernel, "endo", TRUE)) + return (private$.predictShapeAux(data, fdays, today, horizon, local, + best_window_endo, "endo", TRUE)) } if (simtype == "exo") { - return (private$.predictShapeAux(data, - fdays, today, horizon, h_best_exo, kernel, "exo", TRUE)) + return (private$.predictShapeAux(data, fdays, today, horizon, local, + best_window_exo, "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)) + return(private$.predictShapeAux(data, fdays, today, horizon, local, + c(best_window_endo,best_window_exo), "mix", TRUE)) } } ), private = list( # Precondition: "today" is full (no NAs) - .predictShapeAux = function(data, fdays, today, horizon, h, kernel, simtype, final_call) + .predictShapeAux = function(data, fdays, today, horizon, local, window, simtype, + final_call) { - fdays = fdays[ fdays < today ] - # TODO: 3 = magic number - if (length(fdays) < 3) + fdays_cut = fdays[ fdays < today ] + if (length(fdays_cut) <= 1) return (NA) + if (local) + { + # Neighbors: days in "same season"; TODO: 60 == magic number... + fdays = getSimilarDaysIndices(today, data, limit=60, same_season=TRUE, + days_in=fdays_cut) + if (length(fdays) <= 1) + return (NA) + levelToday = data$getLevel(today) + distances = sapply(fdays, function(i) abs(data$getLevel(i)-levelToday)) + #TODO: 2, 3, 5, 10 magic numbers here... + dist_thresh = 2 + min_neighbs = min(3,length(fdays)) + repeat + { + same_pollution = (distances <= dist_thresh) + nb_neighbs = sum(same_pollution) + if (nb_neighbs >= min_neighbs) #will eventually happen + break + dist_thresh = dist_thresh + 3 + } + fdays = fdays[same_pollution] + max_neighbs = 10 + if (nb_neighbs > max_neighbs) + { + # Keep only max_neighbs closest neighbors + fdays = fdays[ + sort(distances[same_pollution],index.return=TRUE)$ix[1:max_neighbs] ] + } + if (length(fdays) == 1) #the other extreme... + { + if (final_call) + { + private$.params$weights <- 1 + private$.params$indices <- fdays + private$.params$window <- 1 + } + return ( data$getSerie(fdays[1])[1:horizon] ) #what else?! + } + } + else + fdays = fdays_cut #no conditioning + if (simtype != "exo") { - h_endo = ifelse(simtype=="mix", h[1], h) + # Compute endogen similarities using given window + window_endo = ifelse(simtype=="mix", window[1], window) # Distances from last observed day to days in the past serieToday = data$getSerie(today) @@ -103,26 +146,18 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", }) sd_dist = sd(distances2) - if (sd_dist < .Machine$double.eps) + 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_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 - } + simils_endo = exp(-distances2/(sd_dist*window_endo^2)) } if (simtype != "endo") { - h_exo = ifelse(simtype=="mix", h[2], h) + # Compute exogen similarities using given window + h_exo = ifelse(simtype=="mix", window[2], window) M = matrix( nrow=1+length(fdays), ncol=1+length(data$getExo(today)) ) M[1,] = c( data$getLevel(today), as.double(data$getExo(today)) ) @@ -149,16 +184,7 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", # 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 - } + simils_exo = exp(-distances2/(sd_dist*window_exo^2)) } similarities = @@ -172,7 +198,7 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", prediction = rep(0, horizon) for (i in seq_along(fdays)) - prediction = prediction + similarities[i] * data$getCenteredSerie(fdays[i]+1)[1:horizon] + prediction = prediction + similarities[i] * data$getSerie(fdays[i]+1)[1:horizon] if (final_call) { @@ -180,11 +206,11 @@ NeighborsForecaster = R6::R6Class("NeighborsForecaster", private$.params$indices <- fdays private$.params$window <- if (simtype=="endo") - h_endo + window_endo else if (simtype=="exo") - h_exo + window_exo else #mix - c(h_endo,h_exo) + c(window_endo,window_exo) } return (prediction) diff --git a/pkg/R/F_Neighbors2.R b/pkg/R/F_Neighbors2.R deleted file mode 100644 index ee40f61..0000000 --- a/pkg/R/F_Neighbors2.R +++ /dev/null @@ -1,228 +0,0 @@ -#' @include Forecaster.R -#' -#' Neighbors2 Forecaster -#' -#' Predict tomorrow as a weighted combination of "futures of the past" days. -#' Inherits \code{\link{Forecaster}} -#' -Neighbors2Forecaster = R6::R6Class("Neighbors2Forecaster", - inherit = Forecaster, - - public = list( - predictShape = function(data, today, memory, horizon, ...) - { - # (re)initialize computed parameters - private$.params <- list("weights"=NA, "indices"=NA, "window"=NA) - - # Do not forecast on days with NAs (TODO: softer condition...) - if (any(is.na(data$getCenteredSerie(today)))) - return (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 ( private$.predictShapeAux(data, - fdays, today, horizon, list(...)$h_window, kernel, simtype, TRUE) ) - } - - # Indices of similar days for cross-validation; TODO: 20 = magic number - cv_days = getSimilarDaysIndices(today, data, limit=20, same_season=FALSE, days_in=fdays) - - # 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 seq_along(cv_days)) - { - # mix_strategy is never used here (simtype != "mix"), therefore left blank - 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(cv_days[i]+1)[1:horizon] - prediction)^2) - } - } - return (error / nb_jours) - } - - if (simtype != "endo") - { - h_best_exo = optimize( - errorOnLastNdays, c(0,7), kernel=kernel, simtype="exo")$minimum - } - if (simtype != "exo") - { - h_best_endo = optimize( - errorOnLastNdays, c(0,7), kernel=kernel, simtype="endo")$minimum - } - - 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, simtype, final_call) - { - fdays_cut = fdays[ fdays < today ] - # TODO: 3 = magic number - if (length(fdays_cut) < 3) - return (NA) - - # Neighbors: days in "same season"; TODO: 60 == magic number... - fdays = getSimilarDaysIndices(today, data, limit=60, same_season=TRUE, days_in=fdays_cut) - if (length(fdays) <= 1) - return (NA) - levelToday = data$getLevel(today) - distances = sapply(fdays, function(i) abs(data$getLevel(i)-levelToday)) - #TODO: 2, 3, 5, 10 magic numbers here... - dist_thresh = 2 - min_neighbs = min(3,length(fdays)) - repeat - { - same_pollution = (distances <= dist_thresh) - nb_neighbs = sum(same_pollution) - if (nb_neighbs >= min_neighbs) #will eventually happen - break - dist_thresh = dist_thresh + 3 - } - fdays = fdays[same_pollution] - max_neighbs = 10 - if (nb_neighbs > max_neighbs) - { - # Keep only max_neighbs closest neighbors - fdays = fdays[ sort(distances[same_pollution],index.return=TRUE)$ix[1:max_neighbs] ] - } - if (length(fdays) == 1) #the other extreme... - { - if (final_call) - { - private$.params$weights <- 1 - private$.params$indices <- fdays - private$.params$window <- 1 - } - return ( data$getSerie(fdays[1])[1:horizon] ) #what else?! - } - - if (simtype != "exo") - { - h_endo = ifelse(simtype=="mix", h[1], h) - - # Distances from last observed day to days in the past - serieToday = data$getSerie(today) - distances2 = sapply(fdays, function(i) { - delta = serieToday - data$getSerie(i) - mean(delta^2) - }) - - 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 - } - } - - if (simtype != "endo") - { - h_exo = ifelse(simtype=="mix", h[2], h) - - 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])) ) - - sigma = cov(M) #NOTE: robust covariance is way too slow - # TODO: 10 == magic number; more robust way == det, or always ginv() - sigma_inv = - if (length(fdays) > 10) - solve(sigma) - else - MASS::ginv(sigma) - - # Distances from last observed day to days in the past - distances2 = sapply(seq_along(fdays), 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 - } - } - - similarities = - if (simtype == "exo") - simils_exo - else if (simtype == "endo") - simils_endo - else #mix - simils_endo * simils_exo - similarities = similarities / sum(similarities) - - prediction = rep(0, horizon) - for (i in seq_along(fdays)) - prediction = prediction + similarities[i] * data$getSerie(fdays[i]+1)[1:horizon] - - if (final_call) - { - prediction = prediction - mean(prediction) #predict centered serie (artificial...) - private$.params$weights <- similarities - private$.params$indices <- fdays - private$.params$window <- - if (simtype=="endo") - h_endo - else if (simtype=="exo") - h_exo - else #mix - c(h_endo,h_exo) - } - - return (prediction) - } - ) -) diff --git a/reports/report.gj b/reports/report.gj index 1a6b8d9..5e57660 100644 --- a/reports/report.gj +++ b/reports/report.gj @@ -1,11 +1,14 @@ -----

Introduction

-J'ai fait quelques essais dans différentes configurations pour la méthode "Neighbors" -(la seule dont on a parlé) et sa variante récente appelée pour l'instant "Neighbors2", -avec simtype="mix" : deux types de similarités prises en compte, puis multiplication des poids. -Pour Neighbors on prédit le saut (par la moyenne pondérée des sauts passés), et pour Neighbors2 -on n'effectue aucun raccordement (prévision directe). +J'ai fait quelques essais dans deux configurations pour la méthode "Neighbors" +(la seule dont on a parlé, incorporant désormais la "variante Bruno/Michel"). + + * avec simtype="mix" et raccordement simple ("Zero") dans le cas "non local", i.e. on va + chercher des voisins n'importe où du moment qu'ils correspondent à deux jours consécutifs sans + valeurs manquantes. + * avec simtype="endo" et raccordement "Neighbor" dans le cas "local" : voisins de même niveau de + pollution et même saison. J'ai systématiquement comparé à une approche naïve : la moyenne des lendemains des jours "similaires" dans tout le passé, ainsi qu'à la persistence -- reproduisant le jour courant ou @@ -40,74 +43,77 @@ indices_np = seq(as.Date("2015-04-26"),as.Date("2015-05-02"),"days") -----

${list_titles[i]}

-----r -p_nn = computeForecast(data, ${list_indices[i]}, "Neighbors", "Neighbors", horizon=H) -p_nn2 = computeForecast(data, ${list_indices[i]}, "Neighbors2", "Zero", horizon=H) -p_az = computeForecast(data, ${list_indices[i]}, "Average", "Zero", horizon=H) -p_pz = computeForecast(data, ${list_indices[i]}, "Persistence", "Zero", horizon=H, same_day=${'TRUE' if loop.index < 2 else 'FALSE'}) +p_n = computeForecast(data, ${list_indices[i]}, "Neighbors", "Zero", horizon=H, + simtype="mix", local=FALSE) +p_l = computeForecast(data, ${list_indices[i]}, "Neighbors", "Neighbors", horizon=H, + simtype="endo", local=TRUE) +p_a = computeForecast(data, ${list_indices[i]}, "Average", "Zero", horizon=H) +p_p = computeForecast(data, ${list_indices[i]}, "Persistence", "Zero", horizon=H, + same_day=${'TRUE' if loop.index < 2 else 'FALSE'}) -----r -e_nn = computeError(data, p_nn, H) -e_nn2 = computeError(data, p_nn2, H) -e_az = computeError(data, p_az, H) -e_pz = computeError(data, p_pz, H) +e_n = computeError(data, p_n, H) +e_l = computeError(data, p_nl, H) +e_a = computeError(data, p_a, H) +e_p = computeError(data, p_p, H) options(repr.plot.width=9, repr.plot.height=7) -plotError(list(e_nn, e_pz, e_az, e_nn2), cols=c(1,2,colors()[258], 4)) +plotError(list(e_n, e_p, e_a, e_l), cols=c(1,2,colors()[258], 4)) -# Noir: Neighbors, bleu: Neighbors2, vert: moyenne, rouge: persistence +# Noir: Neighbors non-local, bleu: Neighbors local, vert: moyenne, rouge: persistence -i_np = which.min(e_nn$abs$indices) -i_p = which.max(e_nn$abs$indices) +i_np = which.min(e_n$abs$indices) +i_p = which.max(e_n$abs$indices) -----r options(repr.plot.width=9, repr.plot.height=4) par(mfrow=c(1,2)) -plotPredReal(data, p_nn, i_np); title(paste("PredReal nn day",i_np)) -plotPredReal(data, p_nn2, i_p); title(paste("PredReal nn day",i_p)) +plotPredReal(data, p_n, i_np); title(paste("PredReal non-loc day",i_np)) +plotPredReal(data, p_n, i_p); title(paste("PredReal non-loc day",i_p)) -plotPredReal(data, p_nn2, i_np); title(paste("PredReal nn2 day",i_np)) -plotPredReal(data, p_nn2, i_p); title(paste("PredReal nn2 day",i_p)) +plotPredReal(data, p_l, i_np); title(paste("PredReal loc day",i_np)) +plotPredReal(data, p_l, i_p); title(paste("PredReal loc day",i_p)) -plotPredReal(data, p_az, i_np); title(paste("PredReal az day",i_np)) -plotPredReal(data, p_az, i_p); title(paste("PredReal az day",i_p)) +plotPredReal(data, p_a, i_np); title(paste("PredReal avg day",i_np)) +plotPredReal(data, p_a, i_p); title(paste("PredReal avg day",i_p)) # Bleu: prévue, noir: réalisée -----r par(mfrow=c(1,2)) -f_np = computeFilaments(data, p_nn, i_np, plot=TRUE); title(paste("Filaments nn day",i_np)) -f_p = computeFilaments(data, p_nn, i_p, plot=TRUE); title(paste("Filaments nn day",i_p)) +f_np_n = computeFilaments(data, p_n, i_np, plot=TRUE); title(paste("Filaments non-loc day",i_np)) +f_p_n = computeFilaments(data, p_n, i_p, plot=TRUE); title(paste("Filaments non-loc day",i_p)) -f_np2 = computeFilaments(data, p_nn2, i_np, plot=TRUE); title(paste("Filaments nn2 day",i_np)) -f_p2 = computeFilaments(data, p_nn2, i_p, plot=TRUE); title(paste("Filaments nn2 day",i_p)) +f_np_l = computeFilaments(data, p_l, i_np, plot=TRUE); title(paste("Filaments loc day",i_np)) +f_p_l = computeFilaments(data, p_l, i_p, plot=TRUE); title(paste("Filaments loc day",i_p)) -----r par(mfrow=c(1,2)) -plotFilamentsBox(data, f_np); title(paste("FilBox nn day",i_np)) -plotFilamentsBox(data, f_p); title(paste("FilBox nn day",i_p)) +plotFilamentsBox(data, f_np_n); title(paste("FilBox non-loc day",i_np)) +plotFilamentsBox(data, f_p_n); title(paste("FilBox non-loc day",i_p)) # Generally too few neighbors: -#plotFilamentsBox(data, f_np2); title(paste("FilBox nn2 day",i_np)) -#plotFilamentsBox(data, f_p2); title(paste("FilBox nn2 day",i_p)) +#plotFilamentsBox(data, f_np_l); title(paste("FilBox loc day",i_np)) +#plotFilamentsBox(data, f_p_l); title(paste("FilBox loc day",i_p)) -----r par(mfrow=c(1,2)) -plotRelVar(data, f_np); title(paste("StdDev nn day",i_np)) -plotRelVar(data, f_p); title(paste("StdDev nn day",i_p)) +plotRelVar(data, f_np_n); title(paste("StdDev non-loc day",i_np)) +plotRelVar(data, f_p_n); title(paste("StdDev non-loc day",i_p)) -plotRelVar(data, f_np2); title(paste("StdDev nn2 day",i_np)) -plotRelVar(data, f_p2); title(paste("StdDev nn2 day",i_p)) +plotRelVar(data, f_np_l); title(paste("StdDev loc day",i_np)) +plotRelVar(data, f_p_l); title(paste("StdDev loc day",i_p)) # Variabilité globale en rouge ; sur les 60 voisins (+ lendemains) en noir -----r par(mfrow=c(1,2)) -plotSimils(p_nn, i_np); title(paste("Weights nn day",i_np)) -plotSimils(p_nn, i_p); title(paste("Weights nn day",i_p)) +plotSimils(p_n, i_np); title(paste("Weights non-loc day",i_np)) +plotSimils(p_n, i_p); title(paste("Weights non-loc day",i_p)) -plotSimils(p_nn2, i_np); title(paste("Weights nn2 day",i_np)) -plotSimils(p_nn2, i_p); title(paste("Weights nn2 day",i_p)) +plotSimils(p_l, i_np); title(paste("Weights loc day",i_np)) +plotSimils(p_l, i_p); title(paste("Weights loc day",i_p)) # - pollué à gauche, + pollué à droite -----r -# Fenêtres sélectionnées dans ]0,7] / nn à gauche, nn2 à droite -p_nn$getParams(i_np)$window -p_nn$getParams(i_p)$window +# Fenêtres sélectionnées dans ]0,7] / non-loc à gauche, loc à droite +p_n$getParams(i_np)$window +p_n$getParams(i_p)$window -p_nn2$getParams(i_np)$window -p_nn2$getParams(i_p)$window +p_l$getParams(i_np)$window +p_l$getParams(i_p)$window % endfor