forecaster.py 57 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540
  1. # Copyright (c) 2017-present, Facebook, Inc.
  2. # All rights reserved.
  3. #
  4. # This source code is licensed under the BSD-style license found in the
  5. # LICENSE file in the root directory of this source tree. An additional grant
  6. # of patent rights can be found in the PATENTS file in the same directory.
  7. from __future__ import absolute_import
  8. from __future__ import division
  9. from __future__ import print_function
  10. from __future__ import unicode_literals
  11. from collections import defaultdict
  12. from datetime import timedelta
  13. import logging
  14. logger = logging.getLogger(__name__)
  15. from matplotlib import pyplot as plt
  16. from matplotlib.dates import MonthLocator, num2date
  17. from matplotlib.ticker import FuncFormatter
  18. import numpy as np
  19. import pandas as pd
  20. # fb-block 1 start
  21. from fbprophet.models import prophet_stan_models
  22. # fb-block 1 end
  23. try:
  24. import pystan # noqa F401
  25. except ImportError:
  26. logger.error('You cannot run prophet without pystan installed')
  27. raise
  28. # fb-block 2
  29. class Prophet(object):
  30. """Prophet forecaster.
  31. Parameters
  32. ----------
  33. growth: String 'linear' or 'logistic' to specify a linear or logistic
  34. trend.
  35. changepoints: List of dates at which to include potential changepoints. If
  36. not specified, potential changepoints are selected automatically.
  37. n_changepoints: Number of potential changepoints to include. Not used
  38. if input `changepoints` is supplied. If `changepoints` is not supplied,
  39. then n_changepoints potential changepoints are selected uniformly from
  40. the first 80 percent of the history.
  41. yearly_seasonality: Fit yearly seasonality.
  42. Can be 'auto', True, False, or a number of Fourier terms to generate.
  43. weekly_seasonality: Fit weekly seasonality.
  44. Can be 'auto', True, False, or a number of Fourier terms to generate.
  45. daily_seasonality: Fit daily seasonality.
  46. Can be 'auto', True, False, or a number of Fourier terms to generate.
  47. holidays: pd.DataFrame with columns holiday (string) and ds (date type)
  48. and optionally columns lower_window and upper_window which specify a
  49. range of days around the date to be included as holidays.
  50. lower_window=-2 will include 2 days prior to the date as holidays. Also
  51. optionally can have a column prior_scale specifying the prior scale for
  52. that holiday.
  53. seasonality_prior_scale: Parameter modulating the strength of the
  54. seasonality model. Larger values allow the model to fit larger seasonal
  55. fluctuations, smaller values dampen the seasonality. Can be specified
  56. for individual seasonalities using add_seasonality.
  57. holidays_prior_scale: Parameter modulating the strength of the holiday
  58. components model, unless overridden in the holidays input.
  59. changepoint_prior_scale: Parameter modulating the flexibility of the
  60. automatic changepoint selection. Large values will allow many
  61. changepoints, small values will allow few changepoints.
  62. mcmc_samples: Integer, if greater than 0, will do full Bayesian inference
  63. with the specified number of MCMC samples. If 0, will do MAP
  64. estimation.
  65. interval_width: Float, width of the uncertainty intervals provided
  66. for the forecast. If mcmc_samples=0, this will be only the uncertainty
  67. in the trend using the MAP estimate of the extrapolated generative
  68. model. If mcmc.samples>0, this will be integrated over all model
  69. parameters, which will include uncertainty in seasonality.
  70. uncertainty_samples: Number of simulated draws used to estimate
  71. uncertainty intervals.
  72. """
  73. def __init__(
  74. self,
  75. growth='linear',
  76. changepoints=None,
  77. n_changepoints=25,
  78. yearly_seasonality='auto',
  79. weekly_seasonality='auto',
  80. daily_seasonality='auto',
  81. holidays=None,
  82. seasonality_prior_scale=10.0,
  83. holidays_prior_scale=10.0,
  84. changepoint_prior_scale=0.05,
  85. mcmc_samples=0,
  86. interval_width=0.80,
  87. uncertainty_samples=1000,
  88. ):
  89. self.growth = growth
  90. self.changepoints = pd.to_datetime(changepoints)
  91. if self.changepoints is not None:
  92. self.n_changepoints = len(self.changepoints)
  93. self.specified_changepoints = True
  94. else:
  95. self.n_changepoints = n_changepoints
  96. self.specified_changepoints = False
  97. self.yearly_seasonality = yearly_seasonality
  98. self.weekly_seasonality = weekly_seasonality
  99. self.daily_seasonality = daily_seasonality
  100. if holidays is not None:
  101. if not (
  102. isinstance(holidays, pd.DataFrame)
  103. and 'ds' in holidays # noqa W503
  104. and 'holiday' in holidays # noqa W503
  105. ):
  106. raise ValueError("holidays must be a DataFrame with 'ds' and "
  107. "'holiday' columns.")
  108. holidays['ds'] = pd.to_datetime(holidays['ds'])
  109. self.holidays = holidays
  110. self.seasonality_prior_scale = float(seasonality_prior_scale)
  111. self.changepoint_prior_scale = float(changepoint_prior_scale)
  112. self.holidays_prior_scale = float(holidays_prior_scale)
  113. self.mcmc_samples = mcmc_samples
  114. self.interval_width = interval_width
  115. self.uncertainty_samples = uncertainty_samples
  116. # Set during fitting
  117. self.start = None
  118. self.y_scale = None
  119. self.logistic_floor = False
  120. self.t_scale = None
  121. self.changepoints_t = None
  122. self.seasonalities = {}
  123. self.extra_regressors = {}
  124. self.stan_fit = None
  125. self.params = {}
  126. self.history = None
  127. self.history_dates = None
  128. self.validate_inputs()
  129. def validate_inputs(self):
  130. """Validates the inputs to Prophet."""
  131. if self.growth not in ('linear', 'logistic'):
  132. raise ValueError(
  133. "Parameter 'growth' should be 'linear' or 'logistic'.")
  134. if self.holidays is not None:
  135. has_lower = 'lower_window' in self.holidays
  136. has_upper = 'upper_window' in self.holidays
  137. if has_lower + has_upper == 1:
  138. raise ValueError('Holidays must have both lower_window and ' +
  139. 'upper_window, or neither')
  140. if has_lower:
  141. if max(self.holidays['lower_window']) > 0:
  142. raise ValueError('Holiday lower_window should be <= 0')
  143. if min(self.holidays['upper_window']) < 0:
  144. raise ValueError('Holiday upper_window should be >= 0')
  145. for h in self.holidays['holiday'].unique():
  146. self.validate_column_name(h, check_holidays=False)
  147. def validate_column_name(self, name, check_holidays=True,
  148. check_seasonalities=True, check_regressors=True):
  149. """Validates the name of a seasonality, holiday, or regressor.
  150. Parameters
  151. ----------
  152. name: string
  153. check_holidays: bool check if name already used for holiday
  154. check_seasonalities: bool check if name already used for seasonality
  155. check_regressors: bool check if name already used for regressor
  156. """
  157. if '_delim_' in name:
  158. raise ValueError('Name cannot contain "_delim_"')
  159. reserved_names = [
  160. 'trend', 'seasonal', 'seasonalities', 'daily', 'weekly', 'yearly',
  161. 'holidays', 'zeros', 'extra_regressors', 'yhat'
  162. ]
  163. rn_l = [n + '_lower' for n in reserved_names]
  164. rn_u = [n + '_upper' for n in reserved_names]
  165. reserved_names.extend(rn_l)
  166. reserved_names.extend(rn_u)
  167. reserved_names.extend([
  168. 'ds', 'y', 'cap', 'floor', 'y_scaled', 'cap_scaled'])
  169. if name in reserved_names:
  170. raise ValueError('Name "{}" is reserved.'.format(name))
  171. if (check_holidays and self.holidays is not None and
  172. name in self.holidays['holiday'].unique()):
  173. raise ValueError(
  174. 'Name "{}" already used for a holiday.'.format(name))
  175. if check_seasonalities and name in self.seasonalities:
  176. raise ValueError(
  177. 'Name "{}" already used for a seasonality.'.format(name))
  178. if check_regressors and name in self.extra_regressors:
  179. raise ValueError(
  180. 'Name "{}" already used for an added regressor.'.format(name))
  181. def setup_dataframe(self, df, initialize_scales=False):
  182. """Prepare dataframe for fitting or predicting.
  183. Adds a time index and scales y. Creates auxiliary columns 't', 't_ix',
  184. 'y_scaled', and 'cap_scaled'. These columns are used during both
  185. fitting and predicting.
  186. Parameters
  187. ----------
  188. df: pd.DataFrame with columns ds, y, and cap if logistic growth. Any
  189. specified additional regressors must also be present.
  190. initialize_scales: Boolean set scaling factors in self from df.
  191. Returns
  192. -------
  193. pd.DataFrame prepared for fitting or predicting.
  194. """
  195. if 'y' in df:
  196. df['y'] = pd.to_numeric(df['y'])
  197. if np.isinf(df['y'].values).any():
  198. raise ValueError('Found infinity in column y.')
  199. df['ds'] = pd.to_datetime(df['ds'])
  200. if df['ds'].isnull().any():
  201. raise ValueError('Found NaN in column ds.')
  202. for name in self.extra_regressors:
  203. if name not in df:
  204. raise ValueError(
  205. 'Regressor "{}" missing from dataframe'.format(name))
  206. df = df.sort_values('ds')
  207. df.reset_index(inplace=True, drop=True)
  208. self.initialize_scales(initialize_scales, df)
  209. if self.logistic_floor:
  210. if 'floor' not in df:
  211. raise ValueError("Expected column 'floor'.")
  212. else:
  213. df['floor'] = 0
  214. if self.growth == 'logistic':
  215. assert 'cap' in df
  216. df['cap_scaled'] = (df['cap'] - df['floor']) / self.y_scale
  217. df['t'] = (df['ds'] - self.start) / self.t_scale
  218. if 'y' in df:
  219. df['y_scaled'] = (df['y'] - df['floor']) / self.y_scale
  220. for name, props in self.extra_regressors.items():
  221. df[name] = pd.to_numeric(df[name])
  222. df[name] = ((df[name] - props['mu']) / props['std'])
  223. if df[name].isnull().any():
  224. raise ValueError('Found NaN in column ' + name)
  225. return df
  226. def initialize_scales(self, initialize_scales, df):
  227. """Initialize model scales.
  228. Sets model scaling factors using df.
  229. Parameters
  230. ----------
  231. initialize_scales: Boolean set the scales or not.
  232. df: pd.DataFrame for setting scales.
  233. """
  234. if not initialize_scales:
  235. return
  236. if self.growth == 'logistic' and 'floor' in df:
  237. self.logistic_floor = True
  238. floor = df['floor']
  239. else:
  240. floor = 0.
  241. self.y_scale = (df['y'] - floor).abs().max()
  242. if self.y_scale == 0:
  243. self.y_scale = 1
  244. self.start = df['ds'].min()
  245. self.t_scale = df['ds'].max() - self.start
  246. for name, props in self.extra_regressors.items():
  247. standardize = props['standardize']
  248. if standardize == 'auto':
  249. if set(df[name].unique()) == set([1, 0]):
  250. # Don't standardize binary variables.
  251. standardize = False
  252. else:
  253. standardize = True
  254. if standardize:
  255. mu = df[name].mean()
  256. std = df[name].std()
  257. if std == 0:
  258. std = mu
  259. self.extra_regressors[name]['mu'] = mu
  260. self.extra_regressors[name]['std'] = std
  261. def set_changepoints(self):
  262. """Set changepoints
  263. Sets m$changepoints to the dates of changepoints. Either:
  264. 1) The changepoints were passed in explicitly.
  265. A) They are empty.
  266. B) They are not empty, and need validation.
  267. 2) We are generating a grid of them.
  268. 3) The user prefers no changepoints be used.
  269. """
  270. if self.changepoints is not None:
  271. if len(self.changepoints) == 0:
  272. pass
  273. else:
  274. too_low = min(self.changepoints) < self.history['ds'].min()
  275. too_high = max(self.changepoints) > self.history['ds'].max()
  276. if too_low or too_high:
  277. raise ValueError(
  278. 'Changepoints must fall within training data.')
  279. else:
  280. # Place potential changepoints evenly through first 80% of history
  281. hist_size = np.floor(self.history.shape[0] * 0.8)
  282. if self.n_changepoints + 1 > hist_size:
  283. self.n_changepoints = hist_size - 1
  284. logger.info(
  285. 'n_changepoints greater than number of observations.'
  286. 'Using {}.'.format(self.n_changepoints)
  287. )
  288. if self.n_changepoints > 0:
  289. cp_indexes = (
  290. np.linspace(0, hist_size, self.n_changepoints + 1)
  291. .round()
  292. .astype(np.int)
  293. )
  294. self.changepoints = (
  295. self.history.iloc[cp_indexes]['ds'].tail(-1)
  296. )
  297. else:
  298. # set empty changepoints
  299. self.changepoints = []
  300. if len(self.changepoints) > 0:
  301. self.changepoints_t = np.sort(np.array(
  302. (self.changepoints - self.start) / self.t_scale))
  303. else:
  304. self.changepoints_t = np.array([0]) # dummy changepoint
  305. def get_changepoint_matrix(self):
  306. """Gets changepoint matrix for history dataframe."""
  307. A = np.zeros((self.history.shape[0], len(self.changepoints_t)))
  308. for i, t_i in enumerate(self.changepoints_t):
  309. A[self.history['t'].values >= t_i, i] = 1
  310. return A
  311. @staticmethod
  312. def fourier_series(dates, period, series_order):
  313. """Provides Fourier series components with the specified frequency
  314. and order.
  315. Parameters
  316. ----------
  317. dates: pd.Series containing timestamps.
  318. period: Number of days of the period.
  319. series_order: Number of components.
  320. Returns
  321. -------
  322. Matrix with seasonality features.
  323. """
  324. # convert to days since epoch
  325. t = np.array(
  326. (dates - pd.datetime(1970, 1, 1))
  327. .dt.total_seconds()
  328. .astype(np.float)
  329. ) / (3600 * 24.)
  330. return np.column_stack([
  331. fun((2.0 * (i + 1) * np.pi * t / period))
  332. for i in range(series_order)
  333. for fun in (np.sin, np.cos)
  334. ])
  335. @classmethod
  336. def make_seasonality_features(cls, dates, period, series_order, prefix):
  337. """Data frame with seasonality features.
  338. Parameters
  339. ----------
  340. cls: Prophet class.
  341. dates: pd.Series containing timestamps.
  342. period: Number of days of the period.
  343. series_order: Number of components.
  344. prefix: Column name prefix.
  345. Returns
  346. -------
  347. pd.DataFrame with seasonality features.
  348. """
  349. features = cls.fourier_series(dates, period, series_order)
  350. columns = [
  351. '{}_delim_{}'.format(prefix, i + 1)
  352. for i in range(features.shape[1])
  353. ]
  354. return pd.DataFrame(features, columns=columns)
  355. def make_holiday_features(self, dates):
  356. """Construct a dataframe of holiday features.
  357. Parameters
  358. ----------
  359. dates: pd.Series containing timestamps used for computing seasonality.
  360. Returns
  361. -------
  362. holiday_features: pd.DataFrame with a column for each holiday.
  363. prior_scale_list: List of prior scales for each holiday column.
  364. """
  365. # Holds columns of our future matrix.
  366. expanded_holidays = defaultdict(lambda: np.zeros(dates.shape[0]))
  367. prior_scales = {}
  368. # Makes an index so we can perform `get_loc` below.
  369. # Strip to just dates.
  370. row_index = pd.DatetimeIndex(dates.apply(lambda x: x.date()))
  371. for _ix, row in self.holidays.iterrows():
  372. dt = row.ds.date()
  373. try:
  374. lw = int(row.get('lower_window', 0))
  375. uw = int(row.get('upper_window', 0))
  376. except ValueError:
  377. lw = 0
  378. uw = 0
  379. ps = float(row.get('prior_scale', self.holidays_prior_scale))
  380. if np.isnan(ps):
  381. ps = float(self.holidays_prior_scale)
  382. if (
  383. row.holiday in prior_scales and prior_scales[row.holiday] != ps
  384. ):
  385. raise ValueError(
  386. 'Holiday {} does not have consistent prior scale '
  387. 'specification.'.format(row.holiday))
  388. if ps <= 0:
  389. raise ValueError('Prior scale must be > 0')
  390. prior_scales[row.holiday] = ps
  391. for offset in range(lw, uw + 1):
  392. occurrence = dt + timedelta(days=offset)
  393. try:
  394. loc = row_index.get_loc(occurrence)
  395. except KeyError:
  396. loc = None
  397. key = '{}_delim_{}{}'.format(
  398. row.holiday,
  399. '+' if offset >= 0 else '-',
  400. abs(offset)
  401. )
  402. if loc is not None:
  403. expanded_holidays[key][loc] = 1.
  404. else:
  405. # Access key to generate value
  406. expanded_holidays[key]
  407. holiday_features = pd.DataFrame(expanded_holidays)
  408. prior_scale_list = [
  409. prior_scales[h.split('_delim_')[0]]
  410. for h in holiday_features.columns
  411. ]
  412. return holiday_features, prior_scale_list
  413. def add_regressor(self, name, prior_scale=None, standardize='auto'):
  414. """Add an additional regressor to be used for fitting and predicting.
  415. The dataframe passed to `fit` and `predict` will have a column with the
  416. specified name to be used as a regressor. When standardize='auto', the
  417. regressor will be standardized unless it is binary. The regression
  418. coefficient is given a prior with the specified scale parameter.
  419. Decreasing the prior scale will add additional regularization. If no
  420. prior scale is provided, self.holidays_prior_scale will be used.
  421. Parameters
  422. ----------
  423. name: string name of the regressor.
  424. prior_scale: optional float scale for the normal prior. If not
  425. provided, self.holidays_prior_scale will be used.
  426. standardize: optional, specify whether this regressor will be
  427. standardized prior to fitting. Can be 'auto' (standardize if not
  428. binary), True, or False.
  429. Returns
  430. -------
  431. The prophet object.
  432. """
  433. if self.history is not None:
  434. raise Exception(
  435. "Regressors must be added prior to model fitting.")
  436. self.validate_column_name(name, check_regressors=False)
  437. if prior_scale is None:
  438. prior_scale = float(self.holidays_prior_scale)
  439. assert prior_scale > 0
  440. self.extra_regressors[name] = {
  441. 'prior_scale': prior_scale,
  442. 'standardize': standardize,
  443. 'mu': 0.,
  444. 'std': 1.,
  445. }
  446. return self
  447. def add_seasonality(self, name, period, fourier_order, prior_scale=None):
  448. """Add a seasonal component with specified period, number of Fourier
  449. components, and prior scale.
  450. Increasing the number of Fourier components allows the seasonality to
  451. change more quickly (at risk of overfitting). Default values for yearly
  452. and weekly seasonalities are 10 and 3 respectively.
  453. Increasing prior scale will allow this seasonality component more
  454. flexibility, decreasing will dampen it. If not provided, will use the
  455. seasonality_prior_scale provided on Prophet initialization (defaults
  456. to 10).
  457. Parameters
  458. ----------
  459. name: string name of the seasonality component.
  460. period: float number of days in one period.
  461. fourier_order: int number of Fourier components to use.
  462. prior_scale: float prior scale for this component.
  463. Returns
  464. -------
  465. The prophet object.
  466. """
  467. if self.history is not None:
  468. raise Exception(
  469. "Seasonality must be added prior to model fitting.")
  470. if name not in ['daily', 'weekly', 'yearly']:
  471. # Allow overwriting built-in seasonalities
  472. self.validate_column_name(name, check_seasonalities=False)
  473. if prior_scale is None:
  474. ps = self.seasonality_prior_scale
  475. else:
  476. ps = float(prior_scale)
  477. if ps <= 0:
  478. raise ValueError('Prior scale must be > 0')
  479. self.seasonalities[name] = {
  480. 'period': period,
  481. 'fourier_order': fourier_order,
  482. 'prior_scale': ps,
  483. }
  484. return self
  485. def make_all_seasonality_features(self, df):
  486. """Dataframe with seasonality features.
  487. Includes seasonality features, holiday features, and added regressors.
  488. Parameters
  489. ----------
  490. df: pd.DataFrame with dates for computing seasonality features and any
  491. added regressors.
  492. Returns
  493. -------
  494. pd.DataFrame with regression features.
  495. list of prior scales for each column of the features dataframe.
  496. """
  497. seasonal_features = []
  498. prior_scales = []
  499. # Seasonality features
  500. for name, props in self.seasonalities.items():
  501. features = self.make_seasonality_features(
  502. df['ds'],
  503. props['period'],
  504. props['fourier_order'],
  505. name,
  506. )
  507. seasonal_features.append(features)
  508. prior_scales.extend(
  509. [props['prior_scale']] * features.shape[1])
  510. # Holiday features
  511. if self.holidays is not None:
  512. features, holiday_priors = self.make_holiday_features(df['ds'])
  513. seasonal_features.append(features)
  514. prior_scales.extend(holiday_priors)
  515. # Additional regressors
  516. for name, props in self.extra_regressors.items():
  517. seasonal_features.append(pd.DataFrame(df[name]))
  518. prior_scales.append(props['prior_scale'])
  519. if len(seasonal_features) == 0:
  520. seasonal_features.append(
  521. pd.DataFrame({'zeros': np.zeros(df.shape[0])}))
  522. prior_scales.append(1.)
  523. return pd.concat(seasonal_features, axis=1), prior_scales
  524. def parse_seasonality_args(self, name, arg, auto_disable, default_order):
  525. """Get number of fourier components for built-in seasonalities.
  526. Parameters
  527. ----------
  528. name: string name of the seasonality component.
  529. arg: 'auto', True, False, or number of fourier components as provided.
  530. auto_disable: bool if seasonality should be disabled when 'auto'.
  531. default_order: int default fourier order
  532. Returns
  533. -------
  534. Number of fourier components, or 0 for disabled.
  535. """
  536. if arg == 'auto':
  537. fourier_order = 0
  538. if name in self.seasonalities:
  539. logger.info(
  540. 'Found custom seasonality named "{name}", '
  541. 'disabling built-in {name} seasonality.'.format(name=name)
  542. )
  543. elif auto_disable:
  544. logger.info(
  545. 'Disabling {name} seasonality. Run prophet with '
  546. '{name}_seasonality=True to override this.'.format(
  547. name=name)
  548. )
  549. else:
  550. fourier_order = default_order
  551. elif arg is True:
  552. fourier_order = default_order
  553. elif arg is False:
  554. fourier_order = 0
  555. else:
  556. fourier_order = int(arg)
  557. return fourier_order
  558. def set_auto_seasonalities(self):
  559. """Set seasonalities that were left on auto.
  560. Turns on yearly seasonality if there is >=2 years of history.
  561. Turns on weekly seasonality if there is >=2 weeks of history, and the
  562. spacing between dates in the history is <7 days.
  563. Turns on daily seasonality if there is >=2 days of history, and the
  564. spacing between dates in the history is <1 day.
  565. """
  566. first = self.history['ds'].min()
  567. last = self.history['ds'].max()
  568. dt = self.history['ds'].diff()
  569. min_dt = dt.iloc[dt.nonzero()[0]].min()
  570. # Yearly seasonality
  571. yearly_disable = last - first < pd.Timedelta(days=730)
  572. fourier_order = self.parse_seasonality_args(
  573. 'yearly', self.yearly_seasonality, yearly_disable, 10)
  574. if fourier_order > 0:
  575. self.seasonalities['yearly'] = {
  576. 'period': 365.25,
  577. 'fourier_order': fourier_order,
  578. 'prior_scale': self.seasonality_prior_scale,
  579. }
  580. # Weekly seasonality
  581. weekly_disable = ((last - first < pd.Timedelta(weeks=2)) or
  582. (min_dt >= pd.Timedelta(weeks=1)))
  583. fourier_order = self.parse_seasonality_args(
  584. 'weekly', self.weekly_seasonality, weekly_disable, 3)
  585. if fourier_order > 0:
  586. self.seasonalities['weekly'] = {
  587. 'period': 7,
  588. 'fourier_order': fourier_order,
  589. 'prior_scale': self.seasonality_prior_scale,
  590. }
  591. # Daily seasonality
  592. daily_disable = ((last - first < pd.Timedelta(days=2)) or
  593. (min_dt >= pd.Timedelta(days=1)))
  594. fourier_order = self.parse_seasonality_args(
  595. 'daily', self.daily_seasonality, daily_disable, 4)
  596. if fourier_order > 0:
  597. self.seasonalities['daily'] = {
  598. 'period': 1,
  599. 'fourier_order': fourier_order,
  600. 'prior_scale': self.seasonality_prior_scale,
  601. }
  602. @staticmethod
  603. def linear_growth_init(df):
  604. """Initialize linear growth.
  605. Provides a strong initialization for linear growth by calculating the
  606. growth and offset parameters that pass the function through the first
  607. and last points in the time series.
  608. Parameters
  609. ----------
  610. df: pd.DataFrame with columns ds (date), y_scaled (scaled time series),
  611. and t (scaled time).
  612. Returns
  613. -------
  614. A tuple (k, m) with the rate (k) and offset (m) of the linear growth
  615. function.
  616. """
  617. i0, i1 = df['ds'].idxmin(), df['ds'].idxmax()
  618. T = df['t'].iloc[i1] - df['t'].iloc[i0]
  619. k = (df['y_scaled'].iloc[i1] - df['y_scaled'].iloc[i0]) / T
  620. m = df['y_scaled'].iloc[i0] - k * df['t'].iloc[i0]
  621. return (k, m)
  622. @staticmethod
  623. def logistic_growth_init(df):
  624. """Initialize logistic growth.
  625. Provides a strong initialization for logistic growth by calculating the
  626. growth and offset parameters that pass the function through the first
  627. and last points in the time series.
  628. Parameters
  629. ----------
  630. df: pd.DataFrame with columns ds (date), cap_scaled (scaled capacity),
  631. y_scaled (scaled time series), and t (scaled time).
  632. Returns
  633. -------
  634. A tuple (k, m) with the rate (k) and offset (m) of the logistic growth
  635. function.
  636. """
  637. i0, i1 = df['ds'].idxmin(), df['ds'].idxmax()
  638. T = df['t'].iloc[i1] - df['t'].iloc[i0]
  639. # Force valid values, in case y > cap or y < 0
  640. C0 = df['cap_scaled'].iloc[i0]
  641. C1 = df['cap_scaled'].iloc[i1]
  642. y0 = max(0.01 * C0, min(0.99 * C0, df['y_scaled'].iloc[i0]))
  643. y1 = max(0.01 * C1, min(0.99 * C1, df['y_scaled'].iloc[i1]))
  644. r0 = C0 / y0
  645. r1 = C1 / y1
  646. if abs(r0 - r1) <= 0.01:
  647. r0 = 1.05 * r0
  648. L0 = np.log(r0 - 1)
  649. L1 = np.log(r1 - 1)
  650. # Initialize the offset
  651. m = L0 * T / (L0 - L1)
  652. # And the rate
  653. k = (L0 - L1) / T
  654. return (k, m)
  655. # fb-block 7
  656. def fit(self, df, **kwargs):
  657. """Fit the Prophet model.
  658. This sets self.params to contain the fitted model parameters. It is a
  659. dictionary parameter names as keys and the following items:
  660. k (Mx1 array): M posterior samples of the initial slope.
  661. m (Mx1 array): The initial intercept.
  662. delta (MxN array): The slope change at each of N changepoints.
  663. beta (MxK matrix): Coefficients for K seasonality features.
  664. sigma_obs (Mx1 array): Noise level.
  665. Note that M=1 if MAP estimation.
  666. Parameters
  667. ----------
  668. df: pd.DataFrame containing the history. Must have columns ds (date
  669. type) and y, the time series. If self.growth is 'logistic', then
  670. df must also have a column cap that specifies the capacity at
  671. each ds.
  672. kwargs: Additional arguments passed to the optimizing or sampling
  673. functions in Stan.
  674. Returns
  675. -------
  676. The fitted Prophet object.
  677. """
  678. if self.history is not None:
  679. raise Exception('Prophet object can only be fit once. '
  680. 'Instantiate a new object.')
  681. history = df[df['y'].notnull()].copy()
  682. if history.shape[0] < 2:
  683. raise ValueError('Dataframe has less than 2 non-NaN rows.')
  684. self.history_dates = pd.to_datetime(df['ds']).sort_values()
  685. history = self.setup_dataframe(history, initialize_scales=True)
  686. self.history = history
  687. self.set_auto_seasonalities()
  688. seasonal_features, prior_scales = (
  689. self.make_all_seasonality_features(history))
  690. self.set_changepoints()
  691. A = self.get_changepoint_matrix()
  692. dat = {
  693. 'T': history.shape[0],
  694. 'K': seasonal_features.shape[1],
  695. 'S': len(self.changepoints_t),
  696. 'y': history['y_scaled'],
  697. 't': history['t'],
  698. 'A': A,
  699. 't_change': self.changepoints_t,
  700. 'X': seasonal_features,
  701. 'sigmas': prior_scales,
  702. 'tau': self.changepoint_prior_scale,
  703. }
  704. if self.growth == 'linear':
  705. kinit = self.linear_growth_init(history)
  706. else:
  707. dat['cap'] = history['cap_scaled']
  708. kinit = self.logistic_growth_init(history)
  709. model = prophet_stan_models[self.growth]
  710. def stan_init():
  711. return {
  712. 'k': kinit[0],
  713. 'm': kinit[1],
  714. 'delta': np.zeros(len(self.changepoints_t)),
  715. 'beta': np.zeros(seasonal_features.shape[1]),
  716. 'sigma_obs': 1,
  717. }
  718. if history['y'].min() == history['y'].max():
  719. # Nothing to fit.
  720. self.params = stan_init()
  721. self.params['sigma_obs'] = 0.
  722. for par in self.params:
  723. self.params[par] = np.array([self.params[par]])
  724. elif self.mcmc_samples > 0:
  725. stan_fit = model.sampling(
  726. dat,
  727. init=stan_init,
  728. iter=self.mcmc_samples,
  729. **kwargs
  730. )
  731. for par in stan_fit.model_pars:
  732. self.params[par] = stan_fit[par]
  733. else:
  734. try:
  735. params = model.optimizing(
  736. dat, init=stan_init, iter=1e4, **kwargs)
  737. except RuntimeError:
  738. params = model.optimizing(
  739. dat, init=stan_init, iter=1e4, algorithm='Newton',
  740. **kwargs
  741. )
  742. for par in params:
  743. self.params[par] = params[par].reshape((1, -1))
  744. # If no changepoints were requested, replace delta with 0s
  745. if len(self.changepoints) == 0:
  746. # Fold delta into the base rate k
  747. self.params['k'] = self.params['k'] + self.params['delta']
  748. self.params['delta'] = np.zeros(self.params['delta'].shape)
  749. return self
  750. # fb-block 8
  751. def predict(self, df=None):
  752. """Predict using the prophet model.
  753. Parameters
  754. ----------
  755. df: pd.DataFrame with dates for predictions (column ds), and capacity
  756. (column cap) if logistic growth. If not provided, predictions are
  757. made on the history.
  758. Returns
  759. -------
  760. A pd.DataFrame with the forecast components.
  761. """
  762. if df is None:
  763. df = self.history.copy()
  764. else:
  765. if df.shape[0] == 0:
  766. raise ValueError('Dataframe has no rows.')
  767. df = self.setup_dataframe(df.copy())
  768. df['trend'] = self.predict_trend(df)
  769. seasonal_components = self.predict_seasonal_components(df)
  770. intervals = self.predict_uncertainty(df)
  771. # Drop columns except ds, cap, floor, and trend
  772. cols = ['ds', 'trend']
  773. if 'cap' in df:
  774. cols.append('cap')
  775. if self.logistic_floor:
  776. cols.append('floor')
  777. # Add in forecast components
  778. df2 = pd.concat((df[cols], intervals, seasonal_components), axis=1)
  779. df2['yhat'] = df2['trend'] + df2['seasonal']
  780. return df2
  781. @staticmethod
  782. def piecewise_linear(t, deltas, k, m, changepoint_ts):
  783. """Evaluate the piecewise linear function.
  784. Parameters
  785. ----------
  786. t: np.array of times on which the function is evaluated.
  787. deltas: np.array of rate changes at each changepoint.
  788. k: Float initial rate.
  789. m: Float initial offset.
  790. changepoint_ts: np.array of changepoint times.
  791. Returns
  792. -------
  793. Vector y(t).
  794. """
  795. # Intercept changes
  796. gammas = -changepoint_ts * deltas
  797. # Get cumulative slope and intercept at each t
  798. k_t = k * np.ones_like(t)
  799. m_t = m * np.ones_like(t)
  800. for s, t_s in enumerate(changepoint_ts):
  801. indx = t >= t_s
  802. k_t[indx] += deltas[s]
  803. m_t[indx] += gammas[s]
  804. return k_t * t + m_t
  805. @staticmethod
  806. def piecewise_logistic(t, cap, deltas, k, m, changepoint_ts):
  807. """Evaluate the piecewise logistic function.
  808. Parameters
  809. ----------
  810. t: np.array of times on which the function is evaluated.
  811. cap: np.array of capacities at each t.
  812. deltas: np.array of rate changes at each changepoint.
  813. k: Float initial rate.
  814. m: Float initial offset.
  815. changepoint_ts: np.array of changepoint times.
  816. Returns
  817. -------
  818. Vector y(t).
  819. """
  820. # Compute offset changes
  821. k_cum = np.concatenate((np.atleast_1d(k), np.cumsum(deltas) + k))
  822. gammas = np.zeros(len(changepoint_ts))
  823. for i, t_s in enumerate(changepoint_ts):
  824. gammas[i] = (
  825. (t_s - m - np.sum(gammas))
  826. * (1 - k_cum[i] / k_cum[i + 1]) # noqa W503
  827. )
  828. # Get cumulative rate and offset at each t
  829. k_t = k * np.ones_like(t)
  830. m_t = m * np.ones_like(t)
  831. for s, t_s in enumerate(changepoint_ts):
  832. indx = t >= t_s
  833. k_t[indx] += deltas[s]
  834. m_t[indx] += gammas[s]
  835. return cap / (1 + np.exp(-k_t * (t - m_t)))
  836. def predict_trend(self, df):
  837. """Predict trend using the prophet model.
  838. Parameters
  839. ----------
  840. df: Prediction dataframe.
  841. Returns
  842. -------
  843. Vector with trend on prediction dates.
  844. """
  845. k = np.nanmean(self.params['k'])
  846. m = np.nanmean(self.params['m'])
  847. deltas = np.nanmean(self.params['delta'], axis=0)
  848. t = np.array(df['t'])
  849. if self.growth == 'linear':
  850. trend = self.piecewise_linear(t, deltas, k, m, self.changepoints_t)
  851. else:
  852. cap = df['cap_scaled']
  853. trend = self.piecewise_logistic(
  854. t, cap, deltas, k, m, self.changepoints_t)
  855. return trend * self.y_scale + df['floor']
  856. def predict_seasonal_components(self, df):
  857. """Predict seasonality components, holidays, and added regressors.
  858. Parameters
  859. ----------
  860. df: Prediction dataframe.
  861. Returns
  862. -------
  863. Dataframe with seasonal components.
  864. """
  865. seasonal_features, _ = self.make_all_seasonality_features(df)
  866. lower_p = 100 * (1.0 - self.interval_width) / 2
  867. upper_p = 100 * (1.0 + self.interval_width) / 2
  868. components = pd.DataFrame({
  869. 'col': np.arange(seasonal_features.shape[1]),
  870. 'component': [x.split('_delim_')[0] for x in seasonal_features.columns],
  871. })
  872. # Add total for all regression components
  873. components = components.append(pd.DataFrame({
  874. 'col': np.arange(seasonal_features.shape[1]),
  875. 'component': 'seasonal',
  876. }))
  877. # Add totals for seasonality, holiday, and extra regressors
  878. components = self.add_group_component(
  879. components, 'seasonalities', self.seasonalities.keys())
  880. if self.holidays is not None:
  881. components = self.add_group_component(
  882. components, 'holidays', self.holidays['holiday'].unique())
  883. components = self.add_group_component(
  884. components, 'extra_regressors', self.extra_regressors.keys())
  885. # Remove the placeholder
  886. components = components[components['component'] != 'zeros']
  887. X = seasonal_features.as_matrix()
  888. data = {}
  889. for component, features in components.groupby('component'):
  890. cols = features.col.tolist()
  891. comp_beta = self.params['beta'][:, cols]
  892. comp_features = X[:, cols]
  893. comp = (
  894. np.matmul(comp_features, comp_beta.transpose())
  895. * self.y_scale # noqa W503
  896. )
  897. data[component] = np.nanmean(comp, axis=1)
  898. data[component + '_lower'] = np.nanpercentile(comp, lower_p,
  899. axis=1)
  900. data[component + '_upper'] = np.nanpercentile(comp, upper_p,
  901. axis=1)
  902. return pd.DataFrame(data)
  903. def add_group_component(self, components, name, group):
  904. """Adds a component with given name that contains all of the components
  905. in group.
  906. Parameters
  907. ----------
  908. components: Dataframe with components.
  909. name: Name of new group component.
  910. group: List of components that form the group.
  911. Returns
  912. -------
  913. Dataframe with components.
  914. """
  915. new_comp = components[components['component'].isin(set(group))].copy()
  916. new_comp['component'] = name
  917. components = components.append(new_comp)
  918. return components
  919. def sample_posterior_predictive(self, df):
  920. """Prophet posterior predictive samples.
  921. Parameters
  922. ----------
  923. df: Prediction dataframe.
  924. Returns
  925. -------
  926. Dictionary with posterior predictive samples for each component.
  927. """
  928. n_iterations = self.params['k'].shape[0]
  929. samp_per_iter = max(1, int(np.ceil(
  930. self.uncertainty_samples / float(n_iterations)
  931. )))
  932. # Generate seasonality features once so we can re-use them.
  933. seasonal_features, _ = self.make_all_seasonality_features(df)
  934. sim_values = {'yhat': [], 'trend': [], 'seasonal': []}
  935. for i in range(n_iterations):
  936. for _j in range(samp_per_iter):
  937. sim = self.sample_model(df, seasonal_features, i)
  938. for key in sim_values:
  939. sim_values[key].append(sim[key])
  940. for k, v in sim_values.items():
  941. sim_values[k] = np.column_stack(v)
  942. return sim_values
  943. def predictive_samples(self, df):
  944. """Sample from the posterior predictive distribution.
  945. Parameters
  946. ----------
  947. df: Dataframe with dates for predictions (column ds), and capacity
  948. (column cap) if logistic growth.
  949. Returns
  950. -------
  951. Dictionary with keys "trend", "seasonal", and "yhat" containing
  952. posterior predictive samples for that component. "seasonal" is the sum
  953. of seasonalities, holidays, and added regressors.
  954. """
  955. df = self.setup_dataframe(df.copy())
  956. sim_values = self.sample_posterior_predictive(df)
  957. return sim_values
  958. def predict_uncertainty(self, df):
  959. """Prediction intervals for yhat and trend.
  960. Parameters
  961. ----------
  962. df: Prediction dataframe.
  963. Returns
  964. -------
  965. Dataframe with uncertainty intervals.
  966. """
  967. sim_values = self.sample_posterior_predictive(df)
  968. lower_p = 100 * (1.0 - self.interval_width) / 2
  969. upper_p = 100 * (1.0 + self.interval_width) / 2
  970. series = {}
  971. for key in ['yhat', 'trend']:
  972. series['{}_lower'.format(key)] = np.nanpercentile(
  973. sim_values[key], lower_p, axis=1)
  974. series['{}_upper'.format(key)] = np.nanpercentile(
  975. sim_values[key], upper_p, axis=1)
  976. return pd.DataFrame(series)
  977. def sample_model(self, df, seasonal_features, iteration):
  978. """Simulate observations from the extrapolated generative model.
  979. Parameters
  980. ----------
  981. df: Prediction dataframe.
  982. seasonal_features: pd.DataFrame of seasonal features.
  983. iteration: Int sampling iteration to use parameters from.
  984. Returns
  985. -------
  986. Dataframe with trend, seasonality, and yhat, each like df['t'].
  987. """
  988. trend = self.sample_predictive_trend(df, iteration)
  989. beta = self.params['beta'][iteration]
  990. seasonal = np.matmul(seasonal_features.as_matrix(), beta) * self.y_scale
  991. sigma = self.params['sigma_obs'][iteration]
  992. noise = np.random.normal(0, sigma, df.shape[0]) * self.y_scale
  993. return pd.DataFrame({
  994. 'yhat': trend + seasonal + noise,
  995. 'trend': trend,
  996. 'seasonal': seasonal,
  997. })
  998. def sample_predictive_trend(self, df, iteration):
  999. """Simulate the trend using the extrapolated generative model.
  1000. Parameters
  1001. ----------
  1002. df: Prediction dataframe.
  1003. iteration: Int sampling iteration to use parameters from.
  1004. Returns
  1005. -------
  1006. np.array of simulated trend over df['t'].
  1007. """
  1008. k = self.params['k'][iteration]
  1009. m = self.params['m'][iteration]
  1010. deltas = self.params['delta'][iteration]
  1011. t = np.array(df['t'])
  1012. T = t.max()
  1013. if T > 1:
  1014. # Get the time discretization of the history
  1015. dt = np.diff(self.history['t'])
  1016. dt = np.min(dt[dt > 0])
  1017. # Number of time periods in the future
  1018. N = np.ceil((T - 1) / float(dt))
  1019. S = len(self.changepoints_t)
  1020. prob_change = min(1, (S * (T - 1)) / N)
  1021. n_changes = np.random.binomial(N, prob_change)
  1022. # Sample ts
  1023. changepoint_ts_new = sorted(np.random.uniform(1, T, n_changes))
  1024. else:
  1025. # Case where we're not extrapolating.
  1026. changepoint_ts_new = []
  1027. n_changes = 0
  1028. # Get the empirical scale of the deltas, plus epsilon to avoid NaNs.
  1029. lambda_ = np.mean(np.abs(deltas)) + 1e-8
  1030. # Sample deltas
  1031. deltas_new = np.random.laplace(0, lambda_, n_changes)
  1032. # Prepend the times and deltas from the history
  1033. changepoint_ts = np.concatenate((self.changepoints_t,
  1034. changepoint_ts_new))
  1035. deltas = np.concatenate((deltas, deltas_new))
  1036. if self.growth == 'linear':
  1037. trend = self.piecewise_linear(t, deltas, k, m, changepoint_ts)
  1038. else:
  1039. cap = df['cap_scaled']
  1040. trend = self.piecewise_logistic(t, cap, deltas, k, m,
  1041. changepoint_ts)
  1042. return trend * self.y_scale + df['floor']
  1043. def make_future_dataframe(self, periods, freq='D', include_history=True):
  1044. """Simulate the trend using the extrapolated generative model.
  1045. Parameters
  1046. ----------
  1047. periods: Int number of periods to forecast forward.
  1048. freq: Any valid frequency for pd.date_range, such as 'D' or 'M'.
  1049. include_history: Boolean to include the historical dates in the data
  1050. frame for predictions.
  1051. Returns
  1052. -------
  1053. pd.Dataframe that extends forward from the end of self.history for the
  1054. requested number of periods.
  1055. """
  1056. last_date = self.history_dates.max()
  1057. dates = pd.date_range(
  1058. start=last_date,
  1059. periods=periods + 1, # An extra in case we include start
  1060. freq=freq)
  1061. dates = dates[dates > last_date] # Drop start if equals last_date
  1062. dates = dates[:periods] # Return correct number of periods
  1063. if include_history:
  1064. dates = np.concatenate((np.array(self.history_dates), dates))
  1065. return pd.DataFrame({'ds': dates})
  1066. def plot(self, fcst, ax=None, uncertainty=True, plot_cap=True, xlabel='ds',
  1067. ylabel='y'):
  1068. """Plot the Prophet forecast.
  1069. Parameters
  1070. ----------
  1071. fcst: pd.DataFrame output of self.predict.
  1072. ax: Optional matplotlib axes on which to plot.
  1073. uncertainty: Optional boolean to plot uncertainty intervals.
  1074. plot_cap: Optional boolean indicating if the capacity should be shown
  1075. in the figure, if available.
  1076. xlabel: Optional label name on X-axis
  1077. ylabel: Optional label name on Y-axis
  1078. Returns
  1079. -------
  1080. A matplotlib figure.
  1081. """
  1082. if ax is None:
  1083. fig = plt.figure(facecolor='w', figsize=(10, 6))
  1084. ax = fig.add_subplot(111)
  1085. else:
  1086. fig = ax.get_figure()
  1087. ax.plot(self.history['ds'].values, self.history['y'], 'k.')
  1088. ax.plot(fcst['ds'].values, fcst['yhat'], ls='-', c='#0072B2')
  1089. if 'cap' in fcst and plot_cap:
  1090. ax.plot(fcst['ds'].values, fcst['cap'], ls='--', c='k')
  1091. if self.logistic_floor and 'floor' in fcst and plot_cap:
  1092. ax.plot(fcst['ds'].values, fcst['floor'], ls='--', c='k')
  1093. if uncertainty:
  1094. ax.fill_between(fcst['ds'].values, fcst['yhat_lower'],
  1095. fcst['yhat_upper'], color='#0072B2',
  1096. alpha=0.2)
  1097. ax.grid(True, which='major', c='gray', ls='-', lw=1, alpha=0.2)
  1098. ax.set_xlabel(xlabel)
  1099. ax.set_ylabel(ylabel)
  1100. fig.tight_layout()
  1101. return fig
  1102. def plot_components(self, fcst, uncertainty=True, plot_cap=True,
  1103. weekly_start=0, yearly_start=0):
  1104. """Plot the Prophet forecast components.
  1105. Will plot whichever are available of: trend, holidays, weekly
  1106. seasonality, and yearly seasonality.
  1107. Parameters
  1108. ----------
  1109. fcst: pd.DataFrame output of self.predict.
  1110. uncertainty: Optional boolean to plot uncertainty intervals.
  1111. plot_cap: Optional boolean indicating if the capacity should be shown
  1112. in the figure, if available.
  1113. weekly_start: Optional int specifying the start day of the weekly
  1114. seasonality plot. 0 (default) starts the week on Sunday. 1 shifts
  1115. by 1 day to Monday, and so on.
  1116. yearly_start: Optional int specifying the start day of the yearly
  1117. seasonality plot. 0 (default) starts the year on Jan 1. 1 shifts
  1118. by 1 day to Jan 2, and so on.
  1119. Returns
  1120. -------
  1121. A matplotlib figure.
  1122. """
  1123. # Identify components to be plotted
  1124. components = ['trend']
  1125. if self.holidays is not None and 'holidays' in fcst:
  1126. components.append('holidays')
  1127. components.extend([name for name in self.seasonalities
  1128. if name in fcst])
  1129. if len(self.extra_regressors) > 0 and 'extra_regressors' in fcst:
  1130. components.append('extra_regressors')
  1131. npanel = len(components)
  1132. fig, axes = plt.subplots(npanel, 1, facecolor='w',
  1133. figsize=(9, 3 * npanel))
  1134. for ax, plot in zip(axes, components):
  1135. if plot == 'trend':
  1136. self.plot_forecast_component(
  1137. fcst, 'trend', ax, uncertainty, plot_cap)
  1138. elif plot == 'holidays':
  1139. self.plot_forecast_component(
  1140. fcst, 'holidays', ax, uncertainty, False)
  1141. elif plot == 'weekly':
  1142. self.plot_weekly(
  1143. ax=ax, uncertainty=uncertainty, weekly_start=weekly_start)
  1144. elif plot == 'yearly':
  1145. self.plot_yearly(
  1146. ax=ax, uncertainty=uncertainty, yearly_start=yearly_start)
  1147. elif plot == 'extra_regressors':
  1148. self.plot_forecast_component(
  1149. fcst, 'extra_regressors', ax, uncertainty, False)
  1150. else:
  1151. self.plot_seasonality(
  1152. name=plot, ax=ax, uncertainty=uncertainty)
  1153. fig.tight_layout()
  1154. return fig
  1155. def plot_forecast_component(
  1156. self, fcst, name, ax=None, uncertainty=True, plot_cap=False):
  1157. """Plot a particular component of the forecast.
  1158. Parameters
  1159. ----------
  1160. fcst: pd.DataFrame output of self.predict.
  1161. name: Name of the component to plot.
  1162. ax: Optional matplotlib Axes to plot on.
  1163. uncertainty: Optional boolean to plot uncertainty intervals.
  1164. plot_cap: Optional boolean indicating if the capacity should be shown
  1165. in the figure, if available.
  1166. Returns
  1167. -------
  1168. a list of matplotlib artists
  1169. """
  1170. artists = []
  1171. if not ax:
  1172. fig = plt.figure(facecolor='w', figsize=(10, 6))
  1173. ax = fig.add_subplot(111)
  1174. artists += ax.plot(fcst['ds'].values, fcst[name], ls='-', c='#0072B2')
  1175. if 'cap' in fcst and plot_cap:
  1176. artists += ax.plot(fcst['ds'].values, fcst['cap'], ls='--', c='k')
  1177. if self.logistic_floor and 'floor' in fcst and plot_cap:
  1178. ax.plot(fcst['ds'].values, fcst['floor'], ls='--', c='k')
  1179. if uncertainty:
  1180. artists += [ax.fill_between(
  1181. fcst['ds'].values, fcst[name + '_lower'],
  1182. fcst[name + '_upper'], color='#0072B2', alpha=0.2)]
  1183. ax.grid(True, which='major', c='gray', ls='-', lw=1, alpha=0.2)
  1184. ax.set_xlabel('ds')
  1185. ax.set_ylabel(name)
  1186. return artists
  1187. def seasonality_plot_df(self, ds):
  1188. """Prepare dataframe for plotting seasonal components.
  1189. Parameters
  1190. ----------
  1191. ds: List of dates for column ds.
  1192. Returns
  1193. -------
  1194. A dataframe with seasonal components on ds.
  1195. """
  1196. df_dict = {'ds': ds, 'cap': 1.}
  1197. for name in self.extra_regressors:
  1198. df_dict[name] = 0.
  1199. df = pd.DataFrame(df_dict)
  1200. df = self.setup_dataframe(df)
  1201. return df
  1202. def plot_weekly(self, ax=None, uncertainty=True, weekly_start=0):
  1203. """Plot the weekly component of the forecast.
  1204. Parameters
  1205. ----------
  1206. ax: Optional matplotlib Axes to plot on. One will be created if this
  1207. is not provided.
  1208. uncertainty: Optional boolean to plot uncertainty intervals.
  1209. weekly_start: Optional int specifying the start day of the weekly
  1210. seasonality plot. 0 (default) starts the week on Sunday. 1 shifts
  1211. by 1 day to Monday, and so on.
  1212. Returns
  1213. -------
  1214. a list of matplotlib artists
  1215. """
  1216. artists = []
  1217. if not ax:
  1218. fig = plt.figure(facecolor='w', figsize=(10, 6))
  1219. ax = fig.add_subplot(111)
  1220. # Compute weekly seasonality for a Sun-Sat sequence of dates.
  1221. days = (pd.date_range(start='2017-01-01', periods=7) +
  1222. pd.Timedelta(days=weekly_start))
  1223. df_w = self.seasonality_plot_df(days)
  1224. seas = self.predict_seasonal_components(df_w)
  1225. days = days.weekday_name
  1226. artists += ax.plot(range(len(days)), seas['weekly'], ls='-',
  1227. c='#0072B2')
  1228. if uncertainty:
  1229. artists += [ax.fill_between(range(len(days)),
  1230. seas['weekly_lower'], seas['weekly_upper'],
  1231. color='#0072B2', alpha=0.2)]
  1232. ax.grid(True, which='major', c='gray', ls='-', lw=1, alpha=0.2)
  1233. ax.set_xticks(range(len(days)))
  1234. ax.set_xticklabels(days)
  1235. ax.set_xlabel('Day of week')
  1236. ax.set_ylabel('weekly')
  1237. return artists
  1238. def plot_yearly(self, ax=None, uncertainty=True, yearly_start=0):
  1239. """Plot the yearly component of the forecast.
  1240. Parameters
  1241. ----------
  1242. ax: Optional matplotlib Axes to plot on. One will be created if
  1243. this is not provided.
  1244. uncertainty: Optional boolean to plot uncertainty intervals.
  1245. yearly_start: Optional int specifying the start day of the yearly
  1246. seasonality plot. 0 (default) starts the year on Jan 1. 1 shifts
  1247. by 1 day to Jan 2, and so on.
  1248. Returns
  1249. -------
  1250. a list of matplotlib artists
  1251. """
  1252. artists = []
  1253. if not ax:
  1254. fig = plt.figure(facecolor='w', figsize=(10, 6))
  1255. ax = fig.add_subplot(111)
  1256. # Compute yearly seasonality for a Jan 1 - Dec 31 sequence of dates.
  1257. days = (pd.date_range(start='2017-01-01', periods=365) +
  1258. pd.Timedelta(days=yearly_start))
  1259. df_y = self.seasonality_plot_df(days)
  1260. seas = self.predict_seasonal_components(df_y)
  1261. artists += ax.plot(df_y['ds'], seas['yearly'], ls='-',
  1262. c='#0072B2')
  1263. if uncertainty:
  1264. artists += [ax.fill_between(
  1265. df_y['ds'].values, seas['yearly_lower'],
  1266. seas['yearly_upper'], color='#0072B2', alpha=0.2)]
  1267. ax.grid(True, which='major', c='gray', ls='-', lw=1, alpha=0.2)
  1268. months = MonthLocator(range(1, 13), bymonthday=1, interval=2)
  1269. ax.xaxis.set_major_formatter(FuncFormatter(
  1270. lambda x, pos=None: '{dt:%B} {dt.day}'.format(dt=num2date(x))))
  1271. ax.xaxis.set_major_locator(months)
  1272. ax.set_xlabel('Day of year')
  1273. ax.set_ylabel('yearly')
  1274. return artists
  1275. def plot_seasonality(self, name, ax=None, uncertainty=True):
  1276. """Plot a custom seasonal component.
  1277. Parameters
  1278. ----------
  1279. ax: Optional matplotlib Axes to plot on. One will be created if
  1280. this is not provided.
  1281. uncertainty: Optional boolean to plot uncertainty intervals.
  1282. Returns
  1283. -------
  1284. a list of matplotlib artists
  1285. """
  1286. artists = []
  1287. if not ax:
  1288. fig = plt.figure(facecolor='w', figsize=(10, 6))
  1289. ax = fig.add_subplot(111)
  1290. # Compute seasonality from Jan 1 through a single period.
  1291. start = pd.to_datetime('2017-01-01 0000')
  1292. period = self.seasonalities[name]['period']
  1293. end = start + pd.Timedelta(days=period)
  1294. plot_points = 200
  1295. days = pd.to_datetime(np.linspace(start.value, end.value, plot_points))
  1296. df_y = self.seasonality_plot_df(days)
  1297. seas = self.predict_seasonal_components(df_y)
  1298. artists += ax.plot(df_y['ds'], seas[name], ls='-',
  1299. c='#0072B2')
  1300. if uncertainty:
  1301. artists += [ax.fill_between(
  1302. df_y['ds'].values, seas[name + '_lower'],
  1303. seas[name + '_upper'], color='#0072B2', alpha=0.2)]
  1304. ax.grid(True, which='major', c='gray', ls='-', lw=1, alpha=0.2)
  1305. ax.set_xticks(pd.to_datetime(np.linspace(start.value, end.value, 7)))
  1306. if period <= 2:
  1307. fmt_str = '{dt:%T}'
  1308. elif period < 14:
  1309. fmt_str = '{dt:%m}/{dt:%d} {dt:%R}'
  1310. else:
  1311. fmt_str = '{dt:%m}/{dt:%d}'
  1312. ax.xaxis.set_major_formatter(FuncFormatter(
  1313. lambda x, pos=None: fmt_str.format(dt=num2date(x))))
  1314. ax.set_xlabel('ds')
  1315. ax.set_ylabel(name)
  1316. return artists
  1317. def copy(self, cutoff=None):
  1318. """Copy Prophet object
  1319. Parameters
  1320. ----------
  1321. cutoff: pd.Timestamp or None, default None.
  1322. cuttoff Timestamp for changepoints member variable.
  1323. changepoints are only retained if 'changepoints <= cutoff'
  1324. Returns
  1325. -------
  1326. Prophet class object with the same parameter with model variable
  1327. """
  1328. if self.specified_changepoints:
  1329. changepoints = self.changepoints
  1330. if cutoff is not None:
  1331. # Filter change points '<= cutoff'
  1332. changepoints = changepoints[changepoints <= cutoff]
  1333. else:
  1334. changepoints = None
  1335. return Prophet(
  1336. growth=self.growth,
  1337. n_changepoints=self.n_changepoints,
  1338. changepoints=changepoints,
  1339. yearly_seasonality=self.yearly_seasonality,
  1340. weekly_seasonality=self.weekly_seasonality,
  1341. daily_seasonality=self.daily_seasonality,
  1342. holidays=self.holidays,
  1343. seasonality_prior_scale=self.seasonality_prior_scale,
  1344. changepoint_prior_scale=self.changepoint_prior_scale,
  1345. holidays_prior_scale=self.holidays_prior_scale,
  1346. mcmc_samples=self.mcmc_samples,
  1347. interval_width=self.interval_width,
  1348. uncertainty_samples=self.uncertainty_samples
  1349. )