From 443d475468f7425fecf6011d0bdf35b89c508694 Mon Sep 17 00:00:00 2001
From: Ben Letham <bletham@bletham>
Date: Tue, 28 Feb 2017 10:08:37 +0200
Subject: [PATCH] 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                            | 64 ++++++++----------------
 R/inst/stan/prophet_linear_growth.stan   |  6 +--
 R/inst/stan/prophet_logistic_growth.stan | 17 ++-----
 R/man/get_changepoint_indexes.Rd         | 18 -------
 R/man/get_changepoint_times.Rd           | 18 -------
 R/tests/testthat/test_prophet.R          | 25 +++++++--
 python/fbprophet/forecaster.py           | 62 +++++++++--------------
 python/fbprophet/tests/test_prophet.py   | 23 ++++++++-
 python/stan/prophet_linear_growth.stan   | 13 ++---
 python/stan/prophet_logistic_growth.stan | 24 ++-------
 10 files changed, 98 insertions(+), 172 deletions(-)
 delete mode 100644 R/man/get_changepoint_indexes.Rd
 delete mode 100644 R/man/get_changepoint_times.Rd

diff --git a/R/R/prophet.R b/R/R/prophet.R
index 68b9377..f52e3e0 100644
--- a/R/R/prophet.R
+++ b/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
diff --git a/R/inst/stan/prophet_linear_growth.stan b/R/inst/stan/prophet_linear_growth.stan
index 27c6458..33b01ef 100644
--- a/R/inst/stan/prophet_linear_growth.stan
+++ b/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];
   }
 }
 
diff --git a/R/inst/stan/prophet_logistic_growth.stan b/R/inst/stan/prophet_logistic_growth.stan
index bea299c..3338d13 100644
--- a/R/inst/stan/prophet_logistic_growth.stan
+++ b/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
   }
 }
diff --git a/R/man/get_changepoint_indexes.Rd b/R/man/get_changepoint_indexes.Rd
deleted file mode 100644
index fc207a5..0000000
--- a/R/man/get_changepoint_indexes.Rd
+++ /dev/null
@@ -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.
-}
-
diff --git a/R/man/get_changepoint_times.Rd b/R/man/get_changepoint_times.Rd
deleted file mode 100644
index 36740cd..0000000
--- a/R/man/get_changepoint_times.Rd
+++ /dev/null
@@ -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.
-}
-
diff --git a/R/tests/testthat/test_prophet.R b/R/tests/testthat/test_prophet.R
index 771162c..008688a 100644
--- a/R/tests/testthat/test_prophet.R
+++ b/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))
diff --git a/python/fbprophet/forecaster.py b/python/fbprophet/forecaster.py
index 2860292..aad9880 100644
--- a/python/fbprophet/forecaster.py
+++ b/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)
+            self.changepoints_t = np.array([0])  # dummy changepoint
 
-    def get_changepoint_times(self):
-        cpi = self.get_changepoint_indexes()
-        return np.array(self.history['t'].iloc[cpi])
 
     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
-
+        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':
diff --git a/python/fbprophet/tests/test_prophet.py b/python/fbprophet/tests/test_prophet.py
index dbc50bd..3823026 100644
--- a/python/fbprophet/tests/test_prophet.py
+++ b/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)
 
diff --git a/python/stan/prophet_linear_growth.stan b/python/stan/prophet_linear_growth.stan
index a94ea92..33b01ef 100644
--- a/python/stan/prophet_linear_growth.stan
+++ b/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];
   }
 }
 
diff --git a/python/stan/prophet_logistic_growth.stan b/python/stan/prophet_logistic_growth.stan
index 8d04dfb..3338d13 100644
--- a/python/stan/prophet_logistic_growth.stan
+++ b/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
   }
 }