Category Archives: 时间序列

时间序列的标签

本篇文章是为了介绍一种基于少量样本标记而获得更多样本的方法,论文的原文是《Label-Less: A Semi-Automatic Labeling Tool for KPI Anomalies》,是清华大学与多家公司(必示科技,中国建设银行等)的合作论文。

在时间序列异常检测中,因为标注的成本比较大,于是需要寻找一种较少而高效地标注时间序列异常点的方法。在该论文中,Alibaba,Tencent,Baidu,eBay,Sogou提供了上千条时间序列(每条时间序列大约是2-6个月的时间跨度),作者们进行了 30 条 KPIs 的标注工作。但是其标注成本依旧是很大的,于是作者们想到了一种异常相似搜索(anomaly similarity search)的算法,目标是对已经标注好的时间序列异常模式进行模版搜索。目的就是达到 label-less,也就是较少的标注而获得更多的标注数据。

在本篇论文中,在异常检测的过程中,作者们使用了时间序列的预测模型(time series prediction models)来获得时间序列的特征,使用了孤立森林(Isolation Forest) 来对时间序列的特征来做无监督的异常检测。并且其效果由于 one class svm 算法和 local outlier factor 算法。在搜索的部分,作者使用了加速版的 DTW 算法(accelerated dynamic time warping approach)来做相似度的搜索和模式的匹配。其中也尝试了各种技巧和方法,包括 constrained DTW,LB Keogh 方法,early stopping 算法等工具。

整个 Label-Less 的架构图如下表示:

Fig2.png

其中的 Operators 指的是业务运维人员,面对着无标记的多条时间序列曲线。系统首先会进行无监督的异常检测算法啊,包括时间序列的预处理(归一化等)操作,然后使用差分(Difference),移动平均算法(moving average),带权重的移动平均算法(weighted moving average),指数移动平均(ewma),holt winters,ARIMA 等算法来做特征的提取。此时,对于不同的时间序列预测工具,我们可以得到不同的预测值,然后把预测值减去实际值并且取绝对值,就得到时间序列的误差序列。i.e. |p_{i} - x_{i}| 就作为数据点 x_{i} 的特征。

在这种情况下,由于用了六个时间序列预测算法,因此原始的时间序列 X (n\times 1) 就可以变成特征矩阵 X' (n\times 1)。对于特征矩阵 X' 可以使用 isolation forest 来做无监督的异常检测并且做阈值的设定;如下图所示:

Fig3.png

而另外的一部分的异常相似搜索(anomaly similarity search)是在第一部分的基础上在做的,Unsupervised Anomaly Detection 会输出疑似异常或者候选异常,并且基于已知的异常模板(Anomaly Template)进行相似度的匹配,此时可以使用 accelerated DTW 算法,选择出最相似的 Top-K 异常,然后运维人员进行标注,得到更多的样本。

由于,对于两条长度分别是 mn 的时间序列,DTW 相似度算法的时间复杂度是 O(mn),因此在搜索的时候需要必要的加速工作。在这种地方,作者们使用了 LB-Kim,LB-Keogh,LB-Keogh-Reverse 算法来做搜索的加速工作。而这些的时间复杂度是 O(m+n)。整体的思路是,如果两条时间序列 qc 的 LB-Kim,LB-Keogh,LB-Keogh-Reverse 的下界大于某个阈值,则不计算它们之间的 DTW 距离。否则就开始计算 DTW。并且在计算 DTW 的时候,如果大于下界,则会提前终止(early stopping),不会继续计算下去。如果都没有大于阈值,则把这个候选曲线和 dist 距离放入列表,最后根据列表中的 dist 来做距离的逆序排列。

整体流程如下:

AnomalySimilaritySearch

其运行速度也比直接使用 DTW 快不少:

Table3.png

Label-Less 的交互页面如下所示:

图(a)表示使用无监督算法获得的疑似异常;

图(b)表示使用异常搜索算法获得的异常结果。

Fig10.png

下图则表示模板,m 表示模板的长度,c 表示相似的异常候选集个数;

Fig11

总结:

整体来看,本文提供了一种通过少量人工标注无监督算法相似度算法来获得更多样本的方法。在候选的时间序列条数足够多的时候,是可以进行时间序列的相似度匹配的。这给未来在运维领域提供海量的时间序列标注数据给予了一定的技术支持。

 

Advertisements

时间序列的联动分析

背景介绍

在互联网公司里面,通常都会监控成千上万的时间序列,用于保障整个系统或者平台的稳定性。在这种情况下,如果能够对多条时间序列之间判断其是否相关,则对于监控而言是非常有效的。基于以上的实际情况,清华大学与 Alibaba 集团在2019年一起合作了论文《CoFlux: Robustly Correlating KPIs by Fluctuations for Service Troubleshooting》,并且发表在 IWQos 2019 上。CoFlux 这个方法可以对多条时间序列来做分析,并且主要用途包括以下几点:

  1. 告警压缩和收敛;
  2. 推荐与已知告警相关的 Top N 的告警;
  3. 在已有的业务范围内(例如数据库的实例)构建异常波动传播链;

kpis.png

CoFlux 的整体介绍

从论文的介绍中来看,CoFlux 的输入和输出分别是:

输入:两条时间序列

输出:这两条时间序列的以下信息

  1. 波动相关性:两条时间序列是否存在波动相关性?
  2. 前后顺序:如果两条时间序列相关,那么它们的前后波动顺序是什么?是同时发生异常还是存在固定的前后顺序?
  3. 方向性:如果两条时间序列是波动相关的,那么它们的波动方向是什么?是一致还是相反?

Remark. CoFlux 的关键点就是并没有对时间序列做异常检测算法,而是直接从时间序列的历史数据(历史半个月或者一个月)出发,判断两条时间序列之间的波动相关性,并且进一步的分析先后顺序与波动方向。

从论文的介绍中来看,CoFlux 的流程图如下图所示:

coflux流程图1

如果两条时间序列 XY 存在波动相关性,则需要输出这两条时间序列的波动先后顺序和是否同向波动。如果两条时间序列 XY 并不存在波动相关性的话,则不需要判断波动先后顺序和是否同向波动。

coflux流程图2

CoFlux 的细节阐述

已知一个长度是 n 的时间序列 S=\{s_{1},\cdots,s_{n}\},对于任意一个 detector,可以得到一条关于 S 的预测值曲线 P=\{p_{1},\cdots,p_{n}\}。于是针对某个 detector 可以得到一个波动特征序列 E=\{\epsilon_{1},\cdots,\epsilon_{n}\},其中 \epsilon_{i} = s_{i} - p_{i}1\leq i\leq n。因此,一个detector 可以对应一个波动序列特征,也是一个时间序列。因此,对于 m 个 detector,可以对应 m 条波动特征序列,并且它们的长度都是 n

在 CoFlux 算法的内部,根据不同的参数使用了总共 86 个 detector,大致列举如下:

  • Difference:根据昨天,七天前的数据来做差分;
  • Holt-Winters:\{\alpha,\beta,\gamma\} \in \{0.2,0.4,0.6,0.8\}
  • 历史上的均值 & 历史上的中位数:1,2,3,4 周;
  • TSD & TSD 中位数:1,2,3,4 周;
  • Wavelet:1,3,5,7 天;
  • 移动平均算法:MA,WMA,EWMA。PS:根据作者们的说法,在这里,MA等方法并不适用。

