在机器学习领域，聚类问题一直是一个非常常见的问题。无论是在传统的机器学习（Machine Learning）领域，还是自然语言处理（Natural Language Processing）领域，都可以用聚类算法做很多的事情。例如在数据分析领域，我们可以把某个物品用特征来描述出来，例如该房子的面积，价格，朝向等内容，然后使用聚类算法来把相似的房子聚集到一起；在自然语言处理领域，通常都会寻找一些相似的新闻或者把相似的文本信息聚集到一起，在这种情况下，可以用 Word2Vec 把自然语言处理成向量特征，然后使用 KMeans 等机器学习算法来作聚类。除此之外，另外一种做法是使用 Jaccard 相似度来计算两个文本内容之间的相似性，然后使用层次聚类（Hierarchical Clustering）的方法来作聚类。

本文将会从常见的聚类算法出发，然后介绍时间序列聚类的常见算法。

机器学习的聚类算法

KMeans — 基于距离的机器学习聚类算法

KMeans 算法的目的是把欧氏空间 中的 个节点，基于它们之间的距离公式，把它们划分成 个类别，其中类别 的个数是需要在执行算法之前人为设定的。

从数学语言上来说，假设已知的欧式空间点集为 ，事先设定的类别个数是 ，当然 是必须要满足的，因为类别的数目不能够多于点集的元素个数。算法的目标是寻找到合适的集合 使得 达到最小，其中 表示集合 中的所有点的均值。

上面的 表示欧式空间的欧几里得距离，在这种情况下，除了使用 范数之外，还可以使用 范数和其余的 范数。只要该范数满足距离的三个性质即可，也就是非负数，对称，三角不等式。

层次聚类 — 基于相似性的机器学习聚类算法

层次聚类通常来说有两种方法，一种是凝聚，另外一种是分裂。

所谓凝聚，其大体思想就是在一开始的时候，把点集集合中的每个元素都当做一类，然后计算每两个类之前的相似度，也就是元素与元素之间的距离；然后计算集合与集合之前的距离，把相似的集合放在一起，不相似的集合就不需要合并；不停地重复以上操作，直到达到某个限制条件或者不能够继续合并集合为止。

所谓分裂，正好与聚合方法相反。其大体思想就是在刚开始的时候把所有元素都放在一类里面，然后计算两个元素之间的相似性，把不相似元素或者集合进行划分，直到达到某个限制条件或者不能够继续分裂集合为止。

在层次聚类里面，相似度的计算函数就是关键所在。在这种情况下，可以设置两个元素之间的距离公式，例如欧氏空间中两个点的欧式距离。在这种情况下，距离越小表示两者之间越相似，距离越大则表示两者之间越不相似。除此之外，还可以设置两个元素之间的相似度。例如两个集合中的公共元素的个数就可以作为这两个集合之间的相似性。在文本里面，通常可以计算句子和句子的相似度，简单来看就是计算两个句子之间的公共词语的个数。

时间序列的聚类算法

通过以上的描述，如果要做时间序列的聚类，通常来说也有多种方法来做，可以使用基于距离的聚类算法 KMeans，也可以使用基于相似度计算的层次聚类算法。

时间序列的特征提取

之前写过很多时间序列特征提取的方法，无论是常见的时间序列特征，例如最大值，最小值，均值，中位数，方差，值域等内容之外。还可以计算时间序列的熵以及分桶的情况，其分桶的熵指的是把时间序列的值域进行切分，就像 Lebesgue 积分一样，查看落入那些等分桶的时间序列的概率分布情况，就可以进行时间序列的分类。除了 Binned Entropy 之外，还有 Sample Entropy 等各种各样的特征。除了时域特征之外，也可以对时间序列的频域做特征，例如小波分析，傅里叶分析等等。因此，在这种情况下，其实只要做好了时间序列的特征，使用 KMeans 算法就可以得到时间序列的聚类效果，也就是把相似的曲线放在一起。参考文章：时间序列的表示与信息提取。

