各位老司机晚上好啊，上篇文章主要采用了Frechet Distance进行曲线分类，这篇文章主要采用机器学习的方法来实现曲线分类，基本思路是对训练集先用聚类方法(如Kmeans和Birch等进行聚类，对数据打上标签)，然后在对测试集采用分类方法(决策树，KNN等)进行分类，决定测试集中曲线的类别，实现曲线较为准确的分类。主要内容包括数据处理，聚类分类算法讨论和结果分析。

一. 数据处理

首先是数据采集，数据采集这一部分我在之前的文章”实习生的程序优化:从90min到90s”已经详细描述了数据采集的全过程，这里主要关注的是数据特征的提取和数据的标准化方法。

1.1 数据特征提取

数据特征的提取，可以把把数据特征分成三个部分:统计特征，时域特征，频域特征。当然还有很多中分法，比如在我的下一篇文章，时序特征分析的，可以将数据分为四个部分:趋势特征，周期特征，季节特征和随机特征。

先来说说数据的统计特征，统计特征主要集中在数据的中心趋势和离中趋势，中心趋势例如均值(mean)，中位数(median)，众数(mode)，中列数(midrange)等，离中趋势例如方差(variance)，极差(range)，四分位数(quartiles),四分位数极差(interquartiles range)，变异系数(coefficient of variation, CV)等

在统计特征方面，主要选择了均值，中位数，方差，变异系数。下面来解释下为啥选择这四个数值作为统计特征。

均值(mean): 这里没有直接使用均值的定义，因为均值极易受到极端值的影响。这里使用了截断均值(trimmed mean)，就是去掉极端值后剩余的数据做一个平均，在计算均值前，去掉了极端值正负1%的数据。这样就可以较为客观的反映出数据的中心趋势。

中位数(median): 其实本来想用众数的，不过数据归一化后，众数就不太好找了，也可以用中位数来表达数据的集中趋势，反映出数据的分布情况。

方差(variance): 经常被用来描述数据的离散情况，就是离中趋势和分散度。这里也是用来描述数据的离散程度

变异系数(coefficient of variation, CV): 变异系数定义为标准差与均值的比值，是数据离散程度归一化的一种描述。主要用于均值大于0的情况，优点是不需要参照数据的平均值。这个是为了消除测量尺度差距导致的误差，因为采样曲线有在好几万的数量级，也有在个位数的数量级。当然在使用变异系数时，最好将均值和标准差列出，变异系数的大小，同时受平均数和标准差两个统计量的影响。

其次是数据时域方面的特征。时域方面选择了自相关系数和信息熵作为参考。

自相关系数(auto-correlation): 也称为序列相关，直观的理解，这个系数是反映的是同一事件在两个不同时期的相关程度，就是度量自己过去的行为对现在的影响。这个系数主要在具有周期性的数据分析比较有用，可以有效的区分出具有周期的曲线数据。

信息熵(information-entropy): 这个特征表示曲线的不确定性，换句话说，如果曲线呈现不规则或者随机的变化，则信息熵越大，反之，呈现周期性或者规律性变化，则信息熵越小。这个特征可以有效分辨出不规则的毛刺数据。

最后是频域特征，频域特征我处理的比较粗糙，就是将曲线进行小波变换得到一系列小波系数(低频系数， 高频系数)。

小波变换的实质是：原信号与小波基函数的相似性。小波系数就是小波基函数与原信号相似的系数。（英文文献中是这样解释：The definition of wavelet transform shows that the wavelet analysis is a measure of similarity the basis functions(wavelets)and the original function .The coefficients caculated indicate how close the function is to the danghter wavelet at that particular scale )。

所以，处理的时候我就是简单的将小波系数去平均值，来描述曲线和基函数的相似度，区分不同特征的曲线。这一点确实有问题，关键是我对小波理论没有理解。哪位老司机带带我，怎样更好的使用小波变换描述信号的频域特征?

当然除了这些特征，我还加入了一些自己造出来的特征，例如，为了识别下面的定时任务，我加入了一个”连续为0时间长度”这样一个特征，用来专门区分定时任务。