detectors

根据直觉来看,

  • 对于任何一条时间序列 kpi,总有一个 detector 可以相对准确地提炼到其波动特征;
  • 如果两条时间序列 XY 波动相关,那么 X 的一个波动特征序列与 Y 的一个波动特征序列应该也是相关的;

Remark. 两条时间序列的波动特征可以对齐同一个 detector,也可以不做对齐工作。如果是前者的话,时间复杂度低;后者的话,时间复杂度高。

下图是从时间序列中提取波动特征曲线的案例:

fluxfeatures.png

提炼时间序列的波动曲线特征只是第一步,后续 CoFlux 还有几个关键的步骤:

  • 特征工程的扩大(amplify): 对波动序列特征进行放大,让某些波动序列特征更加明显;
  • Correlation Measurement:用于解决时间序列存在时间前后的漂移,两条时间序列之间存在 lag 的情况,因此需要对其中一条时间序列做平移操作;
  • CoFlux 考虑了历史数据(历史半个月或者一个月)作为参考,并且一个范围内的 kpi 数量不超过 60 条;

下面来一一讲解这些技术方案,对于每一条波动特征曲线(Flux-Features),按照以下几个步骤来进行操作:

Step 1:对波动特征曲线 E=\{\epsilon_{1},\cdots,\epsilon_{n}\} 做 z-score 的归一化,i.e.

\mu = \frac{\sum_{i=1}^{n}\epsilon_{i}}{n},
\delta = \sqrt{\frac{\sum_{i=1}^{n}(\epsilon_{i}-\mu)^{2}}{n}}.

Step 2:对归一化之后的波动特征曲线做特征放大(feature amplification):定义函数 f_{\alpha,\beta}(x) 如下:

f_{\alpha,\beta}(x)= \begin{cases} e^{\alpha\min(x,\beta)} - 1, \text{ when } x\geq 0,\\ -e^{\alpha\min(|x|,\beta)} + 1, \text{ when } x< 0. \end{cases}

E=\{\epsilon_{1},\cdots,\epsilon_{n}\} 放大之后的波动特征曲线(amplified flux feature)就是:\hat{E}=\{f(\epsilon_{1}),\cdots,f(\epsilon_{n})\}.

Step 3:对于两条放大之后的波动特征曲线(amplified flux features)G=\{g_{1},\cdots,g_{\ell}\}H=\{h_{1},\cdots,h_{\ell}\},可以计算它们之间的相关性,先后顺序,是否同向。

G_{s}= \begin{cases} \{0,\cdots,0,g_{1},\cdots, g_{\ell-s}\}, \text{ when } s\geq 0, \\ \{g_{1-s},\cdots,g_{\ell},0,\cdots,0\}, \text{ when } s< 0. \end{cases}

这里的 0 的个数是 |s| 个。其中,-\ell<s<\ell。特别地,当 s=0 时,G_{0}=\{g_{1},\cdots,g_{s}\}=G,那么我们可以定义 G_{s}H 的内积是:R(G_{s},H) = G_{s}\cdot H,

这里的 \cdot 指的是向量之间的内积(inner product)。同时可以定义相关性(Cross Correlation)为:CC(G_{s},H) = \frac{R(G_{s},H)}{\sqrt{R(G_{s},G_{s})\cdot R(H,H)}}.

由于波动有可能是反向的,那么在这里我们不仅要考虑相关性是大于零的情况,也需要考虑小于零的情况。于是,

minCC = \min_{-\ell<s<\ell}CC(G_{s},H),
maxCC = \max_{-\ell<s<\ell}CC(G_{s},H).

则最小值或者最大值的指标分别是

s_{1}=argmin_{-\ell<s<\ell}CC(G_{s},H),
s_{2}=argmax_{-\ell<s<\ell}CC(G_{s},H).


FCC(G,H) = \begin{cases} (minCC, s_{1}), \text{ when } |maxCC|<|minCC|, \\ (maxCC, s_{2}), \text{ when } |maxCC|\geq|minCC|. \end{cases}

从定义中可以看出,FCC(G,H) 是一个元组,里面蕴含着三个信息,分别是相关性,波动方向,前后顺序。FCC(G,H) \in [-1,1],越接近 1 或者 -1 就表示放大之后的波动特征曲线 GH 越相关。正值的 FCC(G,H) 表示 GH 的波动方向相同,是正相关;负值的 FCC(G,H) 表示 GH 的波动方向想法,是负相关。通过对 s<0 或者 s\geq 0 的分析就可以判断先后顺序。因此,CoFlux 方法的是通过对 FCC(G,H) 的分析来得到最终结果的。

在最后的相关性分析里面,其实伪代码正如论文中所示。先考虑是否存在相关性,再考虑基于相关性下的先后顺序和波动方向。

correlationmeasurement

 

CoFlux 的实战效果

从论文中看,CoFlux 的数据集基本上是小于 60 条时间序列曲线。其中包括 CPU,错误率,错误数,内存使用率,成功率等不同的指标。

datasets

datasets2.png

从运行时间上来看,对于一周的时间序列集合(< 60条)而言,CoFlux 基本上能够在 30 分钟内计算完毕,得到最终的运算结果。

executiontime.png

其效果的评价指标基本上就是机器学习中的常见评价指标了,准确率,召回率之类的。

评价指标

从 F1-Score 的评价指标来看,CoFlux 的效果优于其他算法。

experiments.png

告警压缩

如果对时间序列之间进行告警压缩的话,其实可以大量减少运维人员的工作量。在 CoFlux 里面,时间序列曲线被分成了三类,也就是三个颜色最深的模块。因此 21 条时间序列的告警量在实际中有可能只有三条告警。

alarmclustering

告警关联

在实际运维场景中,除了对告警进行压缩之外,也需要对告警进行关联性的分析。例如一条告警发生了,运维人员都希望知道与它相关的其他告警是什么,这样可以方便运维人员定位问题。

alarmcorrelation

构建告警关系链

在一些相对封闭的场景下,例如 mysql 数据库,通过对它里面的时间序列进行分析。不仅可以得到告警之间是否存在相关性,还可以对先后顺序,波动顺序进行分析。

mysql

 

结论

时间序列之间的联动分析是在运维领域场景下的常见技术,不仅可以做告警的压缩,也能够做告警的关联,还能够构建告警的关系链。在未来的工作中,作者们提到将会用深度学习的方法来进行关联和告警的分析,从而进一步加深对时间序列的研究。

时间序列的聚类

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

word2vec1

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

机器学习的聚类算法

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

KMeans 算法的目的是把欧氏空间 \mathbb{R}^{m} 中的 n 个节点,基于它们之间的距离公式,把它们划分成 K 个类别,其中类别 K 的个数是需要在执行算法之前人为设定的。

kmeans1