在提取时间序列的特征之前，通常可以对时间序列进行基线的提取，把时间序列分成基线和误差项。而基线提取的最简单方法就是进行移动平均算法的拟合过程，在这种情况下，可以把原始的时间序列 分成两个部分 和 。i.e. 。有的时候，提取完时间序列的基线之后，其实对时间序列的基线做特征，有的时候分类效果会优于对原始的时间序列做特征。参考文章：两篇关于时间序列的论文。

时间序列的相似度计算

如果要计算时间序列的相似度，通常来说除了欧几里得距离等 距离之外，还可以使用 DTW 等方法。在这种情况下，DTW 是基于动态规划算法来做的，基本想法是根据动态规划原理，来进行时间序列的“扭曲”，从而把时间序列进行必要的错位，计算出最合适的距离。一个简单的例子就是把 和 进行必要的横坐标平移，计算出两条时间序列的最合适距离。但是，从 DTW 的算法描述来看，它的算法复杂度是相对高的，是 量级的，其中 表示时间序列的长度。参考文章：时间序列的搜索。

如果不考虑时间序列的“扭曲”的话，也可以直接使用欧氏距离，无论是 还是 都有它的用武之地。除了距离公式之外，也可以考虑两条时间序列之间的 Pearson 系数，如果两条时间序列相似的话，那么它们之间的 Pearson 系数接近于 1；如果它们之间是负相关的，那么它们之间的 Pearson 系数接近于 -1；如果它们之间没有相关性，Pearson 系数接近于0。除了 Pearson 系数之外，也可以考虑它们之间的线性相关性，毕竟线性相关性与 Pearson 系数是等价的。参考文章：时间序列的相似性。

除此之外，我们也可以用 Auto Encoder 等自编码器技术对时间序列进行特征的编码，也就是说该自编码器的输入层和输出层是恒等的，中间层的神经元个数少于输入层和输出层。在这种情况下，是可以做到对时间序列进行特征的压缩和构造的。除了 Auto Encoder 等无监督方法之外，如果使用其他有监督的神经网络结构的话，例如前馈神经网络，循环神经网络，卷积神经网络等网络结构，可以把归一化之后的时间序列当做输入层，输出层就是时间序列的各种标签，无论是该时间序列的形状种类还是时间序列的异常/正常标签。当该神经网络训练好了之后，中间层的输出都可以作为 Time Series To Vector 的一种模式。i.e. 也就是把时间序列压缩成一个更短一点的向量，然后基于 COSINE 相似度等方法来计算原始时间序列的相似度。参考文章：基于自编码器的时间序列异常检测算法，基于前馈神经网络的时间序列异常检测算法。

总结

如果想对时间序列进行聚类，其方法是非常多的。无论是时间序列的特征构造，还是时间序列的相似度方法，都是需要基于一些人工经验来做的。如果使用深度学习的方法的话，要么就提供大量的标签数据；要么就只能够使用一些无监督的编码器的方法了。本文目前初步介绍了一些时间序列的聚类算法，后续将会基于笔者的学习情况来做进一步的撰写工作。

参考文献

1. 聚类分析：https://en.wikipedia.org/wiki/Cluster_analysis
2. Dynamic Time Warping：https://en.wikipedia.org/wiki/Dynamic_time_warping
3. Pearson Coefficient：https://en.wikipedia.org/wiki/Pearson_correlation_coefficient
4. Auto Encoder：https://en.wikipedia.org/wiki/Autoencoder
5. Word2Vec：https://en.wikipedia.org/wiki/Word2vec，https://samyzaf.com/ML/nlp/nlp.html

在时间序列的众多研究方向上，除了时间序列异常检测，时间序列的相似性，时间序列的趋势预测之外，无论是在量化交易领域还是其余领域，时间序列的单调性都是一个重要课题。本文将会对时间序列的单调性作简单的介绍。

连续函数的单调性

在微积分里面，通常都会研究可微函数的导数，因为导数是反映可微函数单调性的一个重要指标。假设 是定义域 上的可导函数，那么某个点 的导数则定义为：

