fixes and improvements
authorBenjamin Auder <benjamin.auder@somewhere>
Mon, 13 Feb 2017 17:15:20 +0000 (18:15 +0100)
committerBenjamin Auder <benjamin.auder@somewhere>
Mon, 13 Feb 2017 17:15:20 +0000 (18:15 +0100)
15 files changed:
R/D_Neighbors.R
R/D_Persistence.R
R/Data.R
R/Forecast.R
R/S_Average.R
R/S_Neighbors.R
R/S_Persistence.R
R/getData.R
R/getForecast.R
R/plot.R
R/utils.R
reports/report_2017-01-13.rnw [moved from reports/report_13-01-2017.rnw with 100% similarity]
reports/report_2017-02-02.Rnw [moved from reports/report_02-02-2017.Rnw with 100% similarity]
reports/report_2017-02-02.ipynb [moved from reports/report_02-02-2017.ipynb with 98% similarity]
reports/report_2017-03-01.ipynb [new file with mode: 0644]

index f7bb5a8..dba5f37 100644 (file)
@@ -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
 }
index 2f7935a..979bf05 100644 (file)
@@ -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)
+       }
 }
index 311453b..a17e262 100644 (file)
--- 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
                }
        )
index 1a3505b..57817e7 100644 (file)
@@ -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
                }
index 9694693..819531d 100644 (file)
@@ -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
                }
        )
 )
index c44bd9b..b8e32cc 100644 (file)
@@ -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()
index 017402f..08647e3 100644 (file)
@@ -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)
+                       }
                }
        )
 )
index c2e6842..5139a4d 100644 (file)
 #' @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
index c3248a9..e126946 100644 (file)
@@ -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 )
index e5d4753..a7c9395 100644 (file)
--- 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")
index 879a7d7..3777f2d 100644 (file)
--- 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
similarity index 98%
rename from reports/report_02-02-2017.ipynb
rename to reports/report_2017-02-02.ipynb
index 4f49077..00380f7 100644 (file)
   },
   {
    "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 (file)
index 0000000..4eecb50
--- /dev/null
@@ -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
+}