从数学语言上来说,假设已知的欧式空间点集为 \{x_{1},\cdots,x_{n}\},事先设定的类别个数是 K,当然 K\leq n 是必须要满足的,因为类别的数目不能够多于点集的元素个数。算法的目标是寻找到合适的集合 \{S_{i}\}_{1\leq i\leq K} 使得 argmin_{S_{i}}\sum_{x\in S_{i}}||x-\mu_{i}||^{2} 达到最小,其中 \mu_{i} 表示集合 S_{i} 中的所有点的均值。

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

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

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

hierarchicalclustering1

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

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

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

时间序列的聚类算法

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

时间序列的特征提取

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

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

时间序列的相似度计算

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

dtw1

如果不考虑时间序列的“扭曲”的话,也可以直接使用欧氏距离,无论是 L^{1}, L^{2} 还是 L^{p} 都有它的用武之地。除了距离公式之外,也可以考虑两条时间序列之间的 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

在微积分里面,通常都会研究可微函数的导数,因为导数是反映可微函数单调性的一个重要指标。假设 f(x) 是定义域 (a,b) 上的可导函数,那么某个点 x_{0}\in(a,b) 的导数则定义为:

f'(x_{0}) = \lim_{x\rightarrow x_{0}}\frac{f(x)-f(x_{0})}{x-x_{0}}.

对于区间 (a,b) 上的可导函数 f(x) 而言,假设 x_{0}\in (a,b)。如果 f'(x_{0})>0,那么在 x_{0} 的附近,f(x) 是严格单调递增函数;如果 f'(x_{0})<0,那么在 x_{0} 的附近,f(x) 是严格单调递减函数;如果 f'(x_{0})=0,则基于这个事实无法轻易的判断 f(x)x_{0} 附近的单调性。可以参考这两个例子:(1)f(x)=x^{2}x_{0}=0;(2)f(x) = x^{3}x_{0}=0。这两个例子在 x_{0}=0 的导数都是零,并且第一个例子在 x_{0}=0 附近没有单调性,x_{0}=0 就是最小值点;但是第二个例子在 x_{0}=0 处是严格递增的。

平方函数

立方函数

时间序列的单调性

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

考虑时间序列 X_{N} = [x_{1},\cdots,x_{N}] 的一个子序列 [x_{i},x_{i+1},\cdots,x_{j}],其中 i<j。如果存在某个 k\in (i,j] 和一组非负实数 [w_{i}, w_{i+1},\cdots,w_{j}] 使得

\sum_{m=k}^{j}w_{m}x_{m} > \sum_{m=i}^{k-1} w_{m}x_{m}, 其中\sum_{m=k}^{j}w_{m} = \sum_{m=i}^{k-1}w_{m}.

就称时间序列 [x_{i},x_{i+1},\cdots,x_{j}]上涨的趋势。

如果存在某个 k\in (i,j] 和一组非负实数 [w_{i}, w_{i+1},\cdots,w_{j}] 使得

\sum_{m=k}^{j}w_{m}x_{m} < \sum_{m=i}^{k-1} w_{m}x_{m}, 其中 \sum_{m=k}^{j}w_{m} = \sum_{m=i}^{k-1}w_{m}.

就称时间序列 [x_{i},x_{i+1},\cdots,x_{j}]下跌的趋势。

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

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

假设现在有一个时间序列是 X = [x_{1},\cdots,x_{N}],可以考虑第 i 个点 x_{i} 附近的单调性,按照导数的思想来看就是:当 k\geq 1 时,

(x_{i+k}-x_{i})/((i+k)-i) = (x_{i+k}-x_{i})/k,
(x_{i} - x_{i-k})/(i-(i-k)) = (x_{i} -x_{i-k})/k.

考虑特殊的情形,假设 k=1,当第一个公式大于零时,表示 x_{i+1}>x_{i},i.e. 处于单调上升的趋势中。当第一个公式小于零时,表示 x_{i}<x_{i-1},i.e. 处于单调下降的趋势中。

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

双均线1

简单的移动平均算法

在时间序列领域,简单的移动平均算法 (Simple Moving Average) 是最常见的算法之一。假设原始的时间序列是 X=[x_{1},\cdots,x_{N}],如果考虑时间戳 n 的移动平均值,那就是考虑从时间戳 n 开始,历史上某个窗口上面的所有序列的平均值,用数学公式来描述就是:

M_{w}(n) = \frac{x_{n-w+1}+\cdots+x_{n}}{w} = \frac{\sum_{j=n-w+1}^{n}x_{j}}{w},

其中 w\geq 1 指的就是窗口的大小。

命题 1. 假设窗口值 \ell>s\geq 1M_{s}(n) - M_{\ell}(n) >0, 表示短线上穿长线,曲线有上涨的趋势;M_{s}(n) - M_{\ell}(n) <0, 表示短线下穿长线,曲线有下跌的趋势。

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

证明.
根据条件可以得到,n-\ell+1\leq n-s<n-s+1<n。假设 M_{s}(n) > M_{\ell}(n),那么通过数学推导可以得到:

M_{s}(n) > M_{\ell}(n)
\Leftrightarrow \frac{\sum_{j=n-s+1}^{n}x_{j}}{s} > \frac{\sum_{j=n-\ell+1}^{n}x_{j}}{\ell} = \frac{\sum_{j=n-\ell+1}^{n-s}x_{j} + \sum_{j=n-s+1}^{n}x_{j}}{\ell}
\Leftrightarrow M_{s}(n)=\frac{\sum_{j=n-s+1}^{n}x_{j}}{s} > \frac{\sum_{j=n-\ell+1}^{n-s}x_{j}}{\ell-s} = M_{\ell-s}(n-s),

此时说明 x_{n} 历史上的 s 个点的平均值大于 x_{n-s} 历史上的 \ell - s 个点的平均值,该序列有上涨的趋势。反之,如果 M_{s}(n) < M_{\ell}(n),那么该序列有下跌的趋势。

带权重的移动平均算法

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

WMA_{w}(n) = \frac{x_{n-w+1}+2\cdot x_{n-w+2}+\cdots + w\cdot x_{n}}{1+2+\cdots+w} = \frac{\sum_{j=1}^{w}j \cdot x_{n-w+j}}{w\ \cdot (w+1)/2}.

wma

命题 2. 
假设窗口值 \ell > s,那么WMA_{s}(n) - WMA_{\ell}(n) >0, 表示短线上穿长线,曲线有上涨的趋势;WMA_{s}(n) - WMA_{\ell}(n) <0, 表示短线下穿长线,曲线有下跌的趋势。

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

证明.
根据假设条件可以得到:n-\ell + 1 \leq n-s < n-s < n。假设 WMA_{s}(n) > WMA_{\ell}(n),那么

