Prophet 简介

Facebook 去年开源了一个时间序列预测的算法，叫做 fbprophet，它的官方网址与基本介绍来自于以下几个网站：

1. Github：https://github.com/facebook/prophet
2. 官方网址：https://facebook.github.io/prophet/
3. 论文名字与网址：Forecasting at scale，https://peerj.com/preprints/3190/

从官网的介绍来看，Facebook 所提供的 prophet 算法不仅可以处理时间序列存在一些异常值的情况，也可以处理部分缺失值的情形，还能够几乎全自动地预测时间序列未来的走势。从论文上的描述来看，这个 prophet 算法是基于时间序列分解和机器学习的拟合来做的，其中在拟合模型的时候使用了 pyStan 这个开源工具，因此能够在较快的时间内得到需要预测的结果。除此之外，为了方便统计学家，机器学习从业者等人群的使用，prophet 同时提供了 R 语言和 Python 语言的接口。从整体的介绍来看，如果是一般的商业分析或者数据分析的需求，都可以尝试使用这个开源算法来预测未来时间序列的走势。

Prophet 的算法原理

Prophet 数据的输入和输出

首先让我们来看一个常见的时间序列场景，黑色表示原始的时间序列离散点，深蓝色的线表示使用时间序列来拟合所得到的取值，而浅蓝色的线表示时间序列的一个置信区间，也就是所谓的合理的上界和下界。prophet 所做的事情就是：

1. 输入已知的时间序列的时间戳和相应的值；
2. 输入需要预测的时间序列的长度；
3. 输出未来的时间序列走势。
4. 输出结果可以提供必要的统计指标，包括拟合曲线，上界和下界等。

就一般情况而言，时间序列的离线存储格式为时间戳和值这种格式，更多的话可以提供时间序列的 ID，标签等内容。因此，离线存储的时间序列通常都是以下的形式。其中 date 指的是具体的时间戳，category 指的是某条特定的时间序列 id，value 指的是在 date 下这个 category 时间序列的取值，label 指的是人工标记的标签（’0′ 表示异常，’1‘ 表示正常，’unknown’ 表示没有标记或者人工判断不清）。

而 fbprophet 所需要的时间序列也是这种格式的，根据官网的描述，只要用 csv 文件存储两列即可，第一列的名字是 ‘ds’, 第二列的名称是 ‘y’。第一列表示时间序列的时间戳，第二列表示时间序列的取值。通过 prophet 的计算，可以计算出 yhat，yhat_lower，yhat_upper，分别表示时间序列的预测值，预测值的下界，预测值的上界。两份表格如下面的两幅图表示。

Prophet 的算法实现

在时间序列分析领域，有一种常见的分析方法叫做时间序列的分解（Decomposition of Time Series），它把时间序列 分成几个部分，分别是季节项 ，趋势项 ，剩余项 。也就是说对所有的 ，都有

除了加法的形式，还有乘法的形式，也就是：

以上式子等价于 。所以，有的时候在预测模型的时候，会先取对数，然后再进行时间序列的分解，就能得到乘法的形式。在 fbprophet 算法中，作者们基于这种方法进行了必要的改进和优化。

一般来说，在实际生活和生产环节中，除了季节项，趋势项，剩余项之外，通常还有节假日的效应。所以，在 prophet 算法里面，作者同时考虑了以上四项，也就是：

其中 表示趋势项，它表示时间序列在非周期上面的变化趋势； 表示周期项，或者称为季节项，一般来说是以周或者年为单位； 表示节假日项，表示在当天是否存在节假日； 表示误差项或者称为剩余项。Prophet 算法就是通过拟合这几项，然后最后把它们累加起来就得到了时间序列的预测值。

趋势项模型

在 Prophet 算法里面，趋势项有两个重要的函数，一个是基于逻辑回归函数（logistic function）的，另一个是基于分段线性函数（piecewise linear function）的。

首先，我们来介绍一下基于逻辑回归的趋势项是怎么做的。

如果回顾逻辑回归函数的话，一般都会想起这样的形式： 它的导数是 并且 如果增加一些参数的话，那么逻辑回归就可以改写成：

这里的 称为曲线的最大渐近值， 表示曲线的增长率， 表示曲线的中点。当 时，恰好就是大家常见的 sigmoid 函数的形式。从 sigmoid 的函数表达式来看，它满足以下的微分方程：。

