Allow changepoints on dates that aren't in history, and allow for repeated observations on days. Previously we worked with changepoints via their index in the history. Now we work with them using just their value in scaled time.

R/R/prophet.R

@@ -101,8 +101,9 @@ prophet <- function(df = df,
     interval.width = interval.width,
     uncertainty.samples = uncertainty.samples,
     start = NULL,  # This and following attributes are set during fitting
-    end = NULL,
     y.scale = NULL,
+    t.scale = NULL,
+    changepoints.t = NULL,
     stan.fit = NULL,
     params = list(),
     history = NULL
@@ -206,12 +207,10 @@ setup_dataframe <- function(m, df, initialize_scales = FALSE) {
   if (initialize_scales) {
     m$y.scale <- max(df$y)
     m$start <- min(df$ds)
-    m$end <- max(df$ds)
+    m$t.scale <- as.numeric(max(df$ds) - m$start)
   }
 
-  t.scale <- as.numeric(m$end - m$start)
-
-  df$t <- as.numeric(df$ds - m$start) / t.scale
+  df$t <- as.numeric(df$ds - m$start) / m$t.scale
   if (exists('y', where=df)) {
     df$y_scaled <- df$y / m$y.scale
   }
@@ -254,32 +253,13 @@ set_changepoints <- function(m) {
       m$changepoints <- c()
     }
   }
-  return(m)
-}
-
-#' Gets changepoint indexes in history dataframe.
-#'
-#' @param m Prophet object.
-#'
-#' @return array of indexes.
-#'
-get_changepoint_indexes <- function(m) {
-  if (length(m$changepoints) == 0) {
-    return(c(1))
+  if (length(m$changepoints) > 0) {
+    m$changepoints <- zoo::as.Date(m$changepoints)
+    m$changepoints.t <- sort(as.numeric(m$changepoints - m$start) / m$t.scale)
   } else {
-    return(match(zoo::as.Date(m$changepoints), m$history$ds))
+    m$changepoints.t <- c(0)  # dummy changepoint
   }
-}
-
-#' Gets changepoint times, in scaled space.
-#'
-#' @param m Prophet object.
-#'
-#' @return array of times.
-#'
-get_changepoint_times <- function(m) {
-  cpi <- get_changepoint_indexes(m)
-  return(m$history$t[cpi])
+  return(m)
 }
 
 #' Gets changepoint matrix for history dataframe.
@@ -289,11 +269,9 @@ get_changepoint_times <- function(m) {
 #' @return array of indexes.
 #'
 get_changepoint_matrix <- function(m) {
-  changepoint.indexes <- get_changepoint_indexes(m)
-
-  A <- matrix(0, nrow(m$history), length(changepoint.indexes))
-  for (i in 1:length(changepoint.indexes)) {
-    A[changepoint.indexes[i]:nrow(m$history), i] <- 1
+  A <- matrix(0, nrow(m$history), length(m$changepoints.t))
+  for (i in 1:length(m$changepoints.t)) {
+    A[m$history$t >= m$changepoints.t[i], i] <- 1
   }
   return(A)
 }
@@ -470,17 +448,16 @@ fit.prophet <- function(m, df, ...) {
 
   m <- set_changepoints(m)
   A <- get_changepoint_matrix(m)
-  changepoint.indexes <- get_changepoint_indexes(m)
 
   # Construct input to stan
   dat <- list(
     T = nrow(history),
     K = ncol(seasonal.features),
-    S = length(changepoint.indexes),
+    S = length(m$changepoints.t),
     y = history$y_scaled,
     t = history$t,
     A = A,
-    s_indx = array(changepoint.indexes),
+    t_change = array(m$changepoints.t),
     X = as.matrix(seasonal.features),
     sigma = m$seasonality.prior.scale,
     tau = m$changepoint.prior.scale
@@ -499,7 +476,7 @@ fit.prophet <- function(m, df, ...) {
   stan_init <- function() {
     list(k = kinit[1],
          m = kinit[2],
-         delta = array(rep(0, length(changepoint.indexes))),
+         delta = array(rep(0, length(m$changepoints.t))),
          beta = array(rep(0, ncol(seasonal.features))),
          sigma_obs = 1
     )
@@ -649,12 +626,12 @@ predict_trend <- function(model, df) {
   deltas <- colMeans(model$params$delta, na.rm = TRUE)
 
   t <- df$t
-  cpts <- get_changepoint_times(model)
   if (model$growth == 'linear') {
-    trend <- piecewise_linear(t, deltas, k, param.m, cpts)
+    trend <- piecewise_linear(t, deltas, k, param.m, model$changepoints.t)
   } else {
     cap <- df$cap_scaled
-    trend <- piecewise_logistic(t, cap, deltas, k, param.m, cpts)
+    trend <- piecewise_logistic(
+      t, cap, deltas, k, param.m, model$changepoints.t)
   }
   return(trend * model$y.scale)
 }
@@ -785,7 +762,6 @@ sample_predictive_trend <- function(model, df, iteration) {
   deltas <- model$params$delta[iteration,]
 
   t <- df$t
-  changepoint.ts <- get_changepoint_times(model)
   T <- max(t)
 
   if (T > 1) {
@@ -794,7 +770,7 @@ sample_predictive_trend <- function(model, df, iteration) {
     dt <- min(dt[dt > 0])
     # Number of time periods in the future
     N <- ceiling((T - 1) / dt)
-    S <- length(changepoint.ts)
+    S <- length(model$changepoints.t)
     # The history had S split points, over t = [0, 1].
     # The forecast is on [1, T], and should have the same average frequency of
     # rate changes. Thus for N time periods in the future, we want an average
@@ -820,7 +796,7 @@ sample_predictive_trend <- function(model, df, iteration) {
   deltas.new <- extraDistr::rlaplace(n.changes, mu = 0, sigma = lambda)
 
   # Combine with changepoints from the history
-  changepoint.ts <- c(changepoint.ts, changepoint.ts.new)
+  changepoint.ts <- c(model$changepoints.t, changepoint.ts.new)
   deltas <- c(deltas, deltas.new)
 
   # Get the corresponding trend

R/inst/stan/prophet_linear_growth.stan

@@ -3,9 +3,9 @@ data {
   int<lower=1> K;                       // Number of seasonal vectors
   vector[T] t;                            // Day
   vector[T] y;                            // Time-series
-  int S;                                // Number of split points
+  int S;                                // Number of changepoints
   matrix[T, S] A;                   // Split indicators
-  int s_indx[S];                 // Index of split points
+  real t_change[S];                 // Index of changepoints
   matrix[T,K] X;                // season vectors
   real<lower=0> sigma;              // scale on seasonality prior
   real<lower=0> tau;                  // scale on changepoints prior
@@ -23,7 +23,7 @@ transformed parameters {
   vector[S] gamma;                  // adjusted offsets, for piecewise continuity
 
   for (i in 1:S) {
-    gamma[i] = -t[s_indx[i]] * delta[i];
+    gamma[i] = -t_change[i] * delta[i];
   }
 }
 

R/inst/stan/prophet_logistic_growth.stan

@@ -4,24 +4,14 @@ data {
   vector[T] t;                            // Day
   vector[T] cap;                          // Capacities
   vector[T] y;                            // Time-series
-  int S;                                // Number of split points
+  int S;                                // Number of changepoints
   matrix[T, S] A;                   // Split indicators
-  int s_indx[S];                 // Index of split points
+  real t_change[S];                 // Index of changepoints
   matrix[T,K] X;                    // season vectors
   real<lower=0> sigma;              // scale on seasonality prior
   real<lower=0> tau;                  // scale on changepoints prior
 }
 
-
-transformed data {
-  int s_ext[S + 1];  // Segment endpoints
-  for (j in 1:S) {
-    s_ext[j] = s_indx[j];
-  }
-  s_ext[S + 1] = T + 1;  // Used for the m_adj loop below.
-}
-
-
 parameters {
   real k;                            // Base growth rate
   real m;                            // offset
@@ -30,7 +20,6 @@ parameters {
   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
@@ -45,7 +34,7 @@ transformed parameters {
   // Piecewise offsets
   m_pr = m; // The offset in the previous segment
   for (i in 1:S) {
-    gamma[i] = (t[s_indx[i]] - m_pr) * (1 - k_s[i] / k_s[i + 1]);
+    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
   }
 }

R/man/get_changepoint_indexes.Rd

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

R/man/get_changepoint_times.Rd

@@ -1,18 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/prophet.R
-\name{get_changepoint_times}
-\alias{get_changepoint_times}
-\title{Gets changepoint times, in scaled space.}
-\usage{
-get_changepoint_times(m)
-}
-\arguments{
-\item{m}{Prophet object.}
-}
-\value{
-array of times.
-}
-\description{
-Gets changepoint times, in scaled space.
-}
-

R/tests/testthat/test_prophet.R

@@ -30,6 +30,25 @@ test_that("fit_predict_no_changepoints", {
   expect_error(predict(m, future), NA)
 })
 
+test_that("fit_predict_changepoint_not_in_history", {
+  skip_if_not(Sys.getenv('R_ARCH') != '/i386')
+  train_t <- dplyr::mutate(DATA, ds=zoo::as.Date(ds))
+  train_t <- dplyr::filter(train_t, (ds < zoo::as.Date('2013-01-01')) | 
+                                (ds > zoo::as.Date('2014-01-01')))
+  future <- data.frame(ds=DATA$ds)
+  m <- prophet(train_t, changepoints=c('2013-06-06'))
+  expect_error(predict(m, future), NA)
+})
+
+test_that("fit_predict_duplicates", {
+  skip_if_not(Sys.getenv('R_ARCH') != '/i386')
+  train2 <- train
+  train2$y <- train2$y + 10
+  train_t <- rbind(train, train2)
+  m <- prophet(train_t)
+  expect_error(predict(m, future), NA)
+})
+
 test_that("setup_dataframe", {
   history <- train
   m <- prophet(history, fit = FALSE)
@@ -56,7 +75,7 @@ test_that("get_changepoints", {
 
   m <- prophet:::set_changepoints(m)
 
-  cp <- prophet:::get_changepoint_indexes(m)
+  cp <- m$changepoints.t
   expect_equal(length(cp), m$n.changepoints)
   expect_true(min(cp) > 0)
   expect_true(max(cp) < N)
@@ -76,9 +95,9 @@ test_that("get_zero_changepoints", {
   m$history <- history
 
   m <- prophet:::set_changepoints(m)
-  cp <- prophet:::get_changepoint_indexes(m)
+  cp <- m$changepoints.t
   expect_equal(length(cp), 1)
-  expect_equal(cp[1], 1)
+  expect_equal(cp[1], 0)
 
   mat <- prophet:::get_changepoint_matrix(m)
   expect_equal(nrow(mat), floor(N / 2))

python/fbprophet/forecaster.py

@@ -84,8 +84,9 @@ class Prophet(object):
 
         # Set during fitting
         self.start = None
-        self.end = None
         self.y_scale = None
+        self.t_scale = None
+        self.changepoints_t = None
         self.stan_fit = None
         self.params = {}
         self.history = None
@@ -124,14 +125,14 @@ class Prophet(object):
         df['ds'] = pd.to_datetime(df['ds'])
 
         df = df.sort_values('ds')
+        df.reset_index(inplace=True, drop=True)
 
         if initialize_scales:
             self.y_scale = df['y'].max()
-            self.start, self.end = df['ds'].min(), df['ds'].max()
+            self.start = df['ds'].min()
+            self.t_scale = df['ds'].max() - self.start
 
-        t_scale = self.end - self.start
-
-        df['t'] = (df['ds'] - self.start) / t_scale
+        df['t'] = (df['ds'] - self.start) / self.t_scale
         if 'y' in df:
             df['y_scaled'] = df['y'] / self.y_scale
 
@@ -171,31 +172,17 @@ class Prophet(object):
         else:
             # set empty changepoints
             self.changepoints = []
-
-    def get_changepoint_indexes(self):
-        if len(self.changepoints) == 0:
-            return np.array([0])  # a dummy changepoint
+        if len(self.changepoints) > 0:
+            self.changepoints_t = np.sort(np.array(
+                (self.changepoints - self.start) / self.t_scale))
         else:
-            row_index = pd.DatetimeIndex(self.history['ds'])
-            indexes = []
-            for cp in self.changepoints:
-                # In the future this may raise a KeyError, but for now we
-                # should guarantee that all changepoint dates are included in
-                # the historical data.
-                indexes.append(row_index.get_loc(cp))
-            return np.array(indexes).astype(np.int)
-
-    def get_changepoint_times(self):
-        cpi = self.get_changepoint_indexes()
-        return np.array(self.history['t'].iloc[cpi])
+            self.changepoints_t = np.array([0])  # dummy changepoint
 
-    def get_changepoint_matrix(self):
-        changepoint_indexes = self.get_changepoint_indexes()
-
-        A = np.zeros((self.history.shape[0], len(changepoint_indexes)))
-        for i, index in enumerate(changepoint_indexes):
-            A[index:self.history.shape[0], i] = 1
 
+    def get_changepoint_matrix(self):
+        A = np.zeros((self.history.shape[0], len(self.changepoints_t)))
+        for i, t_i in enumerate(self.changepoints_t):
+            A[self.history['t'].values >= t_i, i] = 1
         return A
 
     @staticmethod
@@ -345,7 +332,6 @@ class Prophet(object):
         The fitted Prophet object.
         """
         history = df[df['y'].notnull()].copy()
-        history.reset_index(inplace=True, drop=True)
 
         history = self.setup_dataframe(history, initialize_scales=True)
         self.history = history
@@ -353,17 +339,15 @@ class Prophet(object):
 
         self.set_changepoints()
         A = self.get_changepoint_matrix()
-        changepoint_indexes = self.get_changepoint_indexes()
 
         dat = {
             'T': history.shape[0],
             'K': seasonal_features.shape[1],
-            'S': len(changepoint_indexes),
+            'S': len(self.changepoints_t),
             'y': history['y_scaled'],
             't': history['t'],
             'A': A,
-            # Need to add one because Stan is 1-indexed.
-            's_indx': changepoint_indexes + 1,
+            't_change': self.changepoints_t,
             'X': seasonal_features,
             'sigma': self.seasonality_prior_scale,
             'tau': self.changepoint_prior_scale,
@@ -381,7 +365,7 @@ class Prophet(object):
             return {
                 'k': kinit[0],
                 'm': kinit[1],
-                'delta': np.zeros(len(changepoint_indexes)),
+                'delta': np.zeros(len(self.changepoints_t)),
                 'beta': np.zeros(seasonal_features.shape[1]),
                 'sigma_obs': 1,
             }
@@ -419,7 +403,7 @@ class Prophet(object):
         `df` can be None, in which case we predict only on history.
         """
         if df is None:
-            df = self.history
+            df = self.history.copy()
         else:
             df = self.setup_dataframe(df)
 
@@ -469,12 +453,12 @@ class Prophet(object):
         deltas = np.nanmean(self.params['delta'], axis=0)
 
         t = np.array(df['t'])
-        cpts = self.get_changepoint_times()
         if self.growth == 'linear':
-            trend = self.piecewise_linear(t, deltas, k, m, cpts)
+            trend = self.piecewise_linear(t, deltas, k, m, self.changepoints_t)
         else:
             cap = df['cap_scaled']
-            trend = self.piecewise_logistic(t, cap, deltas, k, m, cpts)
+            trend = self.piecewise_logistic(
+                t, cap, deltas, k, m, self.changepoints_t)
 
         return trend * self.y_scale
 
@@ -565,7 +549,6 @@ class Prophet(object):
         deltas = self.params['delta'][iteration]
 
         t = np.array(df['t'])
-        changepoint_ts = self.get_changepoint_times()
         T = t.max()
 
         if T > 1:
@@ -574,7 +557,7 @@ class Prophet(object):
             dt = np.min(dt[dt > 0])
             # Number of time periods in the future
             N = np.ceil((T - 1) / float(dt))
-            S = len(changepoint_ts)
+            S = len(self.changepoints_t)
 
             prob_change = min(1, (S * (T - 1)) / N)
             n_changes = np.random.binomial(N, prob_change)
@@ -593,7 +576,8 @@ class Prophet(object):
         deltas_new = np.random.laplace(0, lambda_, n_changes)
 
         # Prepend the times and deltas from the history
-        changepoint_ts = np.concatenate((changepoint_ts, changepoint_ts_new))
+        changepoint_ts = np.concatenate((self.changepoints_t,
+                                         changepoint_ts_new))
         deltas = np.concatenate((deltas, deltas_new))
 
         if self.growth == 'linear':

python/fbprophet/tests/test_prophet.py

@@ -57,6 +57,25 @@ class TestProphet(TestCase):
         forecaster.fit(train)
         forecaster.predict(future)
 
+    def test_fit_changepoint_not_in_history(self):
+        train = DATA[(DATA['ds'] < '2013-01-01') | (DATA['ds'] > '2014-01-01')]
+        train[(train['ds'] > '2014-01-01')] += 20
+        future = pd.DataFrame({'ds': DATA['ds']})
+        forecaster = Prophet(changepoints=['2013-06-06'])
+        forecaster.fit(train)
+        forecaster.predict(future)
+
+    def test_fit_predict_duplicates(self):
+        N = DATA.shape[0]
+        train1 = DATA.head(N // 2).copy()
+        train2 = DATA.head(N // 2).copy()
+        train2['y'] += 10
+        train = train1.append(train2)
+        future = pd.DataFrame({'ds': DATA['ds'].tail(N // 2)})
+        forecaster = Prophet()
+        forecaster.fit(train)
+        forecaster.predict(future)
+
     def test_setup_dataframe(self):
         m = Prophet()
         N = DATA.shape[0]
@@ -81,7 +100,7 @@ class TestProphet(TestCase):
 
         m.set_changepoints()
 
-        cp = m.get_changepoint_indexes()
+        cp = m.changepoints_t
         self.assertEqual(cp.shape[0], m.n_changepoints)
         self.assertEqual(len(cp.shape), 1)
         self.assertTrue(cp.min() > 0)
@@ -100,7 +119,7 @@ class TestProphet(TestCase):
         m.history = history
 
         m.set_changepoints()
-        cp = m.get_changepoint_indexes()
+        cp = m.changepoints_t
         self.assertEqual(cp.shape[0], 1)
         self.assertEqual(cp[0], 0)
 

python/stan/prophet_linear_growth.stan

@@ -1,18 +1,11 @@
-# Copyright (c) 2017-present, Facebook, Inc.
-# All rights reserved.
-#
-# This source code is licensed under the BSD-style license found in the
-# LICENSE file in the root directory of this source tree. An additional grant 
-# of patent rights can be found in the PATENTS file in the same directory.
-
 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 split points
+  int S;                                // Number of changepoints
   matrix[T, S] A;                   // Split indicators
-  int s_indx[S];                 // Index of split points
+  real t_change[S];                 // Index of changepoints
   matrix[T,K] X;                // season vectors
   real<lower=0> sigma;              // scale on seasonality prior
   real<lower=0> tau;                  // scale on changepoints prior
@@ -30,7 +23,7 @@ transformed parameters {
   vector[S] gamma;                  // adjusted offsets, for piecewise continuity
 
   for (i in 1:S) {
-    gamma[i] = -t[s_indx[i]] * delta[i];
+    gamma[i] = -t_change[i] * delta[i];
   }
 }
 

python/stan/prophet_logistic_growth.stan

@@ -1,34 +1,17 @@
-# Copyright (c) 2017-present, Facebook, Inc.
-# All rights reserved.
-#
-# This source code is licensed under the BSD-style license found in the
-# LICENSE file in the root directory of this source tree. An additional grant 
-# of patent rights can be found in the PATENTS file in the same directory.
-
 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 split points
+  int S;                                // Number of changepoints
   matrix[T, S] A;                   // Split indicators
-  int s_indx[S];                 // Index of split points
+  real t_change[S];                 // Index of changepoints
   matrix[T,K] X;                    // season vectors
   real<lower=0> sigma;              // scale on seasonality prior
   real<lower=0> tau;                  // scale on changepoints prior
 }
 
-
-transformed data {
-  int s_ext[S + 1];  // Segment endpoints
-  for (j in 1:S) {
-    s_ext[j] = s_indx[j];
-  }
-  s_ext[S + 1] = T + 1;  // Used for the m_adj loop below.
-}
-
-
 parameters {
   real k;                            // Base growth rate
   real m;                            // offset
@@ -37,7 +20,6 @@ parameters {
   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
@@ -52,7 +34,7 @@ transformed parameters {
   // Piecewise offsets
   m_pr = m; // The offset in the previous segment
   for (i in 1:S) {
-    gamma[i] = (t[s_indx[i]] - m_pr) * (1 - k_s[i] / k_s[i + 1]);
+    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
   }
 }