对于区间 上的可导函数 而言，假设 。如果 ，那么在 的附近， 是严格单调递增函数；如果 ，那么在 的附近， 是严格单调递减函数；如果 ，则基于这个事实无法轻易的判断 在 附近的单调性。可以参考这两个例子：（1），；（2），。这两个例子在 的导数都是零，并且第一个例子在 附近没有单调性， 就是最小值点；但是第二个例子在 处是严格递增的。

时间序列的单调性

通常来说，时间序列分成上涨和下跌两种趋势。如果要严格来写的话，当 时，表示时间序列在 这个区间内是严格单调递增的；当 时，表示时间序列在 这个区间内是严格单调下跌的。但是，在现实环境中，较难找到这种严格递增或者严格递减的情况。在大部分情况下，只存在一个上涨或者下跌的趋势，一旦聚焦到某个时间戳附近时间序列是有可能存在抖动性的。所以我们需要给出一个定义，用来描述时间序列在一个区间内的趋势是上升还是下跌。

考虑时间序列 的一个子序列 ，其中 。如果存在某个 和一组非负实数 使得

其中

就称时间序列 有上涨的趋势。

如果存在某个 和一组非负实数 使得

其中

就称时间序列 有下跌的趋势。

时间序列的单调性 — 均线方法

虽然时间序列是离散的，但是却可以把连续函数的思想应用在上面。

假设现在有一个时间序列是 ，可以考虑第 个点 附近的单调性，按照导数的思想来看就是：当 时，

考虑特殊的情形，假设 ，当第一个公式大于零时，表示 ，i.e. 处于单调上升的趋势中。当第一个公式小于零时，表示 ，i.e. 处于单调下降的趋势中。

但是，时间序列有可能有一定的波动性，也就是说时间序列有可能其实看上去是单调上升的，但是有一定的噪声或者毛刺。所以需要想办法处理掉一些噪声和毛刺。于是，就有人提出了以下几种方法。

简单的移动平均算法

在时间序列领域，简单的移动平均算法 (Simple Moving Average) 是最常见的算法之一。假设原始的时间序列是 ，如果考虑时间戳 的移动平均值，那就是考虑从时间戳 开始，历史上某个窗口上面的所有序列的平均值，用数学公式来描述就是：

其中 指的就是窗口的大小。

命题 1. 假设窗口值 ， 表示短线上穿长线，曲线有上涨的趋势； 表示短线下穿长线，曲线有下跌的趋势。

在这里，短线指的是窗口值 所对应的移动平均线，长线指的是窗口值 所对应的移动平均线。

证明.
根据条件可以得到，。假设 ，那么通过数学推导可以得到：

此时说明 历史上的 个点的平均值大于 历史上的 个点的平均值，该序列有上涨的趋势。反之，如果 ，那么该序列有下跌的趋势。

带权重的移动平均算法

如果窗口值是 ，对于简单移动平均算法，那么 每个元素的权重都是 ，它们都是一样的权重。有的时候我们不希望权重都是恒等的，因为近期的点照理来说是比历史悠久的点更加重要，于是有人提出带权重的移动平均算法 (Weighted Moving Average)。从数学上来看，带权重的移动平均算法指的是

命题 2. 
假设窗口值 ，那么 表示短线上穿长线，曲线有上涨的趋势； 表示短线下穿长线，曲线有下跌的趋势。

在这里，短线指的是窗口值 所对应的带权重的移动平均线，长线指的是窗口值 所对应的带权重的移动平均线。

证明.
根据假设条件可以得到：。假设 ，那么

其中 ，这里的 表示 Gauss 取整函数。因为

所以不等式两边的系数都是非负数。而 ，于是距离当前点 的时间序列相比之前的时间序列有上涨的趋势，并且该不等式两边的系数之和是相等的。这是因为

以上等式易得。于是，当 时，表示时间序列有上涨的趋势；当 时，表示时间序列有下跌的趋势。

指数移动平均算法

指数移动平均算法 (Exponentially Weighted Moving Average) 指的也是移动平均算法，但是它的权重并不是线性递减的，而是呈指数形式递减的。具体来说，如果时间序列是 ，那么它的指数移动平均算法就是：

在这里 。

从数学公式可以推导得出：

在这种情况下，假设 ，那么短线和长线则分别是：

