From 09cf9c19b5c6a04bc23c58b7ac8a4bbae2c6827d Mon Sep 17 00:00:00 2001 From: Benjamin Auder Date: Mon, 13 Feb 2017 18:15:20 +0100 Subject: [PATCH] fixes and improvements --- R/D_Neighbors.R | 12 +- R/D_Persistence.R | 17 +- R/Data.R | 6 + R/Forecast.R | 2 +- R/S_Average.R | 24 +- R/S_Neighbors.R | 30 +- R/S_Persistence.R | 23 +- R/getData.R | 39 +-- R/getForecast.R | 7 +- R/plot.R | 19 +- R/utils.R | 36 ++- ...t_13-01-2017.rnw => report_2017-01-13.rnw} | 0 ...t_02-02-2017.Rnw => report_2017-02-02.Rnw} | 0 ...-02-2017.ipynb => report_2017-02-02.ipynb} | 13 +- reports/report_2017-03-01.ipynb | 281 ++++++++++++++++++ 15 files changed, 432 insertions(+), 77 deletions(-) rename reports/{report_13-01-2017.rnw => report_2017-01-13.rnw} (100%) rename reports/{report_02-02-2017.Rnw => report_2017-02-02.Rnw} (100%) rename reports/{report_02-02-2017.ipynb => report_2017-02-02.ipynb} (98%) create mode 100644 reports/report_2017-03-01.ipynb diff --git a/R/D_Neighbors.R b/R/D_Neighbors.R index f7bb5a8..dba5f37 100644 --- a/R/D_Neighbors.R +++ b/R/D_Neighbors.R @@ -4,14 +4,18 @@ #' @inheritParams getZeroDeltaForecast getNeighborsDeltaForecast = function(data, today, memory, horizon, shape_params, ...) { - if (any(is.na(shape_params$weights) | is.na(shape_params$indices))) + first_day = max(1, today-memory) + filter = shape_params$indices >= first_day + indices = shape_params$indices[filter] + weights = shape_params$weights[filter] + if (any(is.na(weights) | is.na(indices))) return (NA) - gaps = sapply(shape_params$indices, function(i) { + gaps = sapply(indices, function(i) { data$getSerie(i+1)[1] - tail(data$getSerie(i), 1) }) - scal_product = shape_params$weights * gaps - norm_fact = sum( shape_params$weights[!is.na(scal_product)] ) + scal_product = weights * gaps + norm_fact = sum( weights[!is.na(scal_product)] ) sum(scal_product, na.rm=TRUE) / norm_fact } diff --git a/R/D_Persistence.R b/R/D_Persistence.R index 2f7935a..979bf05 100644 --- a/R/D_Persistence.R +++ b/R/D_Persistence.R @@ -4,6 +4,19 @@ #' @inheritParams getZeroDeltaForecast getPersistenceDeltaForecast = function(data, today, memory, horizon, shape_params, ...) { - # Return level of last similar "tomorrow" (thus -6) computed on "horizon" - data$getSerie(today-6)[1] - tail(data$getSerie(today-7), 1) + #return gap between end of similar day curve and first day of tomorrow (in the past) + first_day = max(1, today-memory) + index = today-7 + repeat + { + { + last_similar_serie_end = tail( data$getCenteredSerie(index), 1) + last_similar_tomorrow_begin = data$getCenteredSerie(index+1)[1] + index = index - 7 + }; + if (!is.na(last_similar_serie_end) && !is.na(last_similar_tomorrow_begin)) + return (last_similar_tomorrow_begin - last_similar_serie_end); + if (index < first_day) + return (NA) + } } diff --git a/R/Data.R b/R/Data.R index 311453b..a17e262 100644 --- a/R/Data.R +++ b/R/Data.R @@ -42,36 +42,42 @@ Data = setRefClass( { "Get time values at specified index" + index = dateIndexToInteger(index, .self) data[[index]]$time }, getCenteredSerie = function(index) { "Get serie values (centered) at specified index" + index = dateIndexToInteger(index, .self) data[[index]]$serie }, getLevel = function(index) { "Get level for the serie at specified index" + index = dateIndexToInteger(index, .self) data[[index]]$level }, getSerie = function(index) { "Get serie values (centered+level) at specified index" + index = dateIndexToInteger(index, .self) data[[index]]$serie + data[[index]]$level }, getExoHat = function(index) { "Get exogeous predictions at specified index" + index = dateIndexToInteger(index, .self) data[[index]]$exo_hat }, getExoDm1 = function(index) { "Get exogenous measures the day before specified index" + index = dateIndexToInteger(index, .self) data[[index]]$exo_Dm1 } ) diff --git a/R/Forecast.R b/R/Forecast.R index 1a3505b..57817e7 100644 --- a/R/Forecast.R +++ b/R/Forecast.R @@ -49,7 +49,7 @@ Forecast = setRefClass( }, getIndexInData = function(index) { - "Get (day)index at specified index" + "Get (day) index in data where prediction took place" pred[[index]]$index } diff --git a/R/S_Average.R b/R/S_Average.R index 9694693..819531d 100644 --- a/R/S_Average.R +++ b/R/S_Average.R @@ -2,8 +2,8 @@ #' #' @title Average Shape Forecaster #' -#' @description Return the (pointwise) average of the all the centered day curves in the past. -#' Inherits \code{\link{ShapeForecaster}} +#' @description Return the (pointwise) average of the all the (similar) centered day curves +#' in the past. Inherits \code{\link{ShapeForecaster}} AverageShapeForecaster = setRefClass( Class = "AverageShapeForecaster", contains = "ShapeForecaster", @@ -16,12 +16,24 @@ AverageShapeForecaster = setRefClass( predict = function(today, memory, horizon, ...) { avg = rep(0., horizon) - for (i in (today-memory):today) + first_day = max(1, today-memory) + index = today-7 + 1 + nb_no_na_series = 0 + repeat { - if (!any(is.na(data$getSerie(i)))) - avg = avg + data$getCenteredSerie(i) + { + serie_on_horizon = data$getCenteredSerie(index)[1:horizon] + index = index - 7 + }; + if (!any(is.na(serie_on_horizon))) + { + avg = avg + serie_on_horizon + nb_no_na_series = nb_no_na_series + 1 + }; + if (index < first_day) + break } - avg + avg / nb_no_na_series } ) ) diff --git a/R/S_Neighbors.R b/R/S_Neighbors.R index c44bd9b..b8e32cc 100644 --- a/R/S_Neighbors.R +++ b/R/S_Neighbors.R @@ -23,9 +23,33 @@ NeighborsShapeForecaster = setRefClass( if (length(data$getCenteredSerie(1)) < length(data$getCenteredSerie(2))) first_day = 2 - # Predict only on non-NAs days (TODO:) - if (any(is.na(data$getSerie(today)))) - return (NA) + # 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() diff --git a/R/S_Persistence.R b/R/S_Persistence.R index 017402f..08647e3 100644 --- a/R/S_Persistence.R +++ b/R/S_Persistence.R @@ -2,8 +2,8 @@ #' #' @title Persistence Shape Forecaster #' -#' @description Return the last centered last (similar) day curve (may be a lot of NAs, -#' but we cannot do better). Inherits \code{\link{ShapeForecaster}} +#' @description Return the last centered last (similar) day curve. +#' Inherits \code{\link{ShapeForecaster}} PersistenceShapeForecaster = setRefClass( Class = "PersistenceShapeForecaster", contains = "ShapeForecaster", @@ -15,11 +15,20 @@ PersistenceShapeForecaster = setRefClass( }, predict = function(today, memory, horizon, ...) { - #return centered last (similar) day curve (may be a lot of NAs, but we cannot do better) - last_similar_serie = data$getCenteredSerie(today-6)[1:horizon] - if (any(is.na(last_similar_serie))) #TODO: - return (NA) - last_similar_serie + #return centered last (similar) day curve, avoiding NAs until memory is run + first_day = max(1, today-memory) + index = today-7 + 1 + repeat + { + { + last_similar_serie = data$getCenteredSerie(index)[1:horizon] + index = index - 7 + }; + if (!any(is.na(last_similar_serie))) + return (last_similar_serie); + if (index < first_day) + return (NA) + } } ) ) diff --git a/R/getData.R b/R/getData.R index c2e6842..5139a4d 100644 --- a/R/getData.R +++ b/R/getData.R @@ -13,14 +13,13 @@ #' @param date_format How date/time are stored (e.g. year/month/day hour:minutes; #' see \code{strptime}) #' @param working_tz Timezone to work with ("GMT" or e.g. "Europe/Paris") -#' @param predict_at When does the prediction take place ? ab[:cd][:ef] where a,b,c,d,e,f -#' in (0,9) and define an hour[minute[second]]; time must be present in the file +#' @param predict_at When does the prediction take place ? Integer, in hours. Default: 0 #' #' @return An object of class Data #' #' @export getData = function(ts_data, exo_data, - input_tz="GMT", date_format="%d/%m/%Y %H:%M", working_tz="GMT", predict_at="00") + input_tz="GMT", date_format="%d/%m/%Y %H:%M", working_tz="GMT", predict_at=0) { # Sanity checks (not full, but sufficient at this stage) if (!is.character(input_tz) || !is.character(working_tz)) @@ -31,10 +30,9 @@ getData = function(ts_data, exo_data, stop("Bad time-series input (data frame or CSV file)") if (is.character(ts_data)) ts_data = ts_data[1] - pattern_index_in_predict_at = grep("^[0-9]{2}(:[0-9]{2}){0,2}$", predict_at) - if (!is.character(predict_at) || length(pattern_index_in_predict_at) == 0) - stop("Bad predict_at ( ^[0-9]{2}(:[0-9]{2}){0,2}$ )") - predict_at = predict_at[ pattern_index_in_predict_at[1] ] + predict_at = as.integer(predict_at)[1] + if (predict_at<0 || predict_at>23) + stop("Bad predict_at (0-23)") if (!is.character(date_format)) stop("Bad date_format (character)") date_format = date_format[1] @@ -55,11 +53,6 @@ getData = function(ts_data, exo_data, formatted_dates_POSIXlt = strptime(as.character(ts_df[,1]), date_format, tz=input_tz) ts_df[,1] = format(as.POSIXct(formatted_dates_POSIXlt), tz=working_tz, usetz=TRUE) - if (nchar(predict_at) == 2) - predict_at = paste(predict_at,":00",sep="") - if (nchar(predict_at) == 5) - predict_at = paste(predict_at,":00",sep="") - line = 1 #index in PM10 file (24 lines for 1 cell) nb_lines = nrow(ts_df) nb_exos = ( ncol(exo_df) - 1 ) / 2 @@ -69,22 +62,20 @@ getData = function(ts_data, exo_data, { time = c() serie = c() - repeat { - { - time = c(time, ts_df[line,1]) - serie = c(serie, ts_df[line,2]) - line = line + 1 - }; - if (line >= nb_lines + 1 - # NOTE: always second part of date/time, because it has been formatted - || strsplit(as.character(ts_df[line-1,1])," ")[[1]][2] == predict_at) + repeat { - break - }} + { + time = c(time, ts_df[line,1]) + serie = c(serie, ts_df[line,2]) + line = line + 1 + }; + if (line >= nb_lines + 1 || as.POSIXlt(ts_df[line-1,1])$hour == predict_at) + break + } # NOTE: if predict_at does not cut days at midnight, # for the exogenous to be synchronized they need to be shifted - if (predict_at != "00:00:00") + if (predict_at > 0) { exo_hat = as.data.frame(exo_df[max(1,i-1),(1+nb_exos+1):(1+2*nb_exos)]) exo_Dm1 = if (i>=3) as.data.frame(exo_df[i-1,2:(1+nb_exos)]) else NA diff --git a/R/getForecast.R b/R/getForecast.R index c3248a9..e126946 100644 --- a/R/getForecast.R +++ b/R/getForecast.R @@ -42,7 +42,7 @@ getForecast = function(data, indices, memory, horizon, horizon = as.integer(horizon)[1] if (horizon<=0 || horizon>length(data$getCenteredSerie(2))) stop("Horizon too short or too long") - indices = as.integer(indices) + indices = sapply( seq_along(indices), function(i) dateIndexToInteger(indices[i], data) ) if (any(indices<=0 | indices>data$getSize())) stop("Indices out of range") indices = sapply(indices, dateIndexToInteger, data) @@ -56,9 +56,6 @@ getForecast = function(data, indices, memory, horizon, pred = list() for (today in indices) { - #NOTE: To estimate timing... -# print(paste("Predict until index",today)) - #shape always predicted first (on centered series, no scaling taken into account), #with side-effect: optimize some parameters (h, weights, ...) predicted_shape = shape_forecaster$predict(today, memory, horizon, ...) @@ -69,7 +66,7 @@ getForecast = function(data, indices, memory, horizon, #TODO: this way is faster than a call to append(); why ? pred[[length(pred)+1]] = list( # Predict shape and align it on end of current day - serie = predicted_shape + tail( data$getSerie(today), 1 ) - predicted_shape[1], + serie = predicted_shape + tail( data$getSerie(today), 1 ) - predicted_shape[1] + predicted_delta, #add predicted jump params = shape_forecaster$getParameters(), index = today ) diff --git a/R/plot.R b/R/plot.R index e5d4753..a7c9395 100644 --- a/R/plot.R +++ b/R/plot.R @@ -13,9 +13,9 @@ plotPredReal <- function(data, pred, index) par(mar=c(4.7,5,1,1), cex.axis=2, cex.lab=2, lwd=2) measure = data$getSerie(pred$getIndexInData(index)+1)[1:horizon] yrange = range( pred$getSerie(index), measure ) - plot(measure, type="l", ylim=yrange, lwd=3) + plot(measure, type="l", ylim=yrange, lwd=3, xlab="Temps (en heures)", ylab="PM10") par(new=TRUE) - plot(pred$getSerie(index), type="l", col=2, ylim=yrange, lwd=3) + plot(pred$getSerie(index), type="l", col="#0000FF", ylim=yrange, lwd=3, xlab="", ylab="") } #' @title Plot filaments @@ -85,20 +85,23 @@ plotSimils <- function(pred, index) #' @description Draw error graphs, potentially from several runs of \code{getForecast} #' #' @param err Error as returned by \code{getError} +#' @param cols Colors for each error (default: 1,2,3,...) #' #' @seealso \code{\link{plotPredReal}}, \code{\link{plotFilaments}}, \code{\link{plotSimils}} #' \code{\link{plotFbox}} #' #' @export -plotError <- function(err) +plotError <- function(err, cols=seq_along(err)) { + if (!is.null(err$abs)) + err = list(err) par(mfrow=c(2,2), mar=c(4.7,5,1,1), cex.axis=2, cex.lab=2, lwd=2) L = length(err) yrange = range( sapply(1:L, function(index) ( err[[index]]$abs$day ) ), na.rm=TRUE ) for (i in seq_len(L)) { plot(err[[i]]$abs$day, type="l", xlab=ifelse(i==1,"Temps (heures)",""), - ylab=ifelse(i==1,"Moyenne |y - y_hat|",""), ylim=yrange, col=i) + ylab=ifelse(i==1,"Moyenne |y - y_hat|",""), ylim=yrange, col=cols[i]) if (i < L) par(new=TRUE) } @@ -106,7 +109,7 @@ plotError <- function(err) for (i in seq_len(L)) { plot(err[[i]]$abs$indices, type="l", xlab=ifelse(i==1,"Temps (jours)",""), - ylab=ifelse(i==1,"Moyenne |y - y_hat|",""), ylim=yrange, col=i) + ylab=ifelse(i==1,"Moyenne |y - y_hat|",""), ylim=yrange, col=cols[i]) if (i < L) par(new=TRUE) } @@ -114,7 +117,7 @@ plotError <- function(err) for (i in seq_len(L)) { plot(err[[i]]$MAPE$day, type="l", xlab=ifelse(i==1,"Temps (heures)",""), - ylab=ifelse(i==1,"MAPE moyen",""), ylim=yrange, col=i) + ylab=ifelse(i==1,"MAPE moyen",""), ylim=yrange, col=cols[i]) if (i < L) par(new=TRUE) } @@ -122,7 +125,7 @@ plotError <- function(err) for (i in seq_len(L)) { plot(err[[i]]$MAPE$indices, type="l", xlab=ifelse(i==1,"Temps (jours)",""), - ylab=ifelse(i==1,"MAPE moyen",""), ylim=yrange, col=i) + ylab=ifelse(i==1,"MAPE moyen",""), ylim=yrange, col=cols[i]) if (i < L) par(new=TRUE) } @@ -136,7 +139,7 @@ plotError <- function(err) #' @param fiter Optional filter: return TRUE on indices to process #' #' @export -plotFbox <- function(data, filter=function(index) (TRUE)) +plotFbox <- function(data, filter=function(index) TRUE) { if (!requireNamespace("rainbow", quietly=TRUE)) stop("Functional boxplot requires the rainbow package") diff --git a/R/utils.R b/R/utils.R index 879a7d7..3777f2d 100644 --- a/R/utils.R +++ b/R/utils.R @@ -8,23 +8,47 @@ #' @export dateIndexToInteger = function(index, data) { + index = index[1] if (is.numeric(index)) index = as.integer(index) - #Undocumented "private" method to translate index --> date (is needed) if (is.integer(index)) - return (index[1]) - if (is.character(index)) + return (index) + if (inherits(index, "Date") || is.character(index)) { - tryCatch(dt <- as.POSIXct(index[1]), finally=print("Unrecognized index format")) + tryCatch(dt <- as.POSIXct(index), error=function(e) stop("Unrecognized index format")) #TODO: tz arg to difftime ? - integerIndex <- as.integer( difftime(dt, data$getTime(1)) ) + 1 + integerIndex <- round( (as.numeric( difftime(dt, data$getTime(1)) ))[1] ) + 1 if (integerIndex > 0 && integerIndex <= data$getSize()) - return (integerIndex) + { + #WARNING: if series start at date >0h, result must be shifted + date1 = as.POSIXlt(data$getTime(1)[1]) + date2 = as.POSIXlt(data$getTime(2)[1]) + shift = (date1$year==date2$year && date1$mon==date2$mon && date1$mday==date2$mday) + return (integerIndex + ifelse(shift,1,0)) + } stop("Date outside data range") } stop("Unrecognized index format") } +#' @title integerIndexToDate +#' +#' @description Transform an integer index to date index (relative to data) +#' +#' @param index Date (or integer) index +#' @param data Object of class \code{Data} +#' +#' @export +integerIndexToDate = function(index, data) +{ + index = index[1] + if (is.numeric(index)) + index = as.integer(index) + if (!is.integer(index)) + stop("'index' should be an integer") + as.Date( data$getTime(index)[1] ) +} + #' @title getSimilarDaysIndices #' #' @description Find similar days indices in the past diff --git a/reports/report_13-01-2017.rnw b/reports/report_2017-01-13.rnw similarity index 100% rename from reports/report_13-01-2017.rnw rename to reports/report_2017-01-13.rnw diff --git a/reports/report_02-02-2017.Rnw b/reports/report_2017-02-02.Rnw similarity index 100% rename from reports/report_02-02-2017.Rnw rename to reports/report_2017-02-02.Rnw diff --git a/reports/report_02-02-2017.ipynb b/reports/report_2017-02-02.ipynb similarity index 98% rename from reports/report_02-02-2017.ipynb rename to reports/report_2017-02-02.ipynb index 4f49077..00380f7 100644 --- a/reports/report_02-02-2017.ipynb +++ b/reports/report_2017-02-02.ipynb @@ -57,20 +57,11 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": { "collapsed": false }, - "outputs": [ - { - "ename": "ERROR", - "evalue": "Error in eval(expr, envir, enclos): could not find function \"getData\"\n", - "output_type": "error", - "traceback": [ - "Error in eval(expr, envir, enclos): could not find function \"getData\"\nTraceback:\n" - ] - } - ], + "outputs": [], "source": [ "# Voir ?getData pour les arguments\n", "data = getData(ts_data=\"data/pm10_mesures_H_loc.csv\", exo_data=\"data/meteo_extra_noNAs.csv\",\n", diff --git a/reports/report_2017-03-01.ipynb b/reports/report_2017-03-01.ipynb new file mode 100644 index 0000000..4eecb50 --- /dev/null +++ b/reports/report_2017-03-01.ipynb @@ -0,0 +1,281 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "library(talweg)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "data = getData(ts_data=\"../data/pm10_mesures_H_loc.csv\", exo_data=\"../data/meteo_extra_noNAs.csv\",\n", + " input_tz = \"Europe/Paris\", working_tz=\"Europe/Paris\", predict_at=7)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Pollution par chauffage" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "p_ch_nn = getForecast(data, seq(as.Date(\"2015-01-18\"),as.Date(\"2015-01-24\"),\"days\"), Inf, 17,\n", + " \"Neighbors\", \"Neighbors\", simtype=\"mix\", same_season=FALSE, mix_strategy=\"mult\")\n", + "p_ch_pz = getForecast(data, seq(as.Date(\"2015-01-18\"),as.Date(\"2015-01-24\"),\"days\"), Inf, 17,\n", + " \"Persistence\", \"Zero\")\n", + "p_ch_az = getForecast(data, seq(as.Date(\"2015-01-18\"),as.Date(\"2015-01-24\"),\"days\"), Inf, 17,\n", + " \"Average\", \"Zero\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "e_ch_nn = getError(data, p_ch_nn, 17)\n", + "e_ch_pz = getError(data, p_ch_pz, 17)\n", + "e_ch_az = getError(data, p_ch_az, 17)\n", + "options(repr.plot.width=9, repr.plot.height=6)\n", + "plotError(list(e_ch_nn, e_ch_pz, e_ch_az), cols=c(1,2,colors()[258]))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "par(mfrow=c(1,2))\n", + "options(repr.plot.width=9, repr.plot.height=4)\n", + "plotPredReal(data, p_ch_nn, 3)\n", + "plotPredReal(data, p_ch_nn, 4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "par(mfrow=c(1,2))\n", + "plotFilaments(data, p_ch_nn$getIndexInData(3))\n", + "plotFilaments(data, p_ch_nn$getIndexInData(4))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "par(mfrow=c(1,3))\n", + "plotSimils(p_ch_nn, 3)\n", + "plotSimils(p_ch_nn, 4)\n", + "plotSimils(p_ch_nn, 5)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Pollution par épandage" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "p_ep_nn = getForecast(data, seq(as.Date(\"2015-03-15\"),as.Date(\"2015-03-21\"),\"days\"), Inf, 17,\n", + " \"Neighbors\", \"Neighbors\", simtype=\"mix\", same_season=FALSE, mix_strategy=\"mult\")\n", + "p_ep_pz = getForecast(data, seq(as.Date(\"2015-03-15\"),as.Date(\"2015-03-21\"),\"days\"), Inf, 17,\n", + " \"Persistence\", \"Zero\")\n", + "p_ep_az = getForecast(data, seq(as.Date(\"2015-03-15\"),as.Date(\"2015-03-21\"),\"days\"), Inf, 17,\n", + " \"Average\", \"Zero\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "e_ep_nn = getError(data, p_ep_nn, 17)\n", + "e_ep_pz = getError(data, p_ep_pz, 17)\n", + "e_ep_az = getError(data, p_ep_az, 17)\n", + "options(repr.plot.width=9, repr.plot.height=6)\n", + "plotError(list(e_ep_nn, e_ep_pz, e_ep_az), cols=c(1,2,colors()[258]))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "par(mfrow=c(1,2))\n", + "options(repr.plot.width=9, repr.plot.height=4)\n", + "plotPredReal(data, p_ep_nn, 3)\n", + "plotPredReal(data, p_ep_nn, 4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "par(mfrow=c(1,2))\n", + "plotFilaments(data, p_ep_nn$getIndexInData(3))\n", + "plotFilaments(data, p_ep_nn$getIndexInData(4))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "par(mfrow=c(1,3))\n", + "plotSimils(p_ep_nn, 3)\n", + "plotSimils(p_ep_nn, 4)\n", + "plotSimils(p_ep_nn, 5)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Semaine non polluée" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "p_np_nn = getForecast(data, seq(as.Date(\"2015-04-26\"),as.Date(\"2015-05-02\"),\"days\"), Inf, 17,\n", + " \"Neighbors\", \"Neighbors\", simtype=\"mix\", same_season=FALSE, mix_strategy=\"mult\")\n", + "p_np_pz = getForecast(data, seq(as.Date(\"2015-04-26\"),as.Date(\"2015-05-02\"),\"days\"), Inf, 17,\n", + " \"Persistence\", \"Zero\")\n", + "p_np_az = getForecast(data, seq(as.Date(\"2015-04-26\"),as.Date(\"2015-05-02\"),\"days\"), Inf, 17,\n", + " \"Average\", \"Zero\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "e_np_nn = getError(data, p_np_nn, 17)\n", + "e_np_pz = getError(data, p_np_pz, 17)\n", + "e_np_az = getError(data, p_np_az, 17)\n", + "options(repr.plot.width=9, repr.plot.height=6)\n", + "plotError(list(e_np_nn, e_np_pz, e_np_az), cols=c(1,2,colors()[258]))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "par(mfrow=c(1,2))\n", + "options(repr.plot.width=9, repr.plot.height=4)\n", + "plotPredReal(data, p_np_nn, 3)\n", + "plotPredReal(data, p_np_nn, 4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "par(mfrow=c(1,2))\n", + "plotFilaments(data, p_np_nn$getIndexInData(3))\n", + "plotFilaments(data, p_np_nn$getIndexInData(4))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "par(mfrow=c(1,3))\n", + "plotSimils(p_np_nn, 3)\n", + "plotSimils(p_np_nn, 4)\n", + "plotSimils(p_np_nn, 5)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "R", + "language": "R", + "name": "ir" + }, + "language_info": { + "codemirror_mode": "r", + "file_extension": ".r", + "mimetype": "text/x-r-source", + "name": "R", + "pygments_lexer": "r", + "version": "3.3.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} -- 2.44.0