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
   }
 }