diff --git a/.travis.yml b/.travis.yml index 6692785..3c6694e 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,11 +1,18 @@ -language: python -python: - - "2.7" - - "3.6" - -install: - - pip install --upgrade pip - - pip install -U -r python/requirements.txt - -script: - - cd python && python setup.py develop test +matrix: + include: + - language: python + python: + - "2.7" + - "3.6" + + install: + - pip install --upgrade pip + - pip install -U -r python/requirements.txt + + script: + - cd python && python setup.py develop test + - language: R + sudo: false + cache: packages + before_install: + - cd R diff --git a/R/DESCRIPTION b/R/DESCRIPTION index d689bb2..44c6b73 100644 --- a/R/DESCRIPTION +++ b/R/DESCRIPTION @@ -16,13 +16,16 @@ Depends: Rcpp (>= 0.12.0) Imports: dplyr (>= 0.5.0), + dygraphs, extraDistr, ggplot2, grid, rstan (>= 2.14.0), scales, stats, - tidyr (>= 0.6.1) + tidyr (>= 0.6.1), + utils, + xts Suggests: knitr, testthat, diff --git a/R/NAMESPACE b/R/NAMESPACE index 415e5a9..1e460d9 100644 --- a/R/NAMESPACE +++ b/R/NAMESPACE @@ -5,7 +5,9 @@ S3method(predict,prophet) export(add_regressor) export(add_seasonality) export(cross_validation) +export(dyplot.prophet) export(fit.prophet) +export(layer_changepoints) export(make_future_dataframe) export(plot_forecast_component) export(predictive_samples) diff --git a/R/R/diagnostics.R b/R/R/diagnostics.R index 4723298..bb7ee01 100644 --- a/R/R/diagnostics.R +++ b/R/R/diagnostics.R @@ -26,7 +26,7 @@ generate_cutoffs <- function(df, horizon, k, period) { stop('Less data than horizon.') } tzone <- attr(cutoff, "tzone") # Timezone is wiped by putting in array - result <- c(cutoff) + result <- cutoff if (k > 1) { for (i in 2:k) { cutoff <- cutoff - period @@ -74,7 +74,7 @@ simulated_historical_forecasts <- function(model, horizon, units, k, } cutoffs <- generate_cutoffs(df, horizon, k, period) predicts <- data.frame() - for (i in 1:length(cutoffs)) { + for (i in seq_along(cutoffs)) { cutoff <- cutoffs[i] # Copy the model m <- prophet_copy(model, cutoff) @@ -84,7 +84,7 @@ simulated_historical_forecasts <- function(model, horizon, units, k, # Calculate yhat df.predict <- dplyr::filter(df, ds > cutoff, ds <= cutoff + horizon) # Get the columns for the future dataframe - columns <- c('ds') + columns <- 'ds' if (m$growth == 'logistic') { columns <- c(columns, 'cap') if (m$logistic.floor) { diff --git a/R/R/plot.R b/R/R/plot.R new file mode 100644 index 0000000..1264b2b --- /dev/null +++ b/R/R/plot.R @@ -0,0 +1,374 @@ +#' Merge history and forecast for plotting. +#' +#' @param m Prophet object. +#' @param fcst Data frame returned by prophet predict. +#' +#' @importFrom dplyr "%>%" +#' @keywords internal +df_for_plotting <- function(m, fcst) { + # Make sure there is no y in fcst + fcst$y <- NULL + df <- m$history %>% + dplyr::select(ds, y) %>% + dplyr::full_join(fcst, by = "ds") %>% + dplyr::arrange(ds) + return(df) +} + +#' Plot the prophet forecast. +#' +#' @param x Prophet object. +#' @param fcst Data frame returned by predict(m, df). +#' @param uncertainty Boolean indicating if the uncertainty interval for yhat +#' should be plotted. Must be present in fcst as yhat_lower and yhat_upper. +#' @param plot_cap Boolean indicating if the capacity should be shown in the +#' figure, if available. +#' @param xlabel Optional label for x-axis +#' @param ylabel Optional label for y-axis +#' @param ... additional arguments +#' +#' @return A ggplot2 plot. +#' +#' @examples +#' \dontrun{ +#' history <- data.frame(ds = seq(as.Date('2015-01-01'), as.Date('2016-01-01'), by = 'd'), +#' y = sin(1:366/200) + rnorm(366)/10) +#' m <- prophet(history) +#' future <- make_future_dataframe(m, periods = 365) +#' forecast <- predict(m, future) +#' plot(m, forecast) +#' } +#' +#' @export +plot.prophet <- function(x, fcst, uncertainty = TRUE, plot_cap = TRUE, + xlabel = 'ds', ylabel = 'y', ...) { + df <- df_for_plotting(x, fcst) + gg <- ggplot2::ggplot(df, ggplot2::aes(x = ds, y = y)) + + ggplot2::labs(x = xlabel, y = ylabel) + if (exists('cap', where = df) && plot_cap) { + gg <- gg + ggplot2::geom_line( + ggplot2::aes(y = cap), linetype = 'dashed', na.rm = TRUE) + } + if (x$logistic.floor && exists('floor', where = df) && plot_cap) { + gg <- gg + ggplot2::geom_line( + ggplot2::aes(y = floor), linetype = 'dashed', na.rm = TRUE) + } + if (uncertainty && exists('yhat_lower', where = df)) { + gg <- gg + + ggplot2::geom_ribbon(ggplot2::aes(ymin = yhat_lower, ymax = yhat_upper), + alpha = 0.2, + fill = "#0072B2", + na.rm = TRUE) + } + gg <- gg + + ggplot2::geom_point(na.rm=TRUE) + + ggplot2::geom_line(ggplot2::aes(y = yhat), color = "#0072B2", + na.rm = TRUE) + + ggplot2::theme(aspect.ratio = 3 / 5) + return(gg) +} + +#' Plot the components of a prophet forecast. +#' Prints a ggplot2 with panels for trend, weekly and yearly seasonalities if +#' present, and holidays if present. +#' +#' @param m Prophet object. +#' @param fcst Data frame returned by predict(m, df). +#' @param uncertainty Boolean indicating if the uncertainty interval should be +#' plotted for the trend, from fcst columns trend_lower and trend_upper. +#' @param plot_cap Boolean indicating if the capacity should be shown in the +#' figure, if available. +#' @param weekly_start Integer specifying the start day of the weekly +#' seasonality plot. 0 (default) starts the week on Sunday. 1 shifts by 1 day +#' to Monday, and so on. +#' @param yearly_start Integer specifying the start day of the yearly +#' seasonality plot. 0 (default) starts the year on Jan 1. 1 shifts by 1 day +#' to Jan 2, and so on. +#' +#' @return Invisibly return a list containing the plotted ggplot objects +#' +#' @export +#' @importFrom dplyr "%>%" +prophet_plot_components <- function( + m, fcst, uncertainty = TRUE, plot_cap = TRUE, weekly_start = 0, + yearly_start = 0 +) { + # Plot the trend + panels <- list(plot_forecast_component(fcst, 'trend', uncertainty, plot_cap)) + # Plot holiday components, if present. + if (!is.null(m$holidays) && ('holidays' %in% colnames(fcst))) { + panels[[length(panels) + 1]] <- plot_forecast_component( + fcst, 'holidays', uncertainty, FALSE) + } + # Plot weekly seasonality, if present + if ("weekly" %in% colnames(fcst)) { + panels[[length(panels) + 1]] <- plot_weekly(m, uncertainty, weekly_start) + } + # Plot yearly seasonality, if present + if ("yearly" %in% colnames(fcst)) { + panels[[length(panels) + 1]] <- plot_yearly(m, uncertainty, yearly_start) + } + # Plot other seasonalities + for (name in names(m$seasonalities)) { + if (!(name %in% c('weekly', 'yearly')) && + (name %in% colnames(fcst))) { + panels[[length(panels) + 1]] <- plot_seasonality(m, name, uncertainty) + } + } + # Plot extra regressors + if ((length(m$extra_regressors) > 0) + & ('extra_regressors' %in% colnames(fcst))) { + panels[[length(panels) + 1]] <- plot_forecast_component( + fcst, 'extra_regressors', uncertainty, FALSE) + } + + # Make the plot. + grid::grid.newpage() + grid::pushViewport(grid::viewport(layout = grid::grid.layout(length(panels), + 1))) + for (i in seq_along(panels)) { + print(panels[[i]], vp = grid::viewport(layout.pos.row = i, + layout.pos.col = 1)) + } + return(invisible(panels)) +} + +#' Plot a particular component of the forecast. +#' +#' @param fcst Dataframe output of `predict`. +#' @param name String name of the component to plot (column of fcst). +#' @param uncertainty Boolean to plot uncertainty intervals. +#' @param plot_cap Boolean indicating if the capacity should be shown in the +#' figure, if available. +#' +#' @return A ggplot2 plot. +#' +#' @export +plot_forecast_component <- function( + fcst, name, uncertainty = TRUE, plot_cap = FALSE +) { + gg.comp <- ggplot2::ggplot( + fcst, ggplot2::aes_string(x = 'ds', y = name, group = 1)) + + ggplot2::geom_line(color = "#0072B2", na.rm = TRUE) + if (exists('cap', where = fcst) && plot_cap) { + gg.comp <- gg.comp + ggplot2::geom_line( + ggplot2::aes(y = cap), linetype = 'dashed', na.rm = TRUE) + } + if (exists('floor', where = fcst) && plot_cap) { + gg.comp <- gg.comp + ggplot2::geom_line( + ggplot2::aes(y = floor), linetype = 'dashed', na.rm = TRUE) + } + if (uncertainty) { + gg.comp <- gg.comp + + ggplot2::geom_ribbon( + ggplot2::aes_string( + ymin = paste0(name, '_lower'), ymax = paste0(name, '_upper') + ), + alpha = 0.2, + fill = "#0072B2", + na.rm = TRUE) + } + return(gg.comp) +} + +#' Prepare dataframe for plotting seasonal components. +#' +#' @param m Prophet object. +#' @param ds Array of dates for column ds. +#' +#' @return A dataframe with seasonal components on ds. +#' +#' @keywords internal +seasonality_plot_df <- function(m, ds) { + df_list <- list(ds = ds, cap = 1, floor = 0) + for (name in names(m$extra_regressors)) { + df_list[[name]] <- 0 + } + df <- as.data.frame(df_list) + df <- setup_dataframe(m, df)$df + return(df) +} + +#' Plot the weekly component of the forecast. +#' +#' @param m Prophet model object +#' @param uncertainty Boolean to plot uncertainty intervals. +#' @param weekly_start Integer specifying the start day of the weekly +#' seasonality plot. 0 (default) starts the week on Sunday. 1 shifts by 1 day +#' to Monday, and so on. +#' +#' @return A ggplot2 plot. +#' +#' @keywords internal +plot_weekly <- function(m, uncertainty = TRUE, weekly_start = 0) { + # Compute weekly seasonality for a Sun-Sat sequence of dates. + days <- seq(set_date('2017-01-01'), by='d', length.out=7) + as.difftime( + weekly_start, units = "days") + df.w <- seasonality_plot_df(m, days) + seas <- predict_seasonal_components(m, df.w) + seas$dow <- factor(weekdays(df.w$ds), levels=weekdays(df.w$ds)) + + gg.weekly <- ggplot2::ggplot(seas, ggplot2::aes(x = dow, y = weekly, + group = 1)) + + ggplot2::geom_line(color = "#0072B2", na.rm = TRUE) + + ggplot2::labs(x = "Day of week") + if (uncertainty) { + gg.weekly <- gg.weekly + + ggplot2::geom_ribbon(ggplot2::aes(ymin = weekly_lower, + ymax = weekly_upper), + alpha = 0.2, + fill = "#0072B2", + na.rm = TRUE) + } + return(gg.weekly) +} + +#' Plot the yearly component of the forecast. +#' +#' @param m Prophet model object. +#' @param uncertainty Boolean to plot uncertainty intervals. +#' @param yearly_start Integer specifying the start day of the yearly +#' seasonality plot. 0 (default) starts the year on Jan 1. 1 shifts by 1 day +#' to Jan 2, and so on. +#' +#' @return A ggplot2 plot. +#' +#' @keywords internal +plot_yearly <- function(m, uncertainty = TRUE, yearly_start = 0) { + # Compute yearly seasonality for a Jan 1 - Dec 31 sequence of dates. + days <- seq(set_date('2017-01-01'), by='d', length.out=365) + as.difftime( + yearly_start, units = "days") + df.y <- seasonality_plot_df(m, days) + seas <- predict_seasonal_components(m, df.y) + seas$ds <- df.y$ds + + gg.yearly <- ggplot2::ggplot(seas, ggplot2::aes(x = ds, y = yearly, + group = 1)) + + ggplot2::geom_line(color = "#0072B2", na.rm = TRUE) + + ggplot2::labs(x = "Day of year") + + ggplot2::scale_x_datetime(labels = scales::date_format('%B %d')) + if (uncertainty) { + gg.yearly <- gg.yearly + + ggplot2::geom_ribbon(ggplot2::aes(ymin = yearly_lower, + ymax = yearly_upper), + alpha = 0.2, + fill = "#0072B2", + na.rm = TRUE) + } + return(gg.yearly) +} + +#' Plot a custom seasonal component. +#' +#' @param m Prophet model object. +#' @param name String name of the seasonality. +#' @param uncertainty Boolean to plot uncertainty intervals. +#' +#' @return A ggplot2 plot. +#' +#' @keywords internal +plot_seasonality <- function(m, name, uncertainty = TRUE) { + # Compute seasonality from Jan 1 through a single period. + start <- set_date('2017-01-01') + period <- m$seasonalities[[name]]$period + end <- start + period * 24 * 3600 + plot.points <- 200 + days <- seq(from=start, to=end, length.out=plot.points) + df.y <- seasonality_plot_df(m, days) + seas <- predict_seasonal_components(m, df.y) + seas$ds <- df.y$ds + gg.s <- ggplot2::ggplot( + seas, ggplot2::aes_string(x = 'ds', y = name, group = 1)) + + ggplot2::geom_line(color = "#0072B2", na.rm = TRUE) + if (period <= 2) { + fmt.str <- '%T' + } else if (period < 14) { + fmt.str <- '%m/%d %R' + } else { + fmt.str <- '%m/%d' + } + gg.s <- gg.s + + ggplot2::scale_x_datetime(labels = scales::date_format(fmt.str)) + if (uncertainty) { + gg.s <- gg.s + + ggplot2::geom_ribbon( + ggplot2::aes_string( + ymin = paste0(name, '_lower'), ymax = paste0(name, '_upper') + ), + alpha = 0.2, + fill = "#0072B2", + na.rm = TRUE) + } + return(gg.s) +} + +#' Plot the prophet forecast. +#' +#' @param x Prophet object. +#' @param fcst Data frame returned by predict(m, df). +#' @param uncertainty Boolean indicating if the uncertainty interval for yhat +#' should be plotted. Must be present in fcst as yhat_lower and yhat_upper. +#' @param ... additional arguments +#' @importFrom dplyr "%>%" +#' @return A dygraph plot. +#' +#' @examples +#' \dontrun{ +#' history <- data.frame(ds = seq(as.Date('2015-01-01'), as.Date('2016-01-01'), by = 'd'), +#' y = sin(1:366/200) + rnorm(366)/10) +#' m <- prophet(history) +#' future <- make_future_dataframe(m, periods = 365) +#' forecast <- predict(m, future) +#' dyplot.prophet(m, forecast) +#' } +#' +#' @export +dyplot.prophet <- function(x, fcst, uncertainty=TRUE, + ...) +{ + forecast.label='Predicted' + actual.label='Actual' + # create data.frame for plotting + df <- df_for_plotting(x, fcst) + + # build variables to include, or not, the uncertainty data + if(uncertainty && exists("yhat_lower", where = df)) + { + colsToKeep <- c('y', 'yhat', 'yhat_lower', 'yhat_upper') + forecastCols <- c('yhat_lower', 'yhat', 'yhat_upper') + } else + { + colsToKeep <- c('y', 'yhat') + forecastCols <- c('yhat') + } + # convert to xts for easier date handling by dygraph + dfTS <- xts::xts(df %>% dplyr::select_(.dots=colsToKeep), order.by = df$ds) + + # base plot + dyBase <- dygraphs::dygraph(dfTS) + + presAnnotation <- function(dygraph, x, text) { + dygraph %>% + dygraphs::dyAnnotation(x, text, text, attachAtBottom = TRUE) + } + + dyBase <- dyBase %>% + dygraphs::dyOptions(colors = RColorBrewer::brewer.pal(3, "Set1")) %>% + # plot actual values + dygraphs::dySeries('y', label=actual.label) %>% + # plot forecast and ribbon + dygraphs::dySeries(forecastCols, label=forecast.label) %>% + # allow zooming + dygraphs::dyRangeSelector() %>% + # make unzoom button + dygraphs::dyUnzoom() + if (!is.null(x$holidays)) { + for (i in 1:nrow(x$holidays)) { + # make a gray line + dyBase <- dyBase %>% dygraphs::dyEvent(x$holidays$ds[i],color = "rgb(200,200,200)", strokePattern = "solid") + dyBase <- dyBase %>% dygraphs::dyAnnotation(x$holidays$ds[i], x$holidays$holiday[i], x$holidays$holiday[i], attachAtBottom = TRUE) + } + } + return(dyBase) +} + diff --git a/R/R/prophet.R b/R/R/prophet.R index 6d36854..401e24e 100644 --- a/R/R/prophet.R +++ b/R/R/prophet.R @@ -6,7 +6,7 @@ ## of patent rights can be found in the PATENTS file in the same directory. ## Makes R CMD CHECK happy due to dplyr syntax below -globalVariables(c( +utils::globalVariables(c( "ds", "y", "cap", ".", "component", "dow", "doy", "holiday", "holidays", "holidays_lower", "holidays_upper", "ix", "lower", "n", "stat", "trend", "row_number", "extra_regressors", "col", @@ -190,21 +190,21 @@ validate_column_name <- function( 'trend', 'seasonal', 'seasonalities', 'daily', 'weekly', 'yearly', 'holidays', 'zeros', 'extra_regressors', 'yhat' ) - rn_l = paste(reserved_names,"_lower",sep="") - rn_u = paste(reserved_names,"_upper",sep="") + rn_l = paste0(reserved_names,"_lower") + rn_u = paste0(reserved_names,"_upper") reserved_names = c(reserved_names, rn_l, rn_u, c("ds", "y", "cap", "floor", "y_scaled", "cap_scaled")) if(name %in% reserved_names){ stop("Name ", name, " is reserved.") } - if(check_holidays & !is.null(m$holidays) & + if(check_holidays && !is.null(m$holidays) && (name %in% unique(m$holidays$holiday))){ stop("Name ", name, " already used for a holiday.") } - if(check_seasonalities & (!is.null(m$seasonalities[[name]]))){ + if(check_seasonalities && (!is.null(m$seasonalities[[name]]))){ stop("Name ", name, " already used for a seasonality.") } - if(check_regressors & (!is.null(m$seasonalities[[name]]))){ + if(check_regressors && (!is.null(m$seasonalities[[name]]))){ stop("Name ", name, " already used for an added regressor.") } } @@ -274,11 +274,8 @@ set_date <- function(ds = NULL, tz = "GMT") { ds <- as.character(ds) } - if (min(nchar(ds)) < 12) { - ds <- as.POSIXct(ds, format = "%Y-%m-%d", tz = tz) - } else { - ds <- as.POSIXct(ds, format = "%Y-%m-%d %H:%M:%S", tz = tz) - } + fmt <- if (min(nchar(ds)) < 12) "%Y-%m-%d" else "%Y-%m-%d %H:%M:%S" + ds <- as.POSIXct(ds, format = fmt, tz = tz) attr(ds, "tzone") <- tz return(ds) } @@ -462,7 +459,7 @@ set_changepoints <- function(m) { m$changepoints.t <- sort( time_diff(m$changepoints, m$start, "secs")) / m$t.scale } else { - m$changepoints.t <- c(0) # dummy changepoint + m$changepoints.t <- 0 # dummy changepoint } return(m) } @@ -476,7 +473,7 @@ set_changepoints <- function(m) { #' @keywords internal get_changepoint_matrix <- function(m) { A <- matrix(0, nrow(m$history), length(m$changepoints.t)) - for (i in 1:length(m$changepoints.t)) { + for (i in seq_along(m$changepoints.t)) { A[m$history$t >= m$changepoints.t[i], i] <- 1 } return(A) @@ -494,7 +491,7 @@ get_changepoint_matrix <- function(m) { fourier_series <- function(dates, period, series.order) { t <- time_diff(dates, set_date('1970-01-01 00:00:00')) features <- matrix(0, length(t), 2 * series.order) - for (i in 1:series.order) { + for (i in seq_len(series.order)) { x <- as.numeric(2 * i * pi * t / period) features[, i * 2 - 1] <- sin(x) features[, i * 2] <- cos(x) @@ -514,7 +511,7 @@ fourier_series <- function(dates, period, series.order) { #' @keywords internal make_seasonality_features <- function(dates, period, series.order, prefix) { features <- fourier_series(dates, period, series.order) - colnames(features) <- paste(prefix, 1:ncol(features), sep = '_delim_') + colnames(features) <- paste(prefix, seq_len(ncol(features)), sep = '_delim_') return(data.frame(features)) } @@ -541,13 +538,13 @@ make_holiday_features <- function(m, dates) { && !is.na(.$upper_window)) { offsets <- seq(.$lower_window, .$upper_window) } else { - offsets <- c(0) + offsets <- 0 } names <- paste(.$holiday, '_delim_', ifelse(offsets < 0, '-', '+'), abs(offsets), sep = '') dplyr::data_frame(ds = .$ds + offsets * 24 * 3600, holiday = names) }) %>% - dplyr::mutate(x = 1.) %>% + dplyr::mutate(x = 1) %>% tidyr::spread(holiday, x, fill = 0) holiday.features <- data.frame(ds = set_date(dates)) %>% @@ -685,7 +682,7 @@ add_seasonality <- function(m, name, period, fourier.order, prior.scale = NULL) #' #' @keywords internal make_all_seasonality_features <- function(m, df) { - seasonal.features <- data.frame(row.names = 1:nrow(df)) + seasonal.features <- data.frame(row.names = seq_len(nrow(df))) prior.scales <- c() # Seasonality features @@ -713,7 +710,7 @@ make_all_seasonality_features <- function(m, df) { if (ncol(seasonal.features) == 0) { seasonal.features <- data.frame(zeros = rep(0, nrow(df))) - prior.scales <- c(1.) + prior.scales <- 1 } return(list(seasonal.features = seasonal.features, prior.scales = prior.scales)) @@ -1037,9 +1034,7 @@ predict.prophet <- function(object, df = NULL, ...) { cols <- c(cols, 'floor') } df <- df[cols] - df <- df %>% - dplyr::bind_cols(seasonal.components) %>% - dplyr::bind_cols(intervals) + df <- dplyr::bind_cols(df, seasonal.components, intervals) df$yhat <- df$trend + df$seasonal return(df) } @@ -1061,7 +1056,7 @@ piecewise_linear <- function(t, deltas, k, m, changepoint.ts) { # Get cumulative slope and intercept at each t k_t <- rep(k, length(t)) m_t <- rep(m, length(t)) - for (s in 1:length(changepoint.ts)) { + for (s in seq_along(changepoint.ts)) { indx <- t >= changepoint.ts[s] k_t[indx] <- k_t[indx] + deltas[s] m_t[indx] <- m_t[indx] + gammas[s] @@ -1086,14 +1081,14 @@ piecewise_logistic <- function(t, cap, deltas, k, m, changepoint.ts) { # Compute offset changes k.cum <- c(k, cumsum(deltas) + k) gammas <- rep(0, length(changepoint.ts)) - for (i in 1:length(changepoint.ts)) { + for (i in seq_along(changepoint.ts)) { gammas[i] <- ((changepoint.ts[i] - m - sum(gammas)) * (1 - k.cum[i] / k.cum[i + 1])) } # Get cumulative rate and offset at each t k_t <- rep(k, length(t)) m_t <- rep(m, length(t)) - for (s in 1:length(changepoint.ts)) { + for (s in seq_along(changepoint.ts)) { indx <- t >= changepoint.ts[s] k_t[indx] <- k_t[indx] + deltas[s] m_t[indx] <- m_t[indx] + gammas[s] @@ -1140,14 +1135,14 @@ predict_seasonal_components <- function(m, df) { upper.p <- (1 + m$interval.width)/2 components <- dplyr::data_frame(component = colnames(seasonal.features)) %>% - dplyr::mutate(col = 1:n()) %>% + dplyr::mutate(col = seq_len(n())) %>% tidyr::separate(component, c('component', 'part'), sep = "_delim_", extra = "merge", fill = "right") %>% dplyr::select(col, component) # Add total for all regression components components <- rbind( components, - data.frame(col = 1:ncol(seasonal.features), component = 'seasonal')) + data.frame(col = seq_len(ncol(seasonal.features)), component = 'seasonal')) # Add totals for seasonality, holiday, and extra regressors components <- add_group_component( components, 'seasonalities', names(m$seasonalities)) @@ -1164,7 +1159,7 @@ predict_seasonal_components <- function(m, df) { dplyr::group_by(component) %>% dplyr::do({ comp <- (as.matrix(seasonal.features[, .$col]) %*% t(m$params$beta[, .$col, drop = FALSE])) * m$y.scale - dplyr::data_frame(ix = 1:nrow(seasonal.features), + dplyr::data_frame(ix = seq_len(nrow(seasonal.features)), mean = rowMeans(comp, na.rm = TRUE), lower = apply(comp, 1, stats::quantile, lower.p, na.rm = TRUE), @@ -1218,9 +1213,9 @@ sample_posterior_predictive <- function(m, df) { "seasonal" = matrix(, nrow = nrow(df), ncol = nsamp), "yhat" = matrix(, nrow = nrow(df), ncol = nsamp)) - for (i in 1:n.iterations) { + for (i in seq_len(n.iterations)) { # For each set of parameters from MCMC (or just 1 set for MAP), - for (j in 1:samp.per.iter) { + for (j in seq_len(samp.per.iter)) { # Do a simulation with this set of parameters, sim <- sample_model(m, df, seasonal.features, i) # Store the results @@ -1388,310 +1383,6 @@ make_future_dataframe <- function(m, periods, freq = 'day', return(data.frame(ds = dates)) } -#' Merge history and forecast for plotting. -#' -#' @param m Prophet object. -#' @param fcst Data frame returned by prophet predict. -#' -#' @importFrom dplyr "%>%" -#' @keywords internal -df_for_plotting <- function(m, fcst) { - # Make sure there is no y in fcst - fcst$y <- NULL - df <- m$history %>% - dplyr::select(ds, y) %>% - dplyr::full_join(fcst, by = "ds") %>% - dplyr::arrange(ds) - return(df) -} - -#' Plot the prophet forecast. -#' -#' @param x Prophet object. -#' @param fcst Data frame returned by predict(m, df). -#' @param uncertainty Boolean indicating if the uncertainty interval for yhat -#' should be plotted. Must be present in fcst as yhat_lower and yhat_upper. -#' @param plot_cap Boolean indicating if the capacity should be shown in the -#' figure, if available. -#' @param xlabel Optional label for x-axis -#' @param ylabel Optional label for y-axis -#' @param ... additional arguments -#' -#' @return A ggplot2 plot. -#' -#' @examples -#' \dontrun{ -#' history <- data.frame(ds = seq(as.Date('2015-01-01'), as.Date('2016-01-01'), by = 'd'), -#' y = sin(1:366/200) + rnorm(366)/10) -#' m <- prophet(history) -#' future <- make_future_dataframe(m, periods = 365) -#' forecast <- predict(m, future) -#' plot(m, forecast) -#' } -#' -#' @export -plot.prophet <- function(x, fcst, uncertainty = TRUE, plot_cap = TRUE, - xlabel = 'ds', ylabel = 'y', ...) { - df <- df_for_plotting(x, fcst) - gg <- ggplot2::ggplot(df, ggplot2::aes(x = ds, y = y)) + - ggplot2::labs(x = xlabel, y = ylabel) - if (exists('cap', where = df) && plot_cap) { - gg <- gg + ggplot2::geom_line( - ggplot2::aes(y = cap), linetype = 'dashed', na.rm = TRUE) - } - if (x$logistic.floor && exists('floor', where = df) && plot_cap) { - gg <- gg + ggplot2::geom_line( - ggplot2::aes(y = floor), linetype = 'dashed', na.rm = TRUE) - } - if (uncertainty && exists('yhat_lower', where = df)) { - gg <- gg + - ggplot2::geom_ribbon(ggplot2::aes(ymin = yhat_lower, ymax = yhat_upper), - alpha = 0.2, - fill = "#0072B2", - na.rm = TRUE) - } - gg <- gg + - ggplot2::geom_point(na.rm=TRUE) + - ggplot2::geom_line(ggplot2::aes(y = yhat), color = "#0072B2", - na.rm = TRUE) + - ggplot2::theme(aspect.ratio = 3 / 5) - return(gg) -} - -#' Plot the components of a prophet forecast. -#' Prints a ggplot2 with panels for trend, weekly and yearly seasonalities if -#' present, and holidays if present. -#' -#' @param m Prophet object. -#' @param fcst Data frame returned by predict(m, df). -#' @param uncertainty Boolean indicating if the uncertainty interval should be -#' plotted for the trend, from fcst columns trend_lower and trend_upper. -#' @param plot_cap Boolean indicating if the capacity should be shown in the -#' figure, if available. -#' @param weekly_start Integer specifying the start day of the weekly -#' seasonality plot. 0 (default) starts the week on Sunday. 1 shifts by 1 day -#' to Monday, and so on. -#' @param yearly_start Integer specifying the start day of the yearly -#' seasonality plot. 0 (default) starts the year on Jan 1. 1 shifts by 1 day -#' to Jan 2, and so on. -#' -#' @return Invisibly return a list containing the plotted ggplot objects -#' -#' @export -#' @importFrom dplyr "%>%" -prophet_plot_components <- function( - m, fcst, uncertainty = TRUE, plot_cap = TRUE, weekly_start = 0, - yearly_start = 0 -) { - # Plot the trend - panels <- list(plot_forecast_component(fcst, 'trend', uncertainty, plot_cap)) - # Plot holiday components, if present. - if (!is.null(m$holidays) & ('holidays' %in% colnames(fcst))) { - panels[[length(panels) + 1]] <- plot_forecast_component( - fcst, 'holidays', uncertainty, FALSE) - } - # Plot weekly seasonality, if present - if ("weekly" %in% colnames(fcst)) { - panels[[length(panels) + 1]] <- plot_weekly(m, uncertainty, weekly_start) - } - # Plot yearly seasonality, if present - if ("yearly" %in% colnames(fcst)) { - panels[[length(panels) + 1]] <- plot_yearly(m, uncertainty, yearly_start) - } - # Plot other seasonalities - for (name in names(m$seasonalities)) { - if (!(name %in% c('weekly', 'yearly')) && - (name %in% colnames(fcst))) { - panels[[length(panels) + 1]] <- plot_seasonality(m, name, uncertainty) - } - } - # Plot extra regressors - if ((length(m$extra_regressors) > 0) - & ('extra_regressors' %in% colnames(fcst))) { - panels[[length(panels) + 1]] <- plot_forecast_component( - fcst, 'extra_regressors', uncertainty, FALSE) - } - - # Make the plot. - grid::grid.newpage() - grid::pushViewport(grid::viewport(layout = grid::grid.layout(length(panels), - 1))) - for (i in 1:length(panels)) { - print(panels[[i]], vp = grid::viewport(layout.pos.row = i, - layout.pos.col = 1)) - } - return(invisible(panels)) -} - -#' Plot a particular component of the forecast. -#' -#' @param fcst Dataframe output of `predict`. -#' @param name String name of the component to plot (column of fcst). -#' @param uncertainty Boolean to plot uncertainty intervals. -#' @param plot_cap Boolean indicating if the capacity should be shown in the -#' figure, if available. -#' -#' @return A ggplot2 plot. -#' -#' @export -plot_forecast_component <- function( - fcst, name, uncertainty = TRUE, plot_cap = FALSE -) { - gg.comp <- ggplot2::ggplot( - fcst, ggplot2::aes_string(x = 'ds', y = name, group = 1)) + - ggplot2::geom_line(color = "#0072B2", na.rm = TRUE) - if (exists('cap', where = fcst) && plot_cap) { - gg.comp <- gg.comp + ggplot2::geom_line( - ggplot2::aes(y = cap), linetype = 'dashed', na.rm = TRUE) - } - if (exists('floor', where = fcst) && plot_cap) { - gg.comp <- gg.comp + ggplot2::geom_line( - ggplot2::aes(y = floor), linetype = 'dashed', na.rm = TRUE) - } - if (uncertainty) { - gg.comp <- gg.comp + - ggplot2::geom_ribbon( - ggplot2::aes_string( - ymin = paste0(name, '_lower'), ymax = paste0(name, '_upper') - ), - alpha = 0.2, - fill = "#0072B2", - na.rm = TRUE) - } - return(gg.comp) -} - -#' Prepare dataframe for plotting seasonal components. -#' -#' @param m Prophet object. -#' @param ds Array of dates for column ds. -#' -#' @return A dataframe with seasonal components on ds. -#' -#' @keywords internal -seasonality_plot_df <- function(m, ds) { - df_list <- list(ds = ds, cap = 1, floor = 0) - for (name in names(m$extra_regressors)) { - df_list[[name]] <- 0 - } - df <- as.data.frame(df_list) - df <- setup_dataframe(m, df)$df - return(df) -} - -#' Plot the weekly component of the forecast. -#' -#' @param m Prophet model object -#' @param uncertainty Boolean to plot uncertainty intervals. -#' @param weekly_start Integer specifying the start day of the weekly -#' seasonality plot. 0 (default) starts the week on Sunday. 1 shifts by 1 day -#' to Monday, and so on. -#' -#' @return A ggplot2 plot. -#' -#' @keywords internal -plot_weekly <- function(m, uncertainty = TRUE, weekly_start = 0) { - # Compute weekly seasonality for a Sun-Sat sequence of dates. - days <- seq(set_date('2017-01-01'), by='d', length.out=7) + as.difftime( - weekly_start, units = "days") - df.w <- seasonality_plot_df(m, days) - seas <- predict_seasonal_components(m, df.w) - seas$dow <- factor(weekdays(df.w$ds), levels=weekdays(df.w$ds)) - - gg.weekly <- ggplot2::ggplot(seas, ggplot2::aes(x = dow, y = weekly, - group = 1)) + - ggplot2::geom_line(color = "#0072B2", na.rm = TRUE) + - ggplot2::labs(x = "Day of week") - if (uncertainty) { - gg.weekly <- gg.weekly + - ggplot2::geom_ribbon(ggplot2::aes(ymin = weekly_lower, - ymax = weekly_upper), - alpha = 0.2, - fill = "#0072B2", - na.rm = TRUE) - } - return(gg.weekly) -} - -#' Plot the yearly component of the forecast. -#' -#' @param m Prophet model object. -#' @param uncertainty Boolean to plot uncertainty intervals. -#' @param yearly_start Integer specifying the start day of the yearly -#' seasonality plot. 0 (default) starts the year on Jan 1. 1 shifts by 1 day -#' to Jan 2, and so on. -#' -#' @return A ggplot2 plot. -#' -#' @keywords internal -plot_yearly <- function(m, uncertainty = TRUE, yearly_start = 0) { - # Compute yearly seasonality for a Jan 1 - Dec 31 sequence of dates. - days <- seq(set_date('2017-01-01'), by='d', length.out=365) + as.difftime( - yearly_start, units = "days") - df.y <- seasonality_plot_df(m, days) - seas <- predict_seasonal_components(m, df.y) - seas$ds <- df.y$ds - - gg.yearly <- ggplot2::ggplot(seas, ggplot2::aes(x = ds, y = yearly, - group = 1)) + - ggplot2::geom_line(color = "#0072B2", na.rm = TRUE) + - ggplot2::labs(x = "Day of year") + - ggplot2::scale_x_datetime(labels = scales::date_format('%B %d')) - if (uncertainty) { - gg.yearly <- gg.yearly + - ggplot2::geom_ribbon(ggplot2::aes(ymin = yearly_lower, - ymax = yearly_upper), - alpha = 0.2, - fill = "#0072B2", - na.rm = TRUE) - } - return(gg.yearly) -} - -#' Plot a custom seasonal component. -#' -#' @param m Prophet model object. -#' @param name String name of the seasonality. -#' @param uncertainty Boolean to plot uncertainty intervals. -#' -#' @return A ggplot2 plot. -#' -#' @keywords internal -plot_seasonality <- function(m, name, uncertainty = TRUE) { - # Compute seasonality from Jan 1 through a single period. - start <- set_date('2017-01-01') - period <- m$seasonalities[[name]]$period - end <- start + period * 24 * 3600 - plot.points <- 200 - days <- seq(from=start, to=end, length.out=plot.points) - df.y <- seasonality_plot_df(m, days) - seas <- predict_seasonal_components(m, df.y) - seas$ds <- df.y$ds - gg.s <- ggplot2::ggplot( - seas, ggplot2::aes_string(x = 'ds', y = name, group = 1)) + - ggplot2::geom_line(color = "#0072B2", na.rm = TRUE) - if (period <= 2) { - fmt.str <- '%T' - } else if (period < 14) { - fmt.str <- '%m/%d %R' - } else { - fmt.str <- '%m/%d' - } - gg.s <- gg.s + - ggplot2::scale_x_datetime(labels = scales::date_format(fmt.str)) - if (uncertainty) { - gg.s <- gg.s + - ggplot2::geom_ribbon( - ggplot2::aes_string( - ymin = paste0(name, '_lower'), ymax = paste0(name, '_upper') - ), - alpha = 0.2, - fill = "#0072B2", - na.rm = TRUE) - } - return(gg.s) -} - #' Copy Prophet object. #' #' @param m Prophet model object. diff --git a/R/R/utils.R b/R/R/utils.R new file mode 100644 index 0000000..2116e91 --- /dev/null +++ b/R/R/utils.R @@ -0,0 +1,33 @@ +#' Get layers to overlay significant changepoints on prophet forecast plot. +#' +#' @param m Prophet model object. +#' @param threshold Numeric, changepoints where abs(delta) >= threshold are +#' significant. (Default 0.01) +#' @param cp_color Character, line color. (Default "red") +#' @param cp_linetype Character or integer, line type. (Default "dashed") +#' @param trend Logical, if FALSE, do not draw trend line. (Default TRUE) +#' @param ... Other arguments passed on to layers. +#' +#' @return A list of ggplot2 layers. +#' +#' @examples +#' \dontrun{ +#' plot(m, fcst) + layer_changepoints(m) +#' } +#' +#' @export +layer_changepoints <- function(m, threshold = 0.01, cp_color = "red", + cp_linetype = "dashed", trend = TRUE, ...) { + layers <- list() + if (trend) { + trend_layer <- ggplot2::geom_line( + ggplot2::aes_string("ds", "trend"), color = cp_color, ...) + layers <- append(layers, trend_layer) + } + signif_changepoints <- m$changepoints[abs(m$params$delta) >= threshold] + cp_layer <- ggplot2::geom_vline( + xintercept = as.integer(signif_changepoints), color = cp_color, + linetype = cp_linetype, ...) + layers <- append(layers, cp_layer) + return(layers) +} diff --git a/R/man/df_for_plotting.Rd b/R/man/df_for_plotting.Rd index f8d630d..3bbda91 100644 --- a/R/man/df_for_plotting.Rd +++ b/R/man/df_for_plotting.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/prophet.R +% Please edit documentation in R/plot.R \name{df_for_plotting} \alias{df_for_plotting} \title{Merge history and forecast for plotting.} diff --git a/R/man/dyplot.prophet.Rd b/R/man/dyplot.prophet.Rd new file mode 100644 index 0000000..002a7db --- /dev/null +++ b/R/man/dyplot.prophet.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot.R +\name{dyplot.prophet} +\alias{dyplot.prophet} +\title{Plot the prophet forecast.} +\usage{ +dyplot.prophet(x, fcst, uncertainty = TRUE, ...) +} +\arguments{ +\item{x}{Prophet object.} + +\item{fcst}{Data frame returned by predict(m, df).} + +\item{uncertainty}{Boolean indicating if the uncertainty interval for yhat +should be plotted. Must be present in fcst as yhat_lower and yhat_upper.} + +\item{...}{additional arguments} +} +\value{ +A dygraph plot. +} +\description{ +Plot the prophet forecast. +} +\examples{ +\dontrun{ +history <- data.frame(ds = seq(as.Date('2015-01-01'), as.Date('2016-01-01'), by = 'd'), + y = sin(1:366/200) + rnorm(366)/10) +m <- prophet(history) +future <- make_future_dataframe(m, periods = 365) +forecast <- predict(m, future) +dyplot.prophet(m, forecast) +} + +} diff --git a/R/man/layer_changepoints.Rd b/R/man/layer_changepoints.Rd new file mode 100644 index 0000000..674c717 --- /dev/null +++ b/R/man/layer_changepoints.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils.R +\name{layer_changepoints} +\alias{layer_changepoints} +\title{Get layers to overlay significant changepoints on prophet forecast plot.} +\usage{ +layer_changepoints(m, threshold = 0.01, cp_color = "red", + cp_linetype = "dashed", trend = TRUE, ...) +} +\arguments{ +\item{m}{Prophet model object.} + +\item{threshold}{Numeric, changepoints where abs(delta) >= threshold are +significant. (Default 0.01)} + +\item{cp_color}{Character, line color. (Default "red")} + +\item{cp_linetype}{Character or integer, line type. (Default "dashed")} + +\item{trend}{Logical, if FALSE, do not draw trend line. (Default TRUE)} + +\item{...}{Other arguments passed on to layers.} +} +\value{ +A list of ggplot2 layers. +} +\description{ +Get layers to overlay significant changepoints on prophet forecast plot. +} +\examples{ +\dontrun{ +plot(m, fcst) + layer_changepoints(m) +} + +} diff --git a/R/man/plot.prophet.Rd b/R/man/plot.prophet.Rd index 25c520e..38cd604 100644 --- a/R/man/plot.prophet.Rd +++ b/R/man/plot.prophet.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/prophet.R +% Please edit documentation in R/plot.R \name{plot.prophet} \alias{plot.prophet} \title{Plot the prophet forecast.} diff --git a/R/man/plot_forecast_component.Rd b/R/man/plot_forecast_component.Rd index f85f896..06c1fb1 100644 --- a/R/man/plot_forecast_component.Rd +++ b/R/man/plot_forecast_component.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/prophet.R +% Please edit documentation in R/plot.R \name{plot_forecast_component} \alias{plot_forecast_component} \title{Plot a particular component of the forecast.} diff --git a/R/man/plot_seasonality.Rd b/R/man/plot_seasonality.Rd index a92e8b5..e5f2df4 100644 --- a/R/man/plot_seasonality.Rd +++ b/R/man/plot_seasonality.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/prophet.R +% Please edit documentation in R/plot.R \name{plot_seasonality} \alias{plot_seasonality} \title{Plot a custom seasonal component.} diff --git a/R/man/plot_weekly.Rd b/R/man/plot_weekly.Rd index 6476b56..6092345 100644 --- a/R/man/plot_weekly.Rd +++ b/R/man/plot_weekly.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/prophet.R +% Please edit documentation in R/plot.R \name{plot_weekly} \alias{plot_weekly} \title{Plot the weekly component of the forecast.} diff --git a/R/man/plot_yearly.Rd b/R/man/plot_yearly.Rd index 9a7c534..7a227ed 100644 --- a/R/man/plot_yearly.Rd +++ b/R/man/plot_yearly.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/prophet.R +% Please edit documentation in R/plot.R \name{plot_yearly} \alias{plot_yearly} \title{Plot the yearly component of the forecast.} diff --git a/R/man/prophet_plot_components.Rd b/R/man/prophet_plot_components.Rd index eeeea29..4bf0f13 100644 --- a/R/man/prophet_plot_components.Rd +++ b/R/man/prophet_plot_components.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/prophet.R +% Please edit documentation in R/plot.R \name{prophet_plot_components} \alias{prophet_plot_components} \title{Plot the components of a prophet forecast. diff --git a/R/man/seasonality_plot_df.Rd b/R/man/seasonality_plot_df.Rd index 12d2da9..11def86 100644 --- a/R/man/seasonality_plot_df.Rd +++ b/R/man/seasonality_plot_df.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/prophet.R +% Please edit documentation in R/plot.R \name{seasonality_plot_df} \alias{seasonality_plot_df} \title{Prepare dataframe for plotting seasonal components.} diff --git a/README.md b/README.md index 6f3f90a..ec9dd9b 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,7 @@ # Prophet: Automatic Forecasting Procedure +[![Build Status](https://travis-ci.org/facebook/prophet.svg?branch=master)](https://travis-ci.org/facebook/prophet) + Prophet is a procedure for forecasting time series data. It is based on an additive model where non-linear trends are fit with yearly and weekly seasonality, plus holidays. It works best with daily periodicity data with at least one year of historical data. Prophet is robust to missing data, shifts in the trend, and large outliers. Prophet is [open source software](https://code.facebook.com/projects/) released by Facebook's [Core Data Science team](https://research.fb.com/category/data-science/). It is available for download on [CRAN](https://cran.r-project.org/package=prophet) and [PyPI](https://pypi.python.org/pypi/fbprophet/).