【机器学习】Prophet模型训练过程代码浅读+使用pymc复现stan文件中的贝叶斯建模过程

原创

已于 2023-06-22 23:38:31 修改 · 2.5k 阅读

2 ·

CC 4.0 BY-SA版权

文章标签：

#机器学习 #python #数据挖掘

于 2023-06-22 23:37:34 首次发布

本文详细介绍了Facebook开源的时间序列预测模型Prophet的工作原理，包括数据预处理、贝叶斯建模、预测过程。在预处理阶段，涉及了数据校验、标准化和特征提取，如傅里叶变换处理季节性因素。在贝叶斯建模中，使用Stan语言进行模型定义和参数估计。预测阶段则结合趋势项和季节性因素生成预测结果。此外，文章还讨论了使用pymc复现Stan模型的实现过程。

Prophet是Facebook公司开源的一款时间序列预测模型接口，可用R语言或python进行相关时间预测操作，本文将会从它的代码上对其进行模型拟合以及预测的过程进行一个梳理，并使用python中的pymc包进行Stan文件中贝叶斯建模的复现。

一、数据预处理阶段：

（一）setup_dataframe

（二）set_auto_seasonalities

（三）make_all_seasonality_features

一、数据预处理阶段：

此时主要有3个函数进行了数据的预处理以及特征提取操作：

（一）setup_dataframe

这一函数主要是对要进行训练的DataFrame做数据检验及标准化操作：

1、检验DataFrame的时间列以及是否合法

2、检验使用add_regressor函数增加的外部变量是否合法

3、检验使用add_seasonaliy函数增加的季节性因素是否合法

4、重设index并对时间序列数据做标准化

def setup_dataframe(self, df, initialize_scales=False):
    """Prepare dataframe for fitting or predicting.

    Adds a time index and scales y. Creates auxiliary columns 't', 't_ix',
    'y_scaled', and 'cap_scaled'. These columns are used during both
    fitting and predicting.

    Parameters
    ----------
    df: pd.DataFrame with columns ds, y, and cap if logistic growth. Any
        specified additional regressors must also be present.
    initialize_scales: Boolean set scaling factors in self from df.

    Returns
    -------
    pd.DataFrame prepared for fitting or predicting.
    """
    #检验DataFrame的各列是否合法
    if 'y' in df:  # 'y' will be in training data
        df['y'] = pd.to_numeric(df['y'])
        if np.isinf(df['y'].values).any():
            raise ValueError('Found infinity in column y.')
    if df['ds'].dtype == np.int64:
        df['ds'] = df['ds'].astype(str)
    df['ds'] = pd.to_datetime(df['ds'])
    if df['ds'].dt.tz is not None:
        raise ValueError(
            'Column ds has timezone specified, which is not supported. '
            'Remove timezone.'
        )
    if df['ds'].isnull().any():
        raise ValueError('Found NaN in column ds.')
    for name in self.extra_regressors: #add_regressor会增加extra_regressors，表明那一列作为额外的回归量，可以是乘法也可以是加法
        if name not in df:
            raise ValueError(
                'Regressor {name!r} missing from dataframe'
                .format(name=name)
            )
        df[name] = pd.to_numeric(df[name])
        if df[name].isnull().any():
            raise ValueError(
                'Found NaN in column {name!r}'.format(name=name)
            )
    for props in self.seasonalities.values(): #如果没有add_seasonality的话，fit时不会有用，因为这时self.seasonalities没有初始化
        condition_name = props['condition_name']
        if condition_name is not None:
            if condition_name not in df:
                raise ValueError(
                    'Condition {condition_name!r} missing from dataframe'
                    .format(condition_name=condition_name)
                )
            if not df[condition_name].isin([True, False]).all():
                raise ValueError(
                    'Found non-boolean in column {condition_name!r}'
                    .format(condition_name=condition_name)
                )
            df[condition_name] = df[condition_name].astype('bool')

    if df.index.name == 'ds':#重设index
        df.index.name = None
    df = df.sort_values('ds')
    df = df.reset_index(drop=True)

    self.initialize_scales(initialize_scales, df) #对extra_regressor进行标准化做前置准备

    # 进行标准化
    if self.logistic_floor:
        if 'floor' not in df:
            raise ValueError('Expected column "floor".')
    else:
        df['floor'] = 0
    if self.growth == 'logistic':
        if 'cap' not in df:
            raise ValueError(
                'Capacities must be supplied for logistic growth in '
                'column "cap"'
            )
        if (df['cap'] <= df['floor']).any():
            raise ValueError(
                'cap must be greater than floor (which defaults to 0).'
            )
        df['cap_scaled'] = (df['cap'] - df['floor']) / self.y_scale

    df['t'] = (df['ds'] - self.start) / self.t_scale
    if 'y' in df:
        df['y_scaled'] = (df['y'] - df['floor']) / self.y_scale

    for name, props in self.extra_regressors.items():
        df[name] = ((df[name] - props['mu']) / props['std'])
    return df

