Merge of v0.3 into master

R/DESCRIPTION

@@ -1,21 +1,23 @@
 Package: prophet
 Title: Automatic Forecasting Procedure
-Version: 0.2.1.9000
-Date: 2017-11-08
+Version: 0.3
+Date: 2018-06-01
 Authors@R: c(
   person("Sean", "Taylor", email = "sjt@fb.com", role = c("cre", "aut")),
   person("Ben", "Letham", email = "bletham@fb.com", role = "aut")
   )
 Description: Implements a procedure for forecasting time series data 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.
+    an additive model where non-linear trends are fit with yearly, weekly, and
+    daily seasonality, plus holiday effects. It works best with time series
+    that have strong seasonal effects and several seasons of historical data.
+    Prophet is robust to missing data and shifts in the trend, and typically
+    handles outliers well.
 Depends:
     R (>= 3.2.3),
     Rcpp (>= 0.12.0)
 Imports:
     dplyr (>= 0.5.0),
+    dygraphs (>= 1.1.1.4),
     extraDistr,
     ggplot2,
     grid,

R/NAMESPACE

@@ -2,15 +2,17 @@
 
 S3method(plot,prophet)
 S3method(predict,prophet)
+export(add_changepoints_to_plot)
 export(add_regressor)
 export(add_seasonality)
 export(cross_validation)
 export(fit.prophet)
 export(make_future_dataframe)
+export(performance_metrics)
+export(plot_cross_validation_metric)
 export(plot_forecast_component)
 export(predictive_samples)
 export(prophet)
 export(prophet_plot_components)
-export(simulated_historical_forecasts)
 import(Rcpp)
 importFrom(dplyr,"%>%")

R/R/diagnostics.R

