SLS机器学习介绍（01）：时序统计建模-博客-云栖社区-阿里云

安宇雨 - 随手采集
2018-10-31 19:48:31
随手采集
0000-未整理-等待研究

文章系列链接

背景

时序数据是业务监控中最多方法，双十一大盘、业务监控系统、系统性能平台等都可以看到他的身影。根据统计，60%+业务系统监控都是通过时序表达出来。在数量上，以在阿里云神农系统为例，一台机器有3500多个时序数据，以5秒作为采样周期1000台集群计算，每秒钟会产生70W个点。

时序数据从产生到看到一般分为几个阶段：

采集程序通过采样把原始的状态写入队列
流计算程序通过跳动窗口（分钟、15秒等）进行聚合（例如count，avg，p99）
写入NoSQL等存储引擎进行后期可视化与加工
这套架构非常成熟，在集团中也有非常棒的产品例如IoT、LogAgent、Logtail、TT（Hbase）、LogHub、Blink、JStorm、Kepler、TableStore、Hbase、HiTSDB等。

生成时序数据是为了对业务进行监控和报警，但这块目前非常依赖与人的经验，我们来看几个典型的问题：

问题1：业务/交易走势异常发现

例如下图的监控中，在2分钟时间内流量大小下跌30%以上，在2分钟内后又迅速恢复：

undefined

应对这种情况有几种典型做法：

同环比，例如和上一周、上一天同一时刻比较，如果下降或上升一定大小（例如30%），则产生报警
假设系统是平滑的，与上一时刻比较，如果有比较大差距则认为是巨大变化，产生报警
通过统计建模估算大致走势，如果违背则认为是异常报警

这些方法有比较强的假设，并对数据的精度有较为严格的要求，在生产中容易产生误报。

问题2：异常点发现

以线上业务访问为例，从肉眼很难发现隐藏在其中的局部网络流量变化
在分钟监控下，秒级监控中的毛刺被平均

问题3：延时是否是正常

与Pangu团队交流时发现，Latency告警是一个很依赖经验的问题，无法使用例如业务量等传统时序建模手段。主要有几大原因：

Latency 在QPS等请求上升时，一般都是趋于稳定的
Latency 下降并不能说系统出了问题，有可能当前没有访问
Latency 上升可能和系统的压力、请求，机器负载相关，上升30%可能是请求多了排队引擎，并不是软件或系统Hang。

因此Latency是否正常判断，不能以单一要素维度去建模，要考虑到相关因素。

undefined

问题4：在大量实例中发现问题

案例1：在分布式系统下，我们希望所有的实例行为是否都符合预期的，是否有一些实例没有被调度到造成了浪费
案例2：在安全建模中，会对需要对同一类人（例如同一个组织）进行行为建模，其中是否有局部的人行为异常（例如流量在非工作时间段非常大）

这两个案例都会需要从几千个、几十万个实例中去发现不一样的实例。传统的人肉判断方法会显得力不从心。

undefined

我们的思考与方法

作为长期参与飞天操作系统研发，我们不断地在思考这些问题的本质原因和解决方法。在日志服务（SLS）新版本中，我们把机器学习功能在产品推出，希望能替代有经验开发、运维来解决这些问题。

以下系列文章中我们会介绍以上问题的思路及方法，欢迎大家的讨论。

第一篇：时序统计建模

工业级算法设计思路

1. 通过降维减小数据维度

通常，我们的时间序列数据是多维的，多维时间序列可用于描述受监测事物的多个方面状态与情况，然而随着维度的增长，多时间序列异常异常检测的计算时间会快速增加。通过观察发现当时间序列存在与异常发生原因无关的维度时，进行一些相关性分析，去除无关维度不会对异常检测的准确性产生决定性影响。

2. 使用简单的线性判别

在时间序列的研究当中，很多时候少即是多，时间序列受随机性，趋势性影响较大，往往越多的参数，越复杂的模型意味着很难去长时间，更泛化的拟合真实的数据。而模型的复杂通常也意味着复杂度的提升，因此我们提倡在大数据算法的设计当中，可以适当的为了泛化能力以及时间效率，选择较为简单的线性方法进行粗略的异常判别，后续再使用更精确的算法对子区间进行判别。

3. 尽量设计在线的算法

目前大多数异常检测方法均为静态方法，即对历史中特定段落的时间序列进行分析并得出结果。静态的时间序列方法不能应用于实时的时间序列异常检测。然而在许多应用场景中时间序列是不断增长的，因此我们对于实时获得的时间序列中的异常的需求同样迫切。而且在线的算法也意味着效率近似于线性，是一件非常让人愉悦的事情。

4. 通过并行的思想改善耗时部分的性能