WMA_{s}(n) > WMA_{\ell}(n)
\Leftrightarrow \frac{\sum_{j=1}^{s} j \cdot x_{n-s+j}}{s\cdot(s+1)/2} > \frac{\sum_{j=1}^{\ell}j\cdot x_{n-\ell +j}}{\ell\cdot(\ell+1)/2} = \frac{\sum_{j=1}^{\ell-s}j\cdot x_{n-\ell+s} + \sum_{j=\ell -s + 1}^{\ell}j\cdot x_{n-\ell + j}}{\ell\cdot(\ell+1)/2}
\Leftrightarrow \frac{\sum_{j=1}^{s} j \cdot x_{n-s+j}}{s\cdot(s+1)/2} > \frac{\sum_{j=1}^{\ell-s}j\cdot x_{n-\ell+s} + \sum_{j=1}^{s}(j+\ell-s)\cdot x_{n- s + j}}{\ell\cdot(\ell+1)/2}
\Leftrightarrow \sum_{j=1}^{s}\bigg(\frac{j}{s\cdot(s+1)/2} - \frac{j+\ell -s}{\ell\cdot(\ell+1)/2} \bigg) \cdot x_{n-s+j} > \frac{\sum_{j=1}^{\ell-s}j\cdot x_{n-\ell+j}}{\ell\cdot(\ell+1)/2}
\Leftrightarrow \sum_{j=j_{0}}^{s}\bigg(\frac{j}{s\cdot(s+1)/2} - \frac{j+\ell -s}{\ell\cdot(\ell+1)/2} \bigg) \cdot x_{n-s+j} > \frac{\sum_{j=1}^{\ell-s}j\cdot x_{n-\ell+j}}{\ell\cdot(\ell+1)/2}
+ \sum_{j=1}^{j_{0}-1} \bigg(\frac{j+\ell -s}{\ell\cdot(\ell+1)/2}- \frac{j}{s\cdot(s+1)/2}\bigg) \cdot x_{n-s+j},

其中 j_{0}=[s\cdot(s+1)/(\ell + s-1)],这里的 [\cdot] 表示 Gauss 取整函数。因为

\frac{j}{s\cdot(s+1)/2} - \frac{j+\ell -s}{\ell\cdot(\ell+1)/2} \geq 0 \Leftrightarrow j \geq \frac{s\cdot(s+1)}{\ell+s-1},

所以不等式两边的系数都是非负数。而 n-\ell + 1 \leq n - s < n-s+1 < n - s + j_{0} -1 < n - s + j_{0} < n,于是距离当前点 x_{n} 的时间序列相比之前的时间序列有上涨的趋势,并且该不等式两边的系数之和是相等的。这是因为

\sum_{j=j_{0}}^{s}\bigg(\frac{j}{s\cdot(s+1)/2} - \frac{j+\ell -s}{\ell\cdot(\ell+1)/2} \bigg) = \frac{\sum_{j=1}^{\ell-s}j}{\ell\cdot(\ell+1)/2} + \sum_{j=1}^{j_{0}-1} \bigg(\frac{j+\ell -s}{\ell\cdot(\ell+1)/2}- \frac{j}{s\cdot(s+1)/2}\bigg)
\Leftrightarrow \sum_{j=1}^{s}\bigg(\frac{j}{s\cdot(s+1)/2} - \frac{j+\ell -s}{\ell\cdot(\ell+1)/2} \bigg) = \frac{\sum_{j=1}^{\ell-s}j}{\ell\cdot(\ell+1)/2},

以上等式易得。于是,当 WMA_{s}(n) >WMA_{\ell}(n) 时,表示时间序列有上涨的趋势;当 WMA_{s}(n) < WMA_{\ell}(n) 时,表示时间序列有下跌的趋势。

指数移动平均算法

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

\text{EWMA}(\alpha, i) = x_{1}, \text{ when } i = 1,
\text{EWMA}(\alpha, i) = \alpha \cdot x_{i} + (1-\alpha) \cdot \text{EWMA}(\alpha, i-1), \text{ when } i \geq 2,

在这里 \alpha\in (0,1)

ewma

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

\text{EWMA}(\alpha, i) = \alpha x_{i} + \alpha(1-\alpha) x_{i-1} + \cdots \alpha(1-\alpha)^{k}x_{i-k} + (1-\alpha)^{k+1}\text{EWMA}(\alpha, t-(k+1)).

在这种情况下,假设 s<\ell,那么短线和长线则分别是:

\text{EWMA}_{s}(\alpha, n) = \alpha x_{n} + \alpha(1-\alpha) x_{n-1} + \cdots + \alpha(1-\alpha)^{s-2}x_{n-s+2} + (1-\alpha)^{s-1}x_{n-s+1},
\text{EWMA}_{\ell}(\beta, n) = \beta x_{n} + \beta(1-\beta) x_{n-1} + \cdots + \beta(1-\beta)^{\ell-2}x_{n-\ell+2} + (1-\beta)^{\ell-1}x_{n-\ell+1}.

在这里,\alpha 是与 s 相关的值,\beta 是与 \ell 相关的值。

命题 3. 
假设 s<\ell,当 0<\beta<\alpha<\min\{1,1/(s-1)\} 时,\text{EWMA}_{s}(\alpha, n) - \text{EWMA}_{\ell}(\beta, n) > 0, 表示短线上穿长线,曲线有上涨的趋势;\text{EWMA}_{s}(\alpha, n) - \text{EWMA}_{\ell}(\beta, n) <0, 表示短线下穿长线,曲线有下跌的趋势。注:当 s=1 时,1/(s-1) 可以看做 +\infty.

证明.
s=1 时,\text{EWMA}_{s}(\alpha,n) = x_{n}。那么

\text{EWMA}_{s}(\alpha, n) > \text{EWMA}_{\ell}(\beta,n)
\Leftrightarrow x_{n} > \beta x_{n} + \beta(1-\beta) x_{n-1} + \cdots + \beta(1-\beta)^{\ell-2}x_{n-\ell+2} + (1-\beta)^{\ell-1}x_{n-\ell+1}
\Leftrightarrow x_{n} > \beta x_{n-1} + \cdots + \beta(1-\beta)^{\ell-3}x_{n-\ell+2}+ (1-\beta)^{\ell-2}x_{n-\ell+1}.

这表示时间序列有上涨的趋势。反之,当 \text{EWMA}_{s}(\alpha, n) = x_{n} < \text{EWMA}_{\ell}(\beta, n) 时,表示时间序列有下跌的趋势。

s\geq 2 时,根据假设有 0<\beta<\alpha<1/(s-1),并且

\text{EWMA}_{s}(\alpha, n) = \alpha x_{n} + \alpha(1-\alpha) x_{n-1} + \cdots + \alpha(1-\alpha)^{s-2}x_{n-s+2} + (1-\alpha)^{s-1}x_{n-s+1},
\text{EWMA}_{\ell}(\beta, n) = \beta x_{n} + \beta(1-\beta) x_{n-1} + \cdots + \beta(1-\beta)^{\ell-2}x_{n-\ell+2} + (1-\beta)^{\ell-1}x_{n-\ell+1}
= \beta x_{n} + \beta(1-\beta) x_{n-1} + \cdots + \beta(1-\beta)^{s-2}x_{n-s+2} + \beta(1-\beta)^{s-1}x_{n-s+1}
+ \beta(1-\beta)^{s}x_{n-s} + \cdots + (1-\beta)^{\ell-1}x_{n-\ell+1}.

假设 g(x) = x(1-x)^{n},通过计算可以得到 g'(x) = (1-x)^{n-1}(1-(n+1)x),也就是说 g(x)(0, 1/(n+1)) 上是递增函数,在 (1/(n+1), 1) 是递减函数。于是当 0<\beta<\alpha<1/(s-1) 时,

