From 55d7d1e62d6ecd4e5d3ca5849ecc59d1704311c6 Mon Sep 17 00:00:00 2001 From: Ben Letham Date: Fri, 4 May 2018 16:04:29 -0700 Subject: [PATCH] Single stan model with both trends (Py) --- python/MANIFEST.in | 4 +- python/fbprophet/forecaster.py | 15 +- python/fbprophet/models.py | 9 +- python/fbprophet/tests/test_prophet.py | 10 +- python/setup.py | 29 ++-- python/stan/unix/prophet.stan | 123 +++++++++++++++ python/stan/unix/prophet_linear_growth.stan | 40 ----- python/stan/unix/prophet_logistic_growth.stan | 52 ------- python/stan/win/prophet.stan | 142 ++++++++++++++++++ python/stan/win/prophet_linear_growth.stan | 45 ------ python/stan/win/prophet_logistic_growth.stan | 57 ------- 11 files changed, 289 insertions(+), 237 deletions(-) create mode 100644 python/stan/unix/prophet.stan delete mode 100644 python/stan/unix/prophet_linear_growth.stan delete mode 100644 python/stan/unix/prophet_logistic_growth.stan create mode 100644 python/stan/win/prophet.stan delete mode 100644 python/stan/win/prophet_linear_growth.stan delete mode 100644 python/stan/win/prophet_logistic_growth.stan diff --git a/python/MANIFEST.in b/python/MANIFEST.in index 10c5fe3..0b7433d 100644 --- a/python/MANIFEST.in +++ b/python/MANIFEST.in @@ -3,7 +3,7 @@ include stan/win/*.stan include LICENSE # Ensure in-place built models do not get included in the source dist. -prune fbprophet/stan_models +prune fbprophet/stan_model # Necessary for tests to run -include fbprophet/tests/*.csv \ No newline at end of file +include fbprophet/tests/*.csv diff --git a/python/fbprophet/forecaster.py b/python/fbprophet/forecaster.py index 7856229..fd455ee 100644 --- a/python/fbprophet/forecaster.py +++ b/python/fbprophet/forecaster.py @@ -19,7 +19,7 @@ import numpy as np import pandas as pd # fb-block 1 start -from fbprophet.models import prophet_stan_models +from fbprophet.models import prophet_stan_model from fbprophet.plot import ( plot, plot_components, @@ -348,13 +348,6 @@ class Prophet(object): else: self.changepoints_t = np.array([0]) # dummy changepoint - def get_changepoint_matrix(self): - """Gets changepoint matrix for history dataframe.""" - 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 def fourier_series(dates, period, series_order): """Provides Fourier series components with the specified frequency @@ -790,7 +783,6 @@ class Prophet(object): self.make_all_seasonality_features(history)) self.set_changepoints() - A = self.get_changepoint_matrix() dat = { 'T': history.shape[0], @@ -798,20 +790,21 @@ class Prophet(object): 'S': len(self.changepoints_t), 'y': history['y_scaled'], 't': history['t'], - 'A': A, 't_change': self.changepoints_t, 'X': seasonal_features, 'sigmas': prior_scales, 'tau': self.changepoint_prior_scale, + 'trend_indicator': int(self.growth == 'logistic'), } if self.growth == 'linear': + dat['cap'] = np.zeros(self.history.shape[0]) kinit = self.linear_growth_init(history) else: dat['cap'] = history['cap_scaled'] kinit = self.logistic_growth_init(history) - model = prophet_stan_models[self.growth] + model = prophet_stan_model def stan_init(): return { diff --git a/python/fbprophet/models.py b/python/fbprophet/models.py index 1a1fdbe..50c2a35 100644 --- a/python/fbprophet/models.py +++ b/python/fbprophet/models.py @@ -18,20 +18,17 @@ import pkg_resources # fb-block 2 -def get_prophet_stan_model(model): +def get_prophet_stan_model(): """Load compiled Stan model""" # fb-block 3 # fb-block 4 start model_file = pkg_resources.resource_filename( 'fbprophet', - 'stan_models/{}_growth.pkl'.format(model), + 'stan_model/prophet_model.pkl', ) # fb-block 4 end with open(model_file, 'rb') as f: return pickle.load(f) -prophet_stan_models = { - 'linear': get_prophet_stan_model('linear'), - 'logistic': get_prophet_stan_model('logistic'), -} +prophet_stan_model = get_prophet_stan_model() diff --git a/python/fbprophet/tests/test_prophet.py b/python/fbprophet/tests/test_prophet.py index 254701a..c774847 100644 --- a/python/fbprophet/tests/test_prophet.py +++ b/python/fbprophet/tests/test_prophet.py @@ -151,11 +151,7 @@ class TestProphet(TestCase): self.assertEqual(cp.shape[0], m.n_changepoints) self.assertEqual(len(cp.shape), 1) self.assertTrue(cp.min() > 0) - self.assertTrue(cp.max() < N) - - mat = m.get_changepoint_matrix() - self.assertEqual(mat.shape[0], N // 2) - self.assertEqual(mat.shape[1], m.n_changepoints) + self.assertTrue(cp.max() < 1) def test_get_zero_changepoints(self): m = Prophet(n_changepoints=0) @@ -170,10 +166,6 @@ class TestProphet(TestCase): self.assertEqual(cp.shape[0], 1) self.assertEqual(cp[0], 0) - mat = m.get_changepoint_matrix() - self.assertEqual(mat.shape[0], N // 2) - self.assertEqual(mat.shape[1], 1) - def test_override_n_changepoints(self): m = Prophet() history = DATA.head(20).copy() diff --git a/python/setup.py b/python/setup.py index d528ced..66fc591 100644 --- a/python/setup.py +++ b/python/setup.py @@ -20,20 +20,19 @@ if platform.platform().startswith('Win'): PLATFORM = 'win' SETUP_DIR = os.path.dirname(os.path.abspath(__file__)) -MODELS_DIR = os.path.join(SETUP_DIR, 'stan', PLATFORM) -MODELS_TARGET_DIR = os.path.join('fbprophet', 'stan_models') +MODEL_DIR = os.path.join(SETUP_DIR, 'stan', PLATFORM) +MODEL_TARGET_DIR = os.path.join('fbprophet', 'stan_model') -def build_stan_models(target_dir, models_dir=MODELS_DIR): +def build_stan_model(target_dir, model_dir=MODEL_DIR): from pystan import StanModel - for model_type in ['linear', 'logistic']: - model_name = 'prophet_{}_growth.stan'.format(model_type) - target_name = '{}_growth.pkl'.format(model_type) - with open(os.path.join(models_dir, model_name)) as f: - model_code = f.read() - sm = StanModel(model_code=model_code) - with open(os.path.join(target_dir, target_name), 'wb') as f: - pickle.dump(sm, f, protocol=pickle.HIGHEST_PROTOCOL) + model_name = 'prophet.stan' + target_name = 'prophet_model.pkl' + with open(os.path.join(model_dir, model_name)) as f: + model_code = f.read() + sm = StanModel(model_code=model_code) + with open(os.path.join(target_dir, target_name), 'wb') as f: + pickle.dump(sm, f, protocol=pickle.HIGHEST_PROTOCOL) class BuildPyCommand(build_py): @@ -41,9 +40,9 @@ class BuildPyCommand(build_py): def run(self): if not self.dry_run: - target_dir = os.path.join(self.build_lib, MODELS_TARGET_DIR) + target_dir = os.path.join(self.build_lib, MODEL_TARGET_DIR) self.mkpath(target_dir) - build_stan_models(target_dir) + build_stan_model(target_dir) build_py.run(self) @@ -53,9 +52,9 @@ class DevelopCommand(develop): def run(self): if not self.dry_run: - target_dir = os.path.join(self.setup_path, MODELS_TARGET_DIR) + target_dir = os.path.join(self.setup_path, MODEL_TARGET_DIR) self.mkpath(target_dir) - build_stan_models(target_dir) + build_stan_model(target_dir) develop.run(self) diff --git a/python/stan/unix/prophet.stan b/python/stan/unix/prophet.stan new file mode 100644 index 0000000..746b0cb --- /dev/null +++ b/python/stan/unix/prophet.stan @@ -0,0 +1,123 @@ +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 += 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[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 + } + 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 ./ (1 + exp(-(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 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 tau; // Scale on changepoints prior + int trend_indicator; // 0 for linear, 1 for logistic +} + +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 sigma_obs; // Observation noise + vector[K] beta; // Regressor coefficients +} + +transformed parameters { + vector[T] trend; + + if (trend_indicator == 0) { + trend = linear_trend(k, m, delta, t, A, t_change); + } else if (trend_indicator == 1) { + trend = logistic_trend(k, m, delta, t, cap, A, t_change, S); + } +} + +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(trend + X * beta, sigma_obs); +} diff --git a/python/stan/unix/prophet_linear_growth.stan b/python/stan/unix/prophet_linear_growth.stan deleted file mode 100644 index 1639f4a..0000000 --- a/python/stan/unix/prophet_linear_growth.stan +++ /dev/null @@ -1,40 +0,0 @@ -data { - int T; // Sample size - int 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 tau; // scale on changepoints prior -} - -parameters { - real k; // Base growth rate - real m; // offset - vector[S] delta; // Rate adjustments - real 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); -} diff --git a/python/stan/unix/prophet_logistic_growth.stan b/python/stan/unix/prophet_logistic_growth.stan deleted file mode 100644 index 2464ac5..0000000 --- a/python/stan/unix/prophet_logistic_growth.stan +++ /dev/null @@ -1,52 +0,0 @@ -data { - int T; // Sample size - int 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 tau; // scale on changepoints prior -} - -parameters { - real k; // Base growth rate - real m; // offset - vector[S] delta; // Rate adjustments - real 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); -} diff --git a/python/stan/win/prophet.stan b/python/stan/win/prophet.stan new file mode 100644 index 0000000..34e3163 --- /dev/null +++ b/python/stan/win/prophet.stan @@ -0,0 +1,142 @@ +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 += 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[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 + } + return gamma; + } + + vector logistic_trend( + real k, + real m, + vector delta, + vector t, + vector cap, + matrix A, + vector t_change, + int S, + int T + ) { + vector[S] gamma; + vector[T] Y; + + gamma = logistic_gamma(k, m, delta, t_change, S); + for (i in 1:T) { + Y[i] = cap[i] / (1 + exp(-(k + dot_product(A[i], delta)) * (t[i] - (m + dot_product(A[i], gamma))))) + } + return Y; + } + + // Linear trend function + + vector linear_trend( + real k, + real m, + vector delta, + vector t, + matrix A, + vector t_change, + int S, + int T + ) { + vector[S] gamma; + vector[T] Y; + + gamma = (-t_change .* delta); + for (i in 1:T) { + Y[i] = (k + dot_product(A[i], delta)) * t[i] + (m + dot_product(A[i], gamma)) + } + return Y; + } +} + +data { + int T; // Number of time periods + int 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 tau; // Scale on changepoints prior + int trend_indicator; // 0 for linear, 1 for logistic +} + +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 sigma_obs; // Observation noise + vector[K] beta; // Regressor coefficients +} + +transformed parameters { + vector[T] trend; + vector[T] Y; + + if (trend_indicator == 0) { + trend = linear_trend(k, m, delta, t, A, t_change, S, T); + } else if (trend_indicator == 1) { + trend = logistic_trend(k, m, delta, t, cap, A, t_change, S, T); + } + + for (i in 1:T) { + Y[i] = trend[i] + dot_product(X[i], beta); + } +} + +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(Y, sigma_obs); +} diff --git a/python/stan/win/prophet_linear_growth.stan b/python/stan/win/prophet_linear_growth.stan deleted file mode 100644 index 949fab5..0000000 --- a/python/stan/win/prophet_linear_growth.stan +++ /dev/null @@ -1,45 +0,0 @@ -data { - int T; // Sample size - int K; // Number of seasonal vectors - real t[T]; // Day - real y[T]; // Time-series - int S; // Number of changepoints - real A[T, S]; // Split indicators - real t_change[S]; // Index of changepoints - real X[T,K]; // season vectors - vector[K] sigmas; // scale on seasonality prior - real tau; // scale on changepoints prior -} - -parameters { - real k; // Base growth rate - real m; // offset - real delta[S]; // Rate adjustments - real sigma_obs; // Observation noise (incl. seasonal variation) - real beta[K]; // seasonal vector -} - -transformed parameters { - real gamma[S]; // adjusted offsets, for piecewise continuity - - for (i in 1:S) { - gamma[i] = -t_change[i] * delta[i]; - } -} - -model { - real Y[T]; - - //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 - for (i in 1:T) { - Y[i] = (dot_product(A[i], delta) + k) * t[i] + (dot_product(A[i], gamma) + m) + dot_product(X[i], beta); - } - y ~ normal(Y, sigma_obs); -} diff --git a/python/stan/win/prophet_logistic_growth.stan b/python/stan/win/prophet_logistic_growth.stan deleted file mode 100644 index d24b0be..0000000 --- a/python/stan/win/prophet_logistic_growth.stan +++ /dev/null @@ -1,57 +0,0 @@ -data { - int T; // Sample size - int K; // Number of seasonal vectors - real t[T]; // Day - real cap[T]; // Capacities - real y[T]; // Time-series - int S; // Number of changepoints - real A[T, S]; // Split indicators - real t_change[S]; // Index of changepoints - real X[T,K]; // season vectors - vector[K] sigmas; // scale on seasonality prior - real tau; // scale on changepoints prior -} - -parameters { - real k; // Base growth rate - real m; // offset - real delta[S]; // Rate adjustments - real sigma_obs; // Observation noise (incl. seasonal variation) - real beta[K]; // seasonal vector -} - -transformed parameters { - real gamma[S]; // adjusted offsets, for piecewise continuity - real k_s[S + 1]; // 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 { - real Y[T]; - - //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 - for (i in 1:T) { - Y[i] = cap[i] / (1 + exp(-(k + dot_product(A[i], delta)) * (t[i] - (m + dot_product(A[i], gamma))))) + dot_product(X[i], beta); - } - y ~ normal(Y, sigma_obs); -}