图1 定时任务数据类型

1.2 数据标准化方法

还有其他的方法，比如用反函数(正弦，余弦，正切等)将数据映射到[-1,1]区间上，也可以实现规范化，效果也不错。

二．聚类算法实现过程

数据特征提取完之后，把数据集分为训练集和测试集，先用测试集做聚类(无监督学习)打标签，并观察输出聚类结果，调整参数直到聚类结果达到较好的效果。

聚类算法首先选用的是KMeans，这是一种选定初始质心，不断更新质心的值直到聚类结果不在发生变化的算法，Kmeans的基本步骤如下:

1. 从D中随机取k个元素，作为k个簇的各自的质心。

2. 分别计算剩下的元素到k个簇中心的相异度，将这些元素分别划归到相异度最低的簇。

3. 根据聚类结果，重新计算k个簇各自的中心，计算方法是取簇中所有元素各自维度的算术平均数。

4. 将D中全部元素按照新的中心重新聚类。

5. 重复第4步，直到聚类结果不再变化。

6. 将结果输出。

在实际应用时采用scikit-learn中的Kmeans函数，可以指定簇的个数(num_clusters), 选定质心的方式(init是随机还是k-means++, 官网上的例子一般是选的k-means++, 给出的解释是可以提高算法效率),还有质心误差控制tol，就是说在迭代过程中，如果质心的变化小于tol, 则算法结束。

主要调整的参数就是这三个，算法结束后会返回计算出的标签值(lables_)，kmeans聚类结果如下图所示，共3000条数据(分钟粒度)作为训练集(主要展示毛刺，定时任务和周期性变化的数据)。

图2 周期性数据(7天数据，一天不明显)

图3 毛刺(1天数据)

图4 定时任务(1天数据)

关于聚类算法的总结，昨天看到了一篇老司机的总结，写的太好了，这里就不一一总结了(文末有链接，已收藏)。我来说说我用过的方法和一些聚类算法的思考。

Birch聚类算法: Birch算法是一种层次聚类算法(hierarchical cluster)，适用场景主要是大规模数据和较多的簇，基本原理是构建一颗CF树(特征树)，通过将节点不断的加入到CF树中，调整CF树的结构，最后完成聚类。Birch算法是一种增量的俄聚类算法，如果簇不是球形的，Birch不能很好的工作，因为Birch方法用了半径的概念控制聚类的边界。

簇半径表示簇中所有点到簇质心的平均距离。CF中存储的是簇中所有数据点的特性的统计和，所以当我们把一个数据点加入某个簇的时候，那么这个数据点的详细特征，例如属性值，就丢失了，由于这个特征，BIRCH聚类可以在很大程度上对数据集进行压缩。所以，半径和质心的选择会影响CF树的规模。

Scikit-learn中Birch算法主要调整两个参数，一个是n_cluster(簇的个数)，另一个是compute_lables(计算标签)。还有一个branching_factor主要用于调整CF树的规模，这个我不太会调，因为我不太清楚这个参数和样本规模之间的准确关系，所以我是使用的默认值。如果有老司机知道的话，快带我上车。

Mean-Shift聚类算法: 基于质心的一种算法。这个算法不用指定簇的个数，而是给出一个estimate_bandwith这个值，由这个值来决定簇的个数。基本思路是通过一个不断移动的高维圆，规定一个Mean-shift向量，向概率密度高的地方不断移动，将概率密度高(相似的)的数据归为一类。

这个算法结果也不错，在使用过程中estimate_bandwith这个值对最后结果产生主要影响。scikit-learn中estimate_bandwith有两个参数quantile(我理解的是计算mean-shift向量所需要的样本百分比), n_samples每次计算所需要的样本数。对于采集到的样本，quantile=0.2，n_samples=250聚类结果较为合理，将数据集聚类为6种，定时任务和毛刺数据可以区分。

这里就介绍了运行结果比较好的算法，关于其他的聚类算法。老司机们可以参考scikit-learn官网(文末有链接)。

三．分类算法实现过程

