prophet

  36class Prophet:
  37    """Prophet forecaster.
  38
  39    Parameters
  40    ----------
  41    growth: String 'linear', 'logistic' or 'flat' to specify a linear, logistic or
  42        flat trend.
  43    changepoints: List of dates at which to include potential changepoints. If
  44        not specified, potential changepoints are selected automatically.
  45    n_changepoints: Number of potential changepoints to include. Not used
  46        if input `changepoints` is supplied. If `changepoints` is not supplied,
  47        then n_changepoints potential changepoints are selected uniformly from
  48        the first `changepoint_range` proportion of the history.
  49    changepoint_range: Proportion of history in which trend changepoints will
  50        be estimated. Defaults to 0.8 for the first 80%. Not used if
  51        `changepoints` is specified.
  52    yearly_seasonality: Fit yearly seasonality.
  53        Can be 'auto', True, False, or a number of Fourier terms to generate.
  54    weekly_seasonality: Fit weekly seasonality.
  55        Can be 'auto', True, False, or a number of Fourier terms to generate.
  56    daily_seasonality: Fit daily seasonality.
  57        Can be 'auto', True, False, or a number of Fourier terms to generate.
  58    holidays: pd.DataFrame with columns holiday (string) and ds (date type)
  59        and optionally columns lower_window and upper_window which specify a
  60        range of days around the date to be included as holidays.
  61        lower_window=-2 will include 2 days prior to the date as holidays. Also
  62        optionally can have a column prior_scale specifying the prior scale for
  63        that holiday.
  64    seasonality_mode: 'additive' (default) or 'multiplicative'.
  65    seasonality_prior_scale: Parameter modulating the strength of the
  66        seasonality model. Larger values allow the model to fit larger seasonal
  67        fluctuations, smaller values dampen the seasonality. Can be specified
  68        for individual seasonalities using add_seasonality.
  69    holidays_prior_scale: Parameter modulating the strength of the holiday
  70        components model, unless overridden in the holidays input.
  71    changepoint_prior_scale: Parameter modulating the flexibility of the
  72        automatic changepoint selection. Large values will allow many
  73        changepoints, small values will allow few changepoints.
  74    mcmc_samples: Integer, if greater than 0, will do full Bayesian inference
  75        with the specified number of MCMC samples. If 0, will do MAP
  76        estimation.
  77    interval_width: Float, width of the uncertainty intervals provided
  78        for the forecast. If mcmc_samples=0, this will be only the uncertainty
  79        in the trend using the MAP estimate of the extrapolated generative
  80        model. If mcmc.samples>0, this will be integrated over all model
  81        parameters, which will include uncertainty in seasonality.
  82    uncertainty_samples: Number of simulated draws used to estimate
  83        uncertainty intervals. Settings this value to 0 or False will disable
  84        uncertainty estimation and speed up the calculation.
  85    stan_backend: str as defined in StanBackendEnum default: None - will try to
  86        iterate over all available backends and find the working one.
  87    scaling: 'absmax' (default) or 'minmax'.
  88    holidays_mode: 'additive' or 'multiplicative'. Defaults to seasonality_mode.
  89    """
  90
  91    growth: Literal["linear", "logistic", "flat"]
  92    changepoints: pd.Series[pd.Timestamp] | None
  93    n_changepoints: int
  94    specified_changepoints: bool
  95    changepoint_range: float
  96    yearly_seasonality: Literal["auto"] | int
  97    weekly_seasonality: Literal["auto"] | int
  98    daily_seasonality: Literal["auto"] | int
  99    holidays: pd.DataFrame | None
 100    seasonality_mode: _Mode
 101    holidays_mode: _Mode
 102    seasonality_prior_scale: float
 103    holidays_prior_scale: float
 104    changepoint_prior_scale: float
 105    mcmc_samples: int
 106    interval_width: float
 107    uncertainty_samples: int
 108    scaling: Literal["absmax", "minmax"]
 109
 110    start: pd.Timestamp | None
 111    y_min: float | None
 112    y_scale: float | None
 113    logistic_floor: bool
 114    t_scale: pd.Timedelta | None
 115    changepoints_t: npt.NDArray[np.float64] | None
 116    seasonalities: OrderedDict[str, dict[str, Any]]
 117    extra_regressors: OrderedDict[str, dict[str, Any]]
 118    country_holidays: str | None
 119    stan_fit: Any | None
 120    params: dict[str, Any]
 121    history: pd.DataFrame | None
 122    history_dates: pd.Series[pd.Timestamp] | None
 123    train_component_cols: pd.DataFrame | None
 124    component_modes: dict[_Mode, list[str]] | None
 125    train_holiday_names: pd.Series | None
 126    fit_kwargs: dict[str, Any]
 127    stan_backend: IStanBackend | None
 128    # Optional marker set on nested regressor-predictor models for diagnostics/tests.
 129    _regressor_name: str | None
 130
 131    def __init__(
 132        self,
 133        growth: Literal["linear", "logistic", "flat"] = "linear",
 134        changepoints: pd.Series | list[pd.Timestamp] | None = None,
 135        n_changepoints: int = 25,
 136        changepoint_range: float = 0.8,
 137        yearly_seasonality: Literal["auto"] | int = "auto",
 138        weekly_seasonality: Literal["auto"] | int = "auto",
 139        daily_seasonality: Literal["auto"] | int = "auto",
 140        holidays: pd.DataFrame | None = None,
 141        seasonality_mode: _Mode = "additive",
 142        seasonality_prior_scale: SupportsFloat = 10.0,
 143        holidays_prior_scale: SupportsFloat = 10.0,
 144        changepoint_prior_scale: SupportsFloat = 0.05,
 145        mcmc_samples: int = 0,
 146        interval_width: float = 0.80,
 147        uncertainty_samples: int = 1000,
 148        stan_backend: str | None = None,
 149        scaling: Literal["absmax", "minmax"] = "absmax",
 150        holidays_mode: _Mode | None = None,
 151    ) -> None:
 152        self.growth = growth
 153
 154        if changepoints is not None:
 155            self.changepoints = pd.Series(pd.to_datetime(changepoints), name='ds')
 156            self.n_changepoints = len(self.changepoints)
 157            self.specified_changepoints = True
 158        else:
 159            self.changepoints = changepoints
 160            self.n_changepoints = n_changepoints
 161            self.specified_changepoints = False
 162
 163        self.changepoint_range = changepoint_range
 164        self.yearly_seasonality = yearly_seasonality
 165        self.weekly_seasonality = weekly_seasonality
 166        self.daily_seasonality = daily_seasonality
 167        self.holidays = holidays
 168
 169        self.seasonality_mode = seasonality_mode
 170        self.holidays_mode = holidays_mode or self.seasonality_mode
 171
 172        self.seasonality_prior_scale = float(seasonality_prior_scale)
 173        self.changepoint_prior_scale = float(changepoint_prior_scale)
 174        self.holidays_prior_scale = float(holidays_prior_scale)
 175
 176        self.mcmc_samples = mcmc_samples
 177        self.interval_width = interval_width
 178        self.uncertainty_samples = uncertainty_samples
 179        if scaling not in ("absmax", "minmax"):
 180            raise ValueError("scaling must be one of 'absmax' or 'minmax'")
 181        self.scaling = scaling
 182
 183        # Set during fitting or by other methods
 184        self.start = None
 185        self.y_min = None
 186        self.y_scale = None
 187        self.logistic_floor = False
 188        self.t_scale = None
 189        self.changepoints_t = None
 190        self.seasonalities = OrderedDict({})
 191        self.extra_regressors = OrderedDict({})
 192        self.country_holidays = None
 193        self.stan_fit = None
 194        self.params = {}
 195        self.history = None
 196        self.history_dates = None
 197        self.train_component_cols = None
 198        self.component_modes = None
 199        self.train_holiday_names = None
 200        self.fit_kwargs = {}
 201        self._regressor_name = None
 202        self.validate_inputs()
 203        self._load_stan_backend(stan_backend)
 204
 205    def _load_stan_backend(self, stan_backend: str | None) -> None:
 206        if stan_backend is None:
 207            for i in StanBackendEnum:
 208                try:
 209                    logger.debug("Trying to load backend: %s", i.name)
 210                    return self._load_stan_backend(i.name)
 211                except Exception as e:
 212                    logger.debug("Unable to load backend %s (%s), trying the next one", i.name, e)
 213        else:
 214            # pyrefly:ignore[not-callable]
 215            self.stan_backend = StanBackendEnum.get_backend_class(stan_backend)()
 216
 217        assert self.stan_backend
 218        logger.debug("Loaded stan backend: %s", self.stan_backend.get_type())
 219
 220    def validate_inputs(self) -> None:
 221        """Validates the inputs to Prophet."""
 222        if self.growth not in ('linear', 'logistic', 'flat'):
 223            raise ValueError(
 224                'Parameter "growth" should be "linear", "logistic" or "flat".')
 225        if not isinstance(self.changepoint_range, (int, float)):
 226            raise ValueError("changepoint_range must be a number in [0, 1]'")
 227        if ((self.changepoint_range < 0) or (self.changepoint_range > 1)):
 228            raise ValueError('Parameter "changepoint_range" must be in [0, 1]')
 229        if self.holidays is not None:
 230            if not (
 231                isinstance(self.holidays, pd.DataFrame)
 232                and 'ds' in self.holidays  # noqa W503
 233                and 'holiday' in self.holidays  # noqa W503
 234            ):
 235                raise ValueError('holidays must be a DataFrame with "ds" and '
 236                                 '"holiday" columns.')
 237            self.holidays['ds'] = pd.to_datetime(self.holidays['ds'])
 238            if (
 239                self.holidays['ds'].isnull().any()
 240                or self.holidays['holiday'].isnull().any()
 241            ):
 242                raise ValueError('Found a NaN in holidays dataframe.')
 243            has_lower = 'lower_window' in self.holidays
 244            has_upper = 'upper_window' in self.holidays
 245            if has_lower + has_upper == 1:
 246                raise ValueError('Holidays must have both lower_window and ' +
 247                                 'upper_window, or neither')
 248            if has_lower:
 249                if self.holidays['lower_window'].max() > 0:
 250                    raise ValueError('Holiday lower_window should be <= 0')
 251                if self.holidays['upper_window'].min() < 0:
 252                    raise ValueError('Holiday upper_window should be >= 0')
 253            for h in self.holidays['holiday'].unique():
 254                self.validate_column_name(h, check_holidays=False)
 255        if self.seasonality_mode not in ['additive', 'multiplicative']:
 256            raise ValueError(
 257                'seasonality_mode must be "additive" or "multiplicative"'
 258            )
 259        if self.holidays_mode not in ['additive', 'multiplicative']:
 260            raise ValueError(
 261                'holidays_mode must be "additive" or "multiplicative"'
 262            )
 263
 264    def validate_column_name(
 265        self,
 266        name: str,
 267        check_holidays: bool = True,
 268        check_seasonalities: bool = True,
 269        check_regressors: bool = True,
 270    ) -> None:
 271        """Validates the name of a seasonality, holiday, or regressor.
 272
 273        Parameters
 274        ----------
 275        name: string
 276        check_holidays: bool check if name already used for holiday
 277        check_seasonalities: bool check if name already used for seasonality
 278        check_regressors: bool check if name already used for regressor
 279        """
 280        if '_delim_' in name:
 281            raise ValueError('Name cannot contain "_delim_"')
 282        reserved_names = [
 283            'trend', 'additive_terms', 'daily', 'weekly', 'yearly',
 284            'holidays', 'zeros', 'extra_regressors_additive', 'yhat',
 285            'extra_regressors_multiplicative', 'multiplicative_terms',
 286        ]
 287        rn_l = [n + '_lower' for n in reserved_names]
 288        rn_u = [n + '_upper' for n in reserved_names]
 289        reserved_names.extend(rn_l)
 290        reserved_names.extend(rn_u)
 291        reserved_names.extend([
 292            'ds', 'y', 'cap', 'floor', 'y_scaled', 'cap_scaled'])
 293        if name in reserved_names:
 294            raise ValueError(
 295                'Name {name!r} is reserved.'.format(name=name)
 296            )
 297        if (check_holidays and self.holidays is not None and
 298                name in self.holidays['holiday'].unique()):
 299            raise ValueError(
 300                'Name {name!r} already used for a holiday.'.format(name=name)
 301            )
 302        if (check_holidays and self.country_holidays is not None and
 303                name in get_holiday_names(self.country_holidays)):
 304            raise ValueError(
 305                'Name {name!r} is a holiday name in {country_holidays}.'
 306                .format(name=name, country_holidays=self.country_holidays)
 307            )
 308        if check_seasonalities and name in self.seasonalities:
 309            raise ValueError(
 310                'Name {name!r} already used for a seasonality.'
 311                .format(name=name)
 312            )
 313        if check_regressors and name in self.extra_regressors:
 314            raise ValueError(
 315                'Name {name!r} already used for an added regressor.'
 316                .format(name=name)
 317            )
 318
 319    def setup_dataframe(self, df: pd.DataFrame, initialize_scales: bool = False) -> pd.DataFrame:
 320        """Prepare dataframe for fitting or predicting.
 321
 322        Adds a time index and scales y. Creates auxiliary columns 't', 't_ix',
 323        'y_scaled', and 'cap_scaled'. These columns are used during both
 324        fitting and predicting.
 325
 326        Parameters
 327        ----------
 328        df: pd.DataFrame with columns ds, y, and cap if logistic growth. Any
 329            specified additional regressors must also be present.
 330        initialize_scales: Boolean set scaling factors in self from df.
 331
 332        Returns
 333        -------
 334        pd.DataFrame prepared for fitting or predicting.
 335        """
 336        if 'y' in df:  # 'y' will be in training data
 337            df['y'] = pd.to_numeric(df['y'])
 338            if np.isinf(df['y'].values).any():
 339                raise ValueError('Found infinity in column y.')
 340        if df['ds'].dtype == np.int64:
 341            df['ds'] = df['ds'].astype(str)
 342        df['ds'] = pd.to_datetime(df['ds'])
 343        if df['ds'].dt.tz is not None:
 344            raise ValueError(
 345                'Column ds has timezone specified, which is not supported. '
 346                'Remove timezone.'
 347            )
 348        if df['ds'].isnull().any():
 349            raise ValueError('Found NaN in column ds.')
 350        for name in self.extra_regressors:
 351            if name not in df:
 352                raise ValueError(
 353                    'Regressor {name!r} missing from dataframe'
 354                    .format(name=name)
 355                )
 356            df[name] = pd.to_numeric(df[name])
 357            if df[name].isnull().any():
 358                raise ValueError(
 359                    'Found NaN in column {name!r}'.format(name=name)
 360                )
 361        for props in self.seasonalities.values():
 362            condition_name = props['condition_name']
 363            if condition_name is not None:
 364                if condition_name not in df:
 365                    raise ValueError(
 366                        'Condition {condition_name!r} missing from dataframe'
 367                        .format(condition_name=condition_name)
 368                    )
 369                if not df[condition_name].isin([True, False]).all():
 370                    raise ValueError(
 371                        'Found non-boolean in column {condition_name!r}'
 372                        .format(condition_name=condition_name)
 373                    )
 374                df[condition_name] = df[condition_name].astype('bool')
 375
 376        if df.index.name == 'ds':
 377            df.index.name = None
 378        df = df.sort_values('ds', kind='mergesort')
 379        df = df.reset_index(drop=True)
 380
 381        self.initialize_scales(initialize_scales, df)
 382
 383        if self.logistic_floor:
 384            if 'floor' not in df:
 385                raise ValueError('Expected column "floor".')
 386        else:
 387            if self.scaling == "absmax":
 388                df['floor'] = 0.
 389            elif self.scaling == "minmax":
 390                df['floor'] = self.y_min
 391        if self.growth == 'logistic':
 392            if 'cap' not in df:
 393                raise ValueError(
 394                    'Capacities must be supplied for logistic growth in '
 395                    'column "cap"'
 396                )
 397            if (df['cap'] <= df['floor']).any():
 398                raise ValueError(
 399                    'cap must be greater than floor (which defaults to 0).'
 400                )
 401            df['cap_scaled'] = (df['cap'] - df['floor']) / self.y_scale  # pyrefly:ignore[unsupported-operation]
 402
 403        df['t'] = (df['ds'] - self.start) / self.t_scale  # pyrefly:ignore[unsupported-operation]
 404        if 'y' in df:
 405            df['y_scaled'] = (df['y'] - df['floor']) / self.y_scale  # pyrefly:ignore[unsupported-operation]
 406
 407        for name, props in self.extra_regressors.items():
 408            df[name] = ((df[name] - props['mu']) / props['std'])
 409        return df
 410
 411    def initialize_scales(self, initialize_scales: bool, df: pd.DataFrame) -> None:
 412        """Initialize model scales.
 413
 414        Sets model scaling factors using df.
 415
 416        Parameters
 417        ----------
 418        initialize_scales: Boolean set the scales or not.
 419        df: pd.DataFrame for setting scales.
 420        """
 421        if not initialize_scales:
 422            return
 423
 424        if self.growth == 'logistic' and 'floor' in df:
 425            self.logistic_floor = True
 426            if self.scaling == "absmax":
 427                self.y_min = float((df['y'] - df['floor']).abs().min())
 428                self.y_scale = float((df['y'] - df['floor']).abs().max())
 429            elif self.scaling == "minmax":
 430                self.y_min = df['floor'].min()
 431                self.y_scale = float(df['cap'].max() - self.y_min)
 432        else:
 433            if self.scaling == "absmax":
 434                self.y_min = 0.
 435                self.y_scale = float((df['y']).abs().max())
 436            elif self.scaling == "minmax":
 437                self.y_min = df['y'].min()
 438                self.y_scale =  float(df['y'].max() - self.y_min)
 439        if self.y_scale == 0:
 440            self.y_scale = 1.0
 441
 442        self.start = df['ds'].min()
 443        self.t_scale = df['ds'].max() - self.start
 444        for name, props in self.extra_regressors.items():
 445            standardize = props['standardize']
 446            n_vals = len(df[name].unique())
 447            if n_vals < 2:
 448                standardize = False
 449            if standardize == 'auto':
 450                if set(df[name].unique()) == {1, 0}:
 451                    standardize = False #  Don't standardize binary variables.
 452                else:
 453                    standardize = True
 454            if standardize:
 455                mu = float(df[name].mean())
 456                std = float(df[name].std())
 457                self.extra_regressors[name]['mu'] = mu
 458                self.extra_regressors[name]['std'] = std
 459
 460    def set_changepoints(self) -> None:
 461        """Set changepoints
 462
 463        Sets m$changepoints to the dates of changepoints. Either:
 464        1) The changepoints were passed in explicitly.
 465            A) They are empty.
 466            B) They are not empty, and need validation.
 467        2) We are generating a grid of them.
 468        3) The user prefers no changepoints be used.
 469        """
 470        if self.changepoints is not None:
 471            if len(self.changepoints) == 0:
 472                pass
 473            else:
 474                history = cast(pd.DataFrame, self.history)
 475                too_low = min(self.changepoints) < history['ds'].min()
 476                too_high = max(self.changepoints) > history['ds'].max()
 477                if too_low or too_high:
 478                    raise ValueError('Changepoints must fall within training data.')
 479        else:
 480            # Place potential changepoints evenly through first
 481            # `changepoint_range` proportion of the history
 482            history = cast(pd.DataFrame, self.history)
 483            hist_size = int(np.floor(history.shape[0] * self.changepoint_range))
 484            if self.n_changepoints + 1 > hist_size:
 485                self.n_changepoints = hist_size - 1
 486                logger.info(
 487                    'n_changepoints greater than number of observations. '
 488                    'Using {n_changepoints}.'
 489                    .format(n_changepoints=self.n_changepoints)
 490                )
 491            if self.n_changepoints > 0:
 492                cp_indexes = (
 493                    np.linspace(0, hist_size - 1, self.n_changepoints + 1)
 494                        .round()
 495                        .astype(int)
 496                )
 497                self.changepoints = history.iloc[cp_indexes]['ds'].tail(-1)
 498            else:
 499                # set empty changepoints
 500                self.changepoints = pd.Series(pd.to_datetime([]), name='ds')
 501        if len(self.changepoints) > 0:
 502            self.changepoints_t = np.sort(np.array(
 503                (self.changepoints - self.start) / self.t_scale))  # pyrefly:ignore[unsupported-operation]
 504        else:
 505            self.changepoints_t = np.array([0])  # dummy changepoint
 506
 507    @staticmethod
 508    def fourier_series(
 509        dates: pd.Series,
 510        period: float,
 511        series_order: int,
 512    ) -> np.ndarray[tuple[int, int], np.dtype[np.float64]]:
 513        """Provides Fourier series components with the specified frequency
 514        and order.
 515
 516        Parameters
 517        ----------
 518        dates: pd.Series containing timestamps.
 519        period: Number of days of the period.
 520        series_order: Number of components.
 521
 522        Returns
 523        -------
 524        Matrix with seasonality features.
 525        """
 526        if not (series_order >= 1):
 527            raise ValueError("series_order must be >= 1")
 528
 529        epoch = pd.Timestamp("1970-01-01", tz=dates.dt.tz)
 530        t = (dates - epoch).dt.total_seconds() / (24 * 60 * 60)
 531
 532        x_T = np.pi * 2 * t
 533        fourier_components = np.empty((dates.shape[0], 2 * series_order))
 534        for i in range(series_order):
 535            c = (i + 1) / period * x_T
 536            fourier_components[:, 2 * i] = np.sin(c)
 537            fourier_components[:, (2 * i) + 1] = np.cos(c)
 538        return fourier_components
 539
 540    @classmethod
 541    def make_seasonality_features(
 542        cls,
 543        dates: pd.Series[pd.Timestamp],
 544        period: float,
 545        series_order: int,
 546        prefix: str,
 547    ) -> pd.DataFrame:
 548        """Data frame with seasonality features.
 549
 550        Parameters
 551        ----------
 552        cls: Prophet class.
 553        dates: pd.Series containing timestamps.
 554        period: Number of days of the period.
 555        series_order: Number of components.
 556        prefix: Column name prefix.
 557
 558        Returns
 559        -------
 560        pd.DataFrame with seasonality features.
 561        """
 562        features = cls.fourier_series(dates, period, series_order)
 563        columns = [
 564            '{}_delim_{}'.format(prefix, i + 1)
 565            for i in range(features.shape[1])
 566        ]
 567        return pd.DataFrame(features, columns=columns)
 568
 569    def construct_holiday_dataframe(self, dates: pd.Series[pd.Timestamp]) -> pd.DataFrame:
 570        """Construct a dataframe of holiday dates.
 571
 572        Will combine self.holidays with the built-in country holidays
 573        corresponding to input dates, if self.country_holidays is set.
 574
 575        Parameters
 576        ----------
 577        dates: pd.Series containing timestamps used for computing seasonality.
 578
 579        Returns
 580        -------
 581        dataframe of holiday dates, in holiday dataframe format used in
 582        initialization.
 583        """
 584        all_holidays = pd.DataFrame()
 585        if self.holidays is not None:
 586            all_holidays = self.holidays.copy()
 587        if self.country_holidays is not None:
 588            year_list = list({x.year for x in dates})
 589            country_holidays_df = make_holidays_df(
 590                year_list=year_list, country=self.country_holidays
 591            )
 592            all_holidays = pd.concat((all_holidays, country_holidays_df),
 593                                     sort=False)
 594            all_holidays.reset_index(drop=True, inplace=True)
 595        # Drop future holidays not previously seen in training data
 596        if self.train_holiday_names is not None:
 597            # Remove holiday names didn't show up in fit
 598            index_to_drop = all_holidays.index[
 599                np.logical_not(
 600                    all_holidays.holiday.isin(self.train_holiday_names)
 601                )
 602            ]
 603            all_holidays = all_holidays.drop(index_to_drop)
 604            # Add holiday names in fit but not in predict with ds as NA
 605            holidays_to_add = pd.DataFrame({
 606                'holiday': self.train_holiday_names[
 607                    np.logical_not(self.train_holiday_names
 608                                       .isin(all_holidays.holiday))
 609                ]
 610            })
 611            all_holidays = pd.concat((all_holidays, holidays_to_add),
 612                                     sort=False)
 613            all_holidays.reset_index(drop=True, inplace=True)
 614        return all_holidays
 615
 616    def make_holiday_features(
 617        self,
 618        dates: pd.Series[pd.Timestamp],
 619        holidays: pd.DataFrame,
 620    ) -> tuple[pd.DataFrame, list[float], list[str]]:
 621        """Construct a dataframe of holiday features.
 622
 623        Parameters
 624        ----------
 625        dates: pd.Series containing timestamps used for computing seasonality.
 626        holidays: pd.Dataframe containing holidays, as returned by
 627            construct_holiday_dataframe.
 628
 629        Returns
 630        -------
 631        holiday_features: pd.DataFrame with a column for each holiday.
 632        prior_scale_list: List of prior scales for each holiday column.
 633        holiday_names: List of names of holidays
 634        """
 635        # Holds columns of our future matrix.
 636        expanded_holidays = defaultdict(lambda: np.zeros(dates.shape[0]))
 637        prior_scales = {}
 638        # Makes an index so we can perform `get_loc` below.
 639        # Strip to just dates.
 640        row_index = pd.DatetimeIndex(dates.dt.date)
 641
 642        for row in holidays.itertuples():
 643            dt = cast(pd.Timestamp, row.ds).date()
 644            try:
 645                lw = int(getattr(row, 'lower_window', 0))
 646                uw = int(getattr(row, 'upper_window', 0))
 647            except ValueError:
 648                lw = 0
 649                uw = 0
 650            ps = float(getattr(row, 'prior_scale', self.holidays_prior_scale))
 651            if np.isnan(ps):
 652                ps = float(self.holidays_prior_scale)
 653            if row.holiday in prior_scales and prior_scales[row.holiday] != ps:
 654                raise ValueError(
 655                    'Holiday {holiday!r} does not have consistent prior '
 656                    'scale specification.'.format(holiday=row.holiday)
 657                )
 658            if ps <= 0:
 659                raise ValueError('Prior scale must be > 0')
 660            prior_scales[row.holiday] = ps
 661
 662            for offset in range(lw, uw + 1):
 663                occurrence = pd.to_datetime(dt + timedelta(days=offset))
 664                try:
 665                    loc = row_index.get_loc(occurrence)
 666                except KeyError:
 667                    loc = None
 668                key = '{}_delim_{}{}'.format(
 669                    row.holiday,
 670                    '+' if offset >= 0 else '-',
 671                    abs(offset)
 672                )
 673                if loc is not None:
 674                    expanded_holidays[key][loc] = 1.
 675                else:
 676                    expanded_holidays[key]  # Access key to generate value
 677        holiday_features = pd.DataFrame(expanded_holidays)
 678        # Make sure column order is consistent
 679        holiday_features = holiday_features[sorted(holiday_features.columns
 680                                                                   .tolist())]
 681        prior_scale_list = [
 682            prior_scales[h.split('_delim_')[0]]
 683            for h in holiday_features.columns
 684        ]
 685        holiday_names = list(prior_scales.keys())
 686        # Store holiday names used in fit
 687        if self.train_holiday_names is None:
 688            self.train_holiday_names = pd.Series(holiday_names)
 689        return holiday_features, prior_scale_list, holiday_names
 690
 691    def add_regressor(
 692        self,
 693        name: str,
 694        prior_scale: float | None = None,
 695        standardize: Literal['auto'] | bool = 'auto',
 696        mode: _Mode | None = None,
 697        regressor_predictor: bool | dict | None = None,
 698    ) -> Self:
 699        """Add an additional regressor to be used for fitting and predicting.
 700
 701        The dataframe passed to `fit` and `predict` will have a column with the
 702        specified name to be used as a regressor. When standardize='auto', the
 703        regressor will be standardized unless it is binary. The regression
 704        coefficient is given a prior with the specified scale parameter.
 705        Decreasing the prior scale will add additional regularization. If no
 706        prior scale is provided, self.holidays_prior_scale will be used.
 707        Mode can be specified as either 'additive' or 'multiplicative'. If not
 708        specified, self.seasonality_mode will be used. 'additive' means the
 709        effect of the regressor will be added to the trend, 'multiplicative'
 710        means it will multiply the trend.
 711
 712        Parameters
 713        ----------
 714        name: string name of the regressor.
 715        prior_scale: optional float scale for the normal prior. If not
 716            provided, self.holidays_prior_scale will be used.
 717        standardize: optional, specify whether this regressor will be
 718            standardized prior to fitting. Can be 'auto' (standardize if not
 719            binary), True, or False.
 720        mode: optional, 'additive' or 'multiplicative'. Defaults to
 721            self.seasonality_mode.
 722        regressor_predictor: optional. If provided, fits a dedicated Prophet
 723            model to forecast this regressor. Set to True to use default Prophet
 724            parameters, or provide a dict of keyword arguments for the
 725            underlying Prophet constructor. When set, future regressor values
 726            will be generated from this model (with uncertainty if available)
 727            instead of using user-supplied point estimates.
 728
 729        Returns
 730        -------
 731        The prophet object.
 732        """
 733        if self.history is not None:
 734            raise Exception(
 735                "Regressors must be added prior to model fitting.")
 736        self.validate_column_name(name, check_regressors=False)
 737        if prior_scale is None:
 738            prior_scale = float(self.holidays_prior_scale)
 739        if mode is None:
 740            mode = self.seasonality_mode
 741        if prior_scale <= 0:
 742            raise ValueError('Prior scale must be > 0')
 743        if mode not in ['additive', 'multiplicative']:
 744            raise ValueError("mode must be 'additive' or 'multiplicative'")
 745        predictor_spec = None
 746        if regressor_predictor:
 747            if isinstance(regressor_predictor, dict):
 748                predictor_spec = deepcopy(regressor_predictor)
 749            else:
 750                # Truthy non-dict values (typically True) use default Prophet params.
 751                predictor_spec = {}
 752        self.extra_regressors[name] = {
 753            'prior_scale': prior_scale,
 754            'standardize': standardize,
 755            'mu': 0.,
 756            'std': 1.,
 757            'mode': mode,
 758            'predictor_spec': predictor_spec,
 759            'predictor': None,
 760        }
 761        return self
 762
 763    def _fit_regressor_models(self, df: pd.DataFrame):
 764        """Fit regressor-specific Prophet models when configured."""
 765        for name, props in self.extra_regressors.items():
 766            spec = props.get('predictor_spec')
 767            if spec is None:
 768                continue
 769            regressor_df = df[['ds', name]].copy()
 770            regressor_df = regressor_df[regressor_df[name].notnull()]
 771            if regressor_df.shape[0] < 2:
 772                raise ValueError(f"Not enough data to fit regressor model for '{name}'.")
 773            regressor_df = regressor_df.rename(columns={name: 'y'})
 774            logger.info("Fitting regressor model '%s' with %d observations", name, regressor_df.shape[0])
 775            predictor_kwargs = deepcopy(spec)
 776            if 'stan_backend' not in predictor_kwargs and self.stan_backend is not None:
 777                predictor_kwargs['stan_backend'] = self.stan_backend.get_type()
 778            predictor = self.__class__(**predictor_kwargs)
 779            predictor._regressor_name = name  # marker for diagnostics/testing
 780            predictor.fit(regressor_df, **self.fit_kwargs)
 781            props['predictor'] = predictor
 782
 783    def _prepare_regressors_for_predict(
 784        self, df: pd.DataFrame, n_samples: int | None = None
 785    ) -> tuple[pd.DataFrame, dict[str, npt.NDArray[np.float64]] | None]:
 786        """Populate regressor columns using fitted regressor models if present.
 787
 788        Returns updated dataframe and optional regressor sample draws.
 789        """
 790        regressor_samples: dict[str, npt.NDArray[np.float64]] | None = (
 791            {} if n_samples else None
 792        )
 793        last_history_date = self.history['ds'].max() if self.history is not None else None
 794        for name, props in self.extra_regressors.items():
 795            predictor = props.get('predictor')
 796            if props.get('predictor_spec') is None or predictor is None:
 797                continue
 798            if name not in df:
 799                df[name] = np.nan
 800            future_mask = (
 801                df['ds'] > last_history_date
 802                if last_history_date is not None
 803                else df[name].isna()
 804            )
 805            history_mask = ~future_mask
 806            if history_mask.any():
 807                missing_hist = history_mask & df[name].isna()
 808                if missing_hist.any() and self.history is not None and name in self.history:
 809                    hist_lookup = self.history.set_index('ds')[name]
 810                    props = self.extra_regressors[name]
 811                    hist_lookup = hist_lookup * props['std'] + props['mu']
 812                    df.loc[missing_hist, name] = hist_lookup.reindex(
 813                        df.loc[missing_hist, 'ds']
 814                    ).to_numpy()
 815            if future_mask.any():
 816                future_dates = df.loc[future_mask, 'ds']
 817                reg_future = pd.DataFrame({'ds': future_dates})
 818                logger.info(
 819                    "Running regressor model '%s' for %d future dates",
 820                    name,
 821                    len(reg_future),
 822                )
 823                reg_pred = predictor.predict(reg_future)
 824                df.loc[future_mask, name] = reg_pred['yhat'].to_numpy()
 825                if n_samples and regressor_samples is not None:
 826                    sample_dict = predictor.predictive_samples(reg_future)
 827                    reg_samples = sample_dict['yhat']
 828                    samples_full = np.tile(
 829                        df[name].to_numpy()[:, None], (1, reg_samples.shape[1])
 830                    )
 831                    # Use boolean indexing on ndarray, not pandas ExtensionArray 2D slice.
 832                    samples_full[future_mask.to_numpy(), :] = reg_samples
 833                    regressor_samples[name] = samples_full
 834                    logger.info(
 835                        "Collected %d predictive samples for regressor '%s'",
 836                        reg_samples.shape[1],
 837                        name,
 838                    )
 839            elif n_samples and regressor_samples is not None:
 840                n_draws = n_samples if isinstance(n_samples, int) else 1
 841                regressor_samples[name] = np.tile(
 842                    df[name].to_numpy()[:, None], (1, n_draws)
 843                )
 844        if regressor_samples is not None and len(regressor_samples) == 0:
 845            regressor_samples = None
 846        return df, regressor_samples
 847
 848    def add_seasonality(
 849        self,
 850        name: str,
 851        period: float,
 852        fourier_order: int,
 853        prior_scale: float | None = None,
 854        mode: _Mode | None = None,
 855        condition_name: str | None = None,
 856    ) -> Self:
 857        """Add a seasonal component with specified period, number of Fourier
 858        components, and prior scale.
 859
 860        Increasing the number of Fourier components allows the seasonality to
 861        change more quickly (at risk of overfitting). Default values for yearly
 862        and weekly seasonalities are 10 and 3 respectively.
 863
 864        Increasing prior scale will allow this seasonality component more
 865        flexibility, decreasing will dampen it. If not provided, will use the
 866        seasonality_prior_scale provided on Prophet initialization (defaults
 867        to 10).
 868
 869        Mode can be specified as either 'additive' or 'multiplicative'. If not
 870        specified, self.seasonality_mode will be used (defaults to additive).
 871        Additive means the seasonality will be added to the trend,
 872        multiplicative means it will multiply the trend.
 873
 874        If condition_name is provided, the dataframe passed to `fit` and
 875        `predict` should have a column with the specified condition_name
 876        containing booleans which decides when to apply seasonality.
 877
 878        Parameters
 879        ----------
 880        name: string name of the seasonality component.
 881        period: float number of days in one period.
 882        fourier_order: int number of Fourier components to use.
 883        prior_scale: optional float prior scale for this component.
 884        mode: optional 'additive' or 'multiplicative'
 885        condition_name: string name of the seasonality condition.
 886
 887        Returns
 888        -------
 889        The prophet object.
 890        """
 891        if self.history is not None:
 892            raise Exception(
 893                'Seasonality must be added prior to model fitting.')
 894        if name not in ['daily', 'weekly', 'yearly']:
 895            # Allow overwriting built-in seasonalities
 896            self.validate_column_name(name, check_seasonalities=False)
 897        if prior_scale is None:
 898            ps = self.seasonality_prior_scale
 899        else:
 900            ps = float(prior_scale)
 901        if ps <= 0:
 902            raise ValueError('Prior scale must be > 0')
 903        if fourier_order <= 0:
 904            raise ValueError('Fourier Order must be > 0')
 905        if mode is None:
 906            mode = self.seasonality_mode
 907        if mode not in ['additive', 'multiplicative']:
 908            raise ValueError('mode must be "additive" or "multiplicative"')
 909        if condition_name is not None:
 910            self.validate_column_name(condition_name)
 911        self.seasonalities[name] = {
 912            'period': period,
 913            'fourier_order': fourier_order,
 914            'prior_scale': ps,
 915            'mode': mode,
 916            'condition_name': condition_name,
 917        }
 918        return self
 919
 920    def add_country_holidays(self, country_name: str) -> Self:
 921        """Add in built-in holidays for the specified country.
 922
 923        These holidays will be included in addition to any specified on model
 924        initialization.
 925
 926        Holidays will be calculated for arbitrary date ranges in the history
 927        and future. See the online documentation for the list of countries with
 928        built-in holidays.
 929
 930        Built-in country holidays can only be set for a single country.
 931
 932        Parameters
 933        ----------
 934        country_name: Name of the country, like 'UnitedStates' or 'US'
 935
 936        Returns
 937        -------
 938        The prophet object.
 939        """
 940        if self.history is not None:
 941            raise Exception(
 942                "Country holidays must be added prior to model fitting."
 943            )
 944        # Validate names.
 945        for name in get_holiday_names(country_name):
 946            # Allow merging with existing holidays
 947            self.validate_column_name(name, check_holidays=False)
 948        # Set the holidays.
 949        if self.country_holidays is not None:
 950            logger.warning(
 951                'Changing country holidays from {country_holidays!r} to '
 952                '{country_name!r}.'
 953                .format(
 954                    country_holidays=self.country_holidays,
 955                    country_name=country_name,
 956                )
 957            )
 958        self.country_holidays = country_name
 959        return self
 960
 961    def make_all_seasonality_features(self, df: pd.DataFrame) -> tuple[
 962        pd.DataFrame,
 963        list[float],
 964        pd.DataFrame,
 965        dict[_Mode, list[str]],
 966    ]:
 967        """Dataframe with seasonality features.
 968
 969        Includes seasonality features, holiday features, and added regressors.
 970
 971        Parameters
 972        ----------
 973        df: pd.DataFrame with dates for computing seasonality features and any
 974            added regressors.
 975
 976        Returns
 977        -------
 978        pd.DataFrame with regression features.
 979        list of prior scales for each column of the features dataframe.
 980        Dataframe with indicators for which regression components correspond to
 981            which columns.
 982        Dictionary with keys 'additive' and 'multiplicative' listing the
 983            component names for each mode of seasonality.
 984        """
 985        seasonal_features = []
 986        prior_scales = []
 987        modes: dict[_Mode, list[str]] = {'additive': [], 'multiplicative': []}
 988
 989        # Seasonality features
 990        for name, props in self.seasonalities.items():
 991            features = self.make_seasonality_features(
 992                df['ds'],
 993                props['period'],
 994                props['fourier_order'],
 995                name,
 996            )
 997            if props['condition_name'] is not None:
 998                features[~df[props['condition_name']]] = 0
 999            seasonal_features.append(features)