那么，如果使用分离变量法来求解微分方程 就可以得到：

.

但是在现实环境中，函数 的三个参数 不可能都是常数，而很有可能是随着时间的迁移而变化的，因此，在 Prophet 里面，作者考虑把这三个参数全部换成了随着时间而变化的函数，也就是 。

除此之外，在现实的时间序列中，曲线的走势肯定不会一直保持不变，在某些特定的时候或者有着某种潜在的周期曲线会发生变化，这种时候，就有学者会去研究变点检测，也就是所谓 change point detection。例如下面的这幅图的 就是时间序列的两个变点。

在 Prophet 里面，是需要设置变点的位置的，而每一段的趋势和走势也是会根据变点的情况而改变的。在程序里面有两种方法，一种是通过人工指定的方式指定变点的位置；另外一种是通过算法来自动选择。在默认的函数里面，Prophet 会选择 n_changepoints = 25 个变点，然后设置变点的范围是前 80%（changepoint_range），也就是在时间序列的前 80% 的区间内会设置变点。通过 forecaster.py 里面的 set_changepoints 函数可以知道，首先要看一些边界条件是否合理，例如时间序列的点数是否少于 n_changepoints 等内容；其次如果边界条件符合，那变点的位置就是均匀分布的，这一点可以通过 np.linspace 这个函数看出来。

下面假设已经放置了 个变点了，并且变点的位置是在时间戳 上，那么在这些时间戳上，我们就需要给出增长率的变化，也就是在时间戳 上发生的 change in rate。可以假设有这样一个向量： 其中 表示在时间戳 上的增长率的变化量。如果一开始的增长率我们使用 来代替的话，那么在时间戳 上的增长率就是 ，通过一个指示函数 就是

那么在时间戳 上面的增长率就是 一旦变化量 确定了，另外一个参数 也要随之确定。在这里需要把线段的边界处理好，因此通过数学计算可以得到：

所以，分段的逻辑回归增长模型就是：

其中，

在逻辑回归函数里面，有一个参数是需要提前设置的，那就是 Capacity，也就是所谓的 ，在使用 Prophet 的 growth = ‘logistic’ 的时候，需要提前设置好 的取值才行。

再次，我们来介绍一下基于分段线性函数的趋势项是怎么做的。众所周知，线性函数指的是 而分段线性函数指的是在每一个子区间上，函数都是线性函数，但是在整段区间上，函数并不完全是线性的。正如下图所示，分段线性函数就是一个折线的形状。

因此，基于分段线性函数的模型形如：

其中 表示增长率（growth rate）， 表示增长率的变化量， 表示 offset parameter。而这两种方法（分段线性函数与逻辑回归函数）最大的区别就是 的设置不一样，在分段线性函数中， 注意：这与之前逻辑回归函数中的设置是不一样的。

在 prophet 的源代码中，forecast.py 这个函数里面包含了最关键的步骤，其中 piecewise_logistic 函数表示了前面所说的基于逻辑回归的增长函数，它的输入包含了 cap 这个指标，因此需要用户事先指定 capacity。而在 piecewise_linear 这个函数中，是不需要 capacity 这个指标的，因此 m = Prophet() 这个函数默认的使用 growth = ‘linear’ 这个增长函数，也可以写作 m = Prophet(growth = ‘linear’)；如果想用 growth = ‘logistic’，就要这样写：

m = Prophet(growth='logistic')
df['cap'] = 6
m.fit(df)
future = m.make_future_dataframe(periods=prediction_length, freq='min')
future['cap'] = 6

变点的选择（Changepoint Selection）

在介绍变点之前，先要介绍一下 Laplace 分布，它的概率密度函数为：

其中 表示位置参数， 表示尺度参数。Laplace 分布与正态分布有一定的差异。

在 Prophet 算法中，是需要给出变点的位置，个数，以及增长的变化率的。因此，有三个比较重要的指标，那就是

1. changepoint_range，
2. n_changepoint，
3. changepoint_prior_scale。

changepoint_range 指的是百分比，需要在前 changepoint_range 那么长的时间序列中设置变点，在默认的函数中是 changepoint_range = 0.8。n_changepoint 表示变点的个数，在默认的函数中是 n_changepoint = 25。changepoint_prior_scale 表示变点增长率的分布情况，在论文中， ，这里的 就是 change_point_scale。