训练集聚类完成打上标签过后，就可以对测试集进行分类了。分类算法我主要尝试了两种，决策树和KNN。

先来看下决策树，scikit-learn中DecisionTreeClassifier提供了很多参数，详细的解释如下图所示，实际应用的时候我只调整了一个参数criterion就是选择采用信息熵还是基尼系数构建决策树。

决策树明显的优点就是结果易于解释，整个过程也比较好理解。看了一些论文和博客，感觉决策树算法很像人类做决策的过程。计算复杂度还好，泛化能力高，预测出来的结果较为准确。只需要少量的数据就可以做出较为精确的预测。

决策树的缺点就是容易过拟合，设置叶节点所需的最小样本数(min_samples_split)可以避免这个问题。还有就是通过剪枝也可以避免过拟合的问题，关于如何剪枝我没有进行深入的研究，哪位老司机懂得话，萌新求上车，带带我啊。还有就是决策树生成的时候一般采用的是贪心算法，可能会陷入局部最优解，达不到全局最优。

采用决策树方法分类的结果如下图所示:

图5 毛刺数据(1天数据)

图6 定时任务(1天数据)

图7 周期性规律数据(7天数据)

可以看出，决策树分类算法比较有效，可以将毛刺数据和定时任务和周期性任务加以区分。

但是决策树生成的规则有点复杂，我用一天的数据进行训练，决策树一共生成了近2000条分类规则。这些过于复杂的分类规则会不会影响测试集的分类结果呢，这个有待进一步深究，就分类效果看来，复杂的规则没有对分类结果产生重大的影响。

再来看下KNN(K近邻)算法，这种被称作懒惰学习(lazy-learning)的方法。KNN是基于距离的一种分类方法，基本思路是如果一个样本在特征空间中的k个最相似(即特征空间中最邻近)的样本中的大多数属于某一个类别，则该样本也属于这个类别。

KNN的优点就是简单，不需要估计参数和训练，适合对多分类问题进行分类。

KNN的缺点是计算量大，因为对于每一个待分类的数据都要计算距离，开销大。而且当出现样本倾斜时，可能会使分类效果变差。

KNN算法现在出现了一些变种算法，采用权值的方法（和该样本距离小的邻居权值大）来改进，Han等人于2002年尝试利用贪心法，针对文件分类实做可调整权重的k最近邻居法WAkNN (weighted adjusted k nearest neighbor)，以促进分类效果；而Li等人于2004年提出由于不同分类的文件本身有数量上有差异，因此也应该依照训练集合中各种分类的文件数量，选取不同数目的最近邻居，来参与分类。

其他的一些分类算法比如说SVM(这个比较适合二分类)，随机森林(RandomForest,这个我也尝试了下，直白点就是多个决策树做决策，效果不错)。还有一些基于贝叶斯的分类方法，比如scikit-learn中的GaussaianNB(这个我没有尝试，因为没有完全理解算法原理，回学校会关注下贝叶斯分类器的原理)。

四．结尾

这篇文章的思路和做法和之前在KM上的文章很类似，这里只是详细的描述曲线分类过程中自己遇到的问题和解决方法。但是还是有一些缺陷的，比如对于定时任务的分类:

可以看出，定时任务数据中混入了一些其他的数据，仔细观察了下属于毛刺数据。一开始怀疑是聚类或者分类算法没有选好，但是在分析试验多个算法过后，范县结果并没有明显的改善。仔细分析过后，觉得是数据特征提取的问题，之前虽然加入了一个”连续为0的时间长度”这样一个特征，然而还是没有区分出来。或许可以考虑小波变换，从频域的角度考虑下这个问题，因为看起来毛刺数据和定时任务应该在频率上有明显的分别，这个我回去好好学一下小波变换相关的知识。

下一篇文章萌新想从时间序列挖掘方面对监控平台的曲线预测做的工作做一个介绍，时间序列模型我觉得还是挺有意思的。

图8 定时任务混入了一些毛刺数据

参考文献:

1.scikit-learn官方网站: http://scikit-learn.org(搜索框输入clustering可以看到所有聚类算 法的说明和对比)