在这里， 是与 相关的值， 是与 相关的值。

命题 3. 
假设 ，当 时， 表示短线上穿长线，曲线有上涨的趋势； 表示短线下穿长线，曲线有下跌的趋势。注：当 时， 可以看做 .

证明.
当 时，。那么

这表示时间序列有上涨的趋势。反之，当 时，表示时间序列有下跌的趋势。

当 时，根据假设有 ，并且

假设 ，通过计算可以得到 ，也就是说 在 上是递增函数，在 是递减函数。于是当 时，

如果 ，那么 可以写成

说明在这种情况下时间序列有上涨的趋势。如果 ，那么 可以写成

说明在这种情况下，时间序列有上涨的趋势。

反之，当 时，也可以使用同样的方法证明时间序列有下跌的趋势。

时间序列的单调性 — 带状方法

根据时间序列的走势，其实可以按照一定的规则计算出它的置信区间，也就是所谓的上界和下界。当最后一些点超过上界或者低于下界的时候，就可以说明这个时间序列的当前的趋势。

控制图

假设时间序列是 ，为了计算某个时间戳 下 的走势，需要考虑该时间序列历史上的一些点。假设我们考虑 中的所有点，可以计算出均值和方差如下：

那么就可以计算出上界，中间线，下界分别是：

这里的 表示系数，通常选择 。

命题 4. 当 ，那么说明 有上涨的趋势；当 时，那么说明 有下跌的趋势；这里的 UCL 和 LCL 是基于 原理所得到的上下界。

Moving Average 控制图

假设我们考虑的时间序列为 ，那么基于窗口 的移动平均值就是

那么 的方差是

于是，基于移动平均算法的控制图就是：

这里的 表示系数，通常选择 。

命题 5. 当 ，那么说明 有上涨的趋势；当 时，那么说明 有下跌的趋势；这里的 UCL 和 LCL 是基于移动平均算法的控制图所得到的上下界。

EWMA 控制图

假设 ，那么根据指数移动平均算法可以得到：

进一步分析可以得到： 的方差是：

于是，

因此，基于 EWMA 的控制图指的是：

这里的 是系数，通常取 。

命题 6. 当 ，那么说明 有上涨的趋势；当 时，那么说明 有下跌的趋势；这里的 UCL 和 LCL 是基于 EWMA 的控制图所得到的上下界。

时间序列的单调性 — 柱状方法

MACD 方法

MACD 算法是比较常见的用于判断时间序列单调性的方法，它的大致思路分成以下几步：

• 根据长短窗口分别计算两条指数移动平均线(EWMA short, EWMA long)；
• 计算两条指数移动平均线之间的距离，作为离差值(DIF)；
• 计算离差值(DIF)的指数移动平均线，作为DEA；
• 将 (DIF-DEA) * 2 作为 MACD 柱状图。

用数学公式来详细描述就是：令 , , ，基于时间序列 ，可以计算基于指数移动平均的两条线，对于所有的 ，有

其中

进一步可以计算离差值 (DIF) 如下：

令 ，计算 DEA 如下：

最后可以计算 MACD 柱状图，对任意的 ，

命题 7. 关于 MACD 的部分性质如下：

• 当 DIF(n) 与 DEA(n) 都大于零时，表示时间序列有上涨的趋势；
• 当 DIF(n) 与 DEA(n) 都小于零时，表示时间序列有递减的趋势；
• 当 DIF(n) 下穿 DEA(n) 时，此时 MACD(n) 小于零，表示时间序列有下跌的趋势；
• 当 DIF(n) 上穿 DEA(n) 时，此时 MACD(n) 大于零，表示时间序列有上涨的趋势；
• MACD(n) 附近的向上或者向下的面积，可以作为时间序列上涨或者下跌幅度的标志。

PS：算法可以从指数移动平均算法换成移动平均算法或者带权重的移动平均算法，长短线的周期可以不局限于 26 和 12，信号线的周期也不局限于 9。

参考资料