在整个开源框架里面，在默认的场景下，变点的选择是基于时间序列的前 80% 的历史数据，然后通过等分的方法找到 25 个变点（change points），而变点的增长率是满足 Laplace 分布 的。因此，当 趋近于零的时候， 也是趋向于零的，此时的增长函数将变成全段的逻辑回归函数或者线性函数。这一点从 的定义可以轻易地看出。

对未来的预估（Trend Forecast Uncertainty）

从历史上长度为 的数据中，我们可以选择出 个变点，它们所对应的增长率的变化量是 。此时我们需要预测未来，因此也需要设置相应的变点的位置，从代码中看，在 forecaster.py 的 sample_predictive_trend 函数中，通过 Poisson 分布等概率分布方法找到新增的 changepoint_ts_new 的位置，然后与 changepoint_t 拼接在一起就得到了整段序列的 changepoint_ts。

changepoint_ts_new = 1 + np.random.rand(n_changes) * (T - 1)
changepoint_ts = np.concatenate((self.changepoints_t, changepoint_ts_new))

第一行代码的 1 保证了 changepoint_ts_new 里面的元素都大于 change_ts 里面的元素。除了变点的位置之外，也需要考虑 的情况。这里令 ，于是新的增长率的变化量就是按照下面的规则来选择的：当 时，

季节性趋势

几乎所有的时间序列预测模型都会考虑这个因素，因为时间序列通常会随着天，周，月，年等季节性的变化而呈现季节性的变化，也称为周期性的变化。对于周期函数而言，大家能够马上联想到的就是正弦余弦函数。而在数学分析中，区间内的周期性函数是可以通过正弦和余弦的函数来表示的：假设 是以 为周期的函数，那么它的傅立叶级数就是 。

在论文中，作者使用傅立叶级数来模拟时间序列的周期性。假设 表示时间序列的周期， 表示以年为周期， 表示以周为周期。它的傅立叶级数的形式都是：

就作者的经验而言，对于以年为周期的序列（）而言，；对于以周为周期的序列（ ）而言，。这里的参数可以形成列向量：

。

当 时，

当 时，

因此，时间序列的季节项就是： 而 的初始化是 。这里的 是通过 seasonality_prior_scale 来控制的，也就是说 seasonality_prior_scale。这个值越大，表示季节的效应越明显；这个值越小，表示季节的效应越不明显。同时，在代码里面，seasonality_mode 也对应着两种模式，分别是加法和乘法，默认是加法的形式。在开源代码中， 函数是通过 fourier_series 来构建的。

节假日效应（holidays and events）

在现实环境中，除了周末，同样有很多节假日，而且不同的国家有着不同的假期。在 Prophet 里面，通过维基百科里面对各个国家的节假日的描述，hdays.py 收集了各个国家的特殊节假日。除了节假日之外，用户还可以根据自身的情况来设置必要的假期，例如 The Super Bowl，双十一等。

由于每个节假日对时间序列的影响程度不一样，例如春节，国庆节则是七天的假期，对于劳动节等假期来说则假日较短。因此，不同的节假日可以看成相互独立的模型，并且可以为不同的节假日设置不同的前后窗口值，表示该节假日会影响前后一段时间的时间序列。用数学语言来说，对与第 个节假日来说， 表示该节假日的前后一段时间。为了表示节假日效应，我们需要一个相应的指示函数（indicator function），同时需要一个参数 来表示节假日的影响范围。假设我们有 个节假日，那么

其中 和

其中 并且该正态分布是受到 holidays_prior_scale 这个指标影响的。默认值是 10，当值越大时，表示节假日对模型的影响越大；当值越小时，表示节假日对模型的效果越小。用户可以根据自己的情况自行调整。

模型拟合（Model Fitting）

按照以上的解释，我们的时间序列已经可以通过增长项，季节项，节假日项来构建了，i.e.

下一步我们只需要拟合函数就可以了，在 Prophet 里面，作者使用了 pyStan 这个开源工具中的 L-BFGS 方法来进行函数的拟合。具体可以参考 forecast.py 里面的 stan_init 函数。

Prophet 中可以设置的参数

在 Prophet 中，用户一般可以设置以下四种参数：