大数据的异常检测对于检测方法提出了存储能力与计算能力的新要求。单个计算结点的存储能力与计算能力无法满足这些要求，我们需要利用并行化计算的方法改善异常检测方法的检测性能。

如何设计 Master 和 Slaver 结点的算法，以及它们划分，合并的关系，将哪个部分正确的，有效的用于并行计算，是这个方向上的研究重点。

5. 做好时间与检测效果的 Trade Off

针对于大数据所设计出的算法，应该有相应的调整系数，能够平衡效率和检测效果的不同侧重。对于某些不需要高精度检测，但是实时性较强的环境，可以通过调整参数来达到相应的需求。

6. 以框架的方式设计算法

以上所提出的几点，均是这个方向上的算法的设计思路，或者说期望达到的目标，一个算法想要很好的结合这三点是很困难的。因此我们认为或许可以使用框架的方式，将各种手段或者技术进行合理的整合，比如异常序列的粗预警，细预警的分离，多维数据，多指标检测的并行分离，等等手段，都可以通过框架的方式进行有机的整合，最终形成一个良好可用的算法。

时序预测模型概况

时序预测模型概况.png

概况说明

统计模型预测：其中自回归和移动平均模型，主要代表算法是ARMA系列；使用二阶或三阶指数平滑算法，主要代表算法是Holt Winters；对于趋势和周期序列，主要利用序列分解策略，主要代表算法是STL；同时也有利用卡尔曼滤波和随机游走算法进行时序的预测；
机器学习模型预测：消除时序中的噪点，利用鲁棒性回归和相应的随机过程（泊松过程、高斯过程等）配合使用Robust Regression；局部学习算法Local Learning；还有针对历史数据进行建模，利用贝叶斯算法进行时序预测；
深度学习模型预测：利用递归神经网络RNN，以及RNN的变形LSTM模型，进行时序预测；同时还有利用VAE（Variance Auto Encoder）算法将时序数据利用编码器的方式去噪声，并利用网络的强大拟合能力对其中数据缺失等问题得以解决；改进WaveNet模型使用CNN的形式去模拟部分RNN的策略较好的利用曲线的局部特征进行预测；

本文重点

指标趋势预测：通过分析指标历史数据，判断未来一段时间指标趋势及预测值，常见有 Holt-Winters、时序数据分解（STL）、ARMA 系列算法。该算法技术可用于异常检测、容量预测、容量规划等场景。
主要阐述算法基本模型和使用方法。

模型实战

时间序列一般具有如下比较明显的特征：

趋势性：数据呈现某种持续向上或向下的趋势或者规律
季节性：数据呈现季节性，数据按照一定的规则进行往复出现
随机性：一些采集过程中的噪声或一些不规则的波动

时间序列预测法用于短期、中期和长期预测。时间序列的本质特性是承认动态数据之间的相关性或依赖关系，这种相关性表征了系统的“动态”或“记忆”。如果这种相关性可用数学模型描述，则可有系统的过去及现在的取值预测其未来取值。

模型构建流程

统计时序建模.png

数据预处理
- 平滑滤波
- 周期估计
- *

模型原理介绍

1. 数据预处理操作

1.1 指数平滑模型原理

一般常用的指数平滑法为一次指数平滑、二次指数平滑、三次指数平滑。在平滑之前，首先需要有初始值，通常，如果数据项小于20，则初始取值第1、2、3期的平均值；如果数据项大于20，则可以把第1期的值作为初始值，记为S0。

一次指数平滑模型，适用于水平型的数据

St1=α∗xt−1+(1−α)∗St−11

其中，St1为t时刻的平滑值，St−11为t-1时刻的平滑值，xt−1为t-1时刻的实际值，α为常数值。
预测未来时刻的公式为：

St+11=α∗xt+(1−α)∗St1

α选择的经验：当数据波动不大时，选择较小的值；当数据波动较大时，选择较大的值。

二次指数平滑模型，适用于斜坡型时序数据，该算法保留了平滑信息和趋势信息，基本思路就是对平滑之后的数据在平滑一次

St1=α∗xt−1+(1−α)∗St−11

St2=α∗St1+(1−α)∗St−12

其中，St1为t时刻第一期的平滑值，St2为t时刻第二期的平滑值，xt是t时刻的实际值，α为常数值。
预测未来T时刻的计算公式为：

xt+T=AT+BT∗T

其中

AT=2∗St1−St2

BT=11−α(St1−St2)

三次指数平滑（Holt-Winters）该方法适用于含有趋势性和季节性的时序数据

St1=α∗xt−1+(1−α)∗St−11

St2=α∗St1+(1−α)∗St−12

St3=α∗St2+(1−α)∗St−13