2.数据挖掘:概念与技术 第二版(数据特征提取和处理方面的参考书)

3.维基百科关于数据频域时域特征的解释(小波变换，自相关，信息熵等)

4.麦子学院的机器学习课程，关于机器学习的方法我就是看这个课程入门的

5.聚类方法的总结: http://km.oa.com/group/870/articles/show/303685

时域特征

自相关系数与偏自相关系数

自相关（autocorrelation），指时序序列于其自身在不同时间点的互相关性.
离散数据的$l$阶自相关系数定义为

$$autocorrelation_{l} = \frac{1}{(n-l)\sigma ^{2}}\sum_{t=1}^{n-1}(X_{t}-\mu)(X_{t+l}-\mu)$$
$n$是时间序列$X_{i}$的长度, $\sigma^2$为方差，$\mu$为均值

为了更清楚查看自相关系数，这里选择滞后阶数为120。

    from statsmodels.graphics.tsaplots import plot_acf
    plot_acf(df_data['value'],lags=120)
    plt.show()

偏自相关（partial_autocorrelation），滞后k阶偏自相关系数是指在给定中间k-1个随机变量$x_{t-1},x_{t-2} ... x_{t-k+1}$的条件下，或者说在剔除了中间k-1个随机变量的干扰后，$x_{t-k}对x_{t}$影响的相关性度量

$$\alpha_{k} = \frac{Cov(x_{t},x_{t-k}|x_{t-1},...,x_{t-k+1})}{\sqrt{Var(x_{t}|x_{t-1},...,x_{t-k+1})Var(x_{t-k}|x_{t-1},...,x_{t-k+1})}}$$

同样，选取滞后120阶偏自相关系数查看。

    from statsmodels.graphics.tsaplots import plot_pacf
    plot_pacf(df_data['value'],lags=120)
    plt.show()

自相关和偏自相关系数，可以取某一阶的值作为特征值。或者计算Box-Pierce统计值或 Ljung-Box Q*统计值。
Box-Pierce统计：

$$Q = n\sum_{k=1}^{h}r_{k}^{2}$$

The Ljung-Box Q* 统计：

$$Q^{*} = n(n+2)\sum_{k=1}^{h}\frac{r_{k}^{2}}{n-k}$$
其中， $r_{k}$表示k阶自相关系数值，n为观测的时间序列长度，h为最大的滞后阶数。

LB检验是基于一系列滞后阶数，判断序列总体的相关性或者说随机性是否存在。

差分

一阶（前向）差分

$$\Delta f(x_{k}) = f(x_{k+1})-f(x_{k})$$

时间序列1阶差分，即为序列每个值减去前一个值得到的序列。通过差分可以得到平稳的序列。

    plt.plot(np.diff(df_data['value']))
    plt.title("%s" % name)
    plt.show()

特征值可取某阶差分的均值、方差、中位数等聚合函数值作为特征值。

季节和趋势特征

时间序列的几种特征（模式）

• Trend（趋势性），代表了数据长期的增长或者下降的特性。
• Seasonal（季节性），跟日期时间相关的周期性。如一年、一周等，季节性都是已知且固定的频率。
• Cyclic（周期性），由数据上升或下降出现，不固定频率的周期，通常指与经济相关的“商业周期‘”，较长时间（至少两年）的周期性。
• Residual（残差 或者 噪声），剩下的部分视为数据的噪声

一般的时间序列我们只考虑趋势性(T)，季节性(S)和残差(E)。因此可以将时间序列按下式分解为上面的几种模式，分解方式有：
加法模型：

$$y_{t} = S_{t}+T_{t}+E_{t}$$
乘法模型：
$$y_{t} = S_{t}\times T_{t}\times E_{t}$$

时间序列分解有多种方法：Classical decomposition，X11 decomposition，STL decomposition。

Classical时间序列分解 :
该分解方法目前依然有较多应用，但是不建议采用。主要原因：