1. Capacity：在增量函数是逻辑回归函数的时候，需要设置的容量值。
2. Change Points：可以通过 n_changepoints 和 changepoint_range 来进行等距的变点设置，也可以通过人工设置的方式来指定时间序列的变点。
3. 季节性和节假日：可以根据实际的业务需求来指定相应的节假日。
4. 光滑参数： changepoint_prior_scale 可以用来控制趋势的灵活度， seasonality_prior_scale 用来控制季节项的灵活度， holidays prior scale 用来控制节假日的灵活度。

如果不想设置的话，使用 Prophet 默认的参数即可。

Prophet 的实际使用

Prophet 的简单使用

因为 Prophet 所需要的两列名称是 ‘ds’ 和 ‘y’，其中，’ds’ 表示时间戳，’y’ 表示时间序列的值，因此通常来说都需要修改 pd.dataframe 的列名字。如果原来的两列名字是 ‘timestamp’ 和 ‘value’ 的话，只需要这样写：

df = df.rename(columns={'timestamp':'ds', 'value':'y'})

如果 ‘timestamp’ 是使用 unixtime 来记录的，需要修改成 YYYY-MM-DD hh:mm:ss 的形式：

df['ds'] = pd.to_datetime(df['ds'],unit='s')

在一般情况下，时间序列需要进行归一化的操作，而 pd.dataframe 的归一化操作也十分简单：

df['y'] = (df['y'] - df['y'].mean()) / (df['y'].std())

然后就可以初始化模型，然后拟合模型，并且进行时间序列的预测了。

初始化模型：m = Prophet()
拟合模型：m.fit(df)
计算预测值：periods 表示需要预测的点数，freq 表示时间序列的频率。
future = m.make_future_dataframe(periods=30, freq='min')
future.tail()
forecast = m.predict(future)

而 freq 指的是 pd.dataframe 里面的一个指标，’min’ 表示按分钟来收集的时间序列。具体参见文档：http://pandas.pydata.org/pandas-docs/stable/timeseries.html#offset-aliases

在进行了预测操作之后，通常都希望把时间序列的预测趋势画出来：

画出预测图：
m.plot(forecast)
画出时间序列的分量：
m.plot_components(forecast)

如果要画出更详细的指标，例如中间线，上下界，那么可以这样写：

x1 = forecast['ds']
y1 = forecast['yhat']
y2 = forecast['yhat_lower']
y3 = forecast['yhat_upper']
plt.plot(x1,y1)
plt.plot(x1,y2)
plt.plot(x1,y3)
plt.show()

其实 Prophet 预测的结果都放在了变量 forecast 里面，打印结果的话可以这样写：第一行是打印所有时间戳的预测结果，第二行是打印最后五个时间戳的预测结果。

