timelion扩展--Kibana5.4时间序列分析

Posted 2021-04-21 bytebees

分解 DTW ARIMA

引言

在Kibana5.4时序分析一中我们介绍了timelion的.holt()函数，可以用作预测和异常点检测；时序分析二中介绍了eBay的统计学异常检测方法；除了这些分析手段，常用的工具还有时间序列分解、动态时间规整DTW和经典的自回归移动平均模型ARIMA。 

将这些方法提交到timelion的series_function目录下，集成的数据运算函数(chainable functions)可以给kibana提供更丰富数据可视化和分析方法。

时序分解

在第一节介绍Holt-Winters的时候提到，一个时间序列可以分解成baseline、趋势和季节三部分：y(x)=t(x)+s(x)+b(x)，并给出了可加的模型各部分的递推公式。所以时间序列的分解完全可以这样实现，只需要注意一点：就是当时遗留的关于起始值的选取问题。

• 趋势 
  对于二次指数平滑，趋势的初值可以简单的计算两个数据点之差得到；而对于三次指数平滑，数据有了周期性并且给了多个季节的历史数据，所以可以得到更好的趋势初始值。最常用的方法就是计算两个季节之间对应点之差的均值：

• 季节因素
  时间序列减去计算得到的趋势因素，假设一个周期长度为L，得到的数据对季节做平均：

python的statsmodels.api中有时序分解的方法可直接调用，和R的decompose()效果一样，代码和运行结果图如下：

1. import statsmodels.api as sm

2. def tsa_seasonal_decompose(series,freq=None):

3.    result = sm.tsa.seasonal_decompose(series, model="addictive", freq=freq)

4.    #print result.trend, result.seasonal,result.resid,result.observed

5.    return result.plot()

DTW

做数据分析经常需要比较两组数据的相似性，如果两个时间序列的长度不一致，比如语音识别领域中不同人的语速不同，这时传统的采用一一对应方式计算的欧氏距离、曼哈顿距离都会失效，那该怎么计算相似性呢？

• DTW（Dynamic Time Warping）就是这个问题的一种解决方案，是根据两个序列的走势而不是一一对应的点距离来计算相似性。如果两个序列具有相似的形状，但是在时间轴上不对齐，那么在求相似性之前需要先将序列在时间轴上做扭曲(warping)，通过这种时间序列的延伸或者缩短进一步计算。当序列的波形对齐时所做的扭曲就是正确的，这时候两个序列中所有对应点的距离之和应该是最小的。 

  公式:  γ(i,j)=d(qi,cj)+min(γ(i−1,j−1),γ(i−1,j),γ(i,j−1))

• DTW的时间复杂度是O(mn)，m和n分别是两个序列的长度。为了提高效率，对于较长的序列，通常的做法是限制局部搜索做warping，也就是限制一个时间窗w，认为只有w窗口内的序列点才可以做匹配。

def DTWDistance(s1, s2, w):
    DTW={}
    w = max(w, abs(len(s1)-len(s2)))
    for i in range(-1,len(s1)):
        for j in range(-1,len(s2)):
            DTW[(i, j)] = float('inf')
        DTW[(-1, -1)] = 0
    for i in range(len(s1)):
        for j in range(max(0, i-w), min(len(s2), i+w)):
            dist= (s1[i]-s2[j])**2
        DTW[(i, j)] = dist + min(DTW[(i-1, j)],DTW[(i, j-1)], 
 DTW[(i-1, j-1)])
    return np.sqrt(DTW[len(s1)-1, len(s2)-1])

• 有了衡量时间序列相似性的度量之后，就可以对序列集做k-NN聚类。测试过程中，虽然做了局部搜索的warping，每次都要计算两个序列的DTW距离开销仍然很大，所以首先计算LBKeogh距离，因为DTW(s1,s2)>LBKeogh(s1,s2)，如果LBKeogh距离都比目标序列与最近的类中心距离大，就没有继续计算的必要，否则的话再进一步计算DTW距离：