1. Moving Average：https://en.wikipedia.org/wiki/Moving_average
2. Double Exponentially Moving Average：https://www.investopedia.com/articles/trading/10/double-exponential-moving-average.asp
3. Control Chart：https://en.wikipedia.org/wiki/Control_chart
4. MACD：https://www.investopedia.com/terms/m/macd.asp
5. Introduction to Statistical Quality Control 6th edition，Douglas C.Montgomery

引言

在时间序列异常检测中，特征工程往往是非常繁琐而复杂的，怎样才能够减少时间序列的特征工程工作量一直是一个关键问题。在本文中，作者们提出了一个新的思路，使用深度学习的办法来进行端到端的训练，从而减少时间序列的特征工程。

提到深度学习，大家都能够想到卷积神经网络（Convolutional Neural Network ）在图像识别中的优异表现，能够想到循环神经网络（Recurrent Neural Network）在机器翻译和文本挖掘领域中所取得的成绩。而一旦提到时间序列，一般的人都能够想到使用 ARIMA 模型或者 LSTM 模型来拟合周期型的时间序列，或者使用其他算法来进行时间序列的异常检测。在这篇文章中，既不谈 CNN 和 LSTM 等深度学习模型，也不谈如何使用 LSTM 来拟合时间序列，本文将会介绍如何使用前馈神经网络 FNN 来进行时间序列的异常检测。并且将会介绍如何使用前馈神经网络，来拟合各种各样的时间序列特征。本篇论文《Feedforward Neural Network for Time Series Anomaly Detection》目前已经挂在 Arxiv 上，有兴趣的读者可以自行参阅：https://arxiv.org/abs/1812.08389。

时间序列异常检测

时间序列异常检测的目的就是在时间序列中寻找不符合常见规律的异常点，无论是在学术界还是工业界这都是一个非常重要的问题。而时间序列异常检测的算法也是层出不穷，无论是统计学中的控制图理论，还是指数移动平均算法，甚至近些年最火的深度学习，都可以应用在时间序列的异常检测上面。在通常情况下，时间序列的异常点是十分稀少的，正常点是非常多的，因此，通常的套路都是使用统计判别算法和无监督算法作为第一层，把有监督算法作为第二层，形成一个无监督与有监督相结合的框架。使用无监督算法可以过滤掉大量的正常样本，将我们标注的注意力放在少数的候选集上；使用有监督算法可以大量的提升准确率，可以把时间序列异常点精确地挑选出来。这个框架之前也说过多次，因此在这里就不再做赘述。

提到第二层的有监督学习算法，通常来说就包括逻辑回归，随机森林，GBDT，XGBoost，LightGBM 等算法。在使用这些算法的时候，不可避免地就需要构造时间序列的特征，也就是人工撰写特征工程的工作。提到时间序列的特征，一般都会想到各种各样的统计特征，例如最大值，最小值，均值等等。除了统计特征之外，我们还可以使用一些简单的时间序列模型，例如移动平均算法，指数移动平均算法等去拟合现有的时间序列，所得到的拟合值与实际值的差值就可以作为时间序列的拟合特征。除了统计特征和拟合特征之外，我们还可以根据时间序列的走势，例如周期型，毛刺型，定时任务型来构造出时间序列的分类特征，用于时间序列形状的多分类问题。因此，就笔者的个人观点，时间序列的特征大体上可以分成统计特征，拟合特征，周期性特征，分类特征等几大类。

在机器学习领域下，可以使用准确率和召回率来评价一个系统或者一个模型的好坏。在这里，我们可以使用 negative 标签来表示时间序列的异常，使用 positive 标签来表示时间序列的正常。因此模型的召回率，准确率，F1-Score 可以如下表示：

而时间序列异常检测工作也不是一件容易的事情，通常来说它具有以下几个难点：

1. 海量时间序列。通常情况下，时间序列不仅仅是按照天来收集数据的，有可能是按照小时，甚至分钟量级来收集数据。因此，在一些情况下，时间序列的数量和长度都是非常大的。
2. 类别不均衡。一般来说，在时间序列异常检测领域，正常样本是非常多的，异常样本是非常少的。在这种情况下，训练模型的时候通常都会遇到类别不均衡的问题。
3. 样本不完整。通常来说，时间序列异常检测领域，是需要用人工来标注样本的，这与推荐系统是非常不一样的。这种情况下，很难通过人工标注的方式，来获得所有类型的样本数据。
4. 特征工程复杂。时间序列有着自己的特点，通过特征工程的方式，确实可以获得不少的特征，但是随着时间序列种类的变多，特征工程将会越来越复杂。