print(forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']])
print(forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail())

Prophet 的参数设置

Prophet 的默认参数可以在 forecaster.py 中看到：

def __init__(
    self,
    growth='linear',
    changepoints=None,
    n_changepoints=25, 
    changepoint_range=0.8,
    yearly_seasonality='auto',
    weekly_seasonality='auto',
    daily_seasonality='auto',
    holidays=None,
    seasonality_mode='additive',
    seasonality_prior_scale=10.0,
    holidays_prior_scale=10.0,
    changepoint_prior_scale=0.05,
    mcmc_samples=0,
    interval_width=0.80,
    uncertainty_samples=1000,
):

增长函数的设置

在 Prophet 里面，有两个增长函数，分别是分段线性函数（linear）和逻辑回归函数（logistic）。而 m = Prophet() 默认使用的是分段线性函数（linear），并且如果要是用逻辑回归函数的时候，需要设置 capacity 的值，i.e. df[‘cap’] = 100，否则会出错。

m = Prophet()
m = Prophet(growth='linear')
m = Prophet(growth='logistic')

变点的设置

在 Prophet 里面，变点默认的选择方法是前 80% 的点中等距选择 25 个点作为变点，也可以通过以下方法来自行设置变点，甚至可以人为设置某些点。

m = Prophet(n_changepoints=25)
m = Prophet(changepoint_range=0.8)
m = Prophet(changepoint_prior_scale=0.05)
m = Prophet(changepoints=['2014-01-01'])

而变点的作图可以使用：

from fbprophet.plot import add_changepoints_to_plot
fig = m.plot(forecast)
a = add_changepoints_to_plot(fig.gca(), m, forecast)

周期性的设置

通常来说，可以在 Prophet 里面设置周期性，无论是按月还是周其实都是可以设置的，例如：

m = Prophet(weekly_seasonality=False)
m.add_seasonality(name='monthly', period=30.5, fourier_order=5)

m = Prophet(weekly_seasonality=True)
m.add_seasonality(name='weekly', period=7, fourier_order=3, prior_scale=0.1)

节假日的设置

有的时候，由于双十一或者一些特殊节假日，我们可以设置某些天数是节假日，并且设置它的前后影响范围，也就是 lower_window 和 upper_window。

playoffs = pd.DataFrame({
  'holiday': 'playoff',
  'ds': pd.to_datetime(['2008-01-13', '2009-01-03', '2010-01-16',
                        '2010-01-24', '2010-02-07', '2011-01-08',
                        '2013-01-12', '2014-01-12', '2014-01-19',
                        '2014-02-02', '2015-01-11', '2016-01-17',
                        '2016-01-24', '2016-02-07']),
  'lower_window': 0,
  'upper_window': 1,
})
superbowls = pd.DataFrame({
  'holiday': 'superbowl',
  'ds': pd.to_datetime(['2010-02-07', '2014-02-02', '2016-02-07']),
  'lower_window': 0,
  'upper_window': 1,
})
holidays = pd.concat((playoffs, superbowls))

m = Prophet(holidays=holidays, holidays_prior_scale=10.0)

结束语

对于商业分析等领域的时间序列，Prophet 可以进行很好的拟合和预测，但是对于一些周期性或者趋势性不是很强的时间序列，用 Prophet 可能就不合适了。但是，Prophet 提供了一种时序预测的方法，在用户不是很懂时间序列的前提下都可以使用这个工具得到一个能接受的结果。具体是否用 Prophet 则需要根据具体的时间序列来确定。

参考文献：

1. https://otexts.org/fpp2/components.html
2. https://en.wikipedia.org/wiki/Decomposition_of_time_series
3. A review of change point detection methods, CTruong, L. Oudre, N.Vayatis
4. https://github.com/facebook/prophet
5. https://facebook.github.io/prophet/

近期阅读了一篇论文《Rapid Deployment of Anomaly Detection Models for Large Number of Emerging KPI Streams》，这篇文章基于之前的 ROCKA 系统做了一些额外的工作。ROCKA 系统是用来做时间序列的实时聚类的，而这篇文章是在 ROCKA 系统的基础上增加了时间序列异常检测的功能。通常来说，时间序列异常检测可以使用有监督的方法来解决，参考 Opperentice 系统。而本篇文章使用了半监督学习的思路来解决异常检测的问题，下面来详细分析一下这篇文章的细节，本文的作者把这个系统称为 ADS（Anomaly Detection through Self-training）。

数据集的情况

在论文中，作者使用了两份数据集，分别是已经历史上的 70 条时间序列，另外还有新来的 81 条时间序列。在 ADS 系统中，历史上的 70 条时间序列被划分成 5 类，并且已经可以找出每个类的质心位置，并且每条历史上的时间序列通常来说会大于三个星期（3 weeks）。本篇论文的评价指标是 F-Score，也属于机器学习领域里面比较常用的衡量模型效果的指标。整体来看，这篇文章的数据集大约是 200 条时间序列，时间序列的时间间隔通常来说是五分钟（不过其余的运维场景会有一分钟的数据采集粒度），而一般来说都拥有大半年甚至一年的时间跨度。那么时间点的个数预估是 200 * (1440 / 5) * 365。假设异常的数据：正常的数据 = 1：10000（也就是说平均每条时间序列每周至少发生一次异常），于是这批时间序列数据的异常数据量大约是 200 * (1440 / 5) * 365 / 10000 = 2102，也就是说异常的样本大约是 2102 个左右，剩下的都是正常的样本。PS：当然如果异常的数据：正常的数据的比例大于 1：10000 的话，异常的样本还会更多一些。整体来看，时间序列异常检测是一个样本极其不均衡的场景。

ADS 的系统架构

按照作者之前论文的经验，时间序列异常检测通常都是先做聚类，然后再根据每一个类的特点来做一个异常检测模型，之前的技术架构就是 ROCKA + Opperentice。因为 ROCKA 可以根据时间序列的走势和趋势来进行时间序列的实时分类/聚类，然后 Opperentice 就是做时间序列异常检测的模型。在本文的场景下，作者把 70 条时间序列分成了5 类，因此只需要维护 5 个时间序列的异常检测模型就可以了。当然把时间序列切分成更多的类也是可以的，只是需要维护的时间序列异常检测就变多了，人工成本会加大。

如果看到上面两幅图，有心的读者一定会发现其实 ADS 就是基于 ROCKA 所做的工作。ADS 先对时间序列进行了分类，然后进行了特征提取的工作，再通过半监督学习模型，最后进行异常检测。也就是说，ADS 会走下面四个步骤：

1. ADS 先把历史上的时间序列进行聚类；
2. 通过算法获得每一个类的质心，并且标记出质心曲线的异常点和正常点；
3. 对新来的时间序列进行实时聚类，划分到合适的类别；
4. 基于新来的时间序列（没有标记）和历史上的时间序列（有标记）使用无监督算法来重新训练一个新的模型，进行该类别的时间序列异常检测。

ADS 的细节

对于时间序列的聚类框架 ROCKA，之前的一篇 BLOG 里面已经详细介绍过，这里将不会再赘述。而 ADS 的另一个模块就是半监督学习算法 Contrastive Pessimistic Likelihood Estimation（CPLE），详细的论文细节可以参考论文《Contrastive Pessimistic Likelihood Estimation for Semi-Supervised Classification》。CPLE 有几个好处：

1. CPLE 是半监督学习算法中比较健壮的，因为它并没有过多的假设条件，并且也符合这篇文章的业务场景，同时拥有质心曲线（有标记）和新来的曲线（无标记），使用半监督学习也是符合常理的。除了 CPLE，其实在实战过程中也可以多尝试其他的半监督模型，具体可以参考周志华的《机器学习》。
2. CPLE 的复杂度比较低，计算快。
3. CPLE 支持增量学习，因此，当越来越多新的时间序列进入 ADS，这个模型也会随之而调整并提高准确率。

整体来看，ADS = ROCKA + CPLE，而在论文中，它的对比模型就是 ROCKA + Opperentice。而且在 CPLE 中，也使用了与 Opperentice 系统类似的特征，如下图所示。

其实，从本质上来看，就是半监督学习与有监督学习在这份数据集合上面的比较。从这篇论文里面所展示的数据来看，CPLE 有一定的优势。

整体来看，本篇文章介绍了时间序列异常检测的一种方案，也就是把时间序列先进行聚类的操作，然后根据不同的类来进行异常检测。在异常检测的方法中，不仅可以使用 Random Forest，GBDT，XGBoost 等有监督学习方法，也可以使用 CPLE 等半监督算法。具体在业务中如何使用，其实只能够根据具体的数据来进行合理地选择了。

这次整理的就是清华大学裴丹教授所著的两篇与时间序列相关的论文。一篇是关于时间序列聚类的，《Robust and Rapid Clustering of KPIs for Large-Scale Anomaly Detection》；另外一篇文章是关于时间序列异常检测的，重点检测时间序列上下平移的，《Robust and Rapid Adaption for Concept Drift in Software System Anomaly Detection》。本文将会整理一下这两篇文章的关键技术点。

Robust and Rapid Clustering of KPIs for Large-Scale Anomaly Detection

在互联网公司中，通常会拥有海量的的时间序列，而海量的时间序列就有着各种各样的形状和走势。因此，就有学者提出可以先对时间序列进行分类，然后根据不同的类使用不同的检测模型来进行异常检测。如果要做时间序列的分类，就先需要做聚类的操作，无论从 KMeans，DBSCAN，还是层次聚类来说，都会消耗一定的运算时间。所以，如何在较短的时间内进行聚类或者分类的操作则是这个系统的关键之处。于是，这篇文章提出了一个将时间序列快速聚类的方法。

时间序列 -> 时间序列分类

-> 根据每一类时间序列使用不同的异常检测模型

而在做时间序列聚类的时候，也有着不少的挑战。通常挑战来自于以下几点：

1. 形状：通常来说，时间序列随着业务的变化，节假日效应，变更的发布，将会随着时间的迁移而造成形状的变化。
2. 噪声：无论是从数据采集的角度，还是系统处理的角度，甚至服务器的角度，都有可能给时间序列带来一定的噪声数据，而噪声是需要处理掉的。
3. 平移：在定时任务中，有可能由于系统或者人为的原因，时间序列的走势可能会出现一定程度的左右偏移，有可能每天 5:00 起的定时任务由于前序任务的原因而推迟了。
4. 振幅：通常时间序列都存在一条基线，而不同的时间序列有着不同的振幅，振幅决定了这条时间序列的振荡程度，而振幅或者基线其实也是会随着时间的迁移而变化的。

从整篇论文来看，ROCKA 系统是为了做实时的时间序列分类判断的。要想做成实时的分类判断，就需要有离线和在线两个模块。其中离线是为了做模型训练或者聚类的，在线是为了使用离线处理好的模块来做曲线分类的。

从整个系统来看，离线模块需要做以下几件事情：首先需要收集一批时间序列数据，也就是所谓的 Raw Time Series Data（Raw），通过预处理模块，实施基线提取，再进行聚类的操作，获得相应的聚类结果和质心。在线模块同样也要做类似的事情：首先对于每一条新来的时间序列数据，也就是所谓的 New Time Series Data（Raw），通过预处理模块，实施基线提取，然后使用已经聚类好的离线模块来进行实时的分类。

下面，我们来逐一分析每个模块的作用。

预处理模块（Preprocessing）

通常预处理模块包含几个关键点：

1. 缺失值：缺失值指的是在该上报数据的时间戳上并没有相应的数据上报，数据处于缺失的状态。通常的办法就是把数据补齐，而数据补齐的方法有很多种，最简单的就是使用线性插值（Linear Interpolation）的方式来补齐。
2. 标准化：对于一个时间序列而言，有可能它的均值是10万，有可能只有10，但是它们的走势有可能都是一样的。所以在这个时候需要进行归一化的操作。最常见的有两种归一化方法：一种是标准化，另外一种是最大最小值归一化。如果 表示原始的时间序列的话，标准化指的是 ，其中 和 分别表示均值和标准差。最大最小值归一化指的是 ，其中 分别表示这段时间内的最大值与最小值。

基线提取模块（Baseline Extraction）

基线提取指的是把时间序列分成基线和剩余项两个部分，假设时间序列是 ，基线提取就是：

其中 分别指的是 的基线和剩余项。

在处理基线的时候，有几件事情需要注意：

1. 异常值的处理：通常需要移除一些明显异常的值，然后使用线性插值等方法来把这些移除的值补上。
2. 使用简单的移动平均算法加上一个窗口值 来提取基线。假设时间序列是 ，，。也就是说 。
3. 提取基线的方式其实还有很多，使用带权重的移动平均算法（Weighted Moving Average），指数移动平均算法（Exponentially Weighted Moving Average）都可以提取基线，甚至使用深度学习中的 Autoencoder 或者 VAE 算法都能够提取基线。

基于密度的聚类算法（Clustering）

使用预处理和基线提取技术之后，可以得到原始时间序列的基线值，然后根据这些基线值来进行时间序列的聚类操作。一般来说，时间序列的聚类有两种方法：

1. 通过特征工程的方法，从时间序列中提取出必要的时间序列特征，然后使用 KMeans 等算法来进行聚类。
2. 通过相似度计算工具，对比两条时间序列之间的相似度，相似的聚成一类，不相似的就分成两类。

对于第一种方法而言，最重要的就是特征工程；对于第二种方法而言，最重要的是相似度函数的选择。在这篇文章中，作者选择了第二种方法来进行时间序列的聚类。对于两条时间序列 和 而言，为了解决左右平移的问题，需要考虑一个偏移量 ，然后计算它们之间的内积。

通过这个偏移量 就可以计算出最大的相似度，然后计算出两条时间序列之间的距离。i.e.

其中 指的是 Normalized version of Cross-Correlation， 指的是 Shape-based distance。

而聚类的另一个重要指标就是质心的选择，在这里，每个类的质心可以设置为：

实时分类（Assignment）

对于一条新的时间序列（Raw Data），同样需要经过预处理，基线提取等步骤，然后计算它与之前每一个质心的距离。然后进行距离的从小到大排序，最小的那一类可能就是所需要的。当然，当最小距离大于某个 threshold ，就说明这条新的时间序列曲线很可能不属于之前的任何一类。通过人工查看之后，可以考虑新增一类，并且更新之前所做的聚类模型。

聚类与时间序列异常检测

如果要做海量的时间序列异常检测的话，通常有以下两种做法。

1. 先对时间序列进行聚类或者分类，针对不同的时间序列类型来使用不同的模型；
2. 在时间序列异常检测中加入分类特征。

对于第一种方法而言，需要针对不同形状的时间序列维护不同的模型，而且如果第一层的聚类/分类错误了，那么使用的模型也会出现错误。对于第二种方法而言，关键在于样本的积累和分类特征的构建。

Robust and Rapid Adaption for Concept Drift in Software System Anomaly Detection

这篇文章是关于时间序列异常检测的，而清华大学的 Netman 实验室做时间序列异常检测的相关工作比较多。从 2015 年开始的 Opperentice（Random Forest 模型），Funnel（SST 模型，变更相关），到后来的 Donut（VAE 模型），都是时间序列异常检测的相关文章。而这篇问题提到的 StepWise 系统针对的场景是关于指标的迁移的，所谓概念漂移（Concept Drift）指的就是时间序列会随着变更，发布，调度或者其他事件的发生而导致上下漂移或者左右漂移。一个比较典型的例子就是关于网络流量的漂移，通常来说在某个特点的时间点，时间序列会出现猛跌或者猛涨的情况，但是下跌或者上涨之后的走势和历史数据是及其相似的，只是绝对值上有所变化。

正如上图所示，在时间戳 08-02 附近，时间序列出现了一个下跌的情况。但是根据历史数据所计算出来的上界（Upper Bound）和下界（Lower Bound）却会把未来一段时间的序列都判断为异常。但是前后的曲线走势却是一样的，其实只有下跌的那一小段时间应该被判断为异常，其他时间都是正常的。基于以上的业务场景，这篇文章的作者就提出了概念漂移（Concept Drift）的一些方法。

在整个 StepWise 系统背后，有两个比较关键的地方。其中第一个关键之处就是如何判断异常点。在这个业务场景下，需要检测的是上图发生暴跌的点。而暴跌或者暴涨只是业务运维用来描述时间序列的一个词语。在统计学领域，这种点通常称为变点（Change Point）。所以概念漂移的检测可以转化为一个变点检测的问题，正是论文里面写的。

Insight 1：Concept drift detection can be converted into a spike detection problem in the change score stream.

如果是进行变点检测的话，其实可以参考 2016 年的论文 Funnel：Assessing Software Changes in Web-based Services。里面使用了 Singular Spectrum Transform 来进行检测，并且使用 Difference in Difference 来判断时间序列是否真正出现了异常。关于 SST（Singular Spectrum Transform）的具体细节可以参考 Yasser Mohammad, Toyoaki Nishida 所撰写的 Robust Singular Spectrum Transform。SST 算法是基于 SVD 算法的，有一定的时间复杂度。而基于 SST 又有学者提出了 Robust SST，具体可以详见论文。而 StepWise 的其中两步就是 Detection 和 Classification，其中的 Detection 使用了 Improved SST，Classification 使用了Difference in Difference。而 Funnel 系统的其中两步也是 Improved SST 和 Difference in Difference，这与 Funnel 有异曲同工之妙，Funnel 系统详情请见下图。

而突变前后的时间序列可以分别叫做 Old Concept 和 New Concept，由于前后的走势几乎一致，所以本文的第二个关键点就是相似度的判断。

Insight 2：The relationship between new concept and old concept can be almost perfectly fitted by linear regression.

从数学的角度来讲，Insight 2 的陈述是想判断两条时间序列是否相似。假设 和 ，那么以下陈述是等价的。

1. 存在一个线性关系 使得二维点集 能够被很好的拟合好，也就是说此刻的方差较小。
2. 与 的 Pearson 系数很高；
3. 与 标准化（z-score）之后的曲线几乎一致。
4. 分别存在两个值 使得 与 几乎一致。

所以，在 StepWise 中，作者通过线性回归的算法来判断 old concept 与 new concept 是否存在线性关系，也就是说这两者是否只是平移的关系。其实也可以尝试上面所说的其余方法。

整体来看，StepWise 大致分成两个关键步骤。首先使用类似 Funnel 系统的思想，先进行异常检测（SST），再使用判断一下是不是因为软件的变更引起的（DID），最后使用一层过滤逻辑（线性拟合）来判断是否出现了概念漂移的情况。在实际使用过程中，无论是异常检测还是过滤逻辑，都需要根据具体的业务来做，很难找到一个固定的模型来解决所有的难题。