\alpha > \beta,
\alpha(1-\alpha) > \beta(1-\beta),
\cdots
\alpha(1-\alpha)^{s-2} > \beta(1-\beta)^{s-2}.

如果 (1-\alpha)^{s-1} > \beta(1-\beta)^{s-1},那么 \text{EWMA}_{s}(\alpha, n) > \text{EWMA}_{\ell}(\beta, n) 可以写成

(\alpha -\beta)x_{n} +\cdots + (\alpha(1-\alpha)^{s-2}-\beta(1-\beta)^{s-2})x_{n-s+2} + ((1-\alpha)^{s-1}-\beta(1-\beta)^{s-1})x_{n-s+1}
> \beta(1-\beta)^{s}x_{n-s} +\cdots + (1-\beta)^{\ell-1}x_{n-\ell+1},

说明在这种情况下时间序列有上涨的趋势。如果 (1-\alpha)^{s-1} < \beta(1-\beta)^{s-1},那么 \text{EWMA}_{s}(\alpha, n)> \text{EWMA}_{\ell}(\beta, n) 可以写成

(\alpha -\beta)x_{n} + \cdots + (\alpha(1-\alpha)^{s-2}-\beta(1-\beta)^{s-2})x_{n-s+2}
> (\beta(1-\beta)^{s-1} - (1-\alpha)^{s-1})x_{n-s+1} + \beta(1-\beta)^{s}x_{n-s} +\cdots + (1-\beta)^{\ell-1}x_{n-\ell+1},

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

反之,当 \text{EWMA}_{s}(\alpha, n) < \text{EWMA}_{\ell}(\beta, n) 时,也可以使用同样的方法证明时间序列有下跌的趋势。

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

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

控制图1

3-\sigma 控制图

假设时间序列是 X_{N} = [x_{1},\cdots, x_{N}],为了计算某个时间戳 nx_{n} 的走势,需要考虑该时间序列历史上的一些点。假设我们考虑 [x_{1},x_{2},\cdots, x_{n}] 中的所有点,可以计算出均值和方差如下:

\mu = \frac{x_{1}+\cdots+x_{n}}{n},
\sigma^{2} = \frac{(x_{1}-\mu)^{2}+\cdots+(x_{n}-\mu)^{2}}{n}.

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

\text{UCL} = \mu + L \cdot \sigma,
\text{Center Line} = \mu,
\text{LCL} = \mu - L \cdot \sigma,

这里的 L 表示系数,通常选择 L=3

命题 4. x_{n} > \text{UCL},那么说明 x_{n} 有上涨的趋势;当 x_{n} < \text{LCL} 时,那么说明 x_{n} 有下跌的趋势;这里的 UCL 和 LCL 是基于 3-\sigma 原理所得到的上下界。

Moving Average 控制图

假设我们考虑的时间序列为 X_{N} = [x_{1},\cdots, x_{N}],那么基于窗口 w 的移动平均值就是

M_{w}(n) = \frac{x_{n-w+1}+\cdots + x_{n}}{w} = \frac{\sum_{j=n-w+1}^{n}x_{j}}{w}.

那么 M_{w}(n) 的方差是

V(M_{w}) = \frac{1}{w^{2}}\sum_{j=n-w+1}^{n} V(x_{j}) = \frac{1}{w^{2}}\sum_{j=n-w+1}^{n}\sigma^{2} = \frac{\sigma^{2}}{w}.

于是,基于移动平均算法的控制图就是:

\text{UCL} = \mu + L\cdot \frac{\sigma}{\sqrt{w}},
\text{Center Line} = \mu,
\text{LCL} = \mu - L \cdot \frac{\sigma}{\sqrt{w}},

这里的 L 表示系数,通常选择 L=3

命题 5. x_{n} > \text{UCL},那么说明 x_{n} 有上涨的趋势;当 x_{n} < \text{LCL} 时,那么说明 x_{n} 有下跌的趋势;这里的 UCL 和 LCL 是基于移动平均算法的控制图所得到的上下界。

macontrolchart

EWMA 控制图

假设 X_{N} = [x_{1},\cdots, x_{N}],那么根据指数移动平均算法可以得到:

z_{i} = x_{1}, \text{ when } i=1,
z_{i} = \lambda x_{i} + (1-\lambda) z_{i-1}, \text{ when } i\geq 2.

进一步分析可以得到:z_{i} 的方差是:

\sigma_{z_{i}}^{2}= \lambda^{2} \sigma^{2} + (1-\lambda)^{2} \sigma_{z_{i-1}}^{2},

于是,
\sigma_{z_{i}}^{2} = \frac{\lambda^{2}}{1-(1-\lambda)^{2}} \sigma^{2} \Rightarrow \sigma_{z_{i}} = \sqrt{\frac{\lambda}{2-\lambda}}\sigma.

因此,基于 EWMA 的控制图指的是:

\text{UCL} = \mu + L\sigma\sqrt{\frac{\lambda}{2-\lambda}},
\text{Center Line} = \mu,
\text{LCL} = \mu - L\sigma\sqrt{\frac{\lambda}{2-\lambda}},

这里的 L 是系数,通常取 L= 3

命题 6. x_{n} > \text{UCL},那么说明 x_{n} 有上涨的趋势;当 x_{n} < \text{LCL} 时,那么说明 x_{n} 有下跌的趋势;这里的 UCL 和 LCL 是基于 EWMA 的控制图所得到的上下界。

ewmacontrolchart

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

macd1

MACD 方法

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

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

用数学公式来详细描述就是:令 \ell = 26, s = 12, signal = 9,基于时间序列 X_{N} = [x_{1},\cdots,x_{N}],可以计算基于指数移动平均的两条线,对于所有的 1\leq n\leq N,有

\text{EWMA}_{s}(\alpha, n) = (1-\alpha) \cdot \text{EWMA}_{s}(\alpha, n-1) + \alpha \cdot x_{n},
\text{EWMA}_{\ell}(\beta,n) = (1-\beta) \cdot \text{EWMA}_{\ell}(\beta, n-1) + \beta \cdot x_{n},

其中

\alpha = \frac{2}{s+1} = \frac{2}{13},
\beta = \frac{2}{\ell+1} = \frac{2}{27}.

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

\text{DIF}(n) = \text{EWMA}_{s}(\alpha, n) - \text{EWMA}_{\ell}(\beta,n).

\gamma = 2 / (signal + 1),计算 DEA 如下:

\text{DEA}(\gamma, n) = \gamma * \text{DIF}(n) + (1-\gamma) * \text{DEA}(\gamma, n).

最后可以计算 MACD 柱状图,对任意的 \forall \text{ }1\leq n\leq N

\text{MACD}(n) = (\text{DIF}(n) - \text{DEA}(\gamma, n)) * 2.

命题 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

时间序列异常检测

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

异常检测技术框架1

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

时间序列的特征工程1

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

\text{Recall}=\frac{\text{the number of true anomalous points detected}}{\text{the number of true anomalous points}}=\frac{TN}{TN+FP},

\text{Precision}=\frac{\text{the number of true anomalous points detected}}{\text{the number of anomalous points detected}}=\frac{TN}{TN+FN},

\text{F1-Score} = \frac{2 \cdot \text{precision} \cdot \text{recall}}{\text{precision}+\text{recall}}.

