From cdc313584eae01d1ceda895758b6faf3ea9351dc Mon Sep 17 00:00:00 2001
From: Ben Letham <bletham@gmail.com>
Date: Tue, 29 May 2018 21:47:11 -0700
Subject: [PATCH] Switch future changepoint generation to fully continuous
 model

---
 R/R/prophet.R                  | 28 ++++++++--------------------
 python/fbprophet/forecaster.py | 20 +++++++-------------
 2 files changed, 15 insertions(+), 33 deletions(-)

diff --git a/R/R/prophet.R b/R/R/prophet.R
index dabdcf8..8a2759e 100644
--- a/R/R/prophet.R
+++ b/R/R/prophet.R
@@ -1447,30 +1447,18 @@ sample_predictive_trend <- function(model, df, iteration) {
   t <- df$t
   T <- max(t)
 
+  # New changepoints from a Poisson process with rate S on [1, T]
   if (T > 1) {
-    # Get the time discretization of the history
-    dt <- diff(model$history$t)
-    dt <- min(dt[dt > 0])
-    # Number of time periods in the future
-    N <- ceiling((T - 1) / dt)
     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
-    # of S * (T - 1) changepoints in expectation.
-    prob.change <- min(1, (S * (T - 1)) / N)
-    # This calculation works for both history and df not uniformly spaced.
-    n.changes <- stats::rbinom(1, N, prob.change)
-
-    # Sample ts
-    if (n.changes == 0) {
-      changepoint.ts.new <- c()
-    } else {
-      changepoint.ts.new <- sort(stats::runif(n.changes, min = 1, max = T))
-    }
+    n.changes <- stats::rpois(1, S * (T - 1))
+  } else {
+    n.changes <- 0
+  }
+  if (n.changes > 0) {
+    changepoint.ts.new <- 1 + stats::runif(n.changes) * (T - 1)
+    changepoint.ts.new <- sort(changepoint.ts.new)
   } else {
     changepoint.ts.new <- c()
-    n.changes <- 0
   }
 
   # Get the empirical scale of the deltas, plus epsilon to avoid NaNs.
diff --git a/python/fbprophet/forecaster.py b/python/fbprophet/forecaster.py
index acce127..be7ed9e 100644
--- a/python/fbprophet/forecaster.py
+++ b/python/fbprophet/forecaster.py
@@ -1303,23 +1303,17 @@ class Prophet(object):
         t = np.array(df['t'])
         T = t.max()
 
+        # New changepoints from a Poisson process with rate S on [1, T]
         if T > 1:
-            # Get the time discretization of the history
-            dt = np.diff(self.history['t'])
-            dt = np.min(dt[dt > 0])
-            # Number of time periods in the future
-            N = np.ceil((T - 1) / float(dt))
             S = len(self.changepoints_t)
-
-            prob_change = min(1, (S * (T - 1)) / N)
-            n_changes = np.random.binomial(N, prob_change)
-
-            # Sample ts
-            changepoint_ts_new = sorted(np.random.uniform(1, T, n_changes))
+            n_changes = np.random.poisson(S * (T - 1))
         else:
-            # Case where we're not extrapolating.
-            changepoint_ts_new = []
             n_changes = 0
+        if n_changes > 0:
+            changepoint_ts_new = 1 + np.random.rand(n_changes) * (T - 1)
+            changepoint_ts_new.sort()
+        else:
+            changepoint_ts_new = []
 
         # Get the empirical scale of the deltas, plus epsilon to avoid NaNs.
         lambda_ = np.mean(np.abs(deltas)) + 1e-8