基于以上几个难点，本篇论文提出了一种端到端（End to End）的训练方法，可以解决上面的一些问题。

深度学习的简单回顾

其实最简单的深度学习模型还不是 CNN 和 RNN，最简单的深度学习模型应该是前馈神经网络，也就是所谓的 FNN 模型。当隐藏层的层数较少的时候，当前的前馈神经网络可以称为浅层神经网络；当隐藏层的层数达到一定的数量的时候，当前的前馈神经网络就是所谓的深度前馈神经网络。下面就是一个最简单的前馈神经网络的例子，最左侧是输入层，中间有两个隐藏层，最右侧是输出层。

通常来说，前馈神经网络会涉及到必要的矩阵运算，激活函数的设置等。其中，激活函数的选择有很多，有兴趣的读者可以参见 tensorflow 的官网。比较常见的激活函数有 Sigmoid 函数，tanh 函数，relu 函数以及 relu 函数的各种变种形式（Leaky Relu, PreLu, Elu），以及 Softplus 函数等。

详细来说，以上的激活函数的具体函数表达式如下：

深度学习与时间序列的特征工程

通常来说，基于人工的时间序列特征工程会比较复杂，不仅需要包括均值方差等内容，还包括各种各样的特征，如统计特征，拟合特征，分类特征等。在这种情况下，随着时间的迁移，特征工程将会变得越来越复杂，并且在预测的时候，时间复杂度也会大量增加。那么有没有办法来解决这个问题呢？答案是肯定的。时间序列的一部分特征可以按照如下表格 Table 2 来表示：其中包括均值，方差等特征，也包括拟合特征和部分分类特征。

基于 Table 2，本篇论文的主要定理陈述如下：

Main Theorem. 对于任意正整数 ，存在一个前馈神经网络 使得对于所有的时间序列 ，该神经网络的输入和输出分别是 和表格 2 中 的特征层。

下面，我们就来尝试使用深度学习模型来构造出时间序列的统计特征。首先，我们可以从几个简单的统计特征开始构造，那就是加法（add），减法（minus），最大值（max），最小值（min），均值（avg），绝对值（abs）。在构造时间序列 的以上统计特征之前，我们可以先使用神经网络构造出这几种运算方法。

加法 与减法 的构造十分简单，如下图构造即可：

绝对值函数   通过计算可以得到   所以，可以构造如下的神经网络来表示绝对值函数：

最大值函数 通过计算可以得到

所以，只要能够使用前面的神经网络来构造出绝对值模块，然后使用加减法就可以构造出最大值函数。

最小值函数 通过计算可以得到

所以，同样使用前面的神经网络来构造出绝对值模块，然后使用加减法就可以构造出最小值函数。

在这种情况下，只要能够构造出两个元素的最大值，最小值函数，就可以轻易的构造出 个元素的最大最小值函数，因为

平均值函数 指的是

平方函数 这个函数可以使用 Softplus 激活函数来表达。令 Softplus 为

通过计算可以得到：

因此，Softplus 函数的 Taylor Series 是：

因此， 就可以用神经网络来近似表示：

立方函数 这个函数可以使用 Sigmoid 激活函数来表达。因为 Sigmoid 函数的 Taylor Series 是

那么 就可以用神经网络来近似表示：

深度学习与时间序列的统计特征

提到时间序列的统计特征，一般指的都是已知的时间序列 的最大值，最小值等各种各样的统计指标。如果按照上文所描述的，以下特征都可以用神经网络轻松构造出来：

max:

min:

avg:

variance:

skewness:

kurtosis:

difference:

integration:

absolute_sum_of_changes:

mean_change:

mean_second_derivative_central:

除了以上比较容易构造的特征之外，还有一类特征只为了计算个数的，例如 count_above_mean，count_below_mean 分别是为了计算大于均值的元素个数，小于均值的元素个数。那么最重要的就是要构造出计数函数 count。