Table1

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

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

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

深度学习的简单回顾

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

forwardneuralnetworks1

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

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

\sigma(x) = 1/(1+e^{-x}),

\tanh(x) = \sinh(x)/\cosh(x),

ReLU(x) = \max\{0,x\},

Leaky \text{ }ReLu(x) = \mathcal{I}_{\{x<0\}}\cdot(\alpha x) + \mathcal{I}_{\{x\geq 0\}}\cdot(x), \alpha\in \mathbb{R},

ELU(x) = \mathcal{I}_{\{x<0\}}\cdot(\alpha(e^{x}-1)) + \mathcal{I}_{\{x\geq 0\}}\cdot(x),

PreLU(x) = \mathcal{I}_{\{x_{j}<0\}}\cdot(a_{j}x_{j})+\mathcal{I}_{\{x_{j}\geq 0\}}(x_{j}),

selu(x) = \lambda\cdot(\mathcal{I}_{\{x<0\}}\cdot(\alpha e^{x}-\alpha) + \mathcal{I}_{\{x\geq 0\}}\cdot x), \lambda,\alpha\in\mathbb{R},

softplus(x) = \ln(1+e^{x}).

 

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

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

Table2.png

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

Main Theorem. 对于任意正整数 n\geq 1,存在一个前馈神经网络 D 使得对于所有的时间序列 \boldsymbol{X}_{n}=[x_{1},\cdots,x_{n}],该神经网络的输入和输出分别是 \boldsymbol{X}_{n} 和表格 2 中 \boldsymbol{X}_{n} 的特征层。

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

加法 add(x,y) = x+y 与减法 sub(x,y) = x-y 的构造十分简单,如下图构造即可:

绝对值函数 abs(x) = |x|,  通过计算可以得到 abs(x) = relu(x) + relu(-x).  所以,可以构造如下的神经网络来表示绝对值函数:

functionABS

最大值函数 \max(x,y), 通过计算可以得到

\max(x,y) = (|x-y| + x+ y)/2.

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

functionMAX

最小值函数 \min(x,y), 通过计算可以得到

\min(x,y) = (x+y-|x-y|)/2.

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

functionMIN

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