其中，St1为t时刻第一期的平滑值，St2为t时刻第二期的平滑值，St3为t时刻第三期的平滑值，xt是t时刻的实际值，α为常数值。
预测未来T时刻的计算公式为：

xt+T=AT+BT∗T+CT∗T2

其中

AT=3∗St1−3∗St2+St3

BT=α2∗(1−α)2[(6−5α)∗St1−2(5−4α)St2+(4−3α)St3]

CT=α22(1−α)2[St1−2St2+St3]

下图为平台提供的各种滤波操作的结果示意图

默认使用Holt-Winters指数平滑函数进行平滑

* | select ts_smooth_simple(stamp, val) from log limit 100

使用滑动窗口平滑函数进行平滑

* | select ts_smooth_fir(stamp, val, 'rectangle', 5, 1, 'avg') from log limit 100

1.2 时序数据差分

问题假设：任何一个时间序列，可以分解为如下加法形式

Xt=mt+st+Yt

其中，mt表示缓慢变化的函数，st是已知周期为d的函数，称为周期项，Yt是平稳随机噪声项。在估计和抽取时间序列中的确定性成分mt和st时，使残差量Yt是一个平稳过程，进而求得关于Yt的合适的数学模型，分析它的性质，并连同mt和st可达到对Xt预报和控制的目的。

1.2.1 趋势项估计和提取

不含有季节项成分模型

Xt=mt+Yt

不失一般性，假设E(Yt)=0。

对趋势项mt进行最小二乘法进行估计

mt=α0+α1∗t+α2∗t2

最小化∑t(xt−mt)2以确定α0,α1,α2

1.2.2 趋势项和周期项估计和提取

用差分方法产生平稳的数据，使用差分法可消除趋势，介绍下差分算子和延迟算子

一阶差分算子定义为

▽Xt=Xt−Xt−1=(1−B)Xt

其中，B是延迟算子，定义为

BXt=Xt−1

二阶差分算子定义为

▽2Xt=▽(▽Xt)=(1−2B+B2)Xt=Xt−2Xt−1+Xt−2

消除时序中的周期和趋势，再进行预测
时序模型符合如下形式

Xt=mt+st+Yt

其中，E(Yt)=0,st=st+d,t∈T

首先定义d步差分算子

▽dXt=Xt−Xt−d=(1−Bd)Xt

将d步差分算子，作用于时序模型，得到如下表达式

▽dXt=▽d(mt+st+Yt)

利用序列的周期为d，则有如下表达式

▽dXt=mt−mt−d+Yt−Yt−d

利用d步差分算子，使经d步差分后数据不含有季节项，使得差分后的结果进行求解使得残差量满足平稳性。

1.3 时序数据平稳性

时间序列分析的一些主要方法都是假定数据样本是来自平稳（随机）序列，这是进行分析和确立平稳时间序列线性线性模型的前提。检验数据的平稳性，实质上是检验原来产生数据样本的随机序列的平稳性。从本质上讲，时间序列模型是通过历史规律来预测未来的，即挖掘出的历史数据的一些性质在将来保持不变，这样对模型的预测才有意义，而平稳性就是用来刻画时间序列的统计性质关于时间平移的不变性。

由平稳序列的定义知，必须检验一下两个内容：

均值和方差是否为常数；
其自协方差函数是否只与时间间隔相关，而与间隔端点的位置无关；

数据平稳性检验的参数检验方法：
首先根据平稳序列具有不随时间变化的均值函数和方差函数，其自相关（系数）函数仅是时间间隔的一元函数的性质，构造检验统计量；其次确定检验统计量的分布；进一步由所给数据计算相应统计量的值；最后查表作出结论，接受或者拒绝该数据是否是平稳随机序列。

2. 时序建模原理

2.1 自回归移动平均模型

线性回归模型

Yt=β1xt1+β2xt2+β3xt3+...+βrxtr+ϵt

其中，t=1,2,3,...,N,ϵtN(0,σ22)
回归模型表示了因变量的观测值yt对自变量(xt1,xt2,xt3,...,xtr)的相关性，这种模型在数学上处理了很多方法，如最小二乘参数估计，处理含有异常值的稳健回归方法（如M估计），处理多元非参数非线性回归的投影寻踪回归（projection pursuit regresion，PPR）。

ARMA(p,q)混合自回归移动平均过程

Yt=c+ϕ1Yt−1+ϕ2Yt−2+...+ϕpYt−p+ϵt+θ1ϵt−1+θ2ϵt−2+...+θqϵt−q

上述模型描述了p阶自回归模型，q阶移动平均模式。经过证明，ARMA(p, q)过程的平稳性完全取决于回归参数(ϕ1,ϕ2,ϕ3,...,ϕp)，而与移动平均参数(θ1,θ2,...,θq)无关。