@@ -11,68 +11,81 @@ globalVariables(c(
 
 #' Generate cutoff dates
 #'
-#' @param df Dataframe with historical data
-#' @param horizon timediff forecast horizon
-#' @param k integer number of forecast points
+#' @param df Dataframe with historical data.
+#' @param horizon timediff forecast horizon.
+#' @param initial timediff initial window.
 #' @param period timediff Simulated forecasts are done with this period.
 #'
-#' @return Array of datetimes
+#' @return Array of datetimes.
 #'
 #' @keywords internal
-generate_cutoffs <- function(df, horizon, k, period) {
+generate_cutoffs <- function(df, horizon, initial, period) {
   # Last cutoff is (latest date in data) - (horizon).
   cutoff <- max(df$ds) - horizon
-  if (cutoff < min(df$ds)) {
-    stop('Less data than horizon.')
-  }
   tzone <- attr(cutoff, "tzone")  # Timezone is wiped by putting in array
   result <- c(cutoff)
-  if (k > 1) {
-    for (i in 2:k) {
-      cutoff <- cutoff - period
-      # If data does not exist in data range (cutoff, cutoff + horizon]
-      if (!any((df$ds > cutoff) & (df$ds <= cutoff + horizon))) {
+  while (result[length(result)] >= min(df$ds) + initial) {
+    cutoff <- cutoff - period
+    # If data does not exist in data range (cutoff, cutoff + horizon]
+    if (!any((df$ds > cutoff) & (df$ds <= cutoff + horizon))) {
         # Next cutoff point is 'closest date before cutoff in data - horizon'
         closest.date <- max(df$ds[df$ds <= cutoff])
         cutoff <- closest.date - horizon
-      }
-      if (cutoff < min(df$ds)) {
-        warning('Not enough data for requested number of cutoffs! Using ', i)
-        break
-      }
-      result <- c(result, cutoff)
     }
+    result <- c(result, cutoff)
+  }
+  result <- utils::head(result, -1)
+  if (length(result) == 0) {
+    stop(paste(
+      'Less data than horizon after initial window.',
+      'Make horizon or initial shorter.'
+    ))
   }
   # Reset timezones
   attr(result, "tzone") <- tzone
+  message(paste(
+    'Making', length(result), 'forecasts with cutoffs between',
+    result[length(result)], 'and', result[1]
+  ))
   return(rev(result))
 }
 
-#' Simulated historical forecasts.
+#' Cross-validation for time series.
 #'
-#' Make forecasts from k historical cutoff points, working backwards from
-#' (end - horizon) with a spacing of period between each cutoff.
+#' Computes forecasts from historical cutoff points. Beginning from
+#' (end - horizon), works backwards making cutoffs with a spacing of period
+#' until initial is reached.
+#'
+#' When period is equal to the time interval of the data, this is the
+#' technique described in https://robjhyndman.com/hyndsight/tscv/ .
 #'
 #' @param model Fitted Prophet model.
 #' @param horizon Integer size of the horizon
 #' @param units String unit of the horizon, e.g., "days", "secs".
-#' @param k integer number of forecast points
 #' @param period Integer amount of time between cutoff dates. Same units as
-#'  horizon. If not provided, will use 0.5 * horizon.
+#'  horizon. If not provided, 0.5 * horizon is used.
+#' @param initial Integer size of the first training period. If not provided,
+#'  3 * horizon is used. Same units as horizon.
 #'
 #' @return A dataframe with the forecast, actual value, and cutoff date.
 #'
 #' @export
-simulated_historical_forecasts <- function(model, horizon, units, k,
-                                           period = NULL) {
+cross_validation <- function(
+    model, horizon, units, period = NULL, initial = NULL) {
   df <- model$history
-  horizon <- as.difftime(horizon, units = units)
+  te <- max(df$ds)
+  ts <- min(df$ds)
   if (is.null(period)) {
-    period <- horizon / 2
-  } else {
-    period <- as.difftime(period, units = units)
+    period <- 0.5 * horizon
+  }
+  if (is.null(initial)) {
+    initial <- 3 * horizon
   }
-  cutoffs <- generate_cutoffs(df, horizon, k, period)
+  horizon.dt <- as.difftime(horizon, units = units)
+  initial.dt <- as.difftime(initial, units = units)
+  period.dt <- as.difftime(period, units = units)
+
+  cutoffs <- generate_cutoffs(df, horizon.dt, initial.dt, period.dt)
   predicts <- data.frame()
   for (i in 1:length(cutoffs)) {
     cutoff <- cutoffs[i]
@@ -80,11 +93,14 @@ simulated_historical_forecasts <- function(model, horizon, units, k,
     m <- prophet_copy(model, cutoff)
     # Train model
     history.c <- dplyr::filter(df, ds <= cutoff)
+    if (nrow(history.c) < 2) {
+      stop('Less than two datapoints before cutoff. Increase initial window.')
+    }
     m <- fit.prophet(m, history.c)
     # Calculate yhat
-    df.predict <- dplyr::filter(df, ds > cutoff, ds <= cutoff + horizon)
+    df.predict <- dplyr::filter(df, ds > cutoff, ds <= cutoff + horizon.dt)
     # Get the columns for the future dataframe
-    columns <- c('ds')
+    columns <- 'ds'
     if (m$growth == 'logistic') {
       columns <- c(columns, 'cap')
       if (m$logistic.floor) {
@@ -92,7 +108,7 @@ simulated_historical_forecasts <- function(model, horizon, units, k,
       }
     }
     columns <- c(columns, names(m$extra_regressors))
-    future <- df[columns]
+    future <- df.predict[columns]
     yhat <- stats::predict(m, future)
     # Merge yhat, y, and cutoff.
     df.c <- dplyr::inner_join(df.predict, yhat, by = "ds")
@@ -103,44 +119,200 @@ simulated_historical_forecasts <- function(model, horizon, units, k,
   return(predicts)
 }
 
-#' Cross-validation for time series.
+#' Copy Prophet object.
 #'
-#' Computes forecasts from historical cutoff points. Beginning from initial,
-#' makes cutoffs with a spacing of period up to (end - horizon).
+#' @param m Prophet model object.
+#' @param cutoff Date, possibly as string. Changepoints are only retained if
+#'  changepoints <= cutoff.
 #'
-#' When period is equal to the time interval of the data, this is the
-#' technique described in https://robjhyndman.com/hyndsight/tscv/ .
+#' @return An unfitted Prophet model object with the same parameters as the
+#'  input model.
 #'
-#' @param model Fitted Prophet model.
-#' @param horizon Integer size of the horizon
-#' @param units String unit of the horizon, e.g., "days", "secs".
-#' @param period Integer amount of time between cutoff dates. Same units as
-#'  horizon. If not provided, 0.5 * horizon is used.
-#' @param initial Integer size of the first training period. If not provided,
-#'  3 * horizon is used. Same units as horizon.
+#' @keywords internal
+prophet_copy <- function(m, cutoff = NULL) {
+  if (is.null(m$history)) {
+    stop("This is for copying a fitted Prophet object.")
+  }
+
+  if (m$specified.changepoints) {
+    changepoints <- m$changepoints
+    if (!is.null(cutoff)) {
+      cutoff <- set_date(cutoff)
+      changepoints <- changepoints[changepoints <= cutoff]
+    }
+  } else {
+    changepoints <- NULL
+  }
+  # Auto seasonalities are set to FALSE because they are already set in
+  # m$seasonalities.
+  m2 <- prophet(
+    growth = m$growth,
+    changepoints = changepoints,
+    n.changepoints = m$n.changepoints,
+    changepoint.range = m$changepoint.range,
+    yearly.seasonality = FALSE,
+    weekly.seasonality = FALSE,
+    daily.seasonality = FALSE,
+    holidays = m$holidays,
+    seasonality.mode = m$seasonality.mode,
+    seasonality.prior.scale = m$seasonality.prior.scale,
+    changepoint.prior.scale = m$changepoint.prior.scale,
+    holidays.prior.scale = m$holidays.prior.scale,
+    mcmc.samples = m$mcmc.samples,
+    interval.width = m$interval.width,
+    uncertainty.samples = m$uncertainty.samples,
+    fit = FALSE
+  )
+  m2$extra_regressors <- m$extra_regressors
+  m2$seasonalities <- m$seasonalities
+  return(m2)
+}
+
+#' Compute performance metrics from cross-validation results.
 #'
-#' @return A dataframe with the forecast, actual value, and cutoff date.
+#' Computes a suite of performance metrics on the output of cross-validation.
+#' By default the following metrics are included:
+#' 'mse': mean squared error
+#' 'rmse': root mean squared error
+#' 'mae': mean absolute error
+#' 'mape': mean percent error
+#' 'coverage': coverage of the upper and lower intervals
+#'
+#' A subset of these can be specified by passing a list of names as the
+#' `metrics` argument.
+#'
+#' Metrics are calculated over a rolling window of cross validation
+#' predictions, after sorting by horizon. The size of that window (number of
+#' simulated forecast points) is determined by the rolling_window argument,
+#' which specifies a proportion of simulated forecast points to include in
+#' each window. rolling_window=0 will compute it separately for each simulated
+#' forecast point (i.e., 'mse' will actually be squared error with no mean).
+#' The default of rolling_window=0.1 will use 10% of the rows in df in each
+#' window. rolling_window=1 will compute the metric across all simulated
+#' forecast points. The results are set to the right edge of the window.
+#'
+#' The output is a dataframe containing column 'horizon' along with columns
+#' for each of the metrics computed.
+#'
+#' @param df The dataframe returned by cross_validation.
+#' @param metrics An array of performance metrics to compute. If not provided,
+#'  will use c('mse', 'rmse', 'mae', 'mape', 'coverage').
+#' @param rolling_window Proportion of data to use in each rolling window for
+#'  computing the metrics. Should be in [0, 1].
+#'
+#' @return A dataframe with a column for each metric, and column 'horizon'.
 #'
 #' @export
-cross_validation <- function(
-    model, horizon, units, period = NULL, initial = NULL) {
-  te <- max(model$history$ds)
-  ts <- min(model$history$ds)
-  if (is.null(period)) {
-    period <- 0.5 * horizon
+performance_metrics <- function(df, metrics = NULL, rolling_window = 0.1) {
+  valid_metrics <- c('mse', 'rmse', 'mae', 'mape', 'coverage')
+  if (is.null(metrics)) {
+    metrics <- valid_metrics
   }
-  if (is.null(initial)) {
-    initial <- 3 * horizon
+  if (length(metrics) != length(unique(metrics))) {
+    stop('Input metrics must be an array of unique values.')
   }
-  horizon.dt <- as.difftime(horizon, units = units)
-  initial.dt <- as.difftime(initial, units = units)
-  period.dt <- as.difftime(period, units = units)
-  k <- ceiling(
-    as.double((te - horizon.dt) - (ts + initial.dt), units='secs') /
-    as.double(period.dt, units = 'secs')
-  )
-  if (k < 1) {
-    stop('Not enough data for specified horizon, period, and initial.')
+  if (!all(metrics %in% valid_metrics)) {
+    stop(
+      paste('Valid values for metrics are:', paste(metrics, collapse = ", "))
+    )
   }
-  return(simulated_historical_forecasts(model, horizon, units, k, period))
+  df_m <- df
+  df_m$horizon <- df_m$ds - df_m$cutoff
+  df_m <- df_m[order(df_m$horizon),]
+  # Window size
+  w <- as.integer(rolling_window * nrow(df_m))
+  w <- max(w, 1)
+  w <- min(w, nrow(df_m))
+  cols <- c('horizon')
+  for (metric in metrics) {
+    df_m[[metric]] <- get(metric)(df_m, w)
+    cols <- c(cols, metric)
+  }
+  df_m <- df_m[cols]
+  return(stats::na.omit(df_m))
+}
+
+#' Compute a rolling mean of x
+#'
+#' Right-aligned. Padded with NAs on the front so the output is the same
+#' size as x.
+#'
+#' @param x Array.
+#' @param w Integer window size (number of elements).
+#'
+#' @return Rolling mean of x with window size w.
+#'
+#' @keywords internal
+rolling_mean <- function(x, w) {
+  s <- cumsum(c(0, x))
+  prefix <- rep(NA, w - 1)
+  return(c(prefix, (s[(w + 1):length(s)] - s[1:(length(s) - w)]) / w))
+}
+
+# The functions below specify performance metrics for cross-validation results.
+# Each takes as input the output of cross_validation, and returns the statistic
+# as an array, given a window size for rolling aggregation.
+
+#' Mean squared error
+#'
+#' @param df Cross-validation results dataframe.
+#' @param w Aggregation window size.
+#'
+#' @return Array of mean squared errors.
+#'
+#' @keywords internal
+mse <- function(df, w) {
+  se <- (df$y - df$yhat) ** 2
+  return(rolling_mean(se, w))
+}
+
+#' Root mean squared error
+#'
+#' @param df Cross-validation results dataframe.
+#' @param w Aggregation window size.
+#'
+#' @return Array of root mean squared errors.
+#'
+#' @keywords internal
+rmse <- function(df, w) {
+  return(sqrt(mse(df, w)))
+}
+
+#' Mean absolute error
+#'
+#' @param df Cross-validation results dataframe.
+#' @param w Aggregation window size.
+#'
+#' @return Array of mean absolute errors.
+#'
+#' @keywords internal
+mae <- function(df, w) {
+  ae <- abs(df$y - df$yhat)
+  return(rolling_mean(ae, w))
+}
+
+#' Mean absolute percent error
+#'
+#' @param df Cross-validation results dataframe.
+#' @param w Aggregation window size.
+#'
+#' @return Array of mean absolute percent errors.
+#'
+#' @keywords internal
+mape <- function(df, w) {
+  ape <- abs((df$y - df$yhat) / df$y)
+  return(rolling_mean(ape, w))
+}
+
+#' Coverage
+#'
+#' @param df Cross-validation results dataframe.
+#' @param w Aggregation window size.
+#'
+#' @return Array of coverages
+#'
+#' @keywords internal
+coverage <- function(df, w) {
+  is_covered <- (df$y >= df$yhat_lower) & (df$y <= df$yhat_upper)
+  return(rolling_mean(is_covered, w))
 }

R/R/plot.R

@@ -0,0 +1,508 @@
+#' 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 whichever are available of: trend, holidays, weekly
+#' seasonality, yearly seasonality, and additive and multiplicative extra
+#' regressors.
+#'
+#' @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.
+#' @param render_plot Boolean indicating if the plots should be rendered.
+#'  Set to FALSE if you want the function to only return the list of panels.
+#'
+#' @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, render_plot = TRUE
+) {
+  # Plot the trend
+  panels <- list(
+    plot_forecast_component(m, 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(
+      m, 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
+  regressors <- list(additive = FALSE, multiplicative = FALSE)
+  for (name in names(m$extra_regressors)) {
+    regressors[[m$extra_regressors[[name]]$mode]] <- TRUE
+  }
+  for (mode in c('additive', 'multiplicative')) {
+    if ((regressors[[mode]]) &
+        (paste0('extra_regressors_', mode) %in% colnames(fcst))
+    ) {
+      panels[[length(panels) + 1]] <- plot_forecast_component(
+        m, fcst, paste0('extra_regressors_', mode), uncertainty, FALSE)
+    }
+  }
+
+  if (render_plot) {
+    # 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 m Prophet model
+#' @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(
+  m, 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)
+  }
+  if (name %in% m$component.modes$multiplicative) {
+    gg.comp <- gg.comp + ggplot2::scale_y_continuous(labels = scales::percent)
+  }
+  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)
+  }
+  if (m$seasonalities$weekly$mode == 'multiplicative') {
+    gg.weekly <- (
+      gg.weekly + ggplot2::scale_y_continuous(labels = scales::percent)
+    )
+  }
+  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)
+  }
+  if (m$seasonalities$yearly$mode == 'multiplicative') {
+    gg.yearly <- (
+      gg.yearly + ggplot2::scale_y_continuous(labels = scales::percent)
+    )
+  }
+  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)
+  }
+  if (m$seasonalities[[name]]$mode == 'multiplicative') {
+    gg.s <- gg.s + ggplot2::scale_y_continuous(labels = scales::percent)
+  }
+  return(gg.s)
+}
+
+#' 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) + add_changepoints_to_plot(m)
+#' }
+#'
+#' @export
+add_changepoints_to_plot <- 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)
+}
+
+#' 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 %>%
+    # plot actual values
+    dygraphs::dySeries(
+      'y', label=actual.label, color='black', drawPoints=TRUE, strokeWidth=0
+    ) %>%
+    # plot forecast and ribbon
+    dygraphs::dySeries(forecastCols, label=forecast.label, color='blue') %>%
+    # 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)
+}
+
+#' Plot a performance metric vs. forecast horizon from cross validation.
+
+#' Cross validation produces a collection of out-of-sample model predictions
+#' that can be compared to actual values, at a range of different horizons
+#' (distance from the cutoff). This computes a specified performance metric
+#' for each prediction, and aggregated over a rolling window with horizon.
+#'
+#' This uses fbprophet.diagnostics.performance_metrics to compute the metrics.
+#' Valid values of metric are 'mse', 'rmse', 'mae', 'mape', and 'coverage'.
+#'
+#' rolling_window is the proportion of data included in the rolling window of
+#' aggregation. The default value of 0.1 means 10% of data are included in the
+#' aggregation for computing the metric.
+#'
+#' As a concrete example, if metric='mse', then this plot will show the
+#' squared error for each cross validation prediction, along with the MSE
+#' averaged over rolling windows of 10% of the data.
+#'
+#' @param df_cv The output from fbprophet.diagnostics.cross_validation.
+#' @param metric Metric name, one of 'mse', 'rmse', 'mae', 'mape', 'coverage'.
+#' @param rolling_window Proportion of data to use for rolling average of
+#'  metric. In [0, 1]. Defaults to 0.1.
+#'
+#' @return A ggplot2 plot.
+#'
+#' @export
+plot_cross_validation_metric <- function(df_cv, metric, rolling_window=0.1) {
+  df_none <- performance_metrics(df_cv, metrics = metric, rolling_window = 0)
+  df_h <- performance_metrics(
+    df_cv, metrics = metric, rolling_window = rolling_window
+  )
+
+  # Better plotting of difftime
+  # Target ~10 ticks
+  tick_w <- max(as.double(df_none$horizon, units = 'secs')) / 10.
+  # Find the largest time resolution that has <1 unit per bin
+  dts <- c('days', 'hours', 'mins', 'secs')
+  dt_conversions <- c(
+    24 * 60 * 60,
+    60 * 60,
+    60,
+    1
+  )
+  for (i in seq_along(dts)) {
+    if (as.difftime(1, units = dts[i]) < as.difftime(tick_w, units = 'secs')) {
+      break
+    }
+  }
+  df_none$x_plt <- (
+    as.double(df_none$horizon, units = 'secs') / dt_conversions[i]
+  )
+  df_h$x_plt <- as.double(df_h$horizon, units = 'secs') / dt_conversions[i]
+
+  gg <- (
+    ggplot2::ggplot(df_none, ggplot2::aes_string(x = 'x_plt', y = metric)) +
+    ggplot2::labs(x = paste0('Horizon (', dts[i], ')'), y = metric) +
+    ggplot2::geom_point(color = 'gray') +
+    ggplot2::geom_line(
+      data = df_h, ggplot2::aes_string(x = 'x_plt', y = metric), color = 'blue'
+    ) +
+    ggplot2::theme(aspect.ratio = 3 / 5)
+  )
+
+  return(gg)
+}

R/R/zzz.R

@@ -6,9 +6,10 @@
 ## of patent rights can be found in the PATENTS file in the same directory.
 
 .onLoad <- function(libname, pkgname) {
-  .prophet.stan.models <- list(
-    "linear"=get_prophet_stan_model("linear"),
-    "logistic"=get_prophet_stan_model("logistic"))
-  assign(".prophet.stan.models", .prophet.stan.models,
-         envir=parent.env(environment()))
+  .prophet.stan.model <- get_prophet_stan_model()
+  assign(
+    ".prophet.stan.model",
+    .prophet.stan.model,
+    envir=parent.env(environment())
+  )
 }

R/inst/stan/prophet.stan

@@ -0,0 +1,126 @@
+functions {
+  matrix get_changepoint_matrix(vector t, vector t_change, int T, int S) {
+    // Assumes t and t_change are sorted.
+    matrix[T, S] A;
+    row_vector[S] a_row;
+    int cp_idx;
+
+    // Start with an empty matrix.
+    A = rep_matrix(0, T, S);
+    a_row = rep_row_vector(0, S);
+    cp_idx = 1;
+
+    // Fill in each row of A.
+    for (i in 1:T) {
+      while ((cp_idx <= S) && (t[i] >= t_change[cp_idx])) {
+        a_row[cp_idx] = 1;
+        cp_idx = cp_idx + 1;
+      }
+      A[i] = a_row;
+    }
+    return A;
+  }
+
+  // Logistic trend functions
+
+  vector logistic_gamma(real k, real m, vector delta, vector t_change, int S) {
+    vector[S] gamma;  // adjusted offsets, for piecewise continuity
+    vector[S + 1] k_s;  // actual rate in each segment
+    real m_pr;
+
+    // Compute the rate in each segment
+    k_s = append_row(k, k + cumulative_sum(delta));
+
+    // Piecewise offsets
+    m_pr = m; // The offset in the previous segment
+    for (i in 1:S) {
+      gamma[i] = (t_change[i] - m_pr) * (1 - k_s[i] / k_s[i + 1]);
+      m_pr = m_pr + gamma[i];  // update for the next segment
+    }
+    return gamma;
+  }
+
+  vector logistic_trend(
+    real k,
+    real m,
+    vector delta,
+    vector t,
+    vector cap,
+    matrix A,
+    vector t_change,
+    int S
+  ) {
+    vector[S] gamma;
+
+    gamma = logistic_gamma(k, m, delta, t_change, S);
+    return cap .* inv_logit((k + A * delta) .* (t - (m + A * gamma)));
+  }
+
+  // Linear trend function
+
+  vector linear_trend(
+    real k,
+    real m,
+    vector delta,
+    vector t,
+    matrix A,
+    vector t_change
+  ) {
+    return (k + A * delta) .* t + (m + A * (-t_change .* delta));
+  }
+}
+
+data {
+  int T;                // Number of time periods
+  int<lower=1> K;       // Number of regressors
+  vector[T] t;          // Time
+  vector[T] cap;        // Capacities for logistic trend
+  vector[T] y;          // Time series
+  int S;                // Number of changepoints
+  vector[S] t_change;   // Times of trend changepoints
+  matrix[T,K] X;        // Regressors
+  vector[K] sigmas;     // Scale on seasonality prior
+  real<lower=0> tau;    // Scale on changepoints prior
+  int trend_indicator;  // 0 for linear, 1 for logistic
+  vector[K] s_a;        // Indicator of additive features
+  vector[K] s_m;        // Indicator of multiplicative features
+}
+
+transformed data {
+  matrix[T, S] A;
+  A = get_changepoint_matrix(t, t_change, T, S);
+}
+
+parameters {
+  real k;                   // Base trend growth rate
+  real m;                   // Trend offset
+  vector[S] delta;          // Trend rate adjustments
+  real<lower=0> sigma_obs;  // Observation noise
+  vector[K] beta;           // Regressor coefficients
+}
+
+model {
+  //priors
+  k ~ normal(0, 5);
+  m ~ normal(0, 5);
+  delta ~ double_exponential(0, tau);
+  sigma_obs ~ normal(0, 0.5);
+  beta ~ normal(0, sigmas);
+
+  // Likelihood
+  if (trend_indicator == 0) {
+    y ~ normal(
+      linear_trend(k, m, delta, t, A, t_change)
+      .* (1 + X * (beta .* s_m))
+      + X * (beta .* s_a),
+      sigma_obs
+    );
+  } else if (trend_indicator == 1) {
+    y ~ normal(
+      logistic_trend(k, m, delta, t, cap, A, t_change, S)
+      .* (1 + X * (beta .* s_m))
+      + X * (beta .* s_a),
+      sigma_obs
+    );
+  }
+}

R/inst/stan/prophet_linear_growth.stan

@@ -1,40 +0,0 @@
-data {
-  int T;                                // Sample size
-  int<lower=1> K;                       // Number of seasonal vectors
-  vector[T] t;                            // Day
-  vector[T] y;                            // Time-series
-  int S;                                // Number of changepoints
-  matrix[T, S] A;                   // Split indicators
-  real t_change[S];                 // Index of changepoints
-  matrix[T,K] X;                // season vectors
-  vector[K] sigmas;              // scale on seasonality prior
-  real<lower=0> tau;                  // scale on changepoints prior
-}
-
-parameters {
-  real k;                            // Base growth rate
-  real m;                            // offset
-  vector[S] delta;                       // Rate adjustments
-  real<lower=0> sigma_obs;               // Observation noise (incl. seasonal variation)
-  vector[K] beta;                    // seasonal vector
-}
-
-transformed parameters {
-  vector[S] gamma;                  // adjusted offsets, for piecewise continuity
-
-  for (i in 1:S) {
-    gamma[i] = -t_change[i] * delta[i];
-  }
-}
-
-model {
-  //priors
-  k ~ normal(0, 5);
-  m ~ normal(0, 5);
-  delta ~ double_exponential(0, tau);
-  sigma_obs ~ normal(0, 0.5);
-  beta ~ normal(0, sigmas);
-
-  // Likelihood
-  y ~ normal((k + A * delta) .* t + (m + A * gamma) + X * beta, sigma_obs);
-}

R/inst/stan/prophet_logistic_growth.stan

@@ -1,52 +0,0 @@
-data {
-  int T;                                // Sample size
-  int<lower=1> K;                       // Number of seasonal vectors
-  vector[T] t;                            // Day
-  vector[T] cap;                          // Capacities
-  vector[T] y;                            // Time-series
-  int S;                                // Number of changepoints
-  matrix[T, S] A;                   // Split indicators
-  real t_change[S];                 // Index of changepoints
-  matrix[T,K] X;                    // season vectors
-  vector[K] sigmas;               // scale on seasonality prior
-  real<lower=0> tau;                  // scale on changepoints prior
-}
-
-parameters {
-  real k;                            // Base growth rate
-  real m;                            // offset
-  vector[S] delta;                       // Rate adjustments
-  real<lower=0> sigma_obs;               // Observation noise (incl. seasonal variation)
-  vector[K] beta;                    // seasonal vector
-}
-
-transformed parameters {
-  vector[S] gamma;                  // adjusted offsets, for piecewise continuity
-  vector[S + 1] k_s;                 // actual rate in each segment
-  real m_pr;
-
-  // Compute the rate in each segment
-  k_s[1] = k;
-  for (i in 1:S) {
-    k_s[i + 1] = k_s[i] + delta[i];
-  }
-
-  // Piecewise offsets
-  m_pr = m; // The offset in the previous segment
-  for (i in 1:S) {
-    gamma[i] = (t_change[i] - m_pr) * (1 - k_s[i] / k_s[i + 1]);
-    m_pr = m_pr + gamma[i];  // update for the next segment
-  }
-}
-
-model {
-  //priors
-  k ~ normal(0, 5);
-  m ~ normal(0, 5);
-  delta ~ double_exponential(0, tau);
-  sigma_obs ~ normal(0, 0.1);
-  beta ~ normal(0, sigmas);
-
-  // Likelihood
-  y ~ normal(cap ./ (1 + exp(-(k + A * delta) .* (t - (m + A * gamma)))) + X * beta, sigma_obs);
-}

R/man/add_changepoints_to_plot.Rd

@@ -0,0 +1,35 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot.R
+\name{add_changepoints_to_plot}
+\alias{add_changepoints_to_plot}
+\title{Get layers to overlay significant changepoints on prophet forecast plot.}
+\usage{
+add_changepoints_to_plot(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) + add_changepoints_to_plot(m)
+}
+
+}

R/man/add_regressor.Rd

@@ -4,7 +4,8 @@
 \alias{add_regressor}
 \title{Add an additional regressor to be used for fitting and predicting.}
 \usage{
-add_regressor(m, name, prior.scale = NULL, standardize = "auto")
+add_regressor(m, name, prior.scale = NULL, standardize = "auto",
+  mode = NULL)
 }
 \arguments{
 \item{m}{Prophet object.}
@@ -17,6 +18,9 @@ holidays.prior.scale will be used.}
 \item{standardize}{Bool, specify whether this regressor will be standardized
 prior to fitting. Can be 'auto' (standardize if not binary), True, or
 False.}
+
+\item{mode}{Optional, 'additive' or 'multiplicative'. Defaults to
+m$seasonality.mode.}
 }
 \value{
 The prophet model with the regressor added.
@@ -28,4 +32,8 @@ regressor will be standardized unless it is binary. The regression
 coefficient is given a prior with the specified scale parameter.
 Decreasing the prior scale will add additional regularization. If no
 prior scale is provided, holidays.prior.scale will be used.
+Mode can be specified as either 'additive' or 'multiplicative'. If not
+specified, m$seasonality.mode will be used. 'additive' means the effect of
+the regressor will be added to the trend, 'multiplicative' means it will
+multiply the trend.
 }

R/man/add_seasonality.Rd

@@ -5,7 +5,8 @@
 \title{Add a seasonal component with specified period, number of Fourier
 components, and prior scale.}
 \usage{
-add_seasonality(m, name, period, fourier.order, prior.scale = NULL)
+add_seasonality(m, name, period, fourier.order, prior.scale = NULL,
+  mode = NULL)
 }
 \arguments{
 \item{m}{Prophet object.}
@@ -16,7 +17,9 @@ add_seasonality(m, name, period, fourier.order, prior.scale = NULL)
 
 \item{fourier.order}{Int number of Fourier components to use.}
 
-\item{prior.scale}{Float prior scale for this component.}
+\item{prior.scale}{Optional float prior scale for this component.}
+
+\item{mode}{Optional 'additive' or 'multiplicative'.}
 }
 \value{
 The prophet model with the seasonality added.
@@ -30,4 +33,9 @@ seasonalities are 10 and 3 respectively.
 Increasing prior scale will allow this seasonality component more
 flexibility, decreasing will dampen it. If not provided, will use the
 seasonality.prior.scale provided on Prophet initialization (defaults to 10).
+
+Mode can be specified as either 'additive' or 'multiplicative'. If not
+specified, m$seasonality.mode will be used (defaults to 'additive').
+Additive means the seasonality will be added to the trend, multiplicative
+means it will multiply the trend.
 }

R/man/compile_stan_model.Rd

@@ -4,7 +4,7 @@
 \alias{compile_stan_model}
 \title{Compile Stan model}
 \usage{
-compile_stan_model(model)
+compile_stan_model()
 }
 \arguments{
 \item{model}{String 'linear' or 'logistic' to specify a linear or logistic

R/man/coverage.Rd

@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/diagnostics.R
+\name{coverage}
+\alias{coverage}
+\title{Coverage}
+\usage{
+coverage(df, w)
+}
+\arguments{
+\item{df}{Cross-validation results dataframe.}
+
+\item{w}{Aggregation window size.}
+}
+\value{
+Array of coverages
+}
+\description{
+Coverage
+}
+\keyword{internal}

R/man/cross_validation.Rd

@@ -23,8 +23,9 @@ horizon. If not provided, 0.5 * horizon is used.}
 A dataframe with the forecast, actual value, and cutoff date.
 }
 \description{
-Computes forecasts from historical cutoff points. Beginning from initial,
-makes cutoffs with a spacing of period up to (end - horizon).
+Computes forecasts from historical cutoff points. Beginning from
+(end - horizon), works backwards making cutoffs with a spacing of period
+until initial is reached.
 }
 \details{
 When period is equal to the time interval of the data, this is the

R/man/dyplot.prophet.Rd

@@ -0,0 +1,36 @@
+% 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)
+}
+
+}

R/man/generate_cutoffs.Rd

@@ -4,19 +4,19 @@
 \alias{generate_cutoffs}
 \title{Generate cutoff dates}
 \usage{
-generate_cutoffs(df, horizon, k, period)
+generate_cutoffs(df, horizon, initial, period)
 }
 \arguments{
-\item{df}{Dataframe with historical data}
+\item{df}{Dataframe with historical data.}
 
-\item{horizon}{timediff forecast horizon}
+\item{horizon}{timediff forecast horizon.}
 
-\item{k}{integer number of forecast points}
+\item{initial}{timediff initial window.}
 
 \item{period}{timediff Simulated forecasts are done with this period.}
 }
 \value{
-Array of datetimes
+Array of datetimes.
 }
 \description{
 Generate cutoff dates

R/man/get_changepoint_matrix.Rd

@@ -1,18 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/prophet.R
-\name{get_changepoint_matrix}
-\alias{get_changepoint_matrix}
-\title{Gets changepoint matrix for history dataframe.}
-\usage{
-get_changepoint_matrix(m)
-}
-\arguments{
-\item{m}{Prophet object.}
-}
-\value{
-array of indexes.
-}
-\description{
-Gets changepoint matrix for history dataframe.
-}
-\keyword{internal}

R/man/get_prophet_stan_model.Rd

@@ -4,7 +4,7 @@
 \alias{get_prophet_stan_model}
 \title{Load compiled Stan model}
 \usage{
-get_prophet_stan_model(model)
+get_prophet_stan_model()
 }
 \arguments{
 \item{model}{String 'linear' or 'logistic' to specify a linear or logistic

R/man/mae.Rd

@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/diagnostics.R
+\name{mae}
+\alias{mae}
+\title{Mean absolute error}
+\usage{
+mae(df, w)
+}
+\arguments{
+\item{df}{Cross-validation results dataframe.}
+
+\item{w}{Aggregation window size.}
+}
+\value{
+Array of mean absolute errors.
+}
+\description{
+Mean absolute error
+}
+\keyword{internal}

R/man/make_all_seasonality_features.Rd

@@ -18,6 +18,10 @@ List with items
  seasonal.features: Dataframe with regressor features,
  prior.scales: Array of prior scales for each colum of the features
    dataframe.
+ component.cols: Dataframe with indicators for which regression components
+   correspond to which columns.
+ modes: List with keys 'additive' and 'multiplicative' with arrays of
+   component names for each mode of seasonality.
 }
 \description{
 Dataframe with seasonality features.

R/man/make_holiday_features.Rd

@@ -15,6 +15,7 @@ make_holiday_features(m, dates)
 A list with entries
  holiday.features: dataframe with a column for each holiday.
  prior.scales: array of prior scales for each holiday column.
+ holiday.names: array of names of all holidays.
 }
 \description{
 Construct a matrix of holiday features.

R/man/mape.Rd

@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/diagnostics.R
+\name{mape}
+\alias{mape}
+\title{Mean absolute percent error}
+\usage{
+mape(df, w)
+}
+\arguments{
+\item{df}{Cross-validation results dataframe.}
+
+\item{w}{Aggregation window size.}
+}
+\value{
+Array of mean absolute percent errors.
+}
+\description{
+Mean absolute percent error
+}
+\keyword{internal}

R/man/mse.Rd

@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/diagnostics.R
+\name{mse}
+\alias{mse}
+\title{Mean squared error}
+\usage{
+mse(df, w)
+}
+\arguments{
+\item{df}{Cross-validation results dataframe.}
+
+\item{w}{Aggregation window size.}
+}
+\value{
+Array of mean squared errors.
+}
+\description{
+Mean squared error
+}
+\keyword{internal}

R/man/performance_metrics.Rd

@@ -0,0 +1,46 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/diagnostics.R
+\name{performance_metrics}
+\alias{performance_metrics}
+\title{Compute performance metrics from cross-validation results.}
+\usage{
+performance_metrics(df, metrics = NULL, rolling_window = 0.1)
+}
+\arguments{
+\item{df}{The dataframe returned by cross_validation.}
+
+\item{metrics}{An array of performance metrics to compute. If not provided,
+will use c('mse', 'rmse', 'mae', 'mape', 'coverage').}
+
+\item{rolling_window}{Proportion of data to use in each rolling window for
+computing the metrics. Should be in [0, 1].}
+}
+\value{
+A dataframe with a column for each metric, and column 'horizon'.
+}
+\description{
+Computes a suite of performance metrics on the output of cross-validation.
+By default the following metrics are included:
+'mse': mean squared error
+'rmse': root mean squared error
+'mae': mean absolute error
+'mape': mean percent error
+'coverage': coverage of the upper and lower intervals
+}
+\details{
+A subset of these can be specified by passing a list of names as the
+`metrics` argument.
+
+Metrics are calculated over a rolling window of cross validation
+predictions, after sorting by horizon. The size of that window (number of
+simulated forecast points) is determined by the rolling_window argument,
+which specifies a proportion of simulated forecast points to include in
+each window. rolling_window=0 will compute it separately for each simulated
+forecast point (i.e., 'mse' will actually be squared error with no mean).
+The default of rolling_window=0.1 will use 10% of the rows in df in each
+window. rolling_window=1 will compute the metric across all simulated
+forecast points. The results are set to the right edge of the window.
+
+The output is a dataframe containing column 'horizon' along with columns
+for each of the metrics computed.
+}

R/man/plot_cross_validation_metric.Rd

@@ -0,0 +1,36 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot.R
+\name{plot_cross_validation_metric}
+\alias{plot_cross_validation_metric}
+\title{Plot a performance metric vs. forecast horizon from cross validation.
+Cross validation produces a collection of out-of-sample model predictions
+that can be compared to actual values, at a range of different horizons
+(distance from the cutoff). This computes a specified performance metric
+for each prediction, and aggregated over a rolling window with horizon.}
+\usage{
+plot_cross_validation_metric(df_cv, metric, rolling_window = 0.1)
+}
+\arguments{
+\item{df_cv}{The output from fbprophet.diagnostics.cross_validation.}
+
+\item{metric}{Metric name, one of 'mse', 'rmse', 'mae', 'mape', 'coverage'.}
+
+\item{rolling_window}{Proportion of data to use for rolling average of
+metric. In [0, 1]. Defaults to 0.1.}
+}
+\value{
+A ggplot2 plot.
+}
+\description{
+This uses fbprophet.diagnostics.performance_metrics to compute the metrics.
+Valid values of metric are 'mse', 'rmse', 'mae', 'mape', and 'coverage'.
+}
+\details{
+rolling_window is the proportion of data included in the rolling window of
+aggregation. The default value of 0.1 means 10% of data are included in the
+aggregation for computing the metric.
+
+As a concrete example, if metric='mse', then this plot will show the
+squared error for each cross validation prediction, along with the MSE
+averaged over rolling windows of 10% of the data.
+}

R/man/plot_forecast_component.Rd

@@ -4,9 +4,11 @@
 \alias{plot_forecast_component}
 \title{Plot a particular component of the forecast.}
 \usage{
-plot_forecast_component(fcst, name, uncertainty = TRUE, plot_cap = FALSE)
+plot_forecast_component(m, fcst, name, uncertainty = TRUE, plot_cap = FALSE)
 }
 \arguments{
+\item{m}{Prophet model}
+
 \item{fcst}{Dataframe output of `predict`.}
 
 \item{name}{String name of the component to plot (column of fcst).}

R/man/predictive_samples.Rd

@@ -13,9 +13,8 @@ predictive_samples(m, df)
 (column cap) if logistic growth.}
 }
 \value{
-A list with items "trend", "seasonal", and "yhat" containing
- posterior predictive samples for that component. "seasonal" is the sum
- of seasonalities, holidays, and added regressors.
+A list with items "trend" and "yhat" containing
+ posterior predictive samples for that component.
 }
 \description{
 Sample from the posterior predictive distribution.

R/man/prophet.Rd

@@ -5,9 +5,10 @@
 \title{Prophet forecaster.}
 \usage{
 prophet(df = NULL, growth = "linear", changepoints = NULL,
-  n.changepoints = 25, yearly.seasonality = "auto",
-  weekly.seasonality = "auto", daily.seasonality = "auto",
-  holidays = NULL, seasonality.prior.scale = 10,
+  n.changepoints = 25, changepoint.range = 0.8,
+  yearly.seasonality = "auto", weekly.seasonality = "auto",
+  daily.seasonality = "auto", holidays = NULL,
+  seasonality.mode = "additive", seasonality.prior.scale = 10,
   holidays.prior.scale = 10, changepoint.prior.scale = 0.05,
   mcmc.samples = 0, interval.width = 0.8, uncertainty.samples = 1000,
   fit = TRUE, ...)
@@ -29,7 +30,11 @@ automatically.}
 \item{n.changepoints}{Number of potential changepoints to include. Not used
 if input `changepoints` is supplied. If `changepoints` is not supplied,
 then n.changepoints potential changepoints are selected uniformly from the
-first 80 percent of df$ds.}
+first `changepoint.range` proportion of df$ds.}
+
+\item{changepoint.range}{Proportion of history in which trend changepoints
+will be estimated. Defaults to 0.8 for the first 80%. Not used if
+`changepoints` is specified.}
 
 \item{yearly.seasonality}{Fit yearly seasonality. Can be 'auto', TRUE,
 FALSE, or a number of Fourier terms to generate.}
@@ -46,6 +51,8 @@ range of days around the date to be included as holidays. lower_window=-2
 will include 2 days prior to the date as holidays. Also optionally can have
 a column prior_scale specifying the prior scale for each holiday.}
 
+\item{seasonality.mode}{'additive' (default) or 'multiplicative'.}
+
 \item{seasonality.prior.scale}{Parameter modulating the strength of the
 seasonality model. Larger values allow the model to fit larger seasonal
 fluctuations, smaller values dampen the seasonality. Can be specified for

R/man/prophet_copy.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/diagnostics.R
 \name{prophet_copy}
 \alias{prophet_copy}
 \title{Copy Prophet object.}

R/man/prophet_plot_components.Rd

@@ -3,11 +3,12 @@
 \name{prophet_plot_components}
 \alias{prophet_plot_components}
 \title{Plot the components of a prophet forecast.
-Prints a ggplot2 with panels for trend, weekly and yearly seasonalities if
-present, and holidays if present.}
+Prints a ggplot2 with whichever are available of: trend, holidays, weekly
+seasonality, yearly seasonality, and additive and multiplicative extra
+regressors.}
 \usage{
 prophet_plot_components(m, fcst, uncertainty = TRUE, plot_cap = TRUE,
-  weekly_start = 0, yearly_start = 0)
+  weekly_start = 0, yearly_start = 0, render_plot = TRUE)
 }
 \arguments{
 \item{m}{Prophet object.}
@@ -27,12 +28,16 @@ to Monday, and so on.}
 \item{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.}
+
+\item{render_plot}{Boolean indicating if the plots should be rendered.
+Set to FALSE if you want the function to only return the list of panels.}
 }
 \value{
 Invisibly return a list containing the plotted ggplot objects
 }
 \description{
 Plot the components of a prophet forecast.
-Prints a ggplot2 with panels for trend, weekly and yearly seasonalities if
-present, and holidays if present.
+Prints a ggplot2 with whichever are available of: trend, holidays, weekly
+seasonality, yearly seasonality, and additive and multiplicative extra
+regressors.
 }

R/man/regressor_column_matrix.Rd

@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/prophet.R
+\name{regressor_column_matrix}
+\alias{regressor_column_matrix}
+\title{Dataframe indicating which columns of the feature matrix correspond to
+which seasonality/regressor components.}
+\usage{
+regressor_column_matrix(m, seasonal.features, modes)
+}
+\arguments{
+\item{m}{Prophet object.}
+
+\item{seasonal.features}{Constructed seasonal features dataframe.}
+
+\item{modes}{List with keys 'additive' and 'multiplicative' with arrays of
+component names for each mode of seasonality.}
+}
+\value{
+List with items
+ component.cols: A binary indicator dataframe with columns seasonal
+   components and rows columns in seasonal.features. Entry is 1 if that
+   column is used in that component.
+ modes: Updated input with combination components.
+}
+\description{
+Includes combination components, like 'additive_terms'. These combination
+components will be added to the 'modes' input.
+}
+\keyword{internal}

R/man/rmse.Rd

@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/diagnostics.R
+\name{rmse}
+\alias{rmse}
+\title{Root mean squared error}
+\usage{
+rmse(df, w)
+}
+\arguments{
+\item{df}{Cross-validation results dataframe.}
+
+\item{w}{Aggregation window size.}
+}
+\value{
+Array of root mean squared errors.
+}
+\description{
+Root mean squared error
+}
+\keyword{internal}

R/man/rolling_mean.Rd

@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/diagnostics.R
+\name{rolling_mean}
+\alias{rolling_mean}
+\title{Compute a rolling mean of x}
+\usage{
+rolling_mean(x, w)
+}
+\arguments{
+\item{x}{Array.}
+
+\item{w}{Integer window size (number of elements).}
+}
+\value{
+Rolling mean of x with window size w.
+}
+\description{
+Right-aligned. Padded with NAs on the front so the output is the same
+size as x.
+}
+\keyword{internal}

R/man/sample_model.Rd

@@ -4,7 +4,7 @@
 \alias{sample_model}
 \title{Simulate observations from the extrapolated generative model.}
 \usage{
-sample_model(m, df, seasonal.features, iteration)
+sample_model(m, df, seasonal.features, iteration, s_a, s_m)
 }
 \arguments{
 \item{m}{Prophet object.}
@@ -14,9 +14,13 @@ sample_model(m, df, seasonal.features, iteration)
 \item{seasonal.features}{Data frame of seasonal features}
 
 \item{iteration}{Int sampling iteration to use parameters from.}
+
+\item{s_a}{Indicator vector for additive components}
+
+\item{s_m}{Indicator vector for multiplicative components}
 }
 \value{
-List of trend, seasonality, and yhat, each a vector like df$t.
+List of trend and yhat, each a vector like df$t.
 }
 \description{
 Simulate observations from the extrapolated generative model.

R/man/sample_posterior_predictive.Rd

@@ -12,7 +12,8 @@ sample_posterior_predictive(m, df)
 \item{df}{Prediction dataframe.}
 }
 \value{
-List with posterior predictive samples for each component.
+List with posterior predictive samples for the forecast yhat and
+ for the trend component.
 }
 \description{
 Prophet posterior predictive samples.

R/man/simulated_historical_forecasts.Rd

@@ -1,27 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/diagnostics.R
-\name{simulated_historical_forecasts}
-\alias{simulated_historical_forecasts}
-\title{Simulated historical forecasts.}
-\usage{
-simulated_historical_forecasts(model, horizon, units, k, period = NULL)
-}
-\arguments{
-\item{model}{Fitted Prophet model.}
-
-\item{horizon}{Integer size of the horizon}
-
-\item{units}{String unit of the horizon, e.g., "days", "secs".}
-
-\item{k}{integer number of forecast points}
-
-\item{period}{Integer amount of time between cutoff dates. Same units as
-horizon. If not provided, will use 0.5 * horizon.}
-}
-\value{
-A dataframe with the forecast, actual value, and cutoff date.
-}
-\description{
-Make forecasts from k historical cutoff points, working backwards from
-(end - horizon) with a spacing of period between each cutoff.
-}

R/src/install.libs.R

@@ -1,26 +1,21 @@
 
 
-packageStartupMessage('Compiling models (this will take a minute...)')
+packageStartupMessage('Compiling model (this will take a minute...)')
 
 dest <- file.path(R_PACKAGE_DIR, paste0('libs', R_ARCH))
 dir.create(dest, recursive = TRUE, showWarnings = FALSE)
 
-packageStartupMessage(paste('Writing models to:', dest))
+packageStartupMessage(paste('Writing model to:', dest))
 packageStartupMessage(paste('Compiling using binary:', R.home('bin')))
 
-logistic.growth.src <- file.path(R_PACKAGE_SOURCE, 'inst', 'stan', 'prophet_logistic_growth.stan')
-logistic.growth.binary <- file.path(dest, 'prophet_logistic_growth.RData')
-logistic.growth.stanc <- rstan::stanc(logistic.growth.src)
-logistic.growth.stanm <- rstan::stan_model(stanc_ret = logistic.growth.stanc,
-                                           model_name = 'logistic_growth')
-save('logistic.growth.stanm', file = logistic.growth.binary)
+model.src <- file.path(R_PACKAGE_SOURCE, 'inst', 'stan', 'prophet.stan')
+model.binary <- file.path(dest, 'prophet_stan_model.RData')
+model.stanc <- rstan::stanc(model.src)
+model.stanm <- rstan::stan_model(
+  stanc_ret = model.stanc,
+  model_name = 'prophet_model'
+)
+save('model.stanm', file = model.binary)
 
-linear.growth.src <- file.path(R_PACKAGE_SOURCE, 'inst', 'stan', 'prophet_linear_growth.stan')
-linear.growth.binary <- file.path(dest, 'prophet_linear_growth.RData')
-linear.growth.stanc <- rstan::stanc(linear.growth.src)
-linear.growth.stanm <- rstan::stan_model(stanc_ret = linear.growth.stanc,
-                                         model_name = 'linear_growth')
-save('linear.growth.stanm', file = linear.growth.binary)
-
-packageStartupMessage('------ Models successfully compiled!')
+packageStartupMessage('------ Model successfully compiled!')
 packageStartupMessage('You can ignore any compiler warnings above.')

R/tests/testthat/test_diagnostics.R

@@ -4,50 +4,54 @@ context("Prophet diagnostics tests")
 ## Makes R CMD CHECK happy due to dplyr syntax below
 globalVariables(c("y", "yhat"))
 
-DATA <- head(read.csv('data.csv'), 100)
-DATA$ds <- as.Date(DATA$ds)
+DATA_all <- read.csv('data.csv')
+DATA_all$ds <- as.Date(DATA_all$ds)
+DATA <- head(DATA_all, 100)
 
-test_that("simulated_historical_forecasts", {
+test_that("cross_validation", {
   skip_if_not(Sys.getenv('R_ARCH') != '/i386')
   m <- prophet(DATA)
-  k <- 2
-  for (p in c(1, 10)) {
-    for (h in c(1, 3)) {
-      df.shf <- simulated_historical_forecasts(
-        m, horizon = h, units = 'days', k = k, period = p)
-      # All cutoff dates should be less than ds dates
-      expect_true(all(df.shf$cutoff < df.shf$ds))
-      # The unique size of output cutoff should be equal to 'k'
-      expect_equal(length(unique(df.shf$cutoff)), k)
-      expect_equal(max(df.shf$ds - df.shf$cutoff),
-                   as.difftime(h, units = 'days'))
-      dc <- diff(df.shf$cutoff)
-      dc <- min(dc[dc > 0])
-      expect_true(dc >= as.difftime(p, units = 'days'))
-      # Each y in df_shf and DATA with same ds should be equal
-      df.merged <- dplyr::left_join(df.shf, m$history, by="ds")
-      expect_equal(sum((df.merged$y.x - df.merged$y.y) ** 2), 0)
-    }
-  }
+  # Calculate the number of cutoff points
+  te <- max(DATA$ds)
+  ts <- min(DATA$ds)
+  horizon <- as.difftime(4, units = "days")
+  period <- as.difftime(10, units = "days")
+  initial <- as.difftime(115, units = "days")
+  df.cv <- cross_validation(
+    m, horizon = 4, units = "days", period = 10, initial = 115)
+  expect_equal(length(unique(df.cv$cutoff)), 3)
+  expect_equal(max(df.cv$ds - df.cv$cutoff), horizon)
+  expect_true(min(df.cv$cutoff) >= ts + initial)
+  dc <- diff(df.cv$cutoff)
+  dc <- min(dc[dc > 0])
+  expect_true(dc >= period)
+  expect_true(all(df.cv$cutoff < df.cv$ds))
+  # Each y in df.cv and DATA with same ds should be equal
+  df.merged <- dplyr::left_join(df.cv, m$history, by="ds")
+  expect_equal(sum((df.merged$y.x - df.merged$y.y) ** 2), 0)
+  df.cv <- cross_validation(
+    m, horizon = 4, units = "days", period = 10, initial = 135)
+  expect_equal(length(unique(df.cv$cutoff)), 1)
+  expect_error(
+    cross_validation(
+      m, horizon = 10, units = "days", period = 10, initial = 140)
+  )
 })
 
-test_that("simulated_historical_forecasts_logistic", {
+test_that("cross_validation_logistic", {
   skip_if_not(Sys.getenv('R_ARCH') != '/i386')
   df <- DATA
   df$cap <- 40
-  m <- prophet(df, growth='logistic')
-  df.shf <- simulated_historical_forecasts(
-    m, horizon = 3, units = 'days', k = 2, period = 3)
-  # All cutoff dates should be less than ds dates
-  expect_true(all(df.shf$cutoff < df.shf$ds))
-  # The unique size of output cutoff should be equal to 'k'
-  expect_equal(length(unique(df.shf$cutoff)), 2)
-  # Each y in df_shf and DATA with same ds should be equal
-  df.merged <- dplyr::left_join(df.shf, m$history, by="ds")
+  m <- prophet(df, growth = 'logistic')
+  df.cv <- cross_validation(
+    m, horizon = 1, units = "days", period = 1, initial = 140)
+  expect_equal(length(unique(df.cv$cutoff)), 2)
+  expect_true(all(df.cv$cutoff < df.cv$ds))
+  df.merged <- dplyr::left_join(df.cv, m$history, by="ds")
   expect_equal(sum((df.merged$y.x - df.merged$y.y) ** 2), 0)
 })
 
-test_that("simulated_historical_forecasts_extra_regressors", {
+test_that("cross_validation_extra_regressors", {
   skip_if_not(Sys.getenv('R_ARCH') != '/i386')
   df <- DATA
   df$extra <- seq(0, nrow(df) - 1)
@@ -55,43 +59,16 @@ test_that("simulated_historical_forecasts_extra_regressors", {
   m <- add_seasonality(m, name = 'monthly', period = 30.5, fourier.order = 5)
   m <- add_regressor(m, 'extra')
   m <- fit.prophet(m, df)
-  df.shf <- simulated_historical_forecasts(
-    m, horizon = 3, units = 'days', k = 2, period = 3)
-  # All cutoff dates should be less than ds dates
-  expect_true(all(df.shf$cutoff < df.shf$ds))
-  # The unique size of output cutoff should be equal to 'k'
-  expect_equal(length(unique(df.shf$cutoff)), 2)
-  # Each y in df_shf and DATA with same ds should be equal
-  df.merged <- dplyr::left_join(df.shf, m$history, by="ds")
-  expect_equal(sum((df.merged$y.x - df.merged$y.y) ** 2), 0)
-})
-
-test_that("simulated_historical_forecasts_default_value_check", {
-  skip_if_not(Sys.getenv('R_ARCH') != '/i386')
-  m <- prophet(DATA)
-  df.shf1 <- simulated_historical_forecasts(
-    m, horizon = 10, units = 'days', k = 1)
-  df.shf2 <- simulated_historical_forecasts(
-    m, horizon = 10, units = 'days', k = 1, period = 5)
-  expect_equal(sum(dplyr::select(df.shf1 - df.shf2, y, yhat)), 0)
-})
-
-test_that("cross_validation", {
-  skip_if_not(Sys.getenv('R_ARCH') != '/i386')
-  m <- prophet(DATA)
-  # Calculate the number of cutoff points
-  te <- max(DATA$ds)
-  ts <- min(DATA$ds)
-  horizon <- as.difftime(4, units = "days")
-  period <- as.difftime(10, units = "days")
-  k <- 5
   df.cv <- cross_validation(
-    m, horizon = 4, units = "days", period = 10, initial = 90)
-  expect_equal(length(unique(df.cv$cutoff)), k)
-  expect_equal(max(df.cv$ds - df.cv$cutoff), horizon)
+    m, horizon = 4, units = "days", period = 4, initial = 135)
+  expect_equal(length(unique(df.cv$cutoff)), 2)
+  period <- as.difftime(4, units = "days")
   dc <- diff(df.cv$cutoff)
   dc <- min(dc[dc > 0])
   expect_true(dc >= period)
+  expect_true(all(df.cv$cutoff < df.cv$ds))
+  df.merged <- dplyr::left_join(df.cv, m$history, by="ds")
+  expect_equal(sum((df.merged$y.x - df.merged$y.y) ** 2), 0)
 })
 
 test_that("cross_validation_default_value_check", {
@@ -103,3 +80,103 @@ test_that("cross_validation_default_value_check", {
     m, horizon = 32, units = 'days', period = 10, initial = 96)
   expect_equal(sum(dplyr::select(df.cv1 - df.cv2, y, yhat)), 0)
 })
+
+test_that("performance_metrics", {
+  skip_if_not(Sys.getenv('R_ARCH') != '/i386')
+  m <- prophet(DATA)
+  df_cv <- cross_validation(
+    m, horizon = 4, units = "days", period = 10, initial = 90)
+  # Aggregation level none
+  df_none <- performance_metrics(df_cv, rolling_window = 0)
+  expect_true(all(
+    sort(colnames(df_none))
+    == sort(c('horizon', 'coverage', 'mae', 'mape', 'mse', 'rmse'))
+  ))
+  expect_equal(nrow(df_none), 16)
+  # Aggregation level 0.2
+  df_horizon <- performance_metrics(df_cv, rolling_window = 0.2)
+  expect_equal(length(unique(df_horizon$horizon)), 4)
+  expect_equal(nrow(df_horizon), 14)
+  # Aggregation level all
+  df_all <- performance_metrics(df_cv, rolling_window = 1)
+  expect_equal(nrow(df_all), 1)
+  for (metric in c('mse', 'mape', 'mae', 'coverage')) {
+    expect_equal(df_all[[metric]][1], mean(df_none[[metric]]))
+  }
+  # Custom list of metrics
+  df_horizon <- performance_metrics(df_cv, metrics = c('coverage', 'mse'))
+  expect_true(all(
+    sort(colnames(df_horizon)) == sort(c('coverage', 'mse', 'horizon'))
+  ))
+})
+
+test_that("copy", {
+  skip_if_not(Sys.getenv('R_ARCH') != '/i386')
+  df <- DATA_all
+  df$cap <- 200.
+  df$binary_feature <- c(rep(0, 255), rep(1, 255))
+  inputs <- list(
+    growth = c('linear', 'logistic'),
+    yearly.seasonality = c(TRUE, FALSE),
+    weekly.seasonality = c(TRUE, FALSE),
+    daily.seasonality = c(TRUE, FALSE),
+    holidays = c('null', 'insert_dataframe'),
+    seasonality.mode = c('additive', 'multiplicative')
+  )
+  products <- expand.grid(inputs)
+  for (i in 1:length(products)) {
+    if (products$holidays[i] == 'insert_dataframe') {
+      holidays <- data.frame(ds=c('2016-12-25'), holiday=c('x'))
+    } else {
+      holidays <- NULL
+    }
+    m1 <- prophet(
+      growth = as.character(products$growth[i]),
+      changepoints = NULL,
+      n.changepoints = 3,
+      changepoint.range = 0.9,
+      yearly.seasonality = products$yearly.seasonality[i],
+      weekly.seasonality = products$weekly.seasonality[i],
+      daily.seasonality = products$daily.seasonality[i],
+      holidays = holidays,
+      seasonality.prior.scale = 1.1,
+      holidays.prior.scale = 1.1,
+      changepoints.prior.scale = 0.1,
+      mcmc.samples = 100,
+      interval.width = 0.9,
+      uncertainty.samples = 200,
+      fit = FALSE
+    )
+    out <- prophet:::setup_dataframe(m1, df, initialize_scales = TRUE)
+    m1 <- out$m
+    m1$history <- out$df
+    m1 <- prophet:::set_auto_seasonalities(m1)
+    m2 <- prophet:::prophet_copy(m1)
+    # Values should be copied correctly
+    args <- c('growth', 'changepoints', 'n.changepoints', 'holidays',
+              'seasonality.prior.scale', 'holidays.prior.scale',
+              'changepoints.prior.scale', 'mcmc.samples', 'interval.width',
+              'uncertainty.samples', 'seasonality.mode', 'changepoint.range')
+    for (arg in args) {
+      expect_equal(m1[[arg]], m2[[arg]])
+    }
+    expect_equal(FALSE, m2$yearly.seasonality)
+    expect_equal(FALSE, m2$weekly.seasonality)
+    expect_equal(FALSE, m2$daily.seasonality)
+    expect_equal(m1$yearly.seasonality, 'yearly' %in% names(m2$seasonalities))
+    expect_equal(m1$weekly.seasonality, 'weekly' %in% names(m2$seasonalities))
+    expect_equal(m1$daily.seasonality, 'daily' %in% names(m2$seasonalities))
+  }
+  # Check for cutoff and custom seasonality and extra regressors
+  changepoints <- seq.Date(as.Date('2012-06-15'), as.Date('2012-09-15'), by='d')
+  cutoff <- as.Date('2012-07-25')
+  m1 <- prophet(changepoints = changepoints)
+  m1 <- add_seasonality(m1, 'custom', 10, 5)
+  m1 <- add_regressor(m1, 'binary_feature')
+  m1 <- fit.prophet(m1, df)
+  m2 <- prophet:::prophet_copy(m1, cutoff)
+  changepoints <- changepoints[changepoints <= cutoff]
+  expect_equal(prophet:::set_date(changepoints), m2$changepoints)
+  expect_true('custom' %in% names(m2$seasonalities))
+  expect_true('binary_feature' %in% names(m2$extra_regressors))
+})

R/tests/testthat/test_prophet.R

@@ -127,11 +127,26 @@ test_that("get_changepoints", {
   cp <- m$changepoints.t
   expect_equal(length(cp), m$n.changepoints)
   expect_true(min(cp) > 0)
-  expect_true(max(cp) < N)
+  expect_true(max(cp) <= history$t[ceiling(0.8 * length(history$t))])
+})
+
+test_that("set_changepoint_range", {
+  history <- train
+  m <- prophet(history, fit = FALSE, changepoint.range = 0.4)
+
+  out <- prophet:::setup_dataframe(m, history, initialize_scales = TRUE)
+  history <- out$df
+  m <- out$m
+  m$history <- history
 
-  mat <- prophet:::get_changepoint_matrix(m)
-  expect_equal(nrow(mat), floor(N / 2))
-  expect_equal(ncol(mat), m$n.changepoints)
+  m <- prophet:::set_changepoints(m)
+
+  cp <- m$changepoints.t
+  expect_equal(length(cp), m$n.changepoints)
+  expect_true(min(cp) > 0)
+  expect_true(max(cp) <= history$t[ceiling(0.4 * length(history$t))])
+  expect_error(prophet(history, changepoint.range = -0.1))
+  expect_error(prophet(history, changepoint.range = 2))
 })
 
 test_that("get_zero_changepoints", {
@@ -147,10 +162,6 @@ test_that("get_zero_changepoints", {
   cp <- m$changepoints.t
   expect_equal(length(cp), 1)
   expect_equal(cp[1], 0)
-
-  mat <- prophet:::get_changepoint_matrix(m)
-  expect_equal(nrow(mat), floor(N / 2))
-  expect_equal(ncol(mat), 1)
 })
 
 test_that("override_n_changepoints", {
@@ -239,7 +250,7 @@ test_that("piecewise_logistic", {
 })
 
 test_that("holidays", {
-  holidays = data.frame(ds = c('2016-12-25'),
+  holidays <- data.frame(ds = c('2016-12-25'),
                         holiday = c('xmas'),
                         lower_window = c(-1),
                         upper_window = c(0))
@@ -250,12 +261,14 @@ test_that("holidays", {
   out <- prophet:::make_holiday_features(m, df$ds)
   feats <- out$holiday.features
   priors <- out$prior.scales
+  names <- out$holiday.names
   expect_equal(nrow(feats), nrow(df))
   expect_equal(ncol(feats), 2)
   expect_equal(sum(colSums(feats) - c(1, 1)), 0)
   expect_true(all(priors == c(10., 10.)))
+  expect_equal(names, c('xmas'))
 
-  holidays = data.frame(ds = c('2016-12-25'),
+  holidays <- data.frame(ds = c('2016-12-25'),
                         holiday = c('xmas'),
                         lower_window = c(-1),
                         upper_window = c(10))
@@ -263,9 +276,11 @@ test_that("holidays", {
   out <- prophet:::make_holiday_features(m, df$ds)
   feats <- out$holiday.features
   priors <- out$prior.scales
+  names <- out$holiday.names
   expect_equal(nrow(feats), nrow(df))
   expect_equal(ncol(feats), 12)
   expect_true(all(priors == rep(10, 12)))
+  expect_equal(names, c('xmas'))
   # Check prior specifications
   holidays <- data.frame(
     ds = prophet:::set_date(c('2016-12-25', '2017-12-25')),
@@ -277,7 +292,9 @@ test_that("holidays", {
   m <- prophet(holidays = holidays, fit = FALSE)
   out <- prophet:::make_holiday_features(m, df$ds)
   priors <- out$prior.scales
+  names <- out$holiday.names
   expect_true(all(priors == c(5., 5.)))
+  expect_equal(names, c('xmas'))
   # 2 different priors
   holidays2 <- data.frame(
     ds = prophet:::set_date(c('2012-06-06', '2013-06-06')),
@@ -290,7 +307,9 @@ test_that("holidays", {
   m <- prophet(holidays = holidays2, fit = FALSE)
   out <- prophet:::make_holiday_features(m, df$ds)
   priors <- out$prior.scales
+  names <- out$holiday.names
   expect_true(all(priors == c(8, 8, 5, 5)))
+  expect_true(all(sort(names) == c('seans-bday', 'xmas')))
   holidays2 <- data.frame(
     ds = prophet:::set_date(c('2012-06-06', '2013-06-06')),
     holiday = c('seans-bday', 'seans-bday'),
@@ -353,7 +372,8 @@ test_that("auto_weekly_seasonality", {
   expect_equal(m$weekly.seasonality, 'auto')
   m <- fit.prophet(m, train.w)
   expect_true('weekly' %in% names(m$seasonalities))
-  true <- list(period = 7, fourier.order = 3, prior.scale = 10)
+  true <- list(
+    period = 7, fourier.order = 3, prior.scale = 10, mode = 'additive')
   for (name in names(true)) {
     expect_equal(m$seasonalities$weekly[[name]], true[[name]])
   }
@@ -372,7 +392,8 @@ test_that("auto_weekly_seasonality", {
   m <- prophet(train.w)
   expect_false('weekly' %in% names(m$seasonalities))
   m <- prophet(DATA, weekly.seasonality = 2, seasonality.prior.scale = 3)
-  true <- list(period = 7, fourier.order = 2, prior.scale = 3)
+  true <- list(
+    period = 7, fourier.order = 2, prior.scale = 3, mode = 'additive')
   for (name in names(true)) {
     expect_equal(m$seasonalities$weekly[[name]], true[[name]])
   }
@@ -385,7 +406,8 @@ test_that("auto_yearly_seasonality", {
   expect_equal(m$yearly.seasonality, 'auto')
   m <- fit.prophet(m, DATA)
   expect_true('yearly' %in% names(m$seasonalities))
-  true <- list(period = 365.25, fourier.order = 10, prior.scale = 10)
+  true <- list(
+    period = 365.25, fourier.order = 10, prior.scale = 10, mode = 'additive')
   for (name in names(true)) {
     expect_equal(m$seasonalities$yearly[[name]], true[[name]])
   }
@@ -397,7 +419,8 @@ test_that("auto_yearly_seasonality", {
   m <- prophet(train.y, yearly.seasonality = TRUE)
   expect_true('yearly' %in% names(m$seasonalities))
   m <- prophet(DATA, yearly.seasonality = 7, seasonality.prior.scale = 3)
-  true <- list(period = 365.25, fourier.order = 7, prior.scale = 3)
+  true <- list(
+    period = 365.25, fourier.order = 7, prior.scale = 3, mode = 'additive')
   for (name in names(true)) {
     expect_equal(m$seasonalities$yearly[[name]], true[[name]])
   }
@@ -410,7 +433,8 @@ test_that("auto_daily_seasonality", {
   expect_equal(m$daily.seasonality, 'auto')
   m <- fit.prophet(m, DATA2)
   expect_true('daily' %in% names(m$seasonalities))
-  true <- list(period = 1, fourier.order = 4, prior.scale = 10)
+  true <- list(
+    period = 1, fourier.order = 4, prior.scale = 10, mode = 'additive')
   for (name in names(true)) {
     expect_equal(m$seasonalities$daily[[name]], true[[name]])
   }
@@ -422,7 +446,8 @@ test_that("auto_daily_seasonality", {
   m <- prophet(train.y, daily.seasonality = TRUE)
   expect_true('daily' %in% names(m$seasonalities))
   m <- prophet(DATA2, daily.seasonality = 7, seasonality.prior.scale = 3)
-  true <- list(period = 1, fourier.order = 7, prior.scale = 3)
+  true <- list(
+    period = 1, fourier.order = 7, prior.scale = 3, mode = 'additive')
   for (name in names(true)) {
     expect_equal(m$seasonalities$daily[[name]], true[[name]])
   }
@@ -446,7 +471,8 @@ test_that("custom_seasonality", {
                          prior_scale = c(4))
   m <- prophet(holidays=holidays)
   m <- add_seasonality(m, name='monthly', period=30, fourier.order=5)
-  true <- list(period = 30, fourier.order = 5, prior.scale = 10)
+  true <- list(
+    period = 30, fourier.order = 5, prior.scale = 10, mode = 'additive')
   for (name in names(true)) {
     expect_equal(m$seasonalities$monthly[[name]], true[[name]])
   }
@@ -458,12 +484,25 @@ test_that("custom_seasonality", {
   )
   m <- add_seasonality(m, name='weekly', period=30, fourier.order=5)
   # Test priors
-  m <- prophet(holidays = holidays, yearly.seasonality = FALSE)
+  m <- prophet(
+    holidays = holidays, yearly.seasonality = FALSE,
+    seasonality.mode = 'multiplicative')
   m <- add_seasonality(
-    m, name='monthly', period=30, fourier.order=5, prior.scale = 2)
+    m, name='monthly', period=30, fourier.order=5, prior.scale = 2,
+    mode = 'additive')
   m <- fit.prophet(m, DATA)
-  prior.scales <- prophet:::make_all_seasonality_features(
-    m, m$history)$prior.scales
+  expect_equal(m$seasonalities$monthly$mode, 'additive')
+  expect_equal(m$seasonalities$weekly$mode, 'multiplicative')
+  out <- prophet:::make_all_seasonality_features(m, m$history)
+  prior.scales <- out$prior.scales
+  component.cols <- out$component.cols
+  expect_equal(sum(component.cols$monthly), 10)
+  expect_equal(sum(component.cols$special_day), 1)
+  expect_equal(sum(component.cols$weekly), 6)
+  expect_equal(sum(component.cols$additive_terms), 10)
+  expect_equal(sum(component.cols$multiplicative_terms), 7)
+  expect_equal(sum(component.cols$monthly[1:11]), 10)
+  expect_equal(sum(component.cols$weekly[11:17]), 6)
   expect_true(all(prior.scales == c(rep(2, 10), rep(10, 6), 4)))
 })
 
@@ -472,10 +511,13 @@ test_that("added_regressors", {
   m <- prophet()
   m <- add_regressor(m, 'binary_feature', prior.scale=0.2)
   m <- add_regressor(m, 'numeric_feature', prior.scale=0.5)
+  m <- add_regressor(
+    m, 'numeric_feature2', prior.scale=0.5, mode = 'multiplicative')
   m <- add_regressor(m, 'binary_feature2', standardize=TRUE)
   df <- DATA
   df$binary_feature <- c(rep(0, 255), rep(1, 255))
   df$numeric_feature <- 0:509
+  df$numeric_feature2 <- 0:509
   # Require all regressors in df
   expect_error(
     fit.prophet(m, df)
@@ -483,7 +525,9 @@ test_that("added_regressors", {
   df$binary_feature2 <- c(rep(1, 100), rep(0, 410))
   m <- fit.prophet(m, df)
   # Check that standardizations are correctly set
-  true <- list(prior.scale = 0.2, mu = 0, std = 1, standardize = 'auto')
+  true <- list(
+    prior.scale = 0.2, mu = 0, std = 1, standardize = 'auto', mode = 'additive'
+  )
   for (name in names(true)) {
     expect_equal(true[[name]], m$extra_regressors$binary_feature[[name]])
   }
@@ -492,6 +536,7 @@ test_that("added_regressors", {
     expect_equal(true[[name]], m$extra_regressors$numeric_feature[[name]],
                  tolerance = 1e-5)
   }
+  expect_equal(m$extra_regressors$numeric_feature2$mode, 'multiplicative')
   true <- list(prior.scale = 10., mu = 0.1960784, std = 0.3974183)
   for (name in names(true)) {
     expect_equal(true[[name]], m$extra_regressors$binary_feature2[[name]],
@@ -506,97 +551,92 @@ test_that("added_regressors", {
   out <- prophet:::make_all_seasonality_features(m, df2)
   seasonal.features <- out$seasonal.features
   prior.scales <- out$prior.scales
-  expect_true('binary_feature' %in% colnames(seasonal.features))
-  expect_true('numeric_feature' %in% colnames(seasonal.features))
-  expect_true('binary_feature2' %in% colnames(seasonal.features))
-  expect_equal(ncol(seasonal.features), 29)
-  expect_true(all(sort(prior.scales[27:29]) == c(0.2, 0.5, 10.)))
+  component.cols <- out$component.cols
+  modes <- out$modes
+  expect_equal(ncol(seasonal.features), 30)
+  r_names <- c('binary_feature', 'numeric_feature', 'binary_feature2')
+  true.priors <- c(0.2, 0.5, 10.)
+  for (i in seq_along(r_names)) {
+    name <- r_names[i]
+    expect_true(name %in% colnames(seasonal.features))
+    expect_equal(sum(component.cols[[name]]), 1)
+    expect_equal(sum(prior.scales * component.cols[[name]]), true.priors[i])
+  }
   # Check that forecast components are reasonable
   future <- data.frame(
-    ds = c('2014-06-01'), binary_feature = c(0), numeric_feature = c(10))
+    ds = c('2014-06-01'),
+    binary_feature = c(0),
+    numeric_feature = c(10),
+    numeric_feature2 = c(10)
+  )
   expect_error(predict(m, future))
   future$binary_feature2 <- 0.
   fcst <- predict(m, future)
-  expect_equal(ncol(fcst), 31)
+  expect_equal(ncol(fcst), 37)
   expect_equal(fcst$binary_feature[1], 0)
-  expect_equal(fcst$extra_regressors[1],
+  expect_equal(fcst$extra_regressors_additive[1],
                fcst$numeric_feature[1] + fcst$binary_feature2[1])
-  expect_equal(fcst$seasonalities[1], fcst$yearly[1] + fcst$weekly[1])
-  expect_equal(fcst$seasonal[1],
-               fcst$seasonalities[1] + fcst$extra_regressors[1])
-  expect_equal(fcst$yhat[1], fcst$trend[1] + fcst$seasonal[1])
-  # Check fails if constant extra regressor
-  df$constant_feature <- 5
+  expect_equal(fcst$extra_regressors_multiplicative[1],
+               fcst$numeric_feature2[1])
+  expect_equal(fcst$additive_terms[1],
+               fcst$yearly[1] + fcst$weekly[1]
+               + fcst$extra_regressors_additive[1])
+  expect_equal(fcst$multiplicative_terms[1],
+               fcst$extra_regressors_multiplicative[1])
+  expect_equal(
+    fcst$yhat[1],
+    fcst$trend[1] * (1 + fcst$multiplicative_terms[1]) + fcst$additive_terms[1]
+  )
+  # Check works with constant extra regressor of 0
+  df$constant_feature <- 0
+  m <- prophet()
+  m <- add_regressor(m, 'constant_feature', standardize = TRUE)
+  m <- fit.prophet(m, df)
+  expect_equal(m$extra_regressors$constant_feature$std, 1)
+})
+
+test_that("set_seasonality_mode", {
+  skip_if_not(Sys.getenv('R_ARCH') != '/i386')
   m <- prophet()
-  m <- add_regressor(m, 'constant_feature')
-  expect_error(fit.prophet(m, df))
+  expect_equal(m$seasonality.mode, 'additive')
+  m <- prophet(seasonality.mode = 'multiplicative')
+  expect_equal(m$seasonality.mode, 'multiplicative')
+  expect_error(prophet(seasonality.mode = 'batman'))
 })
 
-test_that("copy", {
+test_that("seasonality_modes", {
   skip_if_not(Sys.getenv('R_ARCH') != '/i386')
+  holidays <- data.frame(ds = c('2016-12-25'),
+                        holiday = c('xmas'),
+                        lower_window = c(-1),
+                        upper_window = c(0))
+  m <- prophet(seasonality.mode = 'multiplicative', holidays = holidays)
+  m <- add_seasonality(
+    m, name = 'monthly', period = 30, fourier.order = 3, mode = 'additive')
+  m <- add_regressor(m, name = 'binary_feature', mode = 'additive')
+  m <- add_regressor(m, name = 'numeric_feature')
+  # Construct seasonal features
   df <- DATA
-  df$cap <- 200.
   df$binary_feature <- c(rep(0, 255), rep(1, 255))
-  inputs <- list(
-    growth = c('linear', 'logistic'),
-    yearly.seasonality = c(TRUE, FALSE),
-    weekly.seasonality = c(TRUE, FALSE),
-    daily.seasonality = c(TRUE, FALSE),
-    holidays = c('null', 'insert_dataframe')
+  df$numeric_feature <- 0:509
+  out <- prophet:::setup_dataframe(m, df, initialize_scales = TRUE)
+  df <- out$df
+  m <- out$m
+  m$history <- df
+  m <- prophet:::set_auto_seasonalities(m)
+  out <- prophet:::make_all_seasonality_features(m, df)
+  component.cols <- out$component.cols
+  modes <- out$modes
+  expect_equal(sum(component.cols$additive_terms), 7)
+  expect_equal(sum(component.cols$multiplicative_terms), 29)
+  expect_equal(
+    sort(modes$additive),
+    c('additive_terms', 'binary_feature', 'extra_regressors_additive',
+      'monthly')
+  )
+  expect_equal(
+    sort(modes$multiplicative),
+    c('extra_regressors_multiplicative', 'holidays', 'multiplicative_terms',
+      'numeric_feature', 'weekly', 'xmas', 'yearly')
   )
-  products <- expand.grid(inputs)
-  for (i in 1:length(products)) {
-    if (products$holidays[i] == 'insert_dataframe') {
-      holidays <- data.frame(ds=c('2016-12-25'), holiday=c('x'))
-    } else {
-      holidays <- NULL
-    }
-    m1 <- prophet(
-      growth = as.character(products$growth[i]),
-      changepoints = NULL,
-      n.changepoints = 3,
-      yearly.seasonality = products$yearly.seasonality[i],
-      weekly.seasonality = products$weekly.seasonality[i],
-      daily.seasonality = products$daily.seasonality[i],
-      holidays = holidays,
-      seasonality.prior.scale = 1.1,
-      holidays.prior.scale = 1.1,
-      changepoints.prior.scale = 0.1,
-      mcmc.samples = 100,
-      interval.width = 0.9,
-      uncertainty.samples = 200,
-      fit = FALSE
-    )
-    out <- prophet:::setup_dataframe(m1, df, initialize_scales = TRUE)
-    m1 <- out$m
-    m1$history <- out$df
-    m1 <- prophet:::set_auto_seasonalities(m1)
-    m2 <- prophet:::prophet_copy(m1)
-    # Values should be copied correctly
-    args <- c('growth', 'changepoints', 'n.changepoints', 'holidays',
-              'seasonality.prior.scale', 'holidays.prior.scale',
-              'changepoints.prior.scale', 'mcmc.samples', 'interval.width',
-              'uncertainty.samples')
-    for (arg in args) {
-      expect_equal(m1[[arg]], m2[[arg]])
-    }
-    expect_equal(FALSE, m2$yearly.seasonality)
-    expect_equal(FALSE, m2$weekly.seasonality)
-    expect_equal(FALSE, m2$daily.seasonality)
-    expect_equal(m1$yearly.seasonality, 'yearly' %in% names(m2$seasonalities))
-    expect_equal(m1$weekly.seasonality, 'weekly' %in% names(m2$seasonalities))
-    expect_equal(m1$daily.seasonality, 'daily' %in% names(m2$seasonalities))
-  }
-  # Check for cutoff and custom seasonality and extra regressors
-  changepoints <- seq.Date(as.Date('2012-06-15'), as.Date('2012-09-15'), by='d')
-  cutoff <- as.Date('2012-07-25')
-  m1 <- prophet(changepoints = changepoints)
-  m1 <- add_seasonality(m1, 'custom', 10, 5)
-  m1 <- add_regressor(m1, 'binary_feature')
-  m1 <- fit.prophet(m1, df)
-  m2 <- prophet:::prophet_copy(m1, cutoff)
-  changepoints <- changepoints[changepoints <= cutoff]
-  expect_equal(prophet:::set_date(changepoints), m2$changepoints)
-  expect_true('custom' %in% names(m2$seasonalities))
-  expect_true('binary_feature' %in% names(m2$extra_regressors))
 })

R/tests/testthat/test_stan_functions.R

@@ -0,0 +1,109 @@
+library(prophet)
+context("Prophet stan model tests")
+
+fn <- tryCatch({
+  rstan::expose_stan_functions(
+    rstan::stanc(file="../../inst/stan/prophet.stan")
+  )
+}, error = function(e) {
+  rstan::expose_stan_functions(
+    rstan::stanc(file=system.file("stan/prophet.stan", package="prophet"))
+  )
+})
+
+DATA <- read.csv('data.csv')
+N <- nrow(DATA)
+train <- DATA[1:floor(N / 2), ]
+future <- DATA[(ceiling(N/2) + 1):N, ]
+
+DATA2 <- read.csv('data2.csv')
+
+DATA$ds <- prophet:::set_date(DATA$ds)
+DATA2$ds <- prophet:::set_date(DATA2$ds)
+
+test_that("get_changepoint_matrix", {
+  skip_if_not(Sys.getenv('R_ARCH') != '/i386')
+  history <- train
+  m <- prophet(history, fit = FALSE)
+
+  out <- prophet:::setup_dataframe(m, history, initialize_scales = TRUE)
+  history <- out$df
+  m <- out$m
+  m$history <- history
+
+  m <- prophet:::set_changepoints(m)
+
+  cp <- m$changepoints.t
+
+  mat <- get_changepoint_matrix(history$t, cp, nrow(history), length(cp))
+  expect_equal(nrow(mat), floor(N / 2))
+  expect_equal(ncol(mat), m$n.changepoints)
+  # Compare to the R implementation
+  A <- matrix(0, nrow(history), length(cp))
+  for (i in 1:length(cp)) {
+    A[history$t >= cp[i], i] <- 1
+  }
+  expect_true(all(A == mat))
+})
+
+test_that("get_zero_changepoints", {
+  skip_if_not(Sys.getenv('R_ARCH') != '/i386')
+  history <- train
+  m <- prophet(history, n.changepoints = 0, fit = FALSE)
+  
+  out <- prophet:::setup_dataframe(m, history, initialize_scales = TRUE)
+  m <- out$m
+  history <- out$df
+  m$history <- history
+
+  m <- prophet:::set_changepoints(m)
+  cp <- m$changepoints.t
+  
+  mat <- get_changepoint_matrix(history$t, cp, nrow(history), length(cp))
+  expect_equal(nrow(mat), floor(N / 2))
+  expect_equal(ncol(mat), 1)
+  expect_true(all(mat == 1))
+})
+
+test_that("linear_trend", {
+  skip_if_not(Sys.getenv('R_ARCH') != '/i386')
+  t <- seq(0, 10)
+  m <- 0
+  k <- 1.0
+  deltas <- c(0.5)
+  changepoint.ts <- c(5)
+  A <- get_changepoint_matrix(t, changepoint.ts, length(t), 1)
+
+  y <- linear_trend(k, m, deltas, t, A, changepoint.ts)
+  y.true <- c(0, 1, 2, 3, 4, 5, 6.5, 8, 9.5, 11, 12.5)
+  expect_equal(y, y.true)
+
+  t <- t[8:length(t)]
+  A <- get_changepoint_matrix(t, changepoint.ts, length(t), 1)
+  y.true <- y.true[8:length(y.true)]
+  y <- linear_trend(k, m, deltas, t, A, changepoint.ts)
+  expect_equal(y, y.true)
+})
+
+test_that("piecewise_logistic", {
+  skip_if_not(Sys.getenv('R_ARCH') != '/i386')
+  t <- seq(0, 10)
+  cap <- rep(10, 11)
+  m <- 0
+  k <- 1.0
+  deltas <- c(0.5)
+  changepoint.ts <- c(5)
+  A <- get_changepoint_matrix(t, changepoint.ts, length(t), 1)
+
+  y <- logistic_trend(k, m, deltas, t, cap, A, changepoint.ts, 1)
+  y.true <- c(5.000000, 7.310586, 8.807971, 9.525741, 9.820138, 9.933071,
+              9.984988, 9.996646, 9.999252, 9.999833, 9.999963)
+  expect_equal(y, y.true, tolerance = 1e-6)
+  
+  t <- t[8:length(t)]
+  A <- get_changepoint_matrix(t, changepoint.ts, length(t), 1)
+  y.true <- y.true[8:length(y.true)]
+  cap <- cap[8:length(cap)]
+  y <- logistic_trend(k, m, deltas, t, cap, A, changepoint.ts, 1)
+  expect_equal(y, y.true, tolerance = 1e-6)
+})

R/vignettes/quick_start.Rmd

@@ -16,6 +16,7 @@ library(prophet)
 library(dplyr)
 ```
 
+This document provides a very brief introduction to the Prophet API. For a detailed guide on using Prophet, please visit the main site at [https://facebook.github.io/prophet/](https://facebook.github.io/prophet/).
 
 Prophet uses the normal model fitting API.  We provide a `prophet` function that performs fitting and returns a model object.  You can then call `predict` and `plot` on this model object.
 
@@ -50,7 +51,7 @@ You can use the generic `plot` function to plot the forecast, but you must also
 plot(m, forecast)
 ```
 
-Just as in Python, you can plot the components of the forecast.  In R, you use the `prophet_plot_components` function instead of an instance method:
+You can plot the components of the forecast using the `prophet_plot_components` function:
 
 ```{r}
 prophet_plot_components(m, forecast)

README.md

@@ -1,6 +1,8 @@
 # Prophet: Automatic Forecasting Procedure
 
-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.
+[![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 based on an additive model where non-linear trends are fit with yearly, weekly, and daily seasonality, plus holiday effects. It works best with time series that have strong seasonal effects and several seasons of historical data. Prophet is robust to missing data and shifts in the trend, and typically handles outliers well.
 
 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/).
 
@@ -14,7 +16,7 @@ Prophet is [open source software](https://code.facebook.com/projects/) released
 - Prophet R package: https://cran.r-project.org/package=prophet
 - Prophet Python package: https://pypi.python.org/pypi/fbprophet/
 - Release blogpost: https://research.fb.com/prophet-forecasting-at-scale/
-- Prophet paper, "Forecasting at Scale": https://peerj.com/preprints/3190.pdf
+- Prophet paper: Sean J. Taylor, Benjamin Letham (2018) Forecasting at scale. The American Statistician 72(1):37-45 (https://peerj.com/preprints/3190.pdf).
 
 ## Installation in R
 
@@ -42,17 +44,19 @@ Prophet is on PyPI, so you can use pip to install it:
 $ pip install fbprophet
 ```
 
-The major dependency that Prophet has is `pystan`.   PyStan has its own [installation instructions](http://pystan.readthedocs.io/en/latest/installation_beginner.html).
+The major dependency that Prophet has is `pystan`.   PyStan has its own [installation instructions](http://pystan.readthedocs.io/en/latest/installation_beginner.html). Install pystan with pip before using pip to install fbprophet.
 
 After installation, you can [get started!](https://facebook.github.io/prophet/docs/quick_start.html#python-api)
 
+If you upgrade the version of PyStan installed on your system, you may need to reinstall fbprophet ([see here](https://github.com/facebook/prophet/issues/324)).
+
 ### Windows
 
 On Windows, PyStan requires a compiler so you'll need to [follow the instructions](http://pystan.readthedocs.io/en/latest/windows.html).  The key step is installing a recent [C++ compiler](http://landinghub.visualstudio.com/visual-cpp-build-tools).
 
 ### Linux
 
-Make sure compilers (gcc, g++, build-essential) and Python development tools (python-dev, python3-dev) are installed. If you are using a VM, be aware that you will need at least 4GB of memory to install fbprophet, and at least 2GB of memory to use fbprophet.
+Make sure compilers (gcc, g++, build-essential) and Python development tools (python-dev, python3-dev) are installed. In Red Hat systems, install the packages gcc64 and gcc64-c++. If you are using a VM, be aware that you will need at least 4GB of memory to install fbprophet, and at least 2GB of memory to use fbprophet.
 
 ### Anaconda
 
@@ -60,6 +64,15 @@ Use `conda install gcc` to set up gcc. The easiest way to install Prophet is thr
 
 ## Changelog
 
+### Version 0.3 (2018.06.01)
+
+- Multiplicative seasonality
+- Cross validation error metrics and visualizations
+- Parameter to set range of potential changepoints
+- Unified Stan model for both trend types
+- Improved future trend uncertainty for sub-daily data
+- Bugfixes
+
 ### Version 0.2.1 (2017.11.08)
 
 - Bugfixes

docs/_data/nav_docs.yml

@@ -4,7 +4,8 @@
   - id: quick_start
   - id: saturating_forecasts
   - id: trend_changepoints
-  - id: seasonality_and_holiday_effects
+  - id: seasonality,_holiday_effects,_and_regressors
+  - id: multiplicative_seasonality
   - id: uncertainty_intervals
   - id: outliers
   - id: non-daily_data

docs/_docs/contributing.md

@@ -1,17 +1,17 @@
 ---
 layout: docs
 docid: "contributing"
-title: "How to Contribute"
+title: "Getting Help and Contributing"
 permalink: /docs/contributing.html
 ---
 
-Prophet has an non-fixed release cycle but we will be making bugfixes in response to user feedback and adding features.  Its current state is Beta (v0.2), we expect no obvious bugs. Please let us know if you encounter a bug by [filing an issue](https://github.com/facebook/prophet/issues).
+Prophet has a non-fixed release cycle but we will be making bugfixes in response to user feedback and adding features.  Its current state is Beta (v0.3), we expect no obvious bugs. Please let us know if you encounter a bug by [filing an issue](https://github.com/facebook/prophet/issues). Github issues is also the right place to ask questions about using Prophet.
 
 We appreciate all contributions. If you are planning to contribute back bug-fixes, please do so without any further discussion.
 
-If you plan to contribute new features or extensions to the core, please first open an issue and discuss the feature with us. Sending a pull request is fine, too, but if it is a large change we suggest you run it by us first.
+If you plan to contribute new features or extensions to the core, please first open an issue and discuss the feature with us. Sending a pull request is fine too, but it will likely be merged more quickly if any design decisions are settled on beforehand in an issue.
 
-We require that any API changes or feature additions are made available for both Python and R in parallel.
+The R and Python versions are kept feature identical, but new features can be implemented for each method in separate commits.
 
 ## Documentation
 

docs/_docs/installation.md

@@ -33,7 +33,7 @@ Prophet is on PyPI, so you can use pip to install it:
 $ pip install fbprophet
 ```
 
-The major dependency that Prophet has is `pystan`.   PyStan has its own [installation instructions](http://pystan.readthedocs.io/en/latest/installation_beginner.html).
+The major dependency that Prophet has is `pystan`.   PyStan has its own [installation instructions](http://pystan.readthedocs.io/en/latest/installation_beginner.html). Install pystan with pip before using pip to install fbprophet.
 
 After installation, you can [get started!](quick_start.html#python-api)
 
@@ -43,7 +43,7 @@ On Windows, PyStan requires a compiler so you'll need to [follow the instruction
 
 ### Linux
 
-Make sure compilers (gcc, g++) and Python development tools (python-dev) are installed. If you are using a VM, be aware that you will need at least 4GB of memory to install fbprophet, and at least 2GB of memory to use fbprophet.
+Make sure compilers (gcc, g++, build-essential) and Python development tools (python-dev, python3-dev) are installed. In Red Hat systems, install the packages gcc64 and gcc64-c++. If you are using a VM, be aware that you will need at least 4GB of memory to install fbprophet, and at least 2GB of memory to use fbprophet.
 
 ### Anaconda
 

docs/_docs/multiplicative_seasonality.md

@@ -0,0 +1,81 @@
+---
+layout: docs
+docid: "multiplicative_seasonality"
+title: "Multiplicative Seasonality"
+permalink: /docs/multiplicative_seasonality.html
+---
+By default Prophet fits additive seasonalities, meaning the effect of the seasonality is added to the trend to get the forecast. This time series of the number of air passengers is an example of when additive seasonality does not work:
+
+```R
+# R
+df <- read.csv('../examples/example_air_passengers.csv')
+m <- prophet(df)
+future <- make_future_dataframe(m, 50, freq = 'm')
+forecast <- predict(m, future)
+plot(m, forecast)
+```
+```python
+# Python
+df = pd.read_csv('../examples/example_air_passengers.csv')
+m = Prophet()
+m.fit(df)
+future = m.make_future_dataframe(50, freq='MS')
+forecast = m.predict(future)
+fig = m.plot(forecast)
+```
+ 
+![png](/prophet/static/multiplicative_seasonality_files/multiplicative_seasonality_4_0.png) 
+
+
+This time series has a clear yearly cycle, but the seasonality in the forecast is too large at the start of the time series and too small at the end. In this time series, the seasonality is not a constant additive factor as assumed by Prophet, rather it grows with the trend. This is multiplicative seasonality.
+
+Prophet can model multiplicative seasonality by setting `seasonality_mode='multiplicative'` in the input arguments:
+
+```R
+# R
+m <- prophet(df, seasonality.mode = 'multiplicative')
+forecast <- predict(m, future)
+plot(m, forecast)
+```
+```python
+# Python
+m = Prophet(seasonality_mode='multiplicative')
+m.fit(df)
+forecast = m.predict(future)
+fig = m.plot(forecast)
+```
+ 
+![png](/prophet/static/multiplicative_seasonality_files/multiplicative_seasonality_7_0.png) 
+
+
+The components figure will now show the seasonality as a percent of the trend:
+
+```R
+# R
+prophet_plot_components(m, forecast)
+```
+```python
+# Python
+fig = m.plot_components(forecast)
+```
+ 
+![png](/prophet/static/multiplicative_seasonality_files/multiplicative_seasonality_10_0.png) 
+
+
+With `seasonality_mode='multiplicative'`, holiday effects will also be modeled as multiplicative. Any added seasonalities or extra regressors will by default use whatever `seasonality_mode` is set to, but can be overriden by specifying `mode='additive'` or `mode='multiplicative'` as an argument when adding the seasonality or regressor.
+
+For example, this block sets the built-in seasonalities to multiplicative, but includes an additive quarterly seasonality and an additive regressor:
+
+```R
+# R
+m <- prophet(seasonality.mode = 'multiplicative')
+m <- add_seasonality(m, 'quarterly', period = 91.25, fourier.order = 8, mode = 'additive')
+m <- add_regressor(m, 'regressor', mode = 'additive')
+```
+```python
+# Python
+m = Prophet(seasonality_mode='multiplicative')
+m.add_seasonality('quarterly', period=91.25, fourier_order=8, mode='additive')
+m.add_regressor('regressor', mode='additive')
+```
+Additive and multiplicative extra regressors will show up in separate panels on the components plot.

docs/_docs/outliers.md

@@ -8,22 +8,20 @@ There are two main ways that outliers can affect Prophet forecasts. Here we make
 
 ```R
 # R
-df <- read.csv('../examples/example_wp_R_outliers1.csv')
-df$y <- log(df$y)
+df <- read.csv('../examples/example_wp_log_R_outliers1.csv')
 m <- prophet(df)
 future <- make_future_dataframe(m, periods = 1096)
 forecast <- predict(m, future)
-plot(m, forecast);
+plot(m, forecast)
 ```
 ```python
 # Python
-df = pd.read_csv('../examples/example_wp_R_outliers1.csv')
-df['y'] = np.log(df['y'])
+df = pd.read_csv('../examples/example_wp_log_R_outliers1.csv')
 m = Prophet()
 m.fit(df)
 future = m.make_future_dataframe(periods=1096)
 forecast = m.predict(future)
-m.plot(forecast);
+fig = m.plot(forecast)
 ```
  
 ![png](/prophet/static/outliers_files/outliers_4_0.png) 
@@ -40,13 +38,13 @@ outliers <- (as.Date(df$ds) > as.Date('2010-01-01')
 df$y[outliers] = NA
 m <- prophet(df)
 forecast <- predict(m, future)
-plot(m, forecast);
+plot(m, forecast)
 ```
 ```python
 # Python
 df.loc[(df['ds'] > '2010-01-01') & (df['ds'] < '2011-01-01'), 'y'] = None
 model = Prophet().fit(df)
-model.plot(model.predict(future));
+fig = model.plot(model.predict(future))
 ```
  
 ![png](/prophet/static/outliers_files/outliers_7_0.png) 
@@ -56,22 +54,20 @@ In the above example the outliers messed up the uncertainty estimation but did n
 
 ```R
 # R
-df <- read.csv('../examples/example_wp_R_outliers2.csv')
-df$y = log(df$y)
+df <- read.csv('../examples/example_wp_log_R_outliers2.csv')
 m <- prophet(df)
 future <- make_future_dataframe(m, periods = 1096)
 forecast <- predict(m, future)
-plot(m, forecast);
+plot(m, forecast)
 ```
 ```python
 # Python
-df = pd.read_csv('../examples/example_wp_R_outliers2.csv')
-df['y'] = np.log(df['y'])
+df = pd.read_csv('../examples/example_wp_log_R_outliers2.csv')
 m = Prophet()
 m.fit(df)
 future = m.make_future_dataframe(periods=1096)
 forecast = m.predict(future)
-m.plot(forecast);
+fig = m.plot(forecast)
 ```
  
 ![png](/prophet/static/outliers_files/outliers_10_0.png) 
@@ -86,13 +82,13 @@ outliers <- (as.Date(df$ds) > as.Date('2015-06-01')
 df$y[outliers] = NA
 m <- prophet(df)
 forecast <- predict(m, future)
-plot(m, forecast);
+plot(m, forecast)
 ```
 ```python
 # Python
 df.loc[(df['ds'] > '2015-06-01') & (df['ds'] < '2015-06-30'), 'y'] = None
 m = Prophet().fit(df)
-m.plot(m.predict(future));
+fig = m.plot(m.predict(future))
 ```
  
 ![png](/prophet/static/outliers_files/outliers_13_0.png) 

docs/_docs/saturating_forecasts.md

@@ -10,40 +10,37 @@ By default, Prophet uses a linear model for its forecast. When forecasting growt
 
 Prophet allows you to make forecasts using a [logistic growth](https://en.wikipedia.org/wiki/Logistic_function) trend model, with a specified carrying capacity. We illustrate this with the log number of page visits to the [R (programming language)](https://en.wikipedia.org/wiki/R_%28programming_language%29) page on Wikipedia:
 
-```python
-# Python
-df = pd.read_csv('../examples/example_wp_R.csv')
-import numpy as np
-df['y'] = np.log(df['y'])
-```
 ```R
 # R
-df <- read.csv('../examples/example_wp_R.csv')
-df$y <- log(df$y)
+df <- read.csv('../examples/example_wp_log_R.csv')
 ```
-We must specify the carrying capacity in a column `cap`. Here we will assume a particular value, but this would usually be set using data or expertise about the market size.
-
 ```python
 # Python
-df['cap'] = 8.5
+df = pd.read_csv('../examples/example_wp_log_R.csv')
 ```
+We must specify the carrying capacity in a column `cap`. Here we will assume a particular value, but this would usually be set using data or expertise about the market size.
+
 ```R
 # R
 df$cap <- 8.5
 ```
+```python
+# Python
+df['cap'] = 8.5
+```
 The important things to note are that `cap` must be specified for every row in the dataframe, and that it does not have to be constant. If the market size is growing, then `cap` can be an increasing sequence.
 
 We then fit the model as before, except pass in an additional argument to specify logistic growth:
 
+```R
+# R
+m <- prophet(df, growth = 'logistic')
+```
 ```python
 # Python
 m = Prophet(growth='logistic')
 m.fit(df)
 ```
-```R
-# R
-m <- prophet(df, growth = 'logistic')
-```
 We make a dataframe for future predictions as before, except we must also specify the capacity in the future. Here we keep capacity constant at the same value as in the history, and forecast 3 years into the future:
 
 ```R
@@ -51,14 +48,14 @@ We make a dataframe for future predictions as before, except we must also specif
 future <- make_future_dataframe(m, periods = 1826)
 future$cap <- 8.5
 fcst <- predict(m, future)
-plot(m, fcst);
+plot(m, fcst)
 ```
 ```python
 # Python
 future = m.make_future_dataframe(periods=1826)
 future['cap'] = 8.5
 fcst = m.predict(future)
-m.plot(fcst);
+fig = m.plot(fcst)
 ```
  
 ![png](/prophet/static/saturating_forecasts_files/saturating_forecasts_13_0.png) 
@@ -91,7 +88,7 @@ future['floor'] = 1.5
 m = Prophet(growth='logistic')
 m.fit(df)
 fcst = m.predict(future)
-m.plot(fcst);
+fig = m.plot(fcst)
 ```
  
 ![png](/prophet/static/saturating_forecasts_files/saturating_forecasts_16_0.png) 

docs/_docs/uncertainty_intervals.md

@@ -15,39 +15,39 @@ One property of this way of measuring uncertainty is that allowing higher flexib
 
 The width of the uncertainty intervals (by default 80%) can be set using the parameter `interval_width`:
 
-```python
-# Python
-forecast = Prophet(interval_width=0.95).fit(df).predict(future)
-```
 ```R
 # R
 m <- prophet(df, interval.width = 0.95)
 forecast <- predict(m, future)
 ```
+```python
+# Python
+forecast = Prophet(interval_width=0.95).fit(df).predict(future)
+```
 Again, these intervals assume that the future will see the same frequency and magnitude of rate changes as the past. This assumption is probably not true, so you should not expect to get accurate coverage on these uncertainty intervals.
 
 ### Uncertainty in seasonality
-By default Prophet will only return uncertainty in the trend and observation noise. To get uncertainty in seasonality, you must do full Bayesian sampling. This is done using the parameter `mcmc.samples` (which defaults to 0). We do this here for the Peyton Manning data from the Quickstart:
+By default Prophet will only return uncertainty in the trend and observation noise. To get uncertainty in seasonality, you must do full Bayesian sampling. This is done using the parameter `mcmc.samples` (which defaults to 0). We do this here for the first six months of the Peyton Manning data from the Quickstart:
 
-```python
-# Python
-m = Prophet(mcmc_samples=300)
-forecast = m.fit(df).predict(future)
-```
 ```R
 # R
 m <- prophet(df, mcmc.samples = 300)
 forecast <- predict(m, future)
 ```
-This replaces the typical MAP estimation with MCMC sampling, and takes much longer - think 10 minutes instead of 10 seconds. If you do full sampling, then you will see the uncertainty in seasonal components when you plot them:
-
 ```python
 # Python
-m.plot_components(forecast);
+m = Prophet(mcmc_samples=300)
+forecast = m.fit(df).predict(future)
 ```
+This replaces the typical MAP estimation with MCMC sampling, and can take much longer depending on how many observations there are - expect several minutes instead of several seconds. If you do full sampling, then you will see the uncertainty in seasonal components when you plot them:
+
 ```R
 # R
-prophet_plot_components(m, forecast);
+prophet_plot_components(m, forecast)
+```
+```python
+# Python
+fig = m.plot_components(forecast)
 ```
  
 ![png](/prophet/static/uncertainty_intervals_files/uncertainty_intervals_10_0.png) 

docs/index.md

@@ -4,7 +4,7 @@ title: Prophet
 id: home
 ---
 
-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 a procedure for forecasting time series data based on an additive model where non-linear trends are fit with yearly, weekly, and daily seasonality, plus holiday effects. It works best with time series that have strong seasonal effects and several seasons of historical data. Prophet is robust to missing data and shifts in the trend, and typically handles outliers well.
 
 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/).