\max(x_{1},\cdots,x_{n}) = \max(x_{1},\max(x_{2},\max(x_{3},\cdots,\max(x_{n-1},x_{n}))),

\min(x_{1},\cdots,x_{n}) = \min(x_{1},\min(x_{2},\max(x_{3},\cdots,\min(x_{n-1},x_{n}))).

平均值函数 avg 指的是 avg(x_{1},\cdots, x_{n}) = (x_{1}+\cdots + x_{n})/n.

functionAVG

平方函数 y = x^{2}, 这个函数可以使用 Softplus 激活函数来表达。令 Softplus 为

f(x) = softplus(x) = \ln(1+e^{x}),

通过计算可以得到:

f(0) = \ln(2),

Df(x) = \sigma(x), Df(0) = 1/2,

D^{2}f(x) = \sigma'(x) = \sigma(x)\cdot(1-\sigma(x)), D^{2}f(0) = 1/4,

D^{3}f(x) = \sigma''(x), D^{3}f(0) = 0,

因此,Softplus 函数的 Taylor Series 是:

f(x) = softplus(x) = f(0) + Df(0)x+ \frac{1}{2!}D^{2}f(0)x^{2} + \frac{1}{3!}D^{3}f(0)x^{3}+o(x^{3})

= \ln(2) +\frac{1}{2}x+\frac{1}{8}x^{2}+o(x^{3}),

因此,x^{2} \approx 8\cdot(f(x) - \ln(2)-\frac{1}{2}x) = 8\cdot(\ln(1+e^{x})-\ln(2)-\frac{1}{2}x). y=x^{2} 就可以用神经网络来近似表示:

functionPower2

立方函数 y = x^{3}, 这个函数可以使用 Sigmoid 激活函数来表达。因为 Sigmoid 函数的 Taylor Series 是

\sigma(x) = \frac{1}{2}+\frac{1}{4}x-\frac{1}{48}x^{3}+o(x^{3}),

那么 x^{3} \approx -48\cdot(\sigma(x) - \frac{1}{2} -\frac{1}{4}x). y=x^{3} 就可以用神经网络来近似表示:

functionPower3

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

提到时间序列的统计特征,一般指的都是已知的时间序列 X_{n} =[x_{1},\cdots,x_{n}] 的最大值,最小值等各种各样的统计指标。如果按照上文所描述的,以下特征都可以用神经网络轻松构造出来:

max:

\max_{1\leq i\leq n}\{x_{1},\cdots,x_{n}\},

min:

\min_{1\leq i\leq n}\{x_{1},\cdots,x_{n}\},

avg:

\mu = \sum_{i=1}^{n}x_{i}/n,

variance:

\sigma^{2}= \sum_{i=1}^{n}(x_{i}-\mu)^{2}/n, \text{ where } \mu = \sum_{i=1}^{n}x_{i}/n,

skewness:

\sum_{i=1}^{n}[(x_{i}-\mu)/\sigma]^{3},

kurtosis:

\sum_{i=1}^{n}[(x_{i}-\mu)/\sigma]^{4},

difference:

x_{2}-x_{1}, x_{3}-x_{2},\cdots, x_{n}-x_{n-1},

integration:

\sum_{i=1}^{n}x_{i},

absolute_sum_of_changes:

E=\sum_{i=1}^{n-1}|x_{i+1}-x_{i}|,

mean_change:

\frac{1}{n}\sum_{i=1}^{n-1}(x_{i+1}-x_{i}) = \frac{1}{n}(x_{n}-x_{1}),

mean_second_derivative_central:

\frac{1}{2n}\sum_{i=1}^{n-2}(x_{i+2}-2x_{i+1}+x_{i}),

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

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

1 \rightarrow 0, 0 \rightarrow 1.

这个逻辑门可以使用逻辑回归函数来估计,可以参见 \sigma 函数的图像,当 x>10 的时候,\sigma(x) \approx 1;x<-10 的时候,\sigma(x)\approx 0. 因此,可以使用函数 f(x) =\sigma(-20x+10) 来估计 NOT 逻辑门。

x=1 时,f(x) = f(1) = \sigma(-10) \approx 0;

x=0 时,f(x) = f(0) = \sigma(10)\approx 1.

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

f_{1}(x) = \sigma(-2\cdot 10^{4} \cdot relu(-x+a) + 10), 可以得到

x>a 时,f_{1}(x) = \sigma(10) \approx 1;

x<a-10^{-3} 时,f_{1}(x) = \sigma(-2\cdot 10^{4}\cdot (a-x) + 10)<\sigma(-10) \approx 0.

因此,所构造的函数 f_{1}(x) 近似于判断待测试值 x 是否大于常数 a.

下面,可以构造一个类似的函数来判断待测试值 x 是否小于常数 a.f_{2}(x) = \sigma(-2\cdot 10^{4} \cdot relu(x-a) + 10), 可以得到

x<a 时,f_{2}(x) = \sigma(10)\approx 1;

x>a+10^{-3} 时,f_{2}(x) = \sigma(-2\cdot 10^{4}\cdot (x-a)+10) < \sigma(-10)\approx 0.

因此,所构造的函数 f_{2}(x) 近似于判断待测试值 x 是否小于常数 a.

回到时间序列的特征 count_above_mean 与 count_below_mean,可以先计算出均值 mean,然后计算时间序列 X_{n}=[x_{1},\cdots,x_{n}] 每个点与均值的差值,然后使用前面的神经网络模块计算出大于零的差值个数与小于零的差值个数即可。

functionCountAboveZero

functionCountBelowZero

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

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

移动平均算法

移动平均算法指的是,已知时间序列 X_{n} = [x_{1},\cdots,x_{n}], 我们可以使用一个窗口值 w\geq 1 得到一组光滑后的时间序列,具体来说就是:

SMA_{j}=\sum_{k=1}^{w}x_{j-w+k}/w = (x_{j-w+1}+\cdots+x_{j})/w,

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

SMA_{n} = \sum_{k=1}^{w}x_{n-w+k}/w = (x_{n-w+1}+\cdots+x_{n})/w.

因此,当前的实际值与光滑后所得到的值的差值就可以作为特征,i.e. SMA_{n}-x_{n} 就可以作为一个特征。然后根据不同的窗口长度 w\geq 1 就可以得到不同的特征值。

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

functionSMA

带权重的移动平均算法

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

WMA_{j} = \sum_{k=1}^{w}k \cdot x_{j-w+k}/\sum_{k=1}^{w}k.

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

WMA_{n} = \sum_{k=1}^{w}k \cdot x_{n-w+k}/\sum_{k=1}^{w}k.

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

functionWMA

指数移动平均算法

指数移动平均算法指的是在已知时间序列的基础上进行加权操作,而权重的大小是呈指数衰减的。用公式来描述就是,已知时间序列 X_{n} = [x_{1},\cdots,x_{n}],

EWMA_{1}=x_{1},

EWMA_{j} = \alpha \cdot x_{j-1} + (1-\alpha)\cdot EWMA_{j-1}, \forall j\geq 1.

从定义上可以得到:

EWMA_{n}

= \alpha[x_{n-1}+(1-\alpha)x_{n-2}+\cdots+(1-\alpha)^{k}x_{n-(k+1)}] + (1-\alpha)^{k+1}EWMA_{n-(k+1)}

\approx \alpha[x_{n-1}+(1-\alpha)x_{n-2}+\cdots+(1-\alpha)^{k}x_{n-(k+1)}]

因此,只需要构建一个加权求和,然后计算 EWMA_{n}-x_{n} 的取值就可以得到特征。所以,神经网络可以构建为如下形式:

functionEWMA

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

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

x_{today}[-1] - x_{yesterday}[-1]x_{today}[-1] - x_{week}[-1] 都是可以通过神经网络构造出来。

mean(x_{today}) - mean(x_{yesterday})mean(x_{today}) - mean(x_{week}) 这一类特征也可以构造出来。

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

\max\{x_{today}[-1]-\max\{x_{yesterday}\}, 0\},

如果当前值 x_{today}[-1] 大于昨天的最大值,就返回它高出的幅度;否则就返回0。

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

\min\{x_{today}[-1]-\min\{x_{week}\},0\},

如果当前值 x_{today}[-1] 小于一周前的最低值,就返回它低于的幅度;否则就返回0。

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

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

在时间序列的分类特征里面,有一种特征叫做值分布特征。假设时间序列的值域在 [0,1] 之内,值分布特征的意思是计算出一个时间序列 X_{n} = [x_{1},\cdots,x_{n}] 的取值在 [0,0.1), [0.1,0.2),\cdots,[0.9,1] 这十个桶的个数,进一步得到它们落入这十个桶的概率是多少。这一类特征可以通过之前所构造的 count 函数来生成。因此,分类特征也是可以通过构造神经网络来形成的。

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

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

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

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

3. 如果神经网络的输入是归一化之后的 raw data,输出是标签 1 或者 0。此时的前馈神经网络需要至少两个以上的隐藏层,才能够达到较好地提取特征的目的。

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

通过前面的陈述,我们可以构造一个端到端(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)。

Figure4

Table3
Figure5

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

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

Table4Table5

Table6

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

Figure8Figure9

论文的主要结论

从本文的主要定理和实验效果来看,前馈神经网络是一个非常有效地可以用作时间序列异常检测的工具。本篇论文不仅提供了一个端到端的训练方法,并且不需要对时间序列进行特征工程的操作。从实验数据来看,使用前馈神经网络(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

Facebook 时间序列预测算法 Prophet 的研究

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 数据的输入和输出

prophetexample1

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

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

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

prophetexample2.png

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

prophetexample3

prophetexample4

Prophet 的算法实现

在时间序列分析领域,有一种常见的分析方法叫做时间序列的分解(Decomposition of Time Series),它把时间序列 y_{t} 分成几个部分,分别是季节项 S_{t},趋势项 T_{t},剩余项 R_{t}。也就是说对所有的 t\geq 0,都有

y_{t} = S_{t} + T_{t} + R_{t}.

除了加法的形式,还有乘法的形式,也就是:

y_{t} = S_{t} \times T_{t} \times R_{t}.

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

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

y(t) = g(t) + s(t) + h(t) + \epsilon_{t}.

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

趋势项模型 g(t)

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

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

如果回顾逻辑回归函数的话,一般都会想起这样的形式:\sigma(x) = 1/(1+e^{-x}), 它的导数是 \sigma'(x) = \sigma(x) \cdot(1-\sigma(x)), 并且 \lim_{x\rightarrow +\infty} \sigma(x) = 1, \lim_{x\rightarrow -\infty} \sigma(x) = 0. 如果增加一些参数的话,那么逻辑回归就可以改写成:

f(x) = C / (1 + e^{-k(x-m)}),

这里的 C 称为曲线的最大渐近值,k 表示曲线的增长率,m 表示曲线的中点。当 C=1, k = 1, m =0 时,恰好就是大家常见的 sigmoid 函数的形式。从 sigmoid 的函数表达式来看,它满足以下的微分方程:y'=y(1-y)

那么,如果使用分离变量法来求解微分方程 y'=y(1-y) 就可以得到:

\frac{y'}{y} + \frac{y'}{1-y} = 1 \Rightarrow \ln\frac{y}{1-y} = 1 \Rightarrow y = 1/(1+K e^{-x}).

但是在现实环境中,函数 f(x) = C / (1+e^{-k(x-m)}) 的三个参数 C, k, m 不可能都是常数,而很有可能是随着时间的迁移而变化的,因此,在 Prophet 里面,作者考虑把这三个参数全部换成了随着时间而变化的函数,也就是 C = C(t), k = k(t), m = m(t)

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

prophetchangepoint1

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

下面假设已经放置了 S 个变点了,并且变点的位置是在时间戳 s_{j}, 1\leq j\leq S 上,那么在这些时间戳上,我们就需要给出增长率的变化,也就是在时间戳 s_{j} 上发生的 change in rate。可以假设有这样一个向量:\boldsymbol{\delta}\in\mathbb{R}^{S}, 其中 \delta_{j} 表示在时间戳 s_{j} 上的增长率的变化量。如果一开始的增长率我们使用 k 来代替的话,那么在时间戳 t 上的增长率就是 k + \sum_{j:t>s_{j}} \delta_{j},通过一个指示函数 \mathbf{a}(t)\in \{0,1\}^{S} 就是

a_{j}(t) = \begin{cases} 1, \text{ if } t\geq s_{j},\\ 0, \text{ otherwise.} \end{cases}

那么在时间戳 t 上面的增长率就是 k + \mathbf{a}^{T}\boldsymbol{\delta}. 一旦变化量 k 确定了,另外一个参数 m 也要随之确定。在这里需要把线段的边界处理好,因此通过数学计算可以得到:

\gamma_{j} = \bigg(s_{j} - m - \sum_{\ell <j} \gamma_{\ell} \bigg) \cdot \bigg( 1- \frac{k + \sum_{\ell < j} \delta_{\ell}}{k + \sum_{\ell\leq j}\delta_{\ell}} \bigg).

所以,分段的逻辑回归增长模型就是:

g(t) = \frac{C(t)}{1+exp(-(k+\boldsymbol{a}(t)^{t}\boldsymbol{\delta}) \cdot (t - (m+\boldsymbol{a}(t)^{T}\boldsymbol{\gamma})},

其中,

\boldsymbol{a}(t) = (a_{1}(t),\cdots,a_{S}(t))^{T},  \boldsymbol{\delta} = (\delta_{1},\cdots,\delta_{S})^{T}, \boldsymbol{\gamma} = (\gamma_{1},\cdots,\gamma_{S})^{T}.

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

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

prophetpiecewiselinear1

因此,基于分段线性函数的模型形如:

g(t)=(k+\boldsymbol{a}(t)\boldsymbol{\delta})\cdot t+(m+\boldsymbol{a}(t)^{T}\boldsymbol{\gamma}),

其中 k 表示增长率(growth rate),\boldsymbol{\delta} 表示增长率的变化量,m 表示 offset parameter。而这两种方法(分段线性函数与逻辑回归函数)最大的区别就是 \boldsymbol{\gamma} 的设置不一样,在分段线性函数中,\boldsymbol{\gamma}=(\gamma_{1},\cdots,\gamma_{S})^{T}, \gamma_{j}=-s_{j}\delta_{j}. 注意:这与之前逻辑回归函数中的设置是不一样的。

在 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 分布,它的概率密度函数为:

f(x|\mu, b) = exp\bigg(-|x-\mu|/b\bigg)/2b,

其中 \mu 表示位置参数,b>0 表示尺度参数。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 表示变点增长率的分布情况,在论文中, \delta_{j} \sim Laplace(0,\tau),这里的 \tau 就是 change_point_scale。

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

对未来的预估(Trend Forecast Uncertainty)

从历史上长度为 T 的数据中,我们可以选择出 S 个变点,它们所对应的增长率的变化量是 \delta_{j} \sim Laplace(0,\tau)。此时我们需要预测未来,因此也需要设置相应的变点的位置,从代码中看,在 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 里面的元素。除了变点的位置之外,也需要考虑 \delta 的情况。这里令 \lambda = \sum_{j=1}^{S}|\delta_{j}|/S,于是新的增长率的变化量就是按照下面的规则来选择的:当 j>T 时,

\delta_{j}=\begin{cases} 0 \text{, with probability } (T-S)/T \\ \sim Laplace(0,\lambda) \text{, with probability } S/T \end{cases}.

季节性趋势

几乎所有的时间序列预测模型都会考虑这个因素,因为时间序列通常会随着天,周,月,年等季节性的变化而呈现季节性的变化,也称为周期性的变化。对于周期函数而言,大家能够马上联想到的就是正弦余弦函数。而在数学分析中,区间内的周期性函数是可以通过正弦和余弦的函数来表示的:假设 f(x) 是以 2\pi 为周期的函数,那么它的傅立叶级数就是 a_{0} + \sum_{n=1}^{\infty}(a_{n}\cos(nx) + b_{n}\sin(nx))

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

s(t) = \sum_{n=1}^{N}\bigg( a_{n}\cos\bigg(\frac{2\pi n t}{P}\bigg) + b_{n}\sin\bigg(\frac{2\pi n t}{P}\bigg)\bigg).

就作者的经验而言,对于以年为周期的序列(P = 365.25)而言,N = 10;对于以周为周期的序列(P = 7 )而言,N = 3。这里的参数可以形成列向量:

\boldsymbol{\beta} = (a_{1},b_{1},\cdots,a_{N},b_{N})^{T}

N = 10 时,

X(t) = \bigg[\cos(\frac{2\pi(1)t}{365.25}),\cdots,\sin(\frac{2\pi(10)t}{365.25})\bigg]

N = 3 时,

X(t) = \bigg[\cos(\frac{2\pi(1)t}{7}),\cdots,\sin(\frac{2\pi(3)t}{7})\bigg]

因此,时间序列的季节项就是:s(t) = X(t) \boldsymbol{\beta},\boldsymbol{\beta} 的初始化是 \boldsymbol{\beta} \sim Normal(0,\sigma^{2})。这里的 \sigma 是通过 seasonality_prior_scale 来控制的,也就是说 \sigma= seasonality_prior_scale。这个值越大,表示季节的效应越明显;这个值越小,表示季节的效应越不明显。同时,在代码里面,seasonality_mode 也对应着两种模式,分别是加法和乘法,默认是加法的形式。在开源代码中,X(t) 函数是通过 fourier_series 来构建的。

节假日效应(holidays and events)

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

prophetholiday1.png

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

h(t)=Z(t) \boldsymbol{\kappa}=\sum_{i=1}^{L} \kappa_{i}\cdot 1_{\{t\in D_{i}\}},

其中 Z(t)=(1_{\{t\in D_{1}\}},\cdots,1_{\{t\in D_{L}\}})\boldsymbol{\kappa}=(\kappa_{1},\cdots,\kappa_{L})^{T}.

其中 \boldsymbol{\kappa}\sim Normal(0,v^{2}) 并且该正态分布是受到 v = holidays_prior_scale 这个指标影响的。默认值是 10,当值越大时,表示节假日对模型的影响越大;当值越小时,表示节假日对模型的效果越小。用户可以根据自己的情况自行调整。

模型拟合(Model Fitting)

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

y(t)=g(t)+s(t)+h(t)+\epsilon

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

Prophet 中可以设置的参数

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

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

prophetdataframefrequency

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

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

prophetexample5prophetexample6

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

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()

prophetexample7

其实 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)

prophetexample8

周期性的设置

通常来说,可以在 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)

prophetexample9

节假日的设置

有的时候,由于双十一或者一些特殊节假日,我们可以设置某些天数是节假日,并且设置它的前后影响范围,也就是 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系统架构
ADS 的系统架构
ROCKA系统架构
ROCKA 的系统架构

如果看到上面两幅图,有心的读者一定会发现其实 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 系统类似的特征,如下图所示。

ADS特征
特征工程

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

ADS效果对比1
Average best F-scores of ADS, iForest, Donut, Opperentice, ROCKA + Opperentice
ADS效果对比2
The New KPI Streams

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