回顾一下 NOT 逻辑计算门是：

这个逻辑门可以使用逻辑回归函数来估计，可以参见 函数的图像，当 的时候， 当 的时候， 因此，可以使用函数 来估计 NOT 逻辑门。

当 时，

当 时，

下面，我们来考虑如何构造出一个函数来判断待测试值 是否大于常数

令 可以得到

当 时，

当 时，

因此，所构造的函数 近似于判断待测试值 是否大于常数

下面，可以构造一个类似的函数来判断待测试值 是否小于常数 令 可以得到

当 时，

当 时，

因此，所构造的函数 近似于判断待测试值 是否小于常数

回到时间序列的特征 count_above_mean 与 count_below_mean，可以先计算出均值 mean，然后计算时间序列 每个点与均值的差值，然后使用前面的神经网络模块计算出大于零的差值个数与小于零的差值个数即可。

深度学习与时间序列的拟合特征

时间序列的拟合特征的基本想法是用一些简单的时间序列算法去拟合数据，然后使用拟合数据和真实数据来形成必要的特征。在这里，我们经常使用的算法包括移动平均算法，带权重的移动平均算法，指数移动平均算法等。下面，我们来看一下如何使用神经网络算法来构造出这几个算法。

移动平均算法

移动平均算法指的是，已知时间序列 我们可以使用一个窗口值 得到一组光滑后的时间序列，具体来说就是：

特别地，如果针对时间序列的最后一个点，就可以得到：

因此，当前的实际值与光滑后所得到的值的差值就可以作为特征，i.e. 就可以作为一个特征。然后根据不同的窗口长度 就可以得到不同的特征值。

用和之前类似的方法，我们同样可以构造出一个神经网络算法来得到这个特征。

带权重的移动平均算法

带权重的移动平均算法指的是计算平均值的时候将不同的点带上不同的数值，i.e.

特别地，如果针对时间序列的最后一个点，就可以得到：

用和之前类似的方法，我们同样可以构造出一个神经网络算法来得到这个特征。

指数移动平均算法

指数移动平均算法指的是在已知时间序列的基础上进行加权操作，而权重的大小是呈指数衰减的。用公式来描述就是，已知时间序列 令

从定义上可以得到：

因此，只需要构建一个加权求和，然后计算 的取值就可以得到特征。所以，神经网络可以构建为如下形式：

深度学习与时间序列的周期性特征

在这里，时间序列的周期性特征就是指当前点与昨天同一个时刻，七天前同一个时刻的差值等指标。可以假设时间序列 可以拆分成三个部分 分别是一周前的数据，昨天的数据，今天的数据，假设它们的长度都是 [n/3]，最后一点都表示不同天但是同一个时刻的取值。所以，同环比特征

与 都是可以通过神经网络构造出来。

与 这一类特征也可以构造出来。

有一些特征时用来计算是否高于历史一段时间的最大值，或者低于历史一段时间的最小值，在这里可以先构造 等函数，再计算两者的差值即可。例如，我们可以构造一个特征用于计算当前值是否高过昨天的峰值，以及超出的幅度是多少。用公式来表示那就是：

如果当前值 大于昨天的最大值，就返回它高出的幅度；否则就返回0。

也可以构造一个特征用于计算当前值是否低于一周前的最低值，以及低于的幅度是多少。用公式来表示那就是：

如果当前值 小于一周前的最低值，就返回它低于的幅度；否则就返回0。

这两个特征只需要使用神经网络表示出 激活函数使用 即可。

深度学习与时间序列的分类特征

在时间序列的分类特征里面，有一种特征叫做值分布特征。假设时间序列的值域在 之内，值分布特征的意思是计算出一个时间序列 的取值在 这十个桶的个数，进一步得到它们落入这十个桶的概率是多少。这一类特征可以通过之前所构造的 count 函数来生成。因此，分类特征也是可以通过构造神经网络来形成的。

深度学习与时间序列的特征总结