• 该方法是基于移动平均实现，因此趋势和残差成分，前后会有一段数据(分别为滑动窗口的1/2)缺失.
• 该方法假设季节成分是重复的，虽然适用很多时间序列，但对于时间跨度较长的序列，该经典分解方法没办法得知季节成分随时间的变化.
• 对异常数据不具有鲁棒性，移动平均受到极端值影响.

    from statsmodels.tsa.seasonal import seasonal_decompose
    result = seasonal_decompose(df_data['value'].values, model='additive',freq=1440)#加法模型
    result.plot()

X11 时间序列分解：
X11是基于经典时间序列分解的改进，克服了Classical decomposition的缺点，分解的趋势成分在整个观测时间段内都存在，季节性成分随时间会缓慢变化，对异常值更加鲁棒，甚至增加了节假日效应和已知预测因子的影响，更多的信息可以查看文献3、4。

STL时间序列分解:
STL（Seasonal and Trend decomposition using Loess）是通用的鲁棒性较好的时间序列分解模型。局部加权回归散点平滑法LOESS （locally weighted scatterplot smoothing） 是针对二维变量，评估非线性关系的一种方法。
STL分解方法更灵活，趋势的光滑性、季节成分的变化速率，都可以进行控制。该方法只支持加法模型。

    from stldecompose import decompose
    stl = decompose(df_data['value'].values, period=1440)
    stl.plot()

时间序列的特征
趋势特征衡量：

$$trend = 1-\frac{Var(Y_{detrend})}{Var(Y_{t})}$$
其中, $Y_{t}$表示原始数据, $Y_{detrend}$表示去除趋势后的数据

季节特征衡量：

$$seasonality = 1-\frac{Var(Y_{deseasonal})}{Var(Y_{t})}$$
其中, $Y_{t}$表示原始数据, $Y_{deseasonal}$表示去除趋势后的数据

赫斯特指数

赫斯特指数（The Hurst exponent）是对时间序列长期记忆性（long-term memory）的衡量，表示数据不来自于周期的长期统计上的依赖性。
DFA（detrended fluctuation analysis）与赫斯特指数相似，但是DFA适用于非平稳的时间序列（均值和方差随时间变化）
计算公式和详细内容可见以下链接
Hurst_exponent，DFA，Hurst指数与股票市场指数的关系，赫斯特指数（Hurst）指数及在 Excel 中的实现

熵

熵为不确定性的量度，因为越随机的信源的熵越大。
香农熵定义为：

$$H = -\sum_{i=1}^{n}p_{i}logp_{i}$$
其中， $p_{i}$表示 $i$出现的概率

信息熵表征信源的平均不确定性。

分桶熵（Binned Entropy）

$$binned\_entropy = - \sum_{k=0}^{min(max\_bins, len(x))} p_k log(p_k) \cdot \mathbf{1}_{(p_k > 0)}$$

其中，$p_{k}$为桶k中样本数量所占的比例

近似熵(Approximate entropy)
在统计学中，近似熵(ApEn)是一种用来量化随时间序列数据波动的规律性和不可预测性的技术.
常规的信息熵和分桶熵没有考虑数据的时序性。
例如下面两个时间序列:

series 1: (10,20,10,20,10,20,10,20,10,20,10,20)
series 2: (10,10,20,10,20,20,20,10,10,20,10,20)

两个序列中10,20 出现的概率都是1/2，因此信息熵无法区分这两个序列。
近似熵实际上是在衡量当维数变化时该时间序列中产生新模式的概率的大小，产生新模式的概率越大，序列就越复杂，对应的近似熵也就越大。若时间序列包含许多重复的模式，则会有较小的近似熵，若为不太可预测的数据，则近似熵更大。
短的时序序列的近似熵依赖于算法的参数，当时序序列长度大于2000会比较稳定。

样本熵（sample entropy）
样本熵是一种新的时间序列复杂性测度方法，样本熵是相对于近似熵的改进，旨在降低近似熵的误差，与已知是随机部分有更加紧密的一致性，其精度更好，不依赖数据长度，对丢失数据不敏感。
具体算法思想可以查看文献5.

参考