def knn(train, test, w):
    preds = []
    for ind, i in enumerate(test):
        min_dist = float("inf")
        closest_seq = []
        for j in train:
            if LB_Keogh(i[:-1], j[:-1], 5) < min_dist:
                dist = DTWDistance(i[:-1], j[:-1], w)
                if dist < min_dist:
                    min_dist = dist
                    closest_seq = j
                preds.append(closest_seq[-1])
    return classification_report(test[:,-1], preds)

LBKeogh距离可参考: http://alexminnaar.com/time-series-classification-and-clustering-with-python.html

ARIMA

ARIMA模型是理论上最常用的时间序列预测模型，称为“自回归移动平均模型”。但是要彻底理解ARIMA，需要掌握很多概念：随机游走和Dickey-Fuller检验、AR和MA、acf和pacf…这里我们不做理论的介绍，只讲操作方法。

上图给出了建立ARIMA的过程。拿到一个时间序列需要首先画图展示出来，然后根据图形判断给定的时间序列是否平稳，如果不平稳，需要做平稳性处理——常用的方法有差分、去趋势等等。而差分的次数就是ARIMA的第一个参数值。

• 平稳序列 
  称一个序列是（弱）平稳的如果满足条件：i)均值不随着时间变化；ii)方差齐性，也就是方差不随着时间发生变化；iii)对任意时间h，协方差cov(s,t)=cov(s+h,t+h)

在得到了平稳的时间序列之后，需要确定ARIMA的另外两个参数。这时候需要画出自回归函数(acf)和偏自回归函数(pacf)图像。

• acf - 对于MA(q)模型，acf会在lag=q位置截尾，acf图中的值落入宽带区域；以下图为例，q选择6；

• pacf - 对于AR(p)模型，pacf会在lag=p的位置截尾，同样值会落入宽带区域；以下图为例，p选择1。

此外，如果是自己写代码还可以使用grid search的方法来确定模型的三个参数。通常我们使用R来进行ARIMA建模，python的statsmodels也提供了对应的方法。

from statsmodels.tsa.arima_model import ARIMA
data_log = np.log(data)
model = ARIMA(data_log, order=(p, 1, q))
model_fit = model.fit(disp=0)
yhat = model_fit.forecast()[0]

总结

时间序列的分析方法还有很多，这三期以Kibana5.4中timelion的请求函数为例，简单介绍了其中的时序分析方法，和未来可在timelion中集成扩展的模型。 

回归模型假设输出变量是相互独立的，而时间序列具有自相关性，建模的时候需要考虑时间上的依赖关系（比如RNN）。但是除了用ARIMA或者Holt-Winters这种统计学方法做时间序列预测，也有不少人使用机器学习的方法建模，下面这篇文章有比较详细的介绍：https://gejza.nipax.cz/_media/stochasticke_procesy:1302.6613.pdf 
虽然机器学习的算法多数都可拓展应用于时序分中，但是常用和有效的（私以为）还是统计学的方法。既然聊到了时间序列，下一期我们就以RNN开始，介绍深度学习和机器学习的相关内容，敬请关注~拜拜~

p.s 最后附上我整理的关于时间序列异常值检测的方法~

好消息

易宝 CTO 陈斌翻译的新书《架构真经》正在京东和亚马逊热卖！

《架构真经》：《架构即未来》姊妹篇，硅谷大咖的干货呈现，互联网架构的50条军规。唐彬、向江旭、叶亚明、段念、吴华鹏、张瑞海、韩军、程炳皓、张云泉、余晨、李大学、霍泰稳联袂力荐。点击“阅读原文”即可购买！

易宝 CTO 陈斌翻译的新书《架构真经》正在京东和亚马逊热卖！

《架构真经》：《架构即未来》姊妹篇，硅谷大咖的干货呈现，互联网架构的50条军规。唐彬、向江旭、叶亚明、段念、吴华鹏、张瑞海、韩军、程炳皓、张云泉、余晨、李大学、霍泰稳联袂力荐。点击“阅读原文”即可购买！