至此，我们已经证明，对于任意长度 ，存在一个神经网络，它的输入和输出分别是原始的时间序列与 Table 2 中的时间序列特征层。整体来看，

1. 存在多个前馈神经网络可以生成时间序列的特征；

2. 深度学习+时间序列异常检测可以实现端到端（End to End）的训练过程，也就是说：输入数据是归一化之后的原始数据（normalized raw data），输出的是两个标签（正常&异常），神经网络的权重可以通过大量数据集和目标函数训练出来。

基于前馈神经网络的时间序列异常检测算法

通过前面的陈述，我们可以构造一个端到端（End to End）的前馈神经网络，意思就是说：前馈神经网络的输入层是原始的时间序列（归一化之后的数据），前馈神经网络的输出层是标签。

在这里，我们考虑的是三天数据的子序列，以 20180810 的 10:00am 为例，考虑当天历史三小时的数据（07:00-10:00），昨天 20180809 前后三小时的数据（07:00-13:00），再考虑一周前 20180803 前后三小时的数据（07:00-13:00）。这样就形成了一个子序列，总共有 903 个点。然后我们可以使用最大最小归一化获得神经网络的输入数据，而输出数据指的就是最后一个点是异常点（label = 0）还是正常（label = 1）。

Figure 5 指出了前馈神经网络的结构图，输入层是归一化之后的时间序列原始数据，中间两层是隐藏层，输出层就是异常或者正常的概率值。而中间层的激活函数可以使用 ReLU 或者 Leaky ReLU，在这里我们通过实验发现 Leaky ReLU 的效果略好于 ReLU。而最后一层的激活函数使用的是 Softmax 函数，输出的两个概率值之和永远都是 1。

在这种神经网络结构下，神经网络的参数量级大约是 10 万量级，在这种情况下，使用少量的几百几千个样本几乎是无法训练出来的。在这里，我们使用了大约 10 万 的样本数据，才得到一个还不错的效果。在这里，我们使用 3-Sigma 算法，EWMA 控制图算法，多项式回归算法，孤立森林算法，XGBoost + 特征工程，前馈神经网络来进行算法的对比。通过数据的对比可以得到，XGBoost 与 DNN 其实差不多，都能够达到实际使用的上线标准。

从深度学习的基础知识可以得到，CNN 的中间层可以用来提取图片的特征，因此，这里的前馈神经网络的隐藏层的输出同样可以作为时间序列的特征层。于是，我们通过实验，基于隐藏层的输出可以作为时间序列的隐藏特征，也就是所谓的 Time Series To Vector。通过 Time Series To Vector，我们可以既可以对时间序列进行聚类（KMeans），也可以对时间序列进行 Cosine 相似度的计算，进而得到同一类时间序列和相似的时间序列。

论文的主要结论

从本文的主要定理和实验效果来看，前馈神经网络是一个非常有效地可以用作时间序列异常检测的工具。本篇论文不仅提供了一个端到端的训练方法，并且不需要对时间序列进行特征工程的操作。从实验数据来看，使用前馈神经网络（feedforward neural network）可以得到与 XGBoost 差不多的效果。并且，前馈神经网络隐藏层的输出可以作为时间序列的隐藏特征（Time Series To Vector），使用 Cosine 相似度或者 KMeans 算法就可以对时间序列进行相似度的计算和聚类操作。在时间序列异常检测领域，使用特征工程 + 有监督算法的方法论比较多，而使用端到端的训练方法，也就是前馈神经网络的方法应该还是相对较少的。因此，端到端的前馈神经网络算法应该是本文的方法与其他方法论的最大不同点。

参考文献

1. 《企业级 AIOps 实施建议》白皮书-V0.6 版本
2. 《腾讯运维的AI实践》— 2018年4月 GOPS 全球运维大会
3. 《Feedforward Neural Network for Time Series Anomaly Detection》，Arxiv，2018年12月18日
4. Github：https://github.com/Tencent/Metis

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），最后使用一层过滤逻辑（线性拟合）来判断是否出现了概念漂移的情况。在实际使用过程中，无论是异常检测还是过滤逻辑，都需要根据具体的业务来做，很难找到一个固定的模型来解决所有的难题。