1000            prior_scales.extend(
1001                [props['prior_scale']] * features.shape[1])
1002            modes[props['mode']].append(name)
1003
1004        # Holiday features
1005        holidays = self.construct_holiday_dataframe(df['ds'])
1006        if len(holidays) > 0:
1007            features, holiday_priors, holiday_names = (
1008                self.make_holiday_features(df['ds'], holidays)
1009            )
1010            seasonal_features.append(features)
1011            prior_scales.extend(holiday_priors)
1012            modes[self.holidays_mode].extend(holiday_names)
1013
1014        # Additional regressors
1015        for name, props in self.extra_regressors.items():
1016            seasonal_features.append(pd.DataFrame(df[name]))
1017            prior_scales.append(props['prior_scale'])
1018            modes[props['mode']].append(name)
1019
1020        # Dummy to prevent empty X
1021        if len(seasonal_features) == 0:
1022            seasonal_features.append(
1023                pd.DataFrame({'zeros': np.zeros(df.shape[0])}))
1024            prior_scales.append(1.)
1025
1026        seasonal_features = pd.concat(seasonal_features, axis=1)
1027        component_cols, modes = self.regressor_column_matrix(
1028            seasonal_features, modes
1029        )
1030        return seasonal_features, prior_scales, component_cols, modes
1031
1032    def regressor_column_matrix(
1033        self,
1034        seasonal_features: pd.DataFrame,
1035        modes: dict[_Mode, list[str]],
1036    ) -> tuple[pd.DataFrame, dict[_Mode, list[str]]]:
1037        """Dataframe indicating which columns of the feature matrix correspond
1038        to which seasonality/regressor components.
1039
1040        Includes combination components, like 'additive_terms'. These
1041        combination components will be added to the 'modes' input.
1042
1043        Parameters
1044        ----------
1045        seasonal_features: Constructed seasonal features dataframe
1046        modes: Dictionary with keys 'additive' and 'multiplicative' listing the
1047            component names for each mode of seasonality.
1048
1049        Returns
1050        -------
1051        component_cols: A binary indicator dataframe with columns seasonal
1052            components and rows columns in seasonal_features. Entry is 1 if
1053            that columns is used in that component.
1054        modes: Updated input with combination components.
1055        """
1056        components = pd.DataFrame({
1057            'col': np.arange(seasonal_features.shape[1]),
1058            'component': [
1059                x.split('_delim_')[0] for x in seasonal_features.columns
1060            ],
1061        })
1062        # Add total for holidays
1063        if self.train_holiday_names is not None:
1064            components = self.add_group_component(
1065                components, 'holidays', self.train_holiday_names.unique())
1066        # Add totals additive and multiplicative components, and regressors
1067        for mode in ('additive', 'multiplicative'):
1068            components = self.add_group_component(
1069                components, mode + '_terms', modes[mode]
1070            )
1071            regressors_by_mode = [
1072                r for r, props in self.extra_regressors.items()
1073                if props['mode'] == mode
1074            ]
1075            components = self.add_group_component(
1076                components, 'extra_regressors_' + mode, regressors_by_mode)
1077            # Add combination components to modes
1078            modes[mode].append(mode + '_terms')
1079            modes[mode].append('extra_regressors_' + mode)
1080        # After all of the additive/multiplicative groups have been added,
1081        modes[self.holidays_mode].append('holidays')
1082        # Convert to a binary matrix
1083        component_cols = pd.crosstab(
1084            components['col'], components['component'],
1085        ).sort_index(level='col')
1086        # Add columns for additive and multiplicative terms, if missing
1087        for name in ['additive_terms', 'multiplicative_terms']:
1088            if name not in component_cols:
1089                component_cols[name] = 0
1090        # Remove the placeholder
1091        component_cols.drop('zeros', axis=1, inplace=True, errors='ignore')
1092        # Validation
1093        if (max(component_cols['additive_terms']
1094            + component_cols['multiplicative_terms']) > 1):
1095            raise Exception('A bug occurred in seasonal components.')
1096        # Compare to the training, if set.
1097        if self.train_component_cols is not None:
1098            component_cols = component_cols[self.train_component_cols.columns]
1099            if not component_cols.equals(self.train_component_cols):
1100                raise Exception('A bug occurred in constructing regressors.')
1101        return component_cols, modes
1102
1103    def add_group_component(
1104        self,
1105        components: pd.DataFrame,
1106        name: str,
1107        group: list[str] | np.ndarray,
1108    ) -> pd.DataFrame:
1109        """Adds a component with given name that contains all of the components
1110        in group.
1111
1112        Parameters
1113        ----------
1114        components: Dataframe with components.
1115        name: Name of new group component.
1116        group: List of components that form the group.
1117
1118        Returns
1119        -------
1120        Dataframe with components.
1121        """
1122        new_comp = components[components['component'].isin(set(group))].copy()
1123        group_cols = new_comp['col'].unique()
1124        if len(group_cols) > 0:
1125            new_comp = pd.DataFrame({'col': group_cols, 'component': name})
1126            components = pd.concat([components, new_comp], ignore_index=True)
1127        return components
1128
1129    def parse_seasonality_args(
1130        self,
1131        name: str,
1132        arg: Literal['auto'] | int,
1133        auto_disable: bool,
1134        default_order: int,
1135    ) -> int:
1136        """Get number of fourier components for built-in seasonalities.
1137
1138        Parameters
1139        ----------
1140        name: string name of the seasonality component.
1141        arg: 'auto', True, False, or number of fourier components as provided.
1142        auto_disable: bool if seasonality should be disabled when 'auto'.
1143        default_order: int default fourier order
1144
1145        Returns
1146        -------
1147        Number of fourier components, or 0 for disabled.
1148        """
1149        if arg == 'auto':
1150            fourier_order = 0
1151            if name in self.seasonalities:
1152                logger.info(
1153                    'Found custom seasonality named {name!r}, disabling '
1154                    'built-in {name!r} seasonality.'.format(name=name)
1155                )
1156            elif auto_disable:
1157                logger.info(
1158                    'Disabling {name} seasonality. Run prophet with '
1159                    '{name}_seasonality=True to override this.'
1160                    .format(name=name)
1161                )
1162            else:
1163                fourier_order = default_order
1164        elif arg is True:
1165            fourier_order = default_order
1166        elif arg is False:
1167            fourier_order = 0
1168        else:
1169            fourier_order = int(arg)
1170        return fourier_order
1171
1172    def set_auto_seasonalities(self) -> None:
1173        """Set seasonalities that were left on auto.
1174
1175        Turns on yearly seasonality if there is >=2 years of history.
1176        Turns on weekly seasonality if there is >=2 weeks of history, and the
1177        spacing between dates in the history is <7 days.
1178        Turns on daily seasonality if there is >=2 days of history, and the
1179        spacing between dates in the history is <1 day.
1180        """
1181        history = cast(pd.DataFrame, self.history)
1182        first = history['ds'].min()
1183        last = history['ds'].max()
1184        dt = history['ds'].diff()
1185        min_dt = dt.iloc[cast(np.ndarray, dt.values).nonzero()[0]].min()
1186
1187        # Yearly seasonality
1188        yearly_disable = last - first < pd.Timedelta(days=730)
1189        fourier_order = self.parse_seasonality_args(
1190            'yearly', self.yearly_seasonality, yearly_disable, 10)
1191        if fourier_order > 0 and yearly_disable:
1192            logger.warning(
1193                'Yearly seasonality is enabled with less than 730 days '
1194                '(approximately 2 years) of history. The model may be '
1195                'under-identified, and the trend/seasonality decomposition '
1196                'can be unstable and dependent on the Prophet/Stan version. '
1197                'Consider disabling yearly seasonality or providing more '
1198                'history.'
1199            )
1200        if fourier_order > 0:
1201            self.seasonalities['yearly'] = {
1202                'period': 365.25,
1203                'fourier_order': fourier_order,
1204                'prior_scale': self.seasonality_prior_scale,
1205                'mode': self.seasonality_mode,
1206                'condition_name': None
1207            }
1208
1209        # Weekly seasonality
1210        weekly_disable = ((last - first < pd.Timedelta(weeks=2)) or
1211                          (min_dt >= pd.Timedelta(weeks=1)))  # pyrefly:ignore[unsupported-operation]
1212        fourier_order = self.parse_seasonality_args(
1213            'weekly', self.weekly_seasonality, weekly_disable, 3)
1214        if fourier_order > 0:
1215            self.seasonalities['weekly'] = {
1216                'period': 7,
1217                'fourier_order': fourier_order,
1218                'prior_scale': self.seasonality_prior_scale,
1219                'mode': self.seasonality_mode,
1220                'condition_name': None
1221            }
1222
1223        # Daily seasonality
1224        daily_disable = ((last - first < pd.Timedelta(days=2)) or
1225                         (min_dt >= pd.Timedelta(days=1)))  # pyrefly:ignore[unsupported-operation]
1226        fourier_order = self.parse_seasonality_args(
1227            'daily', self.daily_seasonality, daily_disable, 4)
1228        if fourier_order > 0:
1229            self.seasonalities['daily'] = {
1230                'period': 1,
1231                'fourier_order': fourier_order,
1232                'prior_scale': self.seasonality_prior_scale,
1233                'mode': self.seasonality_mode,
1234                'condition_name': None
1235            }
1236
1237    @staticmethod
1238    def linear_growth_init(df: pd.DataFrame) -> tuple[float, float]:
1239        """Initialize linear growth.
1240
1241        Provides a strong initialization for linear growth by calculating the
1242        growth and offset parameters that pass the function through the first
1243        and last points in the time series.
1244
1245        Parameters
1246        ----------
1247        df: pd.DataFrame with columns ds (date), y_scaled (scaled time series),
1248            and t (scaled time).
1249
1250        Returns
1251        -------
1252        A tuple (k, m) with the rate (k) and offset (m) of the linear growth
1253        function.
1254        """
1255        i0, i1 = cast(int, df['ds'].idxmin()), cast(int, df['ds'].idxmax())
1256        T = df['t'].iloc[i1] - df['t'].iloc[i0]
1257        k = (df['y_scaled'].iloc[i1] - df['y_scaled'].iloc[i0]) / T
1258        m = df['y_scaled'].iloc[i0] - k * df['t'].iloc[i0]
1259        return (k, m)
1260
1261    @staticmethod
1262    def logistic_growth_init(df: pd.DataFrame) -> tuple[float, float]:
1263        """Initialize logistic growth.
1264
1265        Provides a strong initialization for logistic growth by calculating the
1266        growth and offset parameters that pass the function through the first
1267        and last points in the time series.
1268
1269        Parameters
1270        ----------
1271        df: pd.DataFrame with columns ds (date), cap_scaled (scaled capacity),
1272            y_scaled (scaled time series), and t (scaled time).
1273
1274        Returns
1275        -------
1276        A tuple (k, m) with the rate (k) and offset (m) of the logistic growth
1277        function.
1278        """
1279        i0, i1 = cast(int, df['ds'].idxmin()), cast(int, df['ds'].idxmax())
1280        T = df['t'].iloc[i1] - df['t'].iloc[i0]
1281
1282        # Force valid values, in case y > cap or y < 0
1283        C0 = df['cap_scaled'].iloc[i0]
1284        C1 = df['cap_scaled'].iloc[i1]
1285        y0 = max(0.01 * C0, min(0.99 * C0, df['y_scaled'].iloc[i0]))
1286        y1 = max(0.01 * C1, min(0.99 * C1, df['y_scaled'].iloc[i1]))
1287
1288        r0 = C0 / y0
1289        r1 = C1 / y1
1290
1291        if abs(r0 - r1) <= 0.01:
1292            r0 = 1.05 * r0
1293
1294        L0 = np.log(r0 - 1)
1295        L1 = np.log(r1 - 1)
1296
1297        # Initialize the offset
1298        m = L0 * T / (L0 - L1)
1299        # And the rate
1300        k = (L0 - L1) / T
1301        return (k, m)
1302
1303    @staticmethod
1304    def flat_growth_init(df: pd.DataFrame) -> tuple[float, float]:
1305        """Initialize flat growth.
1306
1307        Provides a strong initialization for flat growth. Sets the growth to 0
1308        and offset parameter as mean of history y_scaled values.
1309
1310        Parameters
1311        ----------
1312        df: pd.DataFrame with columns ds (date), y_scaled (scaled time series),
1313            and t (scaled time).
1314
1315        Returns
1316        -------
1317        A tuple (k, m) with the rate (k) and offset (m) of the linear growth
1318        function.
1319        """
1320        k = 0
1321        m = df['y_scaled'].mean()
1322        return k, m
1323
1324    def preprocess(self, df: pd.DataFrame, **kwargs: Any) -> ModelInputData:
1325        """
1326        Reformats historical data, standardizes y and extra regressors, sets seasonalities and changepoints.
1327
1328        Saves the preprocessed data to the instantiated object, and also returns the relevant components
1329        as a ModelInputData object.
1330        """
1331        if ('ds' not in df) or ('y' not in df):
1332            raise ValueError(
1333                'Dataframe must have columns "ds" and "y" with the dates and '
1334                'values respectively.'
1335            )
1336        history = df[df['y'].notnull()].copy()
1337        if history.shape[0] < 2:
1338            raise ValueError('Dataframe has less than 2 non-NaN rows.')
1339        self.history_dates = pd.to_datetime(pd.Series(df['ds'].unique(), name='ds')).sort_values()
1340
1341        self.history = self.setup_dataframe(history, initialize_scales=True)
1342        self.set_auto_seasonalities()
1343        seasonal_features, prior_scales, component_cols, modes = (
1344            self.make_all_seasonality_features(self.history))
1345        self.train_component_cols = component_cols
1346        self.component_modes = modes
1347        self.fit_kwargs = deepcopy(kwargs)
1348
1349        self.set_changepoints()
1350
1351        if self.growth in ['linear', 'flat']:
1352            cap = np.zeros(self.history.shape[0])
1353        else:
1354            cap = self.history['cap_scaled']
1355
1356        return ModelInputData(
1357            T=self.history.shape[0],
1358            S=len(cast(npt.NDArray[np.float64], self.changepoints_t)),
1359            K=seasonal_features.shape[1],
1360            tau=self.changepoint_prior_scale,
1361            trend_indicator=TrendIndicator[self.growth.upper()].value,
1362            y=self.history['y_scaled'],
1363            t=self.history['t'],
1364            t_change=cast(npt.NDArray[np.float64], self.changepoints_t),
1365            X=seasonal_features,
1366            sigmas=prior_scales,
1367            s_a=component_cols['additive_terms'],
1368            s_m=component_cols['multiplicative_terms'],
1369            cap=cap,
1370        )
1371
1372    def calculate_initial_params(self, num_total_regressors: int) -> ModelParams:
1373        """
1374        Calculates initial parameters for the model based on the preprocessed history.
1375
1376        Parameters
1377        ----------
1378        num_total_regressors: the count of seasonality fourier components plus holidays plus extra regressors.
1379        """
1380        history = cast(pd.DataFrame, self.history)
1381        if self.growth == 'linear':
1382            k, m = self.linear_growth_init(history)
1383        elif self.growth == 'flat':
1384            k, m = self.flat_growth_init(history)
1385        else:
1386            assert self.growth == "logistic"
1387            k, m = self.logistic_growth_init(history)
1388        return ModelParams(
1389            k=k,
1390            m=m,
1391            delta=np.zeros_like(self.changepoints_t),
1392            beta=np.zeros(num_total_regressors),
1393            sigma_obs=1.0,
1394        )
1395
1396    def fit(self, df: pd.DataFrame, **kwargs: Any) -> Self:
1397        """Fit the Prophet model.
1398
1399        This sets self.params to contain the fitted model parameters. It is a
1400        dictionary parameter names as keys and the following items:
1401            k (Mx1 array): M posterior samples of the initial slope.
1402            m (Mx1 array): The initial intercept.
1403            delta (MxN array): The slope change at each of N changepoints.
1404            beta (MxK matrix): Coefficients for K seasonality features.
1405            sigma_obs (Mx1 array): Noise level.
1406        Note that M=1 if MAP estimation.
1407
1408        Parameters
1409        ----------
1410        df: pd.DataFrame containing the history. Must have columns ds (date
1411            type) and y, the time series. If self.growth is 'logistic', then
1412            df must also have a column cap that specifies the capacity at
1413            each ds.
1414        kwargs: Additional arguments passed to the optimizing or sampling
1415            functions in Stan.
1416
1417        Returns
1418        -------
1419        The fitted Prophet object.
1420        """
1421        if self.history is not None:
1422            raise Exception('Prophet object can only be fit once. '
1423                            'Instantiate a new object.')
1424
1425        model_inputs = self.preprocess(df, **kwargs)
1426        self._fit_regressor_models(df)
1427        initial_params = self.calculate_initial_params(model_inputs.K)
1428
1429        dat = dataclasses.asdict(model_inputs)
1430        stan_init = dataclasses.asdict(initial_params)
1431        stan_backend = cast(IStanBackend, self.stan_backend)
1432
1433        history = cast(pd.DataFrame, self.history)
1434        if history['y'].min() == history['y'].max() and \
1435                (self.growth == 'linear' or self.growth == 'flat'):
1436            self.params = stan_init
1437            self.params['sigma_obs'] = 1e-9
1438            for par in self.params:
1439                self.params[par] = np.array([self.params[par]])
1440        elif self.mcmc_samples > 0:
1441            self.params = stan_backend.sampling(stan_init, dat, self.mcmc_samples, **kwargs)
1442        else:
1443            self.params = stan_backend.fit(stan_init, dat, **kwargs)
1444
1445        self.stan_fit = stan_backend.stan_fit
1446        # If no changepoints were requested, replace delta with 0s
1447        if len(cast(pd.Series, self.changepoints)) == 0:
1448            # Fold delta into the base rate k
1449            self.params['k'] = (
1450                self.params['k'] + self.params['delta'].reshape(-1)
1451            )
1452            self.params['delta'] = (np.zeros(self.params['delta'].shape)
1453                                      .reshape((-1, 1)))
1454
1455        return self
1456
1457    def predict(self, df: pd.DataFrame | None = None, vectorized: bool = True) -> pd.DataFrame:
1458        """Predict using the prophet model.
1459
1460        Parameters
1461        ----------
1462        df: pd.DataFrame with dates for predictions (column ds), and capacity
1463            (column cap) if logistic growth. If not provided, predictions are
1464            made on the history.
1465        vectorized: Whether to use a vectorized method to compute uncertainty intervals. Suggest using
1466            True (the default) for much faster runtimes in most cases,
1467            except when (growth = 'logistic' and mcmc_samples > 0).
1468
1469        Returns
1470        -------
1471        A pd.DataFrame with the forecast components.
1472        """
1473        if self.history is None:
1474            raise Exception('Model has not been fit.')
1475
1476        regressor_samples = None
1477        if df is None:
1478            df = self.history.copy()
1479        else:
1480            if df.shape[0] == 0:
1481                raise ValueError('Dataframe has no rows.')
1482            df, regressor_samples = self._prepare_regressors_for_predict(
1483                df.copy(),
1484                self.uncertainty_samples if self.uncertainty_samples else None,
1485            )
1486            df = self.setup_dataframe(df)
1487
1488        if regressor_samples:
1489            standardized_samples = {}
1490            for name, samples in regressor_samples.items():
1491                props = self.extra_regressors[name]
1492                standardized_samples[name] = (samples - props['mu']) / props['std']
1493            regressor_samples = standardized_samples
1494
1495        df['trend'] = self.predict_trend(df)
1496        seasonal_components = self.predict_seasonal_components(df)
1497        if self.uncertainty_samples:
1498            intervals = self.predict_uncertainty(df, vectorized, regressor_samples)
1499        else:
1500            intervals = None
1501
1502        # Drop columns except ds, cap, floor, and trend
1503        cols = ['ds', 'trend']
1504        if 'cap' in df:
1505            cols.append('cap')
1506        if self.logistic_floor:
1507            cols.append('floor')
1508        # Add in forecast components
1509        df2 = pd.concat((df[cols], intervals, seasonal_components), axis=1)
1510        df2['yhat'] = (
1511                df2['trend'] * (1 + df2['multiplicative_terms'])
1512                + df2['additive_terms']
1513        )
1514        return df2
1515
1516    @staticmethod
1517    def piecewise_linear(
1518        t: np.ndarray,
1519        deltas: np.ndarray,
1520        k: float,
1521        m: float,
1522        changepoint_ts: np.ndarray,
1523    ) -> np.ndarray[tuple[int], np.dtype[np.float64]]:
1524        """Evaluate the piecewise linear function.
1525
1526        Parameters
1527        ----------
1528        t: np.array of times on which the function is evaluated.
1529        deltas: np.array of rate changes at each changepoint.
1530        k: Float initial rate.
1531        m: Float initial offset.
1532        changepoint_ts: np.array of changepoint times.
1533
1534        Returns
1535        -------
1536        Vector y(t).
1537        """
1538        deltas_t = (changepoint_ts[None, :] <= t[..., None]) * deltas
1539        k_t = deltas_t.sum(axis=1) + k
1540        m_t = (deltas_t * -changepoint_ts).sum(axis=1) + m
1541        return k_t * t + m_t
1542
1543    @staticmethod
1544    def piecewise_logistic(
1545        t: np.ndarray,
1546        cap: np.ndarray | pd.Series,
1547        deltas: np.ndarray,
1548        k: float,
1549        m: float,
1550        changepoint_ts: np.ndarray,
1551    ) -> np.ndarray[tuple[int], np.dtype[np.float64]]:
1552        """Evaluate the piecewise logistic function.
1553
1554        Parameters
1555        ----------
1556        t: np.array of times on which the function is evaluated.
1557        cap: np.array of capacities at each t.
1558        deltas: np.array of rate changes at each changepoint.
1559        k: Float initial rate.
1560        m: Float initial offset.
1561        changepoint_ts: np.array of changepoint times.
1562
1563        Returns
1564        -------
1565        Vector y(t).
1566        """
1567        # Compute offset changes
1568        # Ensure k and m are scalars for numpy 2.x compatibility
1569        k_scalar: float = np.asarray(k).item() if np.asarray(k).size == 1 else k
1570        m_scalar: float = np.asarray(m).item() if np.asarray(m).size == 1 else m
1571        k_cum = np.concatenate((np.atleast_1d(k_scalar), np.cumsum(deltas) + k_scalar))
1572        gammas = np.zeros(len(changepoint_ts))
1573        for i, t_s in enumerate(changepoint_ts):
1574            gammas[i] = (
1575                    (t_s - m_scalar - np.sum(gammas))
1576                    * (1 - k_cum[i] / k_cum[i + 1])  # noqa W503
1577            )
1578        # Get cumulative rate and offset at each t
1579        k_t = k_scalar * np.ones_like(t)
1580        m_t = m_scalar * np.ones_like(t)
1581        for s, t_s in enumerate(changepoint_ts):
1582            indx = t >= t_s
1583            k_t[indx] += deltas[s]
1584            m_t[indx] += gammas[s]
1585        return cap / (1 + np.exp(-k_t * (t - m_t)))
1586
1587    @staticmethod
1588    def flat_trend(
1589        t: np.ndarray,
1590        m: float,
1591    ) -> np.ndarray[tuple[int], np.dtype[np.float64]]:
1592        """Evaluate the flat trend function.
1593
1594        Parameters
1595        ----------
1596        t: np.array of times on which the function is evaluated.
1597        m: Float initial offset.
1598
1599        Returns
1600        -------
1601        Vector y(t).
1602        """
1603        m_t = m * np.ones_like(t)
1604        return m_t
1605
1606    def predict_trend(
1607        self,
1608        df: pd.DataFrame,
1609    ) -> np.ndarray[tuple[int], np.dtype[np.float64]]:
1610        """Predict trend using the prophet model.
1611
1612        Parameters
1613        ----------
1614        df: Prediction dataframe.
1615
1616        Returns
1617        -------
1618        Vector with trend on prediction dates.
1619        """
1620        k = np.nanmean(self.params['k'])
1621        m = np.nanmean(self.params['m'])
1622        deltas = np.nanmean(self.params['delta'], axis=0)
1623
1624        t = np.array(df['t'])
1625        changepoints_t = cast(npt.NDArray[np.float64], self.changepoints_t)
1626        if self.growth == 'linear':
1627            trend = self.piecewise_linear(t, deltas, k, m, changepoints_t)
1628        elif self.growth == 'logistic':
1629            cap = df['cap_scaled']
1630            trend = self.piecewise_logistic(t, cap, deltas, k, m, changepoints_t)
1631        else:
1632            # constant trend
1633            assert self.growth == "flat"
1634            trend = self.flat_trend(t, m)
1635
1636        return trend * cast(float, self.y_scale) + df['floor']
1637
1638    def predict_seasonal_components(self, df: pd.DataFrame) -> pd.DataFrame:
1639        """Predict seasonality components, holidays, and added regressors.
1640
1641        Parameters
1642        ----------
1643        df: Prediction dataframe.
1644
1645        Returns
1646        -------
1647        Dataframe with seasonal components.
1648        """
1649        seasonal_features, _, component_cols, _ = (
1650            self.make_all_seasonality_features(df)
1651        )
1652        if self.uncertainty_samples:
1653            lower_p = 100 * (1.0 - self.interval_width) / 2
1654            upper_p = 100 * (1.0 + self.interval_width) / 2
1655
1656        X = seasonal_features.values
1657        data = {}
1658        for component in component_cols.columns:
1659            beta_c = self.params['beta'] * component_cols[component].values
1660
1661            comp = np.matmul(X, beta_c.transpose())
1662            assert self.component_modes is not None
1663            if component in self.component_modes['additive']:
1664                comp *= self.y_scale
1665            data[component] = np.nanmean(comp, axis=1)
1666            if self.uncertainty_samples:
1667                # pyrefly:ignore[unbound-name]
1668                data[component + '_lower'] = self.percentile(comp, lower_p, axis=1)
1669                # pyrefly:ignore[unbound-name]
1670                data[component + '_upper'] = self.percentile(comp, upper_p, axis=1)
1671        return pd.DataFrame(data)
1672
1673    def predict_uncertainty(
1674        self,
1675        df: pd.DataFrame,
1676        vectorized: bool,
1677        regressor_samples: dict[str, npt.NDArray[np.float64]] | None = None,
1678    ) -> pd.DataFrame:
1679        """Prediction intervals for yhat and trend.
1680
1681        Parameters
1682        ----------
1683        df: Prediction dataframe.
1684        vectorized: Whether to use a vectorized method for generating future draws.
1685        regressor_samples: Optional draws for regressor values (already standardized)
1686            aligned with df.
1687
1688        Returns
1689        -------
1690        Dataframe with uncertainty intervals.
1691        """
1692        sim_values = self.sample_posterior_predictive(df, vectorized, regressor_samples)
1693
1694        lower_p = 100 * (1.0 - self.interval_width) / 2
1695        upper_p = 100 * (1.0 + self.interval_width) / 2
1696
1697        series = {}
1698        for key in ['yhat', 'trend']:
1699            series['{}_lower'.format(key)] = self.percentile(
1700                sim_values[key], lower_p, axis=1)
1701            series['{}_upper'.format(key)] = self.percentile(
1702                sim_values[key], upper_p, axis=1)
1703
1704        return pd.DataFrame(series)
1705
1706    def sample_posterior_predictive(
1707        self,
1708        df: pd.DataFrame,
1709        vectorized: bool,
1710        regressor_samples: dict[str, np.ndarray] | None = None,
1711    ) -> dict[str, npt.NDArray[np.float64]]:
1712        """Prophet posterior predictive samples.
1713
1714        Parameters
1715        ----------
1716        df: Prediction dataframe.
1717        vectorized: Whether to use a vectorized method to generate future draws.
1718        regressor_samples: Optional draws for regressors (standardized) keyed by
1719            regressor name. If provided, regressor uncertainty will be propagated
1720            through forecast draws.
1721
1722        Returns
1723        -------
1724        Dictionary with posterior predictive samples for the forecast yhat and
1725        for the trend component.
1726        """
1727        n_iterations = self.params['k'].shape[0]
1728        samp_per_iter = max(1, int(np.ceil(
1729            self.uncertainty_samples / float(n_iterations)
1730        )))
1731        # Generate seasonality features once so we can re-use them.
1732        seasonal_features, _, component_cols, _ = (
1733            self.make_all_seasonality_features(df)
1734        )
1735        sim_values = {'yhat': [], 'trend': []}
1736        regressor_positions: dict[str, int] = {}
1737        regressor_draw_counts: dict[str, int] = {}
1738        if regressor_samples:
1739            vectorized = False
1740            for name in regressor_samples:
1741                if name not in seasonal_features.columns:
1742                    continue
1743                pos = seasonal_features.columns.get_loc(name)
1744                if not isinstance(pos, (int, np.integer)):
1745                    raise ValueError(
1746                        f"Expected a unique column position for regressor '{name}'"
1747                    )
1748                regressor_positions[name] = int(pos)
1749                regressor_draw_counts[name] = regressor_samples[name].shape[1]
1750        sample_counter = 0
1751        for i in range(n_iterations):
1752            if vectorized:
1753                sims = self.sample_model_vectorized(
1754                    df=df,
1755                    seasonal_features=seasonal_features,
1756                    iteration=i,
1757                    s_a=component_cols['additive_terms'],
1758                    s_m=component_cols['multiplicative_terms'],
1759                    n_samples=samp_per_iter
1760                )
1761            else:
1762                sims = []
1763                for _ in range(samp_per_iter):
1764                    seasonal_features_sample = seasonal_features
1765                    if regressor_samples:
1766                        seasonal_features_sample = seasonal_features.copy()
1767                        for name, pos in regressor_positions.items():
1768                            draw_count = regressor_draw_counts[name]
1769                            use_idx = sample_counter % draw_count
1770                            seasonal_features_sample.iloc[:, int(pos)] = (
1771                                regressor_samples[name][:, use_idx]
1772                            )
1773                    sims.append(
1774                        self.sample_model(
1775                            df=df,
1776                            seasonal_features=seasonal_features_sample,
1777                            iteration=i,
1778                            s_a=component_cols['additive_terms'],
1779                            s_m=component_cols['multiplicative_terms'],
1780                        )
1781                    )
1782                    sample_counter += 1
1783            for key in sim_values:
1784                for sim in sims:
1785                    sim_values[key].append(sim[key])
1786
1787        return {k: np.column_stack(v) for k, v in sim_values.items()}
1788
1789    def sample_model(
1790        self,
1791        df: pd.DataFrame,
1792        seasonal_features: pd.DataFrame,
1793        iteration: int,
1794        s_a: pd.Series,
1795        s_m: pd.Series,
1796    ) -> dict[str, npt.NDArray[np.float64]]:
1797        """Simulate observations from the extrapolated generative model.
1798
1799        Parameters
1800        ----------
1801        df: Prediction dataframe.
1802        seasonal_features: pd.DataFrame of seasonal features.
1803        iteration: Int sampling iteration to use parameters from.
1804        s_a: Indicator vector for additive components
1805        s_m: Indicator vector for multiplicative components
1806
1807        Returns
1808        -------
1809        Dictionary with `yhat` and `trend`, each like df['t'].
1810        """
1811        trend = self.sample_predictive_trend(df, iteration)
1812
1813        y_scale = cast(float, self.y_scale)
1814        beta = self.params['beta'][iteration]
1815        Xb_a = np.matmul(seasonal_features.values,
1816                         beta * s_a.values) * y_scale
1817        Xb_m = np.matmul(seasonal_features.values, beta * s_m.values)
1818
1819        sigma = self.params['sigma_obs'][iteration]
1820        noise = np.random.normal(0, sigma, df.shape[0]) * y_scale
1821
1822        return {
1823            'yhat': trend * (1 + Xb_m) + Xb_a + noise,
1824            'trend': trend
1825        }
1826
1827    def sample_model_vectorized(
1828        self,
1829        df: pd.DataFrame,
1830        seasonal_features: pd.DataFrame,
1831        iteration: int,
1832        s_a: pd.Series,
1833        s_m: pd.Series,
1834        n_samples: int,
1835    ) -> list[dict[str, np.ndarray]]:
1836        """Simulate observations from the extrapolated generative model. Vectorized version of sample_model().
1837
1838        Parameters
1839        ----------
1840        df: Prediction dataframe.
1841        seasonal_features: pd.DataFrame of seasonal features.
1842        iteration: Int sampling iteration to use parameters from.
1843        s_a: Indicator vector for additive components.
1844        s_m: Indicator vector for multiplicative components.
1845        n_samples: Number of future paths of the trend to simulate.
1846
1847        Returns
1848        -------
1849        List (length n_samples) of dictionaries with arrays for trend and yhat, each ordered like df['t'].
1850        """
1851        # Get the seasonality and regressor components, which are deterministic per iteration
1852        beta = self.params['beta'][iteration]
1853        Xb_a = np.matmul(seasonal_features.values,
1854                        beta * s_a.values) * self.y_scale
1855        Xb_m = np.matmul(seasonal_features.values, beta * s_m.values)
1856        # Get the future trend, which is stochastic per iteration
1857        trends = self.sample_predictive_trend_vectorized(df, n_samples, iteration)  # already on the same scale as the actual data
1858        sigma = self.params['sigma_obs'][iteration]
1859        noise_terms = np.random.normal(0, sigma, trends.shape) * cast(float, self.y_scale)
1860
1861        simulations = []
1862        for trend, noise in zip(trends, noise_terms):
1863            simulations.append({
1864                'yhat': trend * (1 + Xb_m) + Xb_a + noise,
1865                'trend': trend
1866            })
1867        return simulations
1868
1869    def sample_predictive_trend(self, df: pd.DataFrame, iteration: int) -> np.ndarray:
1870        """Simulate the trend using the extrapolated generative model.
1871
1872        Parameters
1873        ----------
1874        df: Prediction dataframe.
1875        iteration: Int sampling iteration to use parameters from.
1876
1877        Returns
1878        -------
1879        np.array of simulated trend over df['t'].
1880        """
1881        k = self.params['k'][iteration]
1882        m = self.params['m'][iteration]
1883        deltas = self.params['delta'][iteration]
1884
1885        t = np.array(df['t'])
1886        T = t.max()
1887
1888        # New changepoints from a Poisson process with rate S on [1, T]
1889        changepoints_t = cast(np.ndarray, self.changepoints_t)
1890        if T > 1:
1891            S = len(changepoints_t)
1892            n_changes = np.random.poisson(S * (T - 1))
1893        else:
1894            n_changes = 0
1895        if n_changes > 0:
1896            changepoint_ts_new = 1 + np.random.rand(n_changes) * (T - 1)
1897            changepoint_ts_new.sort()
1898        else:
1899            changepoint_ts_new = []
1900
1901        # Get the empirical scale of the deltas, plus epsilon to avoid NaNs.
1902        lambda_ = np.mean(np.abs(deltas)) + 1e-8
1903
1904        # Sample deltas
1905        deltas_new = np.random.laplace(0, lambda_, n_changes)
1906
1907        # Prepend the times and deltas from the history
1908        changepoint_ts = np.concatenate((changepoints_t,
1909                                         changepoint_ts_new))
1910        deltas = np.concatenate((deltas, deltas_new))
1911
1912        if self.growth == 'linear':
1913            trend = self.piecewise_linear(t, deltas, k, m, changepoint_ts)
1914        elif self.growth == 'logistic':
1915            cap = df['cap_scaled']
1916            trend = self.piecewise_logistic(t, cap, deltas, k, m,
1917                                            changepoint_ts)
1918        else:
1919            assert self.growth == "flat"
1920            trend = self.flat_trend(t, m)
1921
1922        return trend * cast(float, self.y_scale) + df['floor']
1923
1924    def sample_predictive_trend_vectorized(
1925        self,
1926        df: pd.DataFrame,
1927        n_samples: int,
1928        iteration: int = 0,
1929    ) -> np.ndarray[tuple[int, int], np.dtype[np.float64]]:
1930        """Sample draws of the future trend values. Vectorized version of sample_predictive_trend().
1931
1932        Parameters
1933        ----------
1934        df: Prediction dataframe.
1935        iteration: Int sampling iteration to use parameters from.
1936        n_samples: Number of future paths of the trend to simulate.
1937
1938        Returns
1939        -------
1940        Draws of the trend values with shape (n_samples, len(df)). Values are on the scale of the original data.
1941        """
1942        deltas = self.params["delta"][iteration]
1943        m = self.params["m"][iteration]
1944        k = self.params["k"][iteration]
1945        changepoints_t = cast(np.ndarray, self.changepoints_t)
1946        if self.growth == "linear":
1947            expected = self.piecewise_linear(
1948                cast(np.ndarray, df["t"].values), deltas, k, m, changepoints_t
1949            )
1950        elif self.growth == "logistic":
1951            expected = self.piecewise_logistic(
1952                cast(np.ndarray, df["t"].values),
1953                cast(np.ndarray, df["cap_scaled"].values),
1954                deltas,
1955                k,
1956                m,
1957                changepoints_t,
1958            )
1959        elif self.growth == "flat":
1960            expected = self.flat_trend(cast(np.ndarray, df["t"].values), m)
1961        else:
1962            raise NotImplementedError
1963        uncertainty = self._sample_uncertainty(df, n_samples, iteration)
1964        return (
1965            (np.tile(expected, (n_samples, 1)) + uncertainty) * cast(float, self.y_scale) +
1966            np.tile(cast(np.ndarray, df["floor"].values), (n_samples, 1))
1967        )
1968
1969    def _sample_uncertainty(
1970        self,
1971        df: pd.DataFrame,
1972        n_samples: int,
1973        iteration: int = 0,
1974    ) -> np.ndarray[tuple[int, int], np.dtype[np.float64]]:
1975        """Sample draws of future trend changes, vectorizing as much as possible.
1976
1977        Parameters
1978        ----------
1979        df: pd.DataFrame with columns `t` (time scaled to the model context), trend, and cap.
1980        n_samples: Number of future paths of the trend to simulate
1981        iteration: The iteration of the parameter set to use. Default 0, the first iteration.
1982
1983        Returns
1984        -------
1985        Draws of the trend changes with shape (n_samples, len(df)). Values are standardized.
1986        """
1987        # handle only historic data
1988        if df["t"].max() <= 1:
1989            # there is no trend uncertainty in historic trends
1990            uncertainties = np.zeros((n_samples, len(df)))
1991        else:
1992            future_df = df.loc[df["t"] > 1]
1993            n_length = len(future_df)
1994            # handle 1 length futures by using history
1995            if n_length > 1:
1996                single_diff = np.diff(future_df["t"]).mean()
1997            else:
1998                single_diff = np.diff(cast(pd.DataFrame, self.history)["t"]).mean()
1999            change_likelihood = len(cast(np.ndarray, self.changepoints_t)) * single_diff
2000            deltas = self.params["delta"][iteration]
2001            m = self.params["m"][iteration]
2002            k = self.params["k"][iteration]
2003            mean_delta = np.mean(np.abs(deltas)) + 1e-8
2004            if self.growth == "linear":
2005                mat = self._make_trend_shift_matrix(mean_delta, change_likelihood, n_length, n_samples=n_samples)
2006                uncertainties = mat.cumsum(axis=1).cumsum(axis=1)  # from slope changes to actual values
2007                uncertainties *= single_diff  # scaled by the actual meaning of the slope
2008            elif self.growth == "logistic":
2009                mat = self._make_trend_shift_matrix(mean_delta, change_likelihood, n_length, n_samples=n_samples)
2010                uncertainties = self._logistic_uncertainty(
2011                    mat=mat,
2012                    deltas=deltas,
2013                    k=k,
2014                    m=m,
2015                    cap=cast(np.ndarray, future_df["cap_scaled"].values),
2016                    t_time=cast(np.ndarray, future_df["t"].values),
2017                    n_length=n_length,
2018                    single_diff=single_diff,
2019                )
2020            elif self.growth == "flat":
2021                # no trend uncertainty when there is no growth
2022                uncertainties = np.zeros((n_samples, n_length))
2023            else:
2024                raise NotImplementedError
2025            # handle past included in dataframe
2026            if df["t"].min() <= 1:
2027                past_uncertainty = np.zeros((n_samples, np.sum(df["t"] <= 1)))
2028                uncertainties = np.concatenate([past_uncertainty, uncertainties], axis=1)
2029        return uncertainties
2030
2031    @staticmethod
2032    def _make_trend_shift_matrix(
2033        mean_delta: float, likelihood: float, future_length: int, n_samples: int
2034    ) -> np.ndarray[tuple[int, int], np.dtype[np.float64]]:
2035        """
2036        Creates a matrix of random trend shifts based on historical likelihood and size of shifts.
2037        Can be used for either linear or logistic trend shifts.
2038        Each row represents a different sample of a possible future, and each column is a time step into the future.
2039        """
2040        # create a bool matrix of where these trend shifts should go
2041        bool_slope_change = np.random.uniform(size=(n_samples, future_length)) < likelihood
2042        shift_values = np.random.laplace(0, mean_delta, size=bool_slope_change.shape)
2043        mat = shift_values * bool_slope_change
2044        n_mat = np.hstack([np.zeros((len(mat), 1)), mat])[:, :-1]
2045        mat = (n_mat + mat) / 2
2046        return mat
2047
2048    @staticmethod
2049    def _make_historical_mat_time(
2050        deltas: npt.ArrayLike,
2051        changepoints_t: np.ndarray,
2052        t_time: npt.ArrayLike,
2053        n_row: int = 1,
2054        single_diff: float | None = None,
2055    ) -> tuple[
2056        np.ndarray[tuple[int], np.dtype[np.float64]],
2057        np.ndarray[tuple[int], np.dtype[np.float64]],
2058    ]:
2059        """
2060        Creates a matrix of slope-deltas where these changes occured in training data according to the trained prophet obj
2061        """
2062        if single_diff is None:
2063            single_diff = np.diff(t_time).mean()
2064        prev_time = np.arange(0, 1 + single_diff, single_diff)
2065        idxs = []
2066        for changepoint in changepoints_t:
2067            idxs.append(np.where(prev_time > changepoint)[0][0])
2068        prev_deltas = np.zeros(len(prev_time))
2069        prev_deltas[idxs] = deltas
2070        prev_deltas = np.repeat(prev_deltas.reshape(1, -1), n_row, axis=0)
2071        return prev_deltas, prev_time
2072
2073    def _logistic_uncertainty(
2074        self,
2075        mat: np.ndarray,
2076        deltas: np.ndarray,
2077        k: float,
2078        m: float,
2079        cap: np.ndarray,
2080        t_time: np.ndarray,
2081        n_length: int,
2082        single_diff: float | None = None,
2083    ) -> np.ndarray:
2084        """
2085        Vectorizes prophet's logistic uncertainty by creating a matrix of future possible trends.
2086
2087        Parameters
2088        ----------
2089        mat: A trend shift matrix returned by _make_trend_shift_matrix()
2090        deltas: The size of the trend changes at each changepoint, estimated by the model
2091        k: Float initial rate.
2092        m: Float initial offset.
2093        cap: np.array of capacities at each t.
2094        t_time: The values of t in the model context (i.e. scaled so that anything > 1 represents the future)
2095        n_length: For each path, the number of future steps to simulate
2096        single_diff: The difference between each t step in the model context. Default None, inferred
2097            from t_time.
2098
2099        Returns
2100        -------
2101        A numpy array with shape (n_samples, n_length), representing the width of the uncertainty interval
2102        (standardized, not on the same scale as the actual data values) around 0.
2103        """
2104
2105        def ffill(arr):
2106            mask = arr == 0
2107            idx = np.where(~mask, np.arange(mask.shape[1]), 0)
2108            np.maximum.accumulate(idx, axis=1, out=idx)
2109            return arr[np.arange(idx.shape[0])[:, None], idx]
2110
2111        # for logistic growth we need to evaluate the trend all the way from the start of the train item
2112        historical_mat, historical_time = self._make_historical_mat_time(
2113            deltas, cast(np.ndarray, self.changepoints_t), t_time, len(mat), single_diff
2114        )
2115        mat = np.concatenate([historical_mat, mat], axis=1)
2116        full_t_time = np.concatenate([historical_time, t_time])
2117
2118        # apply logistic growth logic on the slope changes
2119        k_cum = np.concatenate((np.ones((mat.shape[0], 1)) * k, np.where(mat, np.cumsum(mat, axis=1) + k, 0)), axis=1)
2120        k_cum_b = ffill(k_cum)
2121        gammas = np.zeros_like(mat)
2122        for i in range(mat.shape[1]):
2123            x = full_t_time[i] - m - np.sum(gammas[:, :i], axis=1)
2124            ks = 1 - k_cum_b[:, i] / k_cum_b[:, i + 1]
2125            gammas[:, i] = x * ks
2126        # the data before the -n_length is the historical values, which are not needed, so cut the last n_length
2127        k_t = (mat.cumsum(axis=1) + k)[:, -n_length:]
2128        m_t = (gammas.cumsum(axis=1) + m)[:, -n_length:]
2129        sample_trends = cap / (1 + np.exp(-k_t * (t_time - m_t)))
2130        # remove the mean because we only need width of the uncertainty centered around 0
2131        # we will add the width to the main forecast - yhat (which is the mean) - later
2132        return sample_trends - sample_trends.mean(axis=0)
2133
2134    def predictive_samples(self, df: pd.DataFrame, vectorized: bool = True) -> dict[str, npt.NDArray[np.float64]]:
2135        """Sample from the posterior predictive distribution. Returns samples
2136        for the main estimate yhat, and for the trend component. The shape of
2137        each output will be (nforecast x nsamples), where nforecast is the
2138        number of points being forecasted (the number of rows in the input
2139        dataframe) and nsamples is the number of posterior samples drawn.
2140        This is the argument `uncertainty_samples` in the Prophet constructor,
2141        which defaults to 1000.
2142
2143        Parameters
2144        ----------
2145        df: Dataframe with dates for predictions (column ds), and capacity
2146            (column cap) if logistic growth.
2147        vectorized: Whether to use a vectorized method to compute possible draws. Suggest using
2148            True (the default) for much faster runtimes in most cases,
2149            except when (growth = 'logistic' and mcmc_samples > 0).
2150
2151        Returns
2152        -------
2153        Dictionary with keys "trend" and "yhat" containing
2154        posterior predictive samples for that component.
2155        """
2156        regressor_samples = None
2157        df, regressor_samples = self._prepare_regressors_for_predict(
2158            df.copy(),
2159            self.uncertainty_samples if self.uncertainty_samples else None,
2160        )
2161        df = self.setup_dataframe(df)
2162        if regressor_samples:
2163            standardized_samples = {}
2164            for name, samples in regressor_samples.items():
2165                props = self.extra_regressors[name]
2166                standardized_samples[name] = (samples - props['mu']) / props['std']
2167            regressor_samples = standardized_samples
2168        return self.sample_posterior_predictive(df, vectorized, regressor_samples)
2169
2170    def percentile(self, a: npt.ArrayLike, *args: Any, **kwargs: Any) -> np.ndarray:
2171        """
2172        We rely on np.nanpercentile in the rare instances where there
2173        are a small number of bad samples with MCMC that contain NaNs.
2174        However, since np.nanpercentile is far slower than np.percentile,
2175        we only fall back to it if the array contains NaNs. See
2176        https://github.com/facebook/prophet/issues/1310 for more details.
2177        """
2178        fn = np.nanpercentile if np.isnan(a).any() else np.percentile
2179        return fn(a, *args, **kwargs)
2180
2181    def make_future_dataframe(
2182        self, periods: int, freq: str | None = 'D', include_history: bool = True
2183    ) -> pd.DataFrame:
2184        """Simulate the trend using the extrapolated generative model.
2185
2186        Parameters
2187        ----------
2188        periods: Int number of periods to forecast forward.
2189        freq: Any valid frequency for pd.date_range, such as 'D' or 'M'.
2190        include_history: Boolean to include the historical dates in the data
2191            frame for predictions.
2192
2193        Returns
2194        -------
2195        pd.Dataframe that extends forward from the end of self.history for the
2196        requested number of periods.
2197        """
2198        if self.history_dates is None:
2199            raise Exception('Model has not been fit.')
2200        if freq is None:
2201            # taking the tail makes freq inference more reliable
2202            freq = pd.infer_freq(self.history_dates.tail(5))
2203            # returns None if inference failed
2204            if freq is None:
2205                raise Exception('Unable to infer `freq`')
2206        last_date = self.history_dates.max()
2207        dates = pd.date_range(
2208            start=last_date,
2209            periods=periods + 1,  # An extra in case we include start
2210            freq=freq)
2211        dates = dates[dates > last_date]  # Drop start if equals last_date
2212        dates = dates[:periods]  # Return correct number of periods
2213
2214        if include_history:
2215            dates = np.concatenate((np.array(self.history_dates), dates))
2216
2217        return pd.DataFrame({'ds': dates})
2218
2219    def plot(
2220        self,
2221        fcst: pd.DataFrame,
2222        ax: plt.Axes | None = None,
2223        uncertainty: bool = True,
2224        plot_cap: bool = True,
2225        xlabel: str = 'ds',
2226        ylabel: str = 'y',
2227        figsize: tuple[int, int] = (10, 6),
2228        include_legend: bool = False,
2229    ) -> plt.Figure:
2230        """Plot the Prophet forecast.
2231
2232        Parameters
2233        ----------
2234        fcst: pd.DataFrame output of self.predict.
2235        ax: Optional matplotlib axes on which to plot.
2236        uncertainty: Optional boolean to plot uncertainty intervals.
2237        plot_cap: Optional boolean indicating if the capacity should be shown
2238            in the figure, if available.
2239        xlabel: Optional label name on X-axis
2240        ylabel: Optional label name on Y-axis
2241        figsize: Optional tuple width, height in inches.
2242        include_legend: Optional boolean to add legend to the plot.
2243
2244        Returns
2245        -------
2246        A matplotlib figure.
2247        """
2248        return plot(
2249            m=self, fcst=fcst, ax=ax, uncertainty=uncertainty,
2250            plot_cap=plot_cap, xlabel=xlabel, ylabel=ylabel,
2251            figsize=figsize, include_legend=include_legend
2252        )
2253
2254    def plot_components(
2255        self,
2256        fcst: pd.DataFrame,
2257        uncertainty: bool = True,
2258        plot_cap: bool = True,
2259        weekly_start: int = 0,
2260        yearly_start: int = 0,
2261        figsize: tuple[int, int] | None = None,
2262    ) -> plt.Figure:
2263        """Plot the Prophet forecast components.
2264
2265        Will plot whichever are available of: trend, holidays, weekly
2266        seasonality, and yearly seasonality.
2267
2268        Parameters
2269        ----------
2270        fcst: pd.DataFrame output of self.predict.
2271        uncertainty: Optional boolean to plot uncertainty intervals.
2272        plot_cap: Optional boolean indicating if the capacity should be shown
2273            in the figure, if available.
2274        weekly_start: Optional int specifying the start day of the weekly
2275            seasonality plot. 0 (default) starts the week on Sunday. 1 shifts
2276            by 1 day to Monday, and so on.
2277        yearly_start: Optional int specifying the start day of the yearly
2278            seasonality plot. 0 (default) starts the year on Jan 1. 1 shifts
2279            by 1 day to Jan 2, and so on.
2280        figsize: Optional tuple width, height in inches.
2281
2282        Returns
2283        -------
2284        A matplotlib figure.
2285        """
2286        return plot_components(
2287            m=self, fcst=fcst, uncertainty=uncertainty, plot_cap=plot_cap,
2288            weekly_start=weekly_start, yearly_start=yearly_start,
2289            figsize=figsize
2290        )