（二）set_auto_seasonalities

这个函数的作用，是根据数据的时间间隔分别判断是否有Yearly Weekly以及Daily的周期效应并给出内置的对应参数为下一个对季节性因素分解的函数make_all_seasonality_features作准备。

def set_auto_seasonalities(self):
	"""Set seasonalities that were left on auto.

	Turns on yearly seasonality if there is >=2 years of history.
	Turns on weekly seasonality if there is >=2 weeks of history, and the
	spacing between dates in the history is <7 days.
	Turns on daily seasonality if there is >=2 days of history, and the
	spacing between dates in the history is <1 day.
	"""
	first = self.history['ds'].min()
	last = self.history['ds'].max()
	dt = self.history['ds'].diff()
	min_dt = dt.iloc[dt.values.nonzero()[0]].min()
	
	# Yearly seasonality
	yearly_disable = last - first < pd.Timedelta(days=730)
	fourier_order = self.parse_seasonality_args(
		'yearly', self.yearly_seasonality, yearly_disable, 10)
	if fourier_order > 0:
		self.seasonalities['yearly'] = {
			'period': 365.25,
			'fourier_order': fourier_order,
			'prior_scale': self.seasonality_prior_scale,
			'mode': self.seasonality_mode,
			'condition_name': None
		}

	# Weekly seasonality
	weekly_disable = ((last - first < pd.Timedelta(weeks=2)) or
					  (min_dt >= pd.Timedelta(weeks=1)))
	fourier_order = self.parse_seasonality_args(
		'weekly', self.weekly_seasonality, weekly_disable, 3)
	if fourier_order > 0:
		self.seasonalities['weekly'] = {
			'period': 7,
			'fourier_order': fourier_order,
			'prior_scale': self.seasonality_prior_scale,
			'mode': self.seasonality_mode,
			'condition_name': None
		}

	# Daily seasonality
	daily_disable = ((last - first < pd.Timedelta(days=2)) or
					 (min_dt >= pd.Timedelta(days=1)))
	fourier_order = self.parse_seasonality_args(
		'daily', self.daily_seasonality, daily_disable, 4)
	if fourier_order > 0:
		self.seasonalities['daily'] = {
			'period': 1,
			'fourier_order': fourier_order,
			'prior_scale': self.seasonality_prior_scale,
			'mode': self.seasonality_mode,
			'condition_name': None
		}

（三）make_all_seasonality_features

这个函数就是主要将set_auto_seasonalities中给出的各个季节项信息以及节假日因素作出分解并提取特征。

对于季节项信息，主要是调用fourier_series这个函数对各个不同季节信息对应的series_order与peirod做傅里叶变换，其中sieries_order代表要分解出来的序列数量，而peirod代表了季节对应的周期。

def fourier_series(
	dates: pd.Series,
	period: Union[int, float],
	series_order: int,
) -> NDArray[np.float_]:
	"""Provides Fourier series components with the specified frequency
	and order.

	Parameters
	----------
	dates: pd.Series containing timestamps.
	period: Number of days of the period.
	series_order: Number of components.

	Returns
	-------
	Matrix with seasonality features.
	"""
	if not (series_order >= 1):
		raise ValueError("series_order must be >= 1")

	# convert to days since epoch
	t = dates.to_numpy(dtype=np.int64) // NANOSECONDS_TO_SECONDS / (3600 * 24.)

	x_T = t * np.pi * 2
	fourier_components = np.empty((dates.shape[0], 2 * series_order))
	for i in range(series_order):
		c = x_T * (i + 1) / period
		fourier_components[:, 2 * i] = np.sin(c)
		fourier_components[:, (2 * i) + 1] = np.cos(c)
	return fourier_components

此处简单举个例子：

Prophet自带的Yearly参数为例，其中sieries_order为10，代表要分解出10个序列；period为365.25，表示一个周期的长度（也就是频率）为365.25。此处使用的傅里叶级数公式为：

$\sum (a_n cos(2\pi nt/p)+b_nsin(2\pi nt/p))$