plot.py 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467
  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. import logging
  12. import numpy as np
  13. import pandas as pd
  14. from fbprophet.diagnostics import performance_metrics
  15. logging.basicConfig()
  16. logger = logging.getLogger(__name__)
  17. try:
  18. from matplotlib import pyplot as plt
  19. from matplotlib.dates import MonthLocator, num2date
  20. from matplotlib.ticker import FuncFormatter
  21. except ImportError:
  22. logger.error('Importing matplotlib failed. Plotting will not work.')
  23. def plot(
  24. m, fcst, ax=None, uncertainty=True, plot_cap=True, xlabel='ds', ylabel='y',
  25. ):
  26. """Plot the Prophet forecast.
  27. Parameters
  28. ----------
  29. m: Prophet model.
  30. fcst: pd.DataFrame output of m.predict.
  31. ax: Optional matplotlib axes on which to plot.
  32. uncertainty: Optional boolean to plot uncertainty intervals.
  33. plot_cap: Optional boolean indicating if the capacity should be shown
  34. in the figure, if available.
  35. xlabel: Optional label name on X-axis
  36. ylabel: Optional label name on Y-axis
  37. Returns
  38. -------
  39. A matplotlib figure.
  40. """
  41. if ax is None:
  42. fig = plt.figure(facecolor='w', figsize=(10, 6))
  43. ax = fig.add_subplot(111)
  44. else:
  45. fig = ax.get_figure()
  46. fcst_t = fcst['ds'].dt.to_pydatetime()
  47. ax.plot(m.history['ds'].dt.to_pydatetime(), m.history['y'], 'k.')
  48. ax.plot(fcst_t, fcst['yhat'], ls='-', c='#0072B2')
  49. if 'cap' in fcst and plot_cap:
  50. ax.plot(fcst_t, fcst['cap'], ls='--', c='k')
  51. if m.logistic_floor and 'floor' in fcst and plot_cap:
  52. ax.plot(fcst_t, fcst['floor'], ls='--', c='k')
  53. if uncertainty:
  54. ax.fill_between(fcst_t, fcst['yhat_lower'], fcst['yhat_upper'],
  55. color='#0072B2', alpha=0.2)
  56. ax.grid(True, which='major', c='gray', ls='-', lw=1, alpha=0.2)
  57. ax.set_xlabel(xlabel)
  58. ax.set_ylabel(ylabel)
  59. fig.tight_layout()
  60. return fig
  61. def plot_components(
  62. m, fcst, uncertainty=True, plot_cap=True, weekly_start=0, yearly_start=0,
  63. ):
  64. """Plot the Prophet forecast components.
  65. Will plot whichever are available of: trend, holidays, weekly
  66. seasonality, yearly seasonality, and additive and multiplicative extra
  67. regressors.
  68. Parameters
  69. ----------
  70. m: Prophet model.
  71. fcst: pd.DataFrame output of m.predict.
  72. uncertainty: Optional boolean to plot uncertainty intervals.
  73. plot_cap: Optional boolean indicating if the capacity should be shown
  74. in the figure, if available.
  75. weekly_start: Optional int specifying the start day of the weekly
  76. seasonality plot. 0 (default) starts the week on Sunday. 1 shifts
  77. by 1 day to Monday, and so on.
  78. yearly_start: Optional int specifying the start day of the yearly
  79. seasonality plot. 0 (default) starts the year on Jan 1. 1 shifts
  80. by 1 day to Jan 2, and so on.
  81. Returns
  82. -------
  83. A matplotlib figure.
  84. """
  85. # Identify components to be plotted
  86. components = ['trend']
  87. if m.holidays is not None and 'holidays' in fcst:
  88. components.append('holidays')
  89. components.extend([name for name in m.seasonalities
  90. if name in fcst])
  91. regressors = {'additive': False, 'multiplicative': False}
  92. for name, props in m.extra_regressors.items():
  93. regressors[props['mode']] = True
  94. for mode in ['additive', 'multiplicative']:
  95. if regressors[mode] and 'extra_regressors_{}'.format(mode) in fcst:
  96. components.append('extra_regressors_{}'.format(mode))
  97. npanel = len(components)
  98. fig, axes = plt.subplots(npanel, 1, facecolor='w',
  99. figsize=(9, 3 * npanel))
  100. if npanel == 1:
  101. axes = [axes]
  102. for ax, plot_name in zip(axes, components):
  103. if plot_name == 'trend':
  104. plot_forecast_component(
  105. m=m, fcst=fcst, name='trend', ax=ax, uncertainty=uncertainty,
  106. plot_cap=plot_cap,
  107. )
  108. elif plot_name == 'weekly':
  109. plot_weekly(
  110. m=m, ax=ax, uncertainty=uncertainty, weekly_start=weekly_start,
  111. )
  112. elif plot_name == 'yearly':
  113. plot_yearly(
  114. m=m, ax=ax, uncertainty=uncertainty, yearly_start=yearly_start,
  115. )
  116. elif plot_name in [
  117. 'holidays',
  118. 'extra_regressors_additive',
  119. 'extra_regressors_multiplicative',
  120. ]:
  121. plot_forecast_component(
  122. m=m, fcst=fcst, name=plot_name, ax=ax, uncertainty=uncertainty,
  123. plot_cap=False,
  124. )
  125. else:
  126. plot_seasonality(
  127. m=m, name=plot_name, ax=ax, uncertainty=uncertainty,
  128. )
  129. fig.tight_layout()
  130. return fig
  131. def plot_forecast_component(
  132. m, fcst, name, ax=None, uncertainty=True, plot_cap=False,
  133. ):
  134. """Plot a particular component of the forecast.
  135. Parameters
  136. ----------
  137. m: Prophet model.
  138. fcst: pd.DataFrame output of m.predict.
  139. name: Name of the component to plot.
  140. ax: Optional matplotlib Axes to plot on.
  141. uncertainty: Optional boolean to plot uncertainty intervals.
  142. plot_cap: Optional boolean indicating if the capacity should be shown
  143. in the figure, if available.
  144. Returns
  145. -------
  146. a list of matplotlib artists
  147. """
  148. artists = []
  149. if not ax:
  150. fig = plt.figure(facecolor='w', figsize=(10, 6))
  151. ax = fig.add_subplot(111)
  152. fcst_t = fcst['ds'].dt.to_pydatetime()
  153. artists += ax.plot(fcst_t, fcst[name], ls='-', c='#0072B2')
  154. if 'cap' in fcst and plot_cap:
  155. artists += ax.plot(fcst_t, fcst['cap'], ls='--', c='k')
  156. if m.logistic_floor and 'floor' in fcst and plot_cap:
  157. ax.plot(fcst_t, fcst['floor'], ls='--', c='k')
  158. if uncertainty:
  159. artists += [ax.fill_between(
  160. fcst_t, fcst[name + '_lower'], fcst[name + '_upper'],
  161. color='#0072B2', alpha=0.2)]
  162. ax.grid(True, which='major', c='gray', ls='-', lw=1, alpha=0.2)
  163. ax.set_xlabel('ds')
  164. ax.set_ylabel(name)
  165. if name in m.component_modes['multiplicative']:
  166. ax = set_y_as_percent(ax)
  167. return artists
  168. def seasonality_plot_df(m, ds):
  169. """Prepare dataframe for plotting seasonal components.
  170. Parameters
  171. ----------
  172. m: Prophet model.
  173. ds: List of dates for column ds.
  174. Returns
  175. -------
  176. A dataframe with seasonal components on ds.
  177. """
  178. df_dict = {'ds': ds, 'cap': 1., 'floor': 0.}
  179. for name in m.extra_regressors:
  180. df_dict[name] = 0.
  181. df = pd.DataFrame(df_dict)
  182. df = m.setup_dataframe(df)
  183. return df
  184. def plot_weekly(m, ax=None, uncertainty=True, weekly_start=0):
  185. """Plot the weekly component of the forecast.
  186. Parameters
  187. ----------
  188. m: Prophet model.
  189. ax: Optional matplotlib Axes to plot on. One will be created if this
  190. is not provided.
  191. uncertainty: Optional boolean to plot uncertainty intervals.
  192. weekly_start: Optional int specifying the start day of the weekly
  193. seasonality plot. 0 (default) starts the week on Sunday. 1 shifts
  194. by 1 day to Monday, and so on.
  195. Returns
  196. -------
  197. a list of matplotlib artists
  198. """
  199. artists = []
  200. if not ax:
  201. fig = plt.figure(facecolor='w', figsize=(10, 6))
  202. ax = fig.add_subplot(111)
  203. # Compute weekly seasonality for a Sun-Sat sequence of dates.
  204. days = (pd.date_range(start='2017-01-01', periods=7) +
  205. pd.Timedelta(days=weekly_start))
  206. df_w = seasonality_plot_df(m, days)
  207. seas = m.predict_seasonal_components(df_w)
  208. days = days.weekday_name
  209. artists += ax.plot(range(len(days)), seas['weekly'], ls='-',
  210. c='#0072B2')
  211. if uncertainty:
  212. artists += [ax.fill_between(range(len(days)),
  213. seas['weekly_lower'], seas['weekly_upper'],
  214. color='#0072B2', alpha=0.2)]
  215. ax.grid(True, which='major', c='gray', ls='-', lw=1, alpha=0.2)
  216. ax.set_xticks(range(len(days)))
  217. ax.set_xticklabels(days)
  218. ax.set_xlabel('Day of week')
  219. ax.set_ylabel('weekly')
  220. if m.seasonalities['weekly']['mode'] == 'multiplicative':
  221. ax = set_y_as_percent(ax)
  222. return artists
  223. def plot_yearly(m, ax=None, uncertainty=True, yearly_start=0):
  224. """Plot the yearly component of the forecast.
  225. Parameters
  226. ----------
  227. m: Prophet model.
  228. ax: Optional matplotlib Axes to plot on. One will be created if
  229. this is not provided.
  230. uncertainty: Optional boolean to plot uncertainty intervals.
  231. yearly_start: Optional int specifying the start day of the yearly
  232. seasonality plot. 0 (default) starts the year on Jan 1. 1 shifts
  233. by 1 day to Jan 2, and so on.
  234. Returns
  235. -------
  236. a list of matplotlib artists
  237. """
  238. artists = []
  239. if not ax:
  240. fig = plt.figure(facecolor='w', figsize=(10, 6))
  241. ax = fig.add_subplot(111)
  242. # Compute yearly seasonality for a Jan 1 - Dec 31 sequence of dates.
  243. days = (pd.date_range(start='2017-01-01', periods=365) +
  244. pd.Timedelta(days=yearly_start))
  245. df_y = seasonality_plot_df(m, days)
  246. seas = m.predict_seasonal_components(df_y)
  247. artists += ax.plot(
  248. df_y['ds'].dt.to_pydatetime(), seas['yearly'], ls='-', c='#0072B2')
  249. if uncertainty:
  250. artists += [ax.fill_between(
  251. df_y['ds'].dt.to_pydatetime(), seas['yearly_lower'],
  252. seas['yearly_upper'], color='#0072B2', alpha=0.2)]
  253. ax.grid(True, which='major', c='gray', ls='-', lw=1, alpha=0.2)
  254. months = MonthLocator(range(1, 13), bymonthday=1, interval=2)
  255. ax.xaxis.set_major_formatter(FuncFormatter(
  256. lambda x, pos=None: '{dt:%B} {dt.day}'.format(dt=num2date(x))))
  257. ax.xaxis.set_major_locator(months)
  258. ax.set_xlabel('Day of year')
  259. ax.set_ylabel('yearly')
  260. if m.seasonalities['yearly']['mode'] == 'multiplicative':
  261. ax = set_y_as_percent(ax)
  262. return artists
  263. def plot_seasonality(m, name, ax=None, uncertainty=True):
  264. """Plot a custom seasonal component.
  265. Parameters
  266. ----------
  267. m: Prophet model.
  268. name: Seasonality name, like 'daily', 'weekly'.
  269. ax: Optional matplotlib Axes to plot on. One will be created if
  270. this is not provided.
  271. uncertainty: Optional boolean to plot uncertainty intervals.
  272. Returns
  273. -------
  274. a list of matplotlib artists
  275. """
  276. artists = []
  277. if not ax:
  278. fig = plt.figure(facecolor='w', figsize=(10, 6))
  279. ax = fig.add_subplot(111)
  280. # Compute seasonality from Jan 1 through a single period.
  281. start = pd.to_datetime('2017-01-01 0000')
  282. period = m.seasonalities[name]['period']
  283. end = start + pd.Timedelta(days=period)
  284. plot_points = 200
  285. days = pd.to_datetime(np.linspace(start.value, end.value, plot_points))
  286. df_y = seasonality_plot_df(m, days)
  287. seas = m.predict_seasonal_components(df_y)
  288. artists += ax.plot(df_y['ds'].dt.to_pydatetime(), seas[name], ls='-',
  289. c='#0072B2')
  290. if uncertainty:
  291. artists += [ax.fill_between(
  292. df_y['ds'].dt.to_pydatetime(), seas[name + '_lower'],
  293. seas[name + '_upper'], color='#0072B2', alpha=0.2)]
  294. ax.grid(True, which='major', c='gray', ls='-', lw=1, alpha=0.2)
  295. xticks = pd.to_datetime(np.linspace(start.value, end.value, 7)
  296. ).to_pydatetime()
  297. ax.set_xticks(xticks)
  298. if period <= 2:
  299. fmt_str = '{dt:%T}'
  300. elif period < 14:
  301. fmt_str = '{dt:%m}/{dt:%d} {dt:%R}'
  302. else:
  303. fmt_str = '{dt:%m}/{dt:%d}'
  304. ax.xaxis.set_major_formatter(FuncFormatter(
  305. lambda x, pos=None: fmt_str.format(dt=num2date(x))))
  306. ax.set_xlabel('ds')
  307. ax.set_ylabel('{}'.format(name))
  308. if m.seasonalities[name]['mode'] == 'multiplicative':
  309. ax = set_y_as_percent(ax)
  310. return artists
  311. def set_y_as_percent(ax):
  312. yticks = 100 * ax.get_yticks()
  313. yticklabels = ['{0:.4g}%'.format(y) for y in yticks]
  314. ax.set_yticklabels(yticklabels)
  315. return ax
  316. def add_changepoints_to_plot(
  317. ax, m, fcst, threshold=0.01, cp_color='r', cp_linestyle='--', trend=True,
  318. ):
  319. """Add markers for significant changepoints to prophet forecast plot.
  320. Example:
  321. fig = m.plot(forecast)
  322. add_changepoints_to_plot(fig.gca(), m, forecast)
  323. Parameters
  324. ----------
  325. ax: axis on which to overlay changepoint markers.
  326. m: Prophet model.
  327. fcst: Forecast output from m.predict.
  328. threshold: Threshold on trend change magnitude for significance.
  329. cp_color: Color of changepoint markers.
  330. cp_linestyle: Linestyle for changepoint markers.
  331. trend: If True, will also overlay the trend.
  332. Returns
  333. -------
  334. a list of matplotlib artists
  335. """
  336. artists = []
  337. if trend:
  338. artists.append(ax.plot(fcst['ds'], fcst['trend'], c=cp_color))
  339. signif_changepoints = m.changepoints[
  340. np.abs(np.nanmean(m.params['delta'], axis=0)) >= threshold
  341. ]
  342. for cp in signif_changepoints:
  343. artists.append(ax.axvline(x=cp, c=cp_color, ls=cp_linestyle))
  344. return artists
  345. def plot_cross_validation_metric(df_cv, metric, rolling_window=0.1, ax=None):
  346. """Plot a performance metric vs. forecast horizon from cross validation.
  347. Cross validation produces a collection of out-of-sample model predictions
  348. that can be compared to actual values, at a range of different horizons
  349. (distance from the cutoff). This computes a specified performance metric
  350. for each prediction, and aggregated over a rolling window with horizon.
  351. This uses fbprophet.diagnostics.performance_metrics to compute the metrics.
  352. Valid values of metric are 'mse', 'rmse', 'mae', 'mape', and 'coverage'.
  353. rolling_window is the proportion of data included in the rolling window of
  354. aggregation. The default value of 0.1 means 10% of data are included in the
  355. aggregation for computing the metric.
  356. As a concrete example, if metric='mse', then this plot will show the
  357. squared error for each cross validation prediction, along with the MSE
  358. averaged over rolling windows of 10% of the data.
  359. Parameters
  360. ----------
  361. df_cv: The output from fbprophet.diagnostics.cross_validation.
  362. metric: Metric name, one of ['mse', 'rmse', 'mae', 'mape', 'coverage'].
  363. rolling_window: Proportion of data to use for rolling average of metric.
  364. In [0, 1]. Defaults to 0.1.
  365. ax: Optional matplotlib axis on which to plot. If not given, a new figure
  366. will be created.
  367. Returns
  368. -------
  369. a matplotlib figure.
  370. """
  371. if ax is None:
  372. fig = plt.figure(facecolor='w', figsize=(10, 6))
  373. ax = fig.add_subplot(111)
  374. else:
  375. fig = ax.get_figure()
  376. # Get the metric at the level of individual predictions, and with the rolling window.
  377. df_none = performance_metrics(df_cv, metrics=[metric], rolling_window=0)
  378. df_h = performance_metrics(df_cv, metrics=[metric], rolling_window=rolling_window)
  379. # Some work because matplotlib does not handle timedelta
  380. # Target ~10 ticks.
  381. tick_w = max(df_none['horizon'].astype('timedelta64[ns]')) / 10.
  382. # Find the largest time resolution that has <1 unit per bin.
  383. dts = ['D', 'h', 'm', 's', 'ms', 'us', 'ns']
  384. dt_names = [
  385. 'days', 'hours', 'minutes', 'seconds', 'milliseconds', 'microseconds',
  386. 'nanoseconds'
  387. ]
  388. dt_conversions = [
  389. 24 * 60 * 60 * 10 ** 9,
  390. 60 * 60 * 10 ** 9,
  391. 60 * 10 ** 9,
  392. 10 ** 9,
  393. 10 ** 6,
  394. 10 ** 3,
  395. 1.,
  396. ]
  397. for i, dt in enumerate(dts):
  398. if np.timedelta64(1, dt) < np.timedelta64(tick_w, 'ns'):
  399. break
  400. x_plt = df_none['horizon'].astype('timedelta64[ns]').astype(int) / float(dt_conversions[i])
  401. x_plt_h = df_h['horizon'].astype('timedelta64[ns]').astype(int) / float(dt_conversions[i])
  402. ax.plot(x_plt, df_none[metric], '.', alpha=0.5, c='gray')
  403. ax.plot(x_plt_h, df_h[metric], '-', c='b')
  404. ax.grid(True)
  405. ax.set_xlabel('Horizon ({})'.format(dt_names[i]))
  406. ax.set_ylabel(metric)
  407. return fig