* | select ts_predicate_simple(stamp, val, 10, 1, 'avg') from log limit 500

这张图显示原始序列和使用ARMA模型序列预测和异常点检测结果示意图
- 蓝色曲线表示原始序列曲线
- 红色曲线表示ARMA模型预测的序列曲线
- 绿色区域表示在置信度为0.85的条件下预测序列的区间
- 红色的空心点表示判定为异常的点（具体的算法说明后续文章说明）

2.2 STL序列分解模型

STL(Seasonal-Trend Decomposition Procedure Based on Loess)为时序分解中一种常见算法，将某时刻数据分解为趋势分量、周期分量、残差分量，具体公式如下：

Yt=Tt+St+Rt,t=1,2,3,...,N

STL分解算法的具体流程

outer loop:
    计算鲁棒性回归参数
    inner loop:
        Step 1: 去趋势
        Step 2: 周期子序列平滑
        Step 3: 周期子序列的低通滤波
        Step 4: 去除平滑周期子序列趋势
        Step 5: 去周期
        Step 6: 趋势平滑

平台实验结果

* | select ts_decompose(stamp, val, 1, 'avg') from log limit 2000

下图展示利用本平台提供的算法对时序数据的分解结果
- 第一张图展示的是原始时序数据结果
- 第二张图展示的是序列分解出来的趋势序列结果
- 第三张图展示的是序列分解出来的周期序列的结果
- 第四张图展示的是序列分解出来的残差序列的结果

模型选择

给定任意一组参数，均可得到一个模型，但是如何衡量模型的好坏？如何在众多模型中选择一个最好的模型？

关键指标

平均绝对差值（Mean Absolute Difference/Error，MAD/MAE）

M=1n∑i=1n|Actual−Forecast|

平均绝对百分偏差（Mean Absolute Percentage Error/Deviation，MAPE/MAPD）

M=1n∑t=1n|Actualt−ForecasttActualt|

平均平方误差（Mean Squared Difference/Error，MSD/MSE）

M=1n∑i=1n(Actual−Forecast)2

绝对误差和（Sum of Absolute Difference/Error，SAD/SAE）

M=∑i=1n|Actual−Forecast|

AIC/BIC信息量准则

AIC/BIC是评估统计模型的复杂度，衡量统计模型拟合的良好程度，一般情况下：

AIC可以表示为:

AIC=2ln(L)+2K

其中K是参数的数量，L是似然函数。假设条件：模型的误差服从独立正太分布，让n为观测数量，RSS为残差平方和，那么AIC的具体计算过程为：AIC=2ln(RSS/n)+2K，其中增加自由参数的数目提高了拟合的优良性，AIC鼓励数据拟合的优良性，但尽量避免出现过拟合的情况，所以有效的模型应该是AIC值最小的一个。

BIC可以表示为:

BIC=2ln(L)+Klog(n)

其中K是参数的数量，L是似然函数。

AIC和BIC的公司前半部分是一样的，后半部分是惩罚项，当n>=8时，kln(n)>=2K，所以，BIC相比AIC在大数据量时对模型参数惩罚得更多，导致BIC更倾向于选择参数少的简单模型。
- *

时序预测

1. 基于条件期望的预测

根据t时刻观测到的一组变量Xt预测变量Yt+1的值。令Yt+1|t∗表示根据Xt作出的对Yt+1的预测。为了评价这个预测有用性，我们需要给出一个损失函数来表示当我们的预测偏离一个特定的量时我们的关心程度。通过假定一个二次损失函数，可以得到一个便利结果。二次损失函数的含义是选择合适的Yt+1|t∗，使得E(Yt+1−Yt+1|t∗)2最小。

2. 基于线性投影的预测

假设待预测的Yt+1|t∗限定为Xt的线性函数

Yt+1|t∗=α′Xt

假设要求解出一个α值，使得预测误差(Yt+1−α′Xt)无关，即下式成立：

E[(Yt+1−α′Xt)Xt]=0

则预测α′Xt称为Yt+1关于Xt的线性投影。

3. 线性投影和普通最小二乘回归

一个将yt+1于xt联系起来的线性回归模型如下：

yt+1=β′xt+μt

给定一个关于y和x的T个观测值的样本，样本残差平方和定义为

∑t=1T(yt+1−β′xt)2

预测的置信区间

如何理解置信区间？

参数95%的置信度在区间A的含义是：采样100次，计算95%置信度的置信区间，有95次计算所得的区间包含真实值。
真实值不变，变的是置信区间。
通常情况下，95%被作为常用的置信度，原理就在于3σ控制，此时已经有很高的置信度了。
随着置信度的上升，置信区间的跨度也就越大，对参数估计的精度必定降低。
点估计就一个值，精度高，但置信度则低。