[1] Autocorrelation and Partial Autocorrelation. https://machinelearningmastery.com/gentle-introduction-autocorrelation-partial-autocorrelation/
[2] Box-Pierce. https://analysights.wordpress.com/tag/box-pierce-test/
[3] Forecasting: Principles and Practice. https://otexts.org/fpp2/
[4] X11 decomposition. https://www.springer.com/gp/book/9783319318202
[5] 时间序列复杂度和熵. http://www.dlnu.edu.cn/xintong/docs/20121128144202794986.pdf

频域特征

傅里叶变换

傅里叶变换（Fourier transform）是一种线性积分变换，用于信号在时域（或空域）和频域之间的变换。

FFT变换是针对一组数值进行运算的，这组数的长度N是2的整数次幂时，效果和速度会更好。

在python中可以通过numpy.fft计算离散傅里叶变换(点击查看大图)

可以取傅里叶变换后频域特征的均值、方差等聚合函数值作为特征值

参考：
傅里叶分析之掐死教程. https://zhuanlan.zhihu.com/p/19763358
A visual introduction. https://www.youtube.com/watch?v=spUNpyF58BY
频域信号处理. http://bigsec.net/b52/scipydoc/frequency_process.html

小波系数

傅里叶变换只能得到一个频谱，做小波变换却可以得到一个时频谱，进行时频分析。

• 傅里叶变换：知道一段时间内，信号的各个频率分量有分别多少。
• 小波变换：知道一段时间内，信号的各个频率分量有分别多少，以及他们都是什么时候出现的。

傅里叶变换局限性

对非平稳过程，傅里叶变换有局限性。

最上边的是频率始终不变的平稳信号，而下边两个则是频率随着时间改变的非平稳信号。 它们同样包含和最上边信号相同频率的四个成分。

做FFT后，可以发现这三个时域上有巨大差异的信号，频谱(幅值谱)却非常一致。尤其是下边两个个非平稳信号，我们从频域上无法区分它们，因为它们包含的四个频率的信号的成分确实是一样的，只是出现的先后顺序不同。

可见，博里叶变换处理非平稳信号有天生缺陷。它只能获取一段信号总体上包含哪些频率的成分。但是对各成分出现的时刻并无所知。因此时域相差很大的两个信号，可能频谱图一样。

小波变换原理

傅里叶变换把无限长的三角函数作为基函数：

这个基函数会伸缩、平移。缩的窄，对应高频；伸得宽，对应低频。然后这个基函数不断和信号做相乘。某一个尺度（宽窄）下乘出来的结果，就可以理解成信号缩包含的当前尺度对应频率成分有多少。

这一步其实是计算信号和三角函数的相关性。

这两种尺度能乘出一个大的值（相关度高），所以信息包含较多的这两个频率成分，在频谱上着两个频率会出现两个峰。以上，就是粗浅意义上傅里叶变换的原理

小波做的改变就在于，将无限长的三角函数基换成了有限长的会衰减的小波基。

从公式可以看出，不用于傅里叶变换，变量只有频率w,小波变换有两个变量：尺度a(scale)和平移量t(translation)。尺度a控制小波函数的伸缩，平移量t控制小波函数的平移。尺度就对应于频率（反比），平移量t就对应于时间。

当伸缩、平移到这么一种重合情况时，也会相乘得到一个大的值。这时候和傅里叶变换不同的是，这不仅可以知道信号有这样频率的成分，而且知道它在时域上存在的具体位置。

时域信号

傅里叶变换-频域信号

小波变换-时频信号

离散小波变换(Discrete Wavelet Transform)
将数据序列进行小波分解，每一层分解的结果是上次分解得到的低频信号再分解成低频和高频两个部分。

2阶离散小波变换

import pywt
a2, d2, d1 = pywt.wavedec(data, 'db4', mode='sym', level=2)

同理，4阶离散小波变换。

可取小波系数固定某阶的均值作为特征值。

参考：
傅立叶分析和小波分析之间的关系. https://www.zhihu.com/question/22864189
WaveletTutorial. http://web.iitd.ac.in/~sumeet/WaveletTutorial.pdf
小波变换教程. http://blog.jobbole.com/101976/