# 时间序列异常检测—节假日效应的应对之道

1. 元旦：1 天；
2. 春节：7 天；
3. 清明节：3 天；
4. 五一劳动节：5 天；
5. 端午节：3 天；
6. 国庆节：8 天。

1. 周期性的多样性：通过实际案例可以看出，对于不同的时间序列，其周期是完全不一样的，而且在不同的周期上也有着完全不同的表现；
2. KPI 数量巨大：这个通常来说都是智能运维领域中的常见问题；
3. 周期的漂移：一般来说，通过时间序列的走势我们只能够看出一个大致的变化，但是具体到细节的话，周期是存在一定的波动的。例如不一定恰好是 7 天，有可能是 7 天加减 5 分钟之类的周期。这个跟业务的具体场景有关系，也跟当时的实际情况有关。

1. 离线周期性检测（offline periodicity detection）；

$X_{(s)}=\begin{cases}(0,\cdots,0,x_{1},\cdots,x_{m-s}), &\text{ if } s\geq 0 \\ (x_{1-s},x_{1-s+1},\cdots,x_{m},0,\cdots,0), &\text{ else } s<0.\end{cases}$

$CC_{s}(X,Y)=\begin{cases}\sum_{i=1}^{m-s}x_{i}\cdot y_{s+i}, &\text{ if } s\geq 0 \\ \sum_{i=1}^{m+s}x_{i-s}\cdot y_{i}, &\text{ else } s<0.\end{cases}$

$NCC(X,Y)=\max_{s\in[-w,w]\cap\mathbb{Z}}\frac{CC_{s}(X,Y)}{\|x\|_{2}\cdot\|y\|}.$

$SBD(X,Y) = 1-NCC(X,Y).$

# FluxRank: 如何快速地进行机器故障定位

### 贝叶斯网络

1. 如果草已经湿润，求降雨的概率是多少？
2. 如果草已经湿润，求没有降雨且洒水器开启的概率是多少？

$P(R=T,W=T)=P(R=T,S=T,W=T)+P(R=T,S=F,W=T)$

$= P(W=T|R=T,S=T)P(S=T|R=T)P(R=T) + P(W=T|R=T,S=F)P(S=F|R=T)P(R=T)$

$= 0.99*0.01*0.2+0.8*0.99*0.2=0.16038$

$P(W=T)=P(W=T,S=T,R=T)+P(W=T,S=F,R=T)+P(W=T,S=T,R=F)+P(W=T,S=F,R=F)$

$=P(W=T|S=T,R=T)P(S=T|R=T)P(R=T) + P(W=T|S=F,R=T)P(S=F|R=T)P(R=T) + P(W=T|S=T,R=F)P(S=T|R=F)P(R=F)+P(W=T|S=F,R=F)P(S=F|R=F)P(R=F)$

$= 0.99*0.01*0.2+0.8*0.99*0.2+0.9*0.4*0.8+0.0*0.6*0.8=0.44838,$

### FluxRank

FluxRank 这一模块的触发需要服务指标（Service KPI）的异常，因此需要对服务指标（Service KPI）进行异常检测。这里的服务指标通常指的是业务指标，包括某块 APP 的在线人数，某个接口的成功率，某个视频网站的卡顿数等指标。当服务指标出现了异常的时候，就启动 FluxRank 模块进行异常机器定位。

1. 异常检测部分：通过设定阈值或者某个简单的规则来进行异常检测，包括服务的 KPI（Service KPI）和机器的 KPI（machine KPIs）；
2. 手工检查异常的时间段，并且查看在异常的时间段内发生了什么情况；
3. 运维人员根据自身的业务经验来对机器的故障程度做人工排序；
4. 运维人员根据自身的业务经验来对故障进行处理，并且人工给出处理方案。

1. 如何衡量海量 KPIs 的变化程度？在这里不仅有服务的 KPIs，还有机器的 KPIs。而机器的 KPIs 包括内存，硬盘，IO，CPU等诸多固定的指标，那么如何对这些海量的 KPI 曲线进行变化程度的衡量，为后续的指标排序做准备就成为了一个难点；
2. 如何对 KPIs 进行异常性或者重要性的聚类，让运维人员能够一眼看出每个聚簇的差异或者异常程度？
3. 如何对 KPIs 聚类的结果进行排序？

1. 基于 Kenel Density Estimation 用于衡量海量 KPIs 在某一个时间段的变化程度和异常程度；
2. 基于上一步生成的异常程度，对诸多机器所形成的特征使用距离公式或者相似度公式，然后使用 DBSCAN 聚类算法来对机器进行聚类；
3. 在排序部分，对上一步的机器聚类结果进行排序；

#### Change Quantification

1. absolute derivative 方法：个人理解就是对时间序列进行一阶差分操作，然后对一阶差分来做时间序列异常检测，例如 3-sigma 等方法，一旦有明显的变化，就说明当前的时间点出现了突增或者突降；与该方法比较类似的一种方法是：MAD（Median Absolute Deviation）。对于一条时间序列 $X=[x_{1},\cdots,x_{n}]$ 而言，MAD 定义为 $MAD = median_{1\leq i\leq n}(|x_{i}-median(X)|)$，而每个点的异常程度可以定义为：$s_{i}=(x_{i}-median(X))/MAD = (x_{i}-median(X))/median_{1\leq i\leq n}(|x_{i}-median(X)|).$$s_{i}$ 较大或者较小的时候，表示上涨或者下降的异常程度。通过设置相应的阈值，同样可以获得时间序列的异常开始时间。
2. CUSUM 算法也是用于时间序列异常检测的。对于一条时间序列 $X=[x_{1},x_{2},\cdots,x_{n}]$，可以预估它的目标值（target value）$\mu_{0}$，通常可以用均值来估计，也需要计算出这条时间序列的标准差 $\sigma$。通常设定 $\mu_{1}=\mu_{0}+\delta\sigma$$K=\delta\sigma/2=|\mu_{1}-\mu_{0}|/2$。而 Tabular CUSUM 指的是迭代公式 $C_{i}^{+}=\max[0,x_{i}-(\mu_{0}+K)+C_{i-1}^{+}]$$C_{i}^{-}=\max[0,(\mu_{0}-K)-x_{i}+C_{i-1}^{-}]$，初始值是 $C_{0}^{+}=C_{0}^{-}=0$。当累计偏差 $C_{i}^{+}$ 或者 $C_{i}^{-}$ 大于 $H=5\sigma$ 的时候，表示 $x_{i}$ 出现了异常，也就是 out of control。通过这个值，可以获得时间序列开始异常的时间。

#### Change Degree

$P_{o}(\{x_{j}\}|\{x_{i}\}) = \prod_{j=1}^{\ell}P(X\geq x_{j}|\{x_{i}\})$,

$P_{u}(\{x_{j}\}|\{x_{i}\}) = \prod_{j=1}^{\ell}P(X\leq x_{j}|\{x_{i}\})$,

$o=-\frac{1}{\ell}\sum_{j=1}^{\ell}\ln P(X\geq x_{j}|\{x_{i}\})$,

$u =-\frac{1}{\ell}\sum_{j=1}^{\ell}\ln P(X\leq x_{j}|\{x_{i}\})$.

Beta 分布的概率密度函数（probabilisty density function）是 $f(x;\alpha,\beta) = x^{\alpha-1}(1-x)^{\beta-1}/B(\alpha,\beta)$，其中 $B(\alpha,\beta)=\Gamma(\alpha)\Gamma(\beta)/\Gamma(\alpha+\beta)$。在机器 KPIs 中，CPU 等指标可以用 Beta 分布；

#### Digest Distillation

Pearson 系数指的是：$\rho_{X,Y}=\sum_{i=1}^{n}(x_{i}-\overline{x})\cdot(y_{i}-\overline{y})/\sqrt{\sum_{i=1}^{n}(x_{i}-\overline{x})^{2}\cdot\sum_{i=1}^{n}(y_{i}-\overline{y})^{2}},$ 其中 $\overline{x}=\sum_{i=1}^{n}x_{i}/n$$\overline{y}=\sum_{i=1}^{n}y_{i}/n$

Kendall tau 系数指的是：如果 ($x_{i}>x_{j}$$y_{i}>y_{j}$) 或者 ($x_{i}$y_{i})，那么称之为 concordant；如果 ($x_{i}$y_{i}>y_{j}$) 或者 ($x_{i}>x_{j}$$y_{i})，称之为 discordant；如果 $x_{i}=x_{j}$ 或者 $y_{i}=y_{j}$，则既不是 concordant，也不是 discordant。那么 Kendall tau 定义为 $[\text{(number of concordant pairs)}-\text{(number of disordant paris)}] / [n(n-1)/2]$

Spearman 系数指的是：通过原始序列变成秩次变量（rank）（从大到小降序排列即可），$x_{i}$ 将会对应到 $x_{i}'$，后者表示 $x_{i}$ 在从大到小排序之后的序列 $\{x_{i}\}_{1\leq i\leq n}$ 的位置，称之为秩次（rank），得到序列 $X'=[x_{1}',\cdots,x_{n}']$。对原始序列 $Y=[y_{1},\cdots,y_{n}]$ 作同样的操作，得到 $Y'=[y_{1}',\cdots,y_{n}']$。一个相同的值在一列数据中必须有相同的秩次，那么在计算中采用的秩次就是数值在按从大到小排列时所在位置的平均值。如果没有相同的 rank，那么使用公式 $r_{s} = 1-6\sum_{i=1}^{n}d_{i}^{2}/(n(n^{2}-1))$ 进行计算，其中 $d_{i}=x_{i}'-y_{i}'$；如果存在相同的秩次，则对 $X'=[x_{1}',\cdots,x_{n}']$$Y'=[y_{1}',\cdots,y_{n}']$ 来做 Pearson 系数即可，也就是 $\rho_{X',Y'}$

#### Digest Ranking

1. 变化开始时间（change start time）$T_{c}$ 会在失败发生时间 $T_{f}$ 之前；
2. 不同的故障机器 KPIs 的 change start time 是非常接近的；
3. 故障机器的一些 KPIs 的 change degree 是非常大的；
4. 故障机器的占比是与故障原因相关的，故障机器越多说明故障越大；

### 参考资料

1. FluxRank: A Widely-Deployable Framework to Automatically Localizing Root Cause Machines for Software Service Failure Mitigation，Ping Liu，Yu Chen，Xiaohui Nie，Jing Zhu，Shenglin Zhang，Kaixin Sui，Ming Zhang，Dan Pei，ISSRE 2019， Berlin, Germany, Oct 28-31, 2019。
2. Introduction to Statistical Quality Control，6th edition，Douglas C.Montgomery。
3. Bayesian Network：https://en.wikipedia.org/wiki/Bayesian_network

# 一篇关于时间序列异常检测的论文

### 数据集的情况

2. 通过算法获得每一个类的质心，并且标记出质心曲线的异常点和正常点；
3. 对新来的时间序列进行实时聚类，划分到合适的类别；
4. 基于新来的时间序列（没有标记）和历史上的时间序列（有标记）使用无监督算法来重新训练一个新的模型，进行该类别的时间序列异常检测。

1. CPLE 是半监督学习算法中比较健壮的，因为它并没有过多的假设条件，并且也符合这篇文章的业务场景，同时拥有质心曲线（有标记）和新来的曲线（无标记），使用半监督学习也是符合常理的。除了 CPLE，其实在实战过程中也可以多尝试其他的半监督模型，具体可以参考周志华的《机器学习》。
2. CPLE 的复杂度比较低，计算快。

# 两篇关于时间序列的论文

### Robust and Rapid Clustering of KPIs for Large-Scale Anomaly Detection

-> 根据每一类时间序列使用不同的异常检测模型

1. 形状：通常来说，时间序列随着业务的变化，节假日效应，变更的发布，将会随着时间的迁移而造成形状的变化。
2. 噪声：无论是从数据采集的角度，还是系统处理的角度，甚至服务器的角度，都有可能给时间序列带来一定的噪声数据，而噪声是需要处理掉的。
3. 平移：在定时任务中，有可能由于系统或者人为的原因，时间序列的走势可能会出现一定程度的左右偏移，有可能每天 5:00 起的定时任务由于前序任务的原因而推迟了。
4. 振幅：通常时间序列都存在一条基线，而不同的时间序列有着不同的振幅，振幅决定了这条时间序列的振荡程度，而振幅或者基线其实也是会随着时间的迁移而变化的。

#### 预处理模块（Preprocessing）

1. 缺失值：缺失值指的是在该上报数据的时间戳上并没有相应的数据上报，数据处于缺失的状态。通常的办法就是把数据补齐，而数据补齐的方法有很多种，最简单的就是使用线性插值（Linear Interpolation）的方式来补齐。
2. 标准化：对于一个时间序列而言，有可能它的均值是10万，有可能只有10，但是它们的走势有可能都是一样的。所以在这个时候需要进行归一化的操作。最常见的有两种归一化方法：一种是标准化，另外一种是最大最小值归一化。如果 $[x_{1},x_{2}, \cdots, x_{n}]$ 表示原始的时间序列的话，标准化指的是 $\hat{x}_{i} = (x_{i} -\mu) /\sigma$，其中 $\mu$$\sigma$ 分别表示均值和标准差。最大最小值归一化指的是 $\hat{x}_{i} = (x_{i} - \min) / (\max - \min)$，其中 $\max, \min$ 分别表示这段时间内的最大值与最小值。

#### 基线提取模块（Baseline Extraction）

$x_{i} = baseline_{i} + residual_{i},$

1. 异常值的处理：通常需要移除一些明显异常的值，然后使用线性插值等方法来把这些移除的值补上。
2. 使用简单的移动平均算法加上一个窗口值 $w$ 来提取基线。假设时间序列是 $[x_{1},\cdots,x_{n}]$$SMA_{i} = \sum_{j=i-w+1}^{i} x_{j} / w$$R_{i} = x_{i} - SMA_{i}$。也就是说 $x_{i} = SMA_{i} + R_{i}$
3. 提取基线的方式其实还有很多，使用带权重的移动平均算法（Weighted Moving Average），指数移动平均算法（Exponentially Weighted Moving Average）都可以提取基线，甚至使用深度学习中的 Autoencoder 或者 VAE 算法都能够提取基线。

#### 基于密度的聚类算法（Clustering）

1. 通过特征工程的方法，从时间序列中提取出必要的时间序列特征，然后使用 KMeans 等算法来进行聚类。
2. 通过相似度计算工具，对比两条时间序列之间的相似度，相似的聚成一类，不相似的就分成两类。

$CC_{s}(X,Y) = \begin{cases} \sum_{i=1}^{m-s}x_{s+i}\cdot y_{i}, \text{ where } s\geq 0, \\ \sum_{i=1}^{m+s}x_{i} \cdot y_{i-s} \text{ where } s<0.\end{cases}$

$NCC(X,Y) = \max_{s}\bigg(\frac{CC_{s}(X,Y)}{||X||_{2}\cdot||Y||_{2}}\bigg),$

$SBD(X,Y) = 1 - NCC(X,Y).$

$\text{Centroid} = argmin_{X \in \text{cluster}_{i}} \sum_{Y \in \text{cluster}_{i}} SBD(X,Y)^{2}.$

#### 聚类与时间序列异常检测

1. 先对时间序列进行聚类或者分类，针对不同的时间序列类型来使用不同的模型；
2. 在时间序列异常检测中加入分类特征。

### Robust and Rapid Adaption for Concept Drift in Software System Anomaly Detection

Insight 1：Concept drift detection can be converted into a spike detection problem in the change score stream.

Insight 2：The relationship between new concept and old concept can be almost perfectly fitted by linear regression.

1. 存在一个线性关系 $y = kx$ 使得二维点集 $\{(x_{i},y_{i}), 1\leq i\leq n\}$ 能够被很好的拟合好，也就是说此刻的方差较小。
2. $[x_{1},\cdots,x_{n}]$$[y_{1},\cdots,y_{n}]$ 的 Pearson 系数很高；
3. $[x_{1},\cdots,x_{n}]$$[y_{1},\cdots,y_{n}]$ 标准化（z-score）之后的曲线几乎一致。
4. 分别存在两个值 $\mu_{1},\mu_{2}$ 使得 $[x_{1}/\mu_{1},\cdots,x_{n}/\mu_{1}]$$[y_{1}/\mu_{2},\cdots,y_{n}/\mu_{2}]$ 几乎一致。

1. 计算机视觉方向
2. 自然语言处理方向
3. 语音识别方向
4. 机器学习方向

1. 质量保障；
2. 效率提升；
3. 成本管理。

# 基于自编码器的时间序列异常检测算法

### 自编码器

$\phi: X\rightarrow F,$

$\psi: F\rightarrow X,$

$\phi,\psi = argmin_{\phi,\psi}||X-(\psi\circ\phi)X||^{2},$

$z = f(Ax+c),$

$x'=g(Bx+d),$

Claim. 对于自编码器而言，其中隐藏层的维度 $m$ 一定是要小于输入层的维度 $n$ 的。

Proof. 如果 $n=m$，那么令 $A=B=I_{n},$ $c=d=0,$ $f=g=id$ 就可以得到一个自编码器，而这个自编码器对于提取特征没有任何的意义。同理，当 $m>n$ 时，$A$ 是一个 $m\times n$ 矩阵，$B$ 是一个 $n\times m$ 矩阵。从线性代数的角度来看，有无数个矩阵 $A, B$ 满足 $BA=I_{n}$。这种情况下对于提取特征也是没有意义的。而当 $m 时，其实无法找到矩阵 $A,B$ 使得 $BA=I_{n}.$ 如果存在 $BA=I_{n},$ 那么

$n = rank(I_{n})=rank(BA) \leq \min\{rank(A),rank(B)\} \leq \min\{m,m\}=m.$

### 时间序列异常检测：

-> Auto Encoder（Encoder 和 Decoder）

-> 重构后的时间序列

-> 通过重构后的时间序列与原始时间序列的整体误差和局部误差来判断异常点

### 总结

1. 从理论上说，它只能对一个时间序列单独训练一个模型，不同类型的时间序列需要使用不同的模型。这样的话，其实维护模型的成本比较高，不太适用于大规模的时间序列异常检测场景；
2. 对周期型的曲线效果比较好，如果是毛刺型的数据，有可能就不太适用；因为长期的毛刺型数据就可以看成正常的数据了。
3. 每次调参需要人为设置一定的阈值，不同的时间序列所需要的阈值是不一样的。

### 参考文献

1. Unsupervised Anomaly Detection via Variational Auto-Encoder for Seasonal KPIs in Web Applications, Haowen XU, etc., 2018
2. Deep Learning, Ian Goodfellow, etc., 2016
3. https://zr9558.com/2016/06/12/replicator-neural-networks/

# 黑盒函数的探索

### 黑盒函数的研究对象

1. 最大值与最小值，i.e. $\max f(x)$$\min f(x).$
2. 根，i.e. $\{x:f(x) = 0\}.$
3. 函数的单调性与凹凸性等。

### 黑盒函数的根

$x = \frac{-b\pm\sqrt{b^{2}-4ac}}{2a}.$

#### 二分法

[介值定理] 如果连续函数 $f(x)$ 的定义域包含 $[a,b],$ 而且通过 $(a,f(a))$$(b,f(b))$ 两点，它也必定通过区间 $[a,b]$ 内的任意一点 $(c,f(c)),$ 其中 $a

1. 先找到一个区间 $[a,b],$ 使得 $f(a)f(b)<0;$
2. 求这个区间的中点 $m=(a+b)/2,$ 并求出 $f(m)$ 的取值；
3. 如果 $f(m)=0,$ 那么 $m$ 就是函数的根；如果 $f(m)f(a)>0,$ 就选择 $[m,b]$ 为新的区间，否则选择 $[a,m]$ 为新的区间；
4. 重复第 2 步和第 3 步直到达到最大迭代次数或者最理想的精度为止。

#### 牛顿法（Newton’s Method）

$x_{1} = x_{0} - f(x_{0})/f'(x_{0}).$

$x_{n+1}= x_{n}-f(x_{n})/f'(x_{n}).$

#### 割线法

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

$f'(x_{0}) \approx \frac{f(x)-f(x_{0})}{x-x_{0}}.$

$x_{n+1} = x_{n} - \frac{x_{n}-x_{n-1}}{f(x_{n})-f(x_{n-1})}f(x_{n}).$

### 黑盒函数的最大值与最小值

#### Weierstrass 逼近定理

1. 闭区间上的连续函数可以用多项式级数一致逼近；
2. 闭区间上的周期为 $2\pi$ 的连续函数可以用三角函数级数一致逼近。

[Weierstrass 逼近定理] 假设 $f(x)$ 是闭区间 $[a,b]$ 连续的实函数。对于任意的 $\epsilon>0$，存在一个多项式 $p(x)$ 使得对于任意的 $x\in[a,b],$$|f(x)-p(x)|<\epsilon.$

#### Lagrange 插值公式

$p(x) = \sum_{i=0}^{n}\bigg(\prod_{0\leq j\leq n, j\neq i}\frac{x-x_{j}}{x_{i}-x_{j}}\bigg)y_{i}.$

#### 粒子群算法（Particle Swarm Optimization）

PSO 最初是为了模拟鸟群等动物的群体运动而形成的一种优化算法。PSO 算法是假设有一个粒子群，根据群体粒子和单个粒子的最优效果，来调整每一个粒子的下一步行动方向。假设粒子的个数是 $N_{p}$，每一个粒子 $\bold{x}_{i}\in \mathbb{R}^{n}$ 都是 $n$ 维欧几里德空间里面的点，同时需要假设粒子的速度 $\bold{v}_{i}\in\mathbb{R}^{n}$。在每一轮迭代中，需要更新两个最值，分别是每一个粒子在历史上的最优值和所有粒子在历史上的最优值，分别记为 $\bold{x}_{i}^{*}$$1\leq i \leq N_{p}$）和 $\bold{x}^{g}$。在第 $t+1$ 次迭代的时候，

$\bold{v}_{i}(t+1) = \bold{v}_{i}(t) + c r_{1}[\bold{x}_{i}^{*}(t) - \bold{x}_{i}(t)] + c r_{2}[\bold{x}^{g}(t) - \bold{x}_{i}(t)],$

$\bold{x}_{i}(t+1) = \bold{x}_{i}(t)+\bold{v}_{i}(t+1), \text{ } 1\leq i\leq N_{p}.$

#### 模拟退火（Simulated Annealing）

Repeat：

a. Repeat:

i. 进行一个随机扰动 $\bold{x} = \bold{x} + \Delta \bold{x}$

ii. 计算 $\Delta E(\bold{x}) = E(\bold{x}+\Delta\bold{x}) - E(\bold{x})$

b. 令 $T = T-\Delta T$

# 时间序列的自回归模型—从线性代数的角度来看

### Fibonacci 序列

$F_{n}=F_{n-1}+F_{n-2}$

0，1，1，2，3，5，8，13，21，34，55，89，144，…

### 求解 Fibonacci 序列的通项公式 －－－ 矩阵对角化

$\left( \begin{array}{c} F_{n+2} \\ F_{n+1} \\ \end{array}\right)= \left( \begin{array}{cc} 1 & 1 \\ 1 & 0 \\ \end{array}\right) \left( \begin{array}{c} F_{n+1} \\ F_{n} \\ \end{array}\right) = A \left( \begin{array}{c} F_{n+1} \\ F_{n} \\ \end{array}\right),$

$\left( \begin{array}{c} F_{n} \\ F_{n-1} \\ \end{array}\right)= A \left( \begin{array}{c} F_{n-1} \\ F_{n-2} \\ \end{array}\right) = \cdots = A^{n-1} \left( \begin{array}{c} F_{1} \\ F_{0} \\ \end{array}\right).$

$P^{-1}AP = diag(\lambda_{1},\cdots,\lambda_{m})$

$\implies AP = P diag(\lambda_{1},\cdots,\lambda_{m})$

$\implies A = P diag(\lambda_{1},\cdots,\lambda_{m})P^{-1}.$

$A^{k} = (PDP^{-1})\cdots(PDP^{-1}) = P D^{k} P^{-1}= P diag(\lambda_{1}^{k},\cdots,\lambda_{m}^{k})P^{-1}.$

$det(\lambda I - A) = 0,$

$A = \left( \begin{array}{cc} 1 & 1 \\ 1 & 0 \\ \end{array}\right)$,

$\vec{\alpha}_{1} = (\phi,1)^{T}, \vec{\alpha}_{2} = (-\phi^{-1},1)$.

$F_{k} = \frac{1}{\sqrt{5}}\bigg(\frac{1+\sqrt{5}}{2}\bigg)^{k} - \frac{1}{\sqrt{5}}\bigg(\frac{1-\sqrt{5}}{2}\bigg)^{k}=\frac{\phi^{k}-(-\phi)^{-k}}{\sqrt{5}}$.

### 时间序列的弱平稳性

1. $E(x_{t})$ 对于所有的 $t\geq 0$ 都是恒定的；
2. $Var(x_{t})$ 对于所有的 $t\geq 0$ 都是恒定的；
3. $x_{t}$$x_{t-h}$ 的协方差对于所有的 $t\geq 0$ 都是恒定的。

$ACF(x_{t},x_{t-h}) = \frac{Covariance(x_{t},x_{t-h})}{\sqrt{Var(x_{t})\cdot Var(x_{t-h})}}$.

$ACF(x_{t},x_{t-h}) = \frac{Covariance(x_{t},x_{t-h})}{Var(x_{t})}$.

### 时间序列的自回归模型（AutoRegression Model）

#### AR(1) 模型

AR(1) 模型指的是时间序列 $\{x_{t}\}_{t\geq 0}$ 在时间戳 $t$ 时刻的取值 $x_{t}$ 与时间戳 $t - 1$ 时刻的取值 $x_{t-1}$ 相关，其公式就是：

$x_{t}=\delta+\phi_{1}x_{t-1}+w_{t}$

1. $w_{t}\sim N(0,\sigma_{w}^{2})$，并且 $w_{t}$ 满足 iid 条件。其中 $N(0,\sigma_{w}^{2})$ 表示 Gauss 正态分布，它的均值是0，方差是 $\sigma_{w}^{2}$
2. $w_{t}$$x_{t}$ 是相互独立的（independent）。
3. $x_{0},x_{1},\cdots$弱平稳的，i.e. 必须满足 $|\phi_{1}|<1$

1. $E(x_{t}) = \delta/(1-\phi_{1})$.
2. $Var(x_{t}) = \sigma_{w}^{2}/(1-\phi_{1}^{2})$.
3. $Covariance(x_{t},x_{t-h}) = \phi_{1}^{h}$.

Proof of 1. 从 AR(1) 的模型出发，可以得到

$E(x_{t}) = E(\delta + \phi_{1}x_{t-1}+w_{t}) = \delta + \phi_{1}E(x_{t-1}) = \delta + \phi_{1}E(x_{t})$,

Proof of 2. 从 AR(1) 的模型出发，可以得到

$Var(x_{t}) = Var(\delta + \phi_{1}x_{t-1}+w_{t})$

$= \phi_{1}^{2}Var(x_{t-1}) +Var(w_{t}) = \phi_{1}^{2}Var(x_{t}) + \sigma_{w}^{2}$,

Proof of 3.$\mu = E(x_{t}), \forall t\geq 0$. 从 $x_{t}$ 的定义出发，可以得到：

$x_{t}-\mu = \phi_{1}(x_{t-1}-\mu)+w_{t}$

$= \phi_{1}^{h}(x_{t-h}-\mu) + \phi_{1}^{h-1}w_{t-h+1}+\cdots+\phi_{1}w_{t-1}+w_{t},$

$\rho_{h} = Covariance(x_{t},x_{t-h}) = \frac{E((x_{t}-\mu)\cdot(x_{t-h}-\mu))}{Var(x_{t})}=\phi_{1}^{h}$.

#### AR(1) 模型与一维动力系统

$f(x) = \phi_{1}x + \delta,$

Method 1.

$f^{n}(x) = \phi_{1}^{n}x+ \frac{1-\phi_{1}^{n}}{1-\phi_{1}}\delta$,

$n\rightarrow \infty$，可以得到 $f^{n}(x)\rightarrow \delta/(1-\phi_{1})$。这与 $E(x_{t}) = \delta/(1-\phi_{1})$ 其实是保持一致的。

Method 2.

$f(x)-\frac{\delta}{1-\phi_{1}} = \phi_{1}(x-\frac{\delta}{1-\phi_{1}})$

$\implies |f(x)-\frac{\delta}{1-\phi_{1}}| <\frac{1+|\phi_{1}|}{2}\cdot|x-\frac{\delta}{1-\phi_{1}}|$

$\implies |f^{n}(x)-\frac{\delta}{1-\phi_{1}}|<\bigg(\frac{1+|\phi_{1}|}{2}\bigg)^{n}\cdot|x-\frac{\delta}{1-\phi_{1}}|$.

$f(x)-\frac{\delta}{1-\phi_{1}} = \phi_{1}(x-\frac{\delta}{1-\phi_{1}})$

$\implies |f(x)-\frac{\delta}{1-\phi_{1}}| >\frac{1+|\phi_{1}|}{2}\cdot|x-\frac{\delta}{1-\phi_{1}}|$

$\implies |f^{n}(x)-\frac{\delta}{1-\phi_{1}}|>\bigg(\frac{1+|\phi_{1}|}{2}\bigg)^{n}\cdot|x-\frac{\delta}{1-\phi_{1}}|$.

#### AR(p) 模型

1. AR(1) 模型形如：

$x_{t}=\delta+\phi_{1}x_{t-1}+w_{t}.$

2. AR(2) 模型形如：

$x_{t} = \delta + \phi_{1}x_{t-1}+\phi_{2}x_{t-2}+w_{t}.$

3. AR(p) 模型形如：

$x_{t} = \delta + \phi_{1}x_{t-1}+\phi_{2}x_{t-2}+\cdots+\phi_{p}x_{t-p}+w_{t}.$

#### AR(p) 模型的稳定性 －－－ 基于线性代数

$x_{t}= \phi_{1}x_{t-1} + \phi_{2}x_{t-2}$.

$\left( \begin{array}{c} x_{t+2} \\ x_{t+1} \\ \end{array}\right)= \left( \begin{array}{cc} \phi_{1} & \phi_{2} \\ 1 & 0 \\ \end{array}\right) \left( \begin{array}{c} x_{t+1} \\ x_{t} \\ \end{array}\right) = A \left( \begin{array}{c} x_{t+1} \\ x_{t} \\ \end{array}\right).$

$x_{t} = \phi_{1}x_{t-1}+\phi_{2}x_{t-2}+\cdots+\phi_{p}x_{t-p}.$

$\left(\begin{array}{c} x_{t+p} \\ x_{t+p-1} \\ \vdots \\ x_{t+1}\\ \end{array}\right) = \left(\begin{array}{ccccc} \phi_{1} & \phi_{2} &\cdots & \phi_{p-1} & \phi_{p} \\ 1 & 0 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 1 & 0 \\ \end{array}\right) \left(\begin{array}{c} x_{t+p-1} \\ x_{t+p-2} \\ \vdots \\ x_{t} \\ \end{array}\right) = A \left(\begin{array}{c} x_{t+p-1} \\ x_{t+p-2} \\ \vdots \\ x_{t} \\ \end{array}\right)$

$\lambda^{p}-\phi_{1}\lambda^{p-1}-\phi_{2}\lambda^{p-2}-\cdots-\phi_{p}=0.$

$x_{t} = \phi_{1}x_{t-1}+\phi_{2}x_{t-2}+\cdots+\phi_{p}x_{t-p}$

# 时间序列的搜索

### 基于动态规划的相似度计算方法

#### DTW 的经典算法

$W^{*} = argmin_{W}(\sqrt{\sum_{k=1}^{K}w_{k}})$.

$DTW(i,j)$

$= dist(q_{i},c_{j}) + \min(DTW(i-1,j-1), DTW(i-1,j), DTW(i,j-1)).$

1. 如果 $dist(q_{i},c_{j}) = (q_{i}-c_{j})^{2}$ 时，则 $\sqrt{DTW(n,n)}$ 表示最后的距离；
2. 如果 $dist(q_{i},c_{j}) = |q_{i}-c_{j}|$ 时，则 $DTW(n,n)$ 表示最后的距离。
3. 如果 $dist(q_{i},c_{j}) = |q_{i}-c_{j}|^{p}$ 时，则 $(DTW(n,n))^{1/p}$ 表示最后的距离。

Remark.

DTW 算法不满足三角不等式。例如：$x={0}, y={1,2}, z={1,2,2}$，则

$DTW(x,z)=5>DTW(x,y)+DTW(y,z) = 3 + 0 =3.$

#### DTW 的加速算法

$DTW(i,j) = dist(q_{i},c_{j}) + \min(DTW(i-1,j-1), DTW(i-1,j), DTW(i,j-1)),$

### 相似时间序列的搜索

Question 1. 给定一条时间序列 $Q$ 和一个时间序列的数据库 $D=\{C_{i}\}_{i=1}^{\infty}$。通过某种相似度或者距离计算方法，计算出给定的时间序列 $Q$ 与时间序列数据库中 $D$ 中最相似的时间序列。

Question 2. 给定一条时间序列 $Q$ 和一个时间序列的数据库 $D=\{C_{i}\}_{i=1}^{\infty}$，以及正整数 $K$。从数据库 $D$ 中寻找与给定的时间序列 $Q$ 最相似的 $K$ 条时间序列。

#### DTW 算法的下界 LB_Kim

$LBKim(Q,C)$

$= \max\{|\max(Q)-\max(C)|,|\min(Q)-\min(C)|,|First(Q)-First(C)|,|Last(Q)-Last(C)| \}.$

#### DTW 算法的下界 LB_Keogh

（1）首先定义时间序列 Q 的上下界。i.e. $Q:\{q_{1},\cdots,q_{n}\}$，给定一个窗口的取值 r，得到 $U_{i}=\max(q_{i-r},q_{i+r})$$L_{i} = \min(q_{i-r},q_{i+r})$

（2）定义公式：

$LBKeogh(Q,C)$

$= \sqrt{\sum_{i=1}^{n}(c_{i}-U_{i})^{2}I_{\{c_{i}>U_{i}\}} + \sum_{i=1}^{n}(c_{i}-L_{i})^{2}I_{\{c_{i}

Lower Bounding Sequential Scan(Q)
best_so_far = infinity
for all sequences in database
LB_dist = lower_bound_distance(C_{i},Q)
if LB_dist < best_so_far
true_dist = DTW(C_{i},Q)
if true_dist < best_so_far
best_so_far = true_dist
index_of_best_match = i
endif
endif
endfor

# 如何理解时间序列？— 从Riemann积分和Lebesgue积分谈起

Riemann 积分和 Lebesgue 积分是数学中两个非常重要的概念。本文将会从 Riemann 积分和 Lebesgue 积分的定义出发，介绍它们各自的性质和联系。

## 积分

### Riemann 积分

Riemann 积分虽然被称为 Riemann 积分，但是在 Riemann 之前就有学者对这类积分进行了详细的研究。早在阿基米德时代，阿基米德为了计算曲线 $x^{2}$ 在 [0,1] 区间上与 X 坐标轴所夹的图形面积，就使用了 Riemann 积分的思想。 他把 [0,1] 区间等长地切割成 n 段，每一段使用一个长方形去逼近 $x^{2}$ 这条曲线的分段面积，再把 n 取得很大，所以得到当 n 趋近于无穷的时候，就知道该面积其实是 1/3。

$\sum_{i=0}^{n-1}f(t_{i})(x_{i+1}-x_{i}).$

$|\sum_{i=0}^{n-1}f(t_{i})(x_{i+1}-x_{i}) - s|<\epsilon.$

### Lebesgue 积分

Riemann 积分是为了计算曲线与 X 轴所围成的面积，而 Lebesgue 积分也是做同样的事情，但是计算面积的方法略有不同。要想直观的解释两种积分的原理，可以参见下图：

Riemann 积分是把一条曲线的底部分成等长的区间，测量每一个区间上的曲线高度，所以总面积就是这些区间与高度所围成的面积和。

Lebesgue 积分是把曲线化成等高线图，每两根相邻等高线的差值是一样的。每根等高线之内含有它所圈着的长度，因此总面积就是这些等高线内的面积之和。

1. Riemann 积分：从一个角落开始一口一口吃，每口都包含所有的配料；
2. Lebesgue 积分：从最上层开始吃，按照“面包-配菜-肉-蛋-面包”的节奏，一层一层来吃。

$\int(\sum_{k}a_{k}1_{S_{k}})d\mu = \sum_{k}a_{k}\int 1_{S_{k}}d\mu = \sum_{k}a_{k}\mu(S_{k}).$

$\int_{E}f d\mu = \sup\{\int_{E}sd\mu: \bold{0}\leq s\leq f\}$,

$\int fd\mu = \int f^{+}d\mu - \int f^{-}d\mu.$.

### Riemann 积分与Lebesgue 积分的关系

$(R)\int_{a}^{b}f(x)dx = (L)\int_{[a,b]}f(x)dx$.

1. $x$ 是有理数时，$D(x) = 1$
2. $x$ 是无理数时，$D(x) = 0$.

Dirichlet 函数是定义在实数轴的函数，并且值域是 $\{0,1\}$，无法画出函数图像，它不是 Riemann 可积的，但是它 Lebesgue 可积。

## 时间序列

### 时间序列的表示 — 基于 Riemann 积分

1. 分段线性逼近（Piecewise Linear Approximation）
2. 分段聚合逼近（Piecewise Aggregate Approximation）
3. 分段常数逼近（Piecewise Constant Approximation）

#### 分段聚合逼近（Piecewise Aggregate Approximation）— 类似 Riemann 积分

$\overline{x}_{i} = \frac{w}{N} \cdot \sum_{j=\frac{N}{w}(i-1)+1}^{\frac{N}{w}i} x_{j}$.

#### 符号特征（Symbolic Approximation）— 类似用简单函数来计算 Lebesgue 积分

SAX 方法的流程如下：

1. 正规化（normalization）：把原始的时间序列映射到一个新的时间序列，新的时间序列满足均值为零，方差为一的条件。
2. 分段表示（PAA）：$\{x_{1},\cdots, x_{N}\} \Rightarrow \{\overline{x}_{1},\cdots,\overline{x}_{w}\}$
3. 符号表示（SAX）：如果 $\overline{x}_{i}，那么 $\hat{X}_{i}=l_{1}$；如果 $z_{(j-1)/\alpha}\leq \overline{x}_{i}，那么 $\hat{X}_{i} = l_{j}$，在这里 $2\leq j\leq \alpha$；如果 $\overline{x}_{i}\geq z_{(\alpha-1)/\alpha}$，那么 $\hat{X}_{i} = l_{\alpha}$

### 时间序列的表示 — 基于 Lebesgue 积分

#### 熵（Entropy）

$\text{entropy}(X) = -\sum_{i=1}^{\infty}P\{x=x_{i}\}\ln(P\{x=x_{i}\})$.

#### 分桶熵（Binned Entropy）

$\text{binned entropy}(X) = -\sum_{k=0}^{\min(maxbin, len(X))} p_{k}\ln(p_{k})\cdot 1_{(p_{k}>0)},$

# 时间序列的相似性

#### 度量空间

1. $d(x,y)\geq 0$，并且 $d(x,y)=0$ 当且仅当 $x=y$;
2. $d(x,y)=d(y,x)$，也就是满足对称性；
3. $d(x,z)\leq d(x,y)+d(y,z)$，也就是三角不等式。

Remark.

#### 内积空间

1. 对于任意的 $x,y\in V$，有 $ =\overline{}.$

2. 共轭双线性形式指的是：

$\forall a\in F, \forall x,y\in V, =a,$

$\forall x,y,z\in F, = + .$

$\forall b\in F, \forall x,y\in V, = \overline{b},$

$\forall x,y,z\in F, = +.$

3. 非负性：$\forall x\in V, \geq 0.$

4. 非退化：从 $V$ 到对偶空间 $V^{*}$ 的映射：$x\mapsto$ 是同构映射。

### 基于欧几里德距离的相似度计算

$d_{L^{1}}(X_{T},Y_{T}) = \sum_{t=1}^{T}|x_{t}-y_{t}|,$

$d_{L^{p}}(X_{T}, Y_{T}) = (\sum_{t=1}^{T}|x_{t}-y_{t}|^{p})^{1/p},$

$d_{L^{2}}(X_{T}, Y_{T}) = (\sum_{t=1}^{T}|x_{t}-y_{t}|^{2})^{1/2},$

$d_{L^{\infty}}(X_{T},Y_{T}) = \max_{1\leq t\leq T}|x_{t}-y_{t}|.$

### 基于相关性的相似度计算方法

#### Pearson 系数（Pearson Coefficient）

$\text{COR}(X_{T},Y_{T}) = \frac{\sum_{t=1}^{T}(x_{t}-\overline{X}_{T})\cdot(y_{t}-\overline{Y}_{T})}{\sqrt{\sum_{t=1}^{T}(x_{t}-\overline{X}_{T})^{2}}\cdot\sqrt{\sum_{t=1}^{T}(y_{t}-\overline{Y}_{T})^{2}}},$

$\overline{X}_{T} = \sum_{t=1}^{T}x_{t}/T$, $\overline{Y}_{T} = \sum_{t=1}^{T}y_{t}/T$

Pearson 系数的性质如下：

1. 如果两条时间序列 $X_{T} = Y_{T}$，则 $\text{COR}(X_{T},Y_{T}) =1$ 表是它们是完全一致的，如果两条时间序列 $X_{T} = -Y_{T}$，则 $\text{COR}(X_{T},Y_{T}) = -1$ 表示它们之间是负相关的。
2. $-1\leq \text{COR}(X_{T},Y_{T})\leq 1$.

$d_{COR,1}(X_{T},Y_{T}) = \sqrt{2\cdot(1-COR(X_{T},Y_{T}))},$

$d_{COR,2}(X_{T},Y_{T}) = \sqrt{\big(\frac{1-COR(X_{T},Y_{T})}{1+COR(X_{T},Y_{T})}\big)^{\beta}},$

#### The First Order Temporal Correlation Coefficient

$\text{CORT}(X_{T},Y_{T}) = \frac{\sum_{t=1}^{T-1}(x_{t+1}-x_{t})\cdot(y_{t+1}-y_{t})}{\sqrt{\sum_{t=1}^{T-1}(x_{t+1}-x_{t})^{2}}\cdot\sqrt{\sum_{t=1}^{T-1}(y_{t+1}-y_{t})^{2}}},$

$\text{CORT}(X_{T},Y_{T})$ 的性质：

1. $-1\leq \text{CORT}(X_{T},Y_{T}) \leq 1$
2. $\text{CORT}(X_{T},Y_{T}) =1$ 表示两条时间序列持有类似的趋势， 它们会同时上涨或者下跌，并且涨幅或者跌幅也是类似的。
3. $\text{CORT}(X_{T},Y_{T})=-1$ 表示两条时间序列的上涨和下跌趋势恰好相反。
4. $\text{CORT}(X_{T},Y_{T})=0$ 表示两条时间序列在单调性方面没有相关性。

$d_{CORT}(X_{T},Y_{T}) = \phi_{k}[CORT(X_{T},Y_{T})]\cdot d(X_{T},Y_{T}),$

$\phi_{k}(u) = 2/(1+\exp(ku)), k\geq 0$

#### 基于自相关系数的距离（Autocorrelation-based distance)

$\hat{\rho}_{k} = \frac{1}{(T-k)\sigma^{2}}\sum_{t=1}^{T-k}(x_{t}-\mu)\cdot(x_{t+k}-\mu)$,

$\hat{\rho}_{X_{T}} = (\hat{\rho}_{1,X_{T}},\cdots,\hat{\rho}_{L,X_{T}})^{T}\in \mathbb{R}^{L},$

$\hat{\rho}_{Y_{T}} = (\hat{\rho}_{1,Y_{T}},\cdots,\hat{\rho}_{L,Y_{T}})^{T}\in\mathbb{R}^{L}.$

$d_{ACF}(X_{T},Y_{T}) = \sqrt{(\hat{\rho}_{X_{T}}-\hat{\rho}_{Y_{T}})^{T}\Omega(\hat{\rho}_{X_{T}}-\hat{\rho}_{Y_{T}})}$.

（1）$\Omega = I_{L}$ 表示单位矩阵。用公式表示就是

$d_{ACFU}(X_{T},Y_{T}) =\sqrt{\sum_{i=1}^{L}(\hat{\rho}_{i,X_{T}}-\hat{\rho}_{i,Y_{T}})^{2}}$.

（2）$\Omega = diag\{p(1-p),p(1-p)^{2},\cdots,p(1-p)^{L}\}$ 表示一个 $L\times L$ 的对角矩阵，其中 $0。此时相当于一个带权重的求和公式。

$d_{ACFU}(X_{T},Y_{T}) =\sqrt{\sum_{i=1}^{L}p(1-p)^{i}(\hat{\rho}_{i,X_{T}}-\hat{\rho}_{i,Y_{T}})^{2}}$.

### 基于周期性的相似度计算方法

$I_{X_{T}}(\lambda_{k}) = T^{-1}|\sum_{t=1}^{T}x_{t}e^{-i\lambda_{k}t}|^{2}$,

$I_{Y_{T}}(\lambda_{k}) = T^{-1}|\sum_{t=1}^{T}y_{t}e^{-i\lambda_{k}t}|^{2}$.

（1）用原始的特征来表示距离：

$d_{P}(X_{T},Y_{T}) = \frac{1}{n}\sqrt{\sum_{k=1}^{n}(I_{X_{T}}(\lambda_{k})-I_{Y_{T}}(\lambda_{k}))^{2}}$.

（2）用正则化之后的特征来描述就是：

$d_{P}(X_{T},Y_{T}) = \frac{1}{n}\sqrt{\sum_{k=1}^{n}(NI_{X_{T}}(\lambda_{k})-NI_{Y_{T}}(\lambda_{k}))^{2}}$,

（3）用取对数之后的特征表示：

$d_{LNP}(X_{T},Y_{T}) = \frac{1}{n}\sqrt{\sum_{k=1}^{n}(\ln NI_{X_{T}}(\lambda_{k})-\ln NI_{Y_{T}}(\lambda_{k}))^{2}}$.

### 基于模型的相似度计算

#### Piccolo 距离

$ARMA(p,q)$ 模型有自己的 AR 表示，因此可以得到相应的一组参数 $(\pi_{1},\pi_{2},\cdots)$，所以，对于每一条时间序列，都可以用一组最优的参数去逼近。如果

$\hat{\prod}_{X_{T}}=(\hat{\pi}_{1,X_{T}},\cdots,\hat{\pi}_{k_{1},X_{T}})^{T},$

$\hat{\prod}_{X_{T}}=(\hat{\pi}_{1,X_{T}},\cdots,\hat{\pi}_{k_{1},X_{T}})^{T}$

$d_{PIC}(X_{T},Y_{T}) =\sqrt{\sum_{j=1}^{k}(\hat{\pi}_{j,X_{T}}'-\hat{\pi}_{j,Y_{T}}')^{2}}$,

#### Maharaj 距离

$d_{MAH}(X_{T},Y_{T}) =\sqrt{T}(\hat{\prod}'_{X_{T}}-\hat{\prod}'_{Y_{T}})^{T}\hat{V}^{-1}(\hat{\prod}'_{X_{T}}-\hat{\prod}'_{Y_{T}}).$

#### 基于 Cepstral 的距离

$\psi_{1}=\phi_{1}$

$\psi_{h}=\phi_{h}+\sum_{m=1}^{h-1}(\phi_{m}-\psi_{h-m})$$1

$\psi_{h}=\sum_{m=1}^{p}(1-\frac{m}{h})\phi_{m}\psi_{h-m}$$p

$d_{LPC, Cep}(X_{T},Y_{T}) =\sqrt{\sum_{i=1}^{T}(\psi_{i,X_{T}}-\psi_{i,Y_{T}})^{2}}$.

# 时间序列的表示与信息提取

### 时间序列的统计特征

$\mu=\frac{1}{T}\sum_{i=1}^{T}x_{i},$

$\sigma^{2} = \sum_{i=1}^{T}\frac{1}{T}(x_{i}-\mu)^{2},$

$\text{skewness}(X) = E[(\frac{X-\mu}{\sigma})^{3}]=\frac{1}{T}\sum_{i=1}^{T}\frac{(x_{i}-\mu)^{3}}{\sigma^{3}},$

$\text{kurtosis}(X) = E[(\frac{X-\mu}{\sigma})^{4}]=\frac{1}{T}\sum_{i=1}^{T}\frac{(x_{i}-\mu)^{4}}{\sigma^{4}} .$

### 时间序列的熵特征

$\text{entropy}(X) = -\sum_{i=1}^{\infty}P\{x=x_{i}\}\ln(P\{x=x_{i}\}).$

#### Binned Entropy

$\text{binned entropy}(X) = -\sum_{k=0}^{\min(maxbin, len(X))} p_{k}\ln(p_{k})\cdot 1_{(p_{k}>0)},$

#### Approximate Entropy

Step 1. 固定两个参数，正整数 $m$ 和正数 $r$，正整数 $m$ 是为了把时间序列进行一个片段的提取，正数 $r$ 是表示时间序列距离的某个参数。i.e. 需要构造新的 $m$ 维向量如下：

$X_{1}(m) = (x_{1},\cdots, x_{m})\in\mathbb{R}^{m},$

$X_{i}(m) = (x_{i},\cdots, x_{m+i-1})\in\mathbb{R}^{m},$

$X_{N-m+1}(m) = (x_{N-m+1},\cdots, x_{N})\in\mathbb{R}^{m}.$

Step 2. 通过新的向量 $X_{1}(m),\cdots, X_{N-m+1}(m)$，可以计算出哪些向量与 $X_{i}$ 较为相似。i.e.

$C_{i}^{m}(r) = (\text{number of }X_{j}(m)\text{ such that } d(X_{i}(m), X_{j}(m))\leq r)/(N-m+1),$

Step 3. 考虑函数

$\Phi^{m}(r) = (N-m+1)^{-1}\cdot \sum_{i=1}^{N-m+1} \ln(C_{i}^{m}(r)),$

Step 4. Approximate Entropy 可以定义为：

$\text{ApEn}(m,r) = \Phi^{m}(r)-\Phi^{m+1}(r).$

Remark.

1. 正整数 $m$ 一般可以取值为 2 或者 3，$r>0$ 会基于具体的时间序列具体调整；
2. 如果某条时间序列具有很多重复的片段（repetitive pattern）或者自相似性（self-similarity pattern），那么它的 Approximate Entropy 就会相对小；反之，如果某条时间序列几乎是随机出现的，那么它的 Approximate Entropy 就会相对较大。

#### Sample Entropy

$A = \#\{\text{vector pairs having } d(X_{i}(m+1),X_{j}(m+1))

$B = \#\{ \text{vector pairs having } d(X_{i}(m), X_{j}(m))

$\text{SampEn} = -\ln(A/B) \geq 0.$

Remark.

1. Sample Entropy 总是非负数；
2. Sample Entropy 越小表示该时间序列具有越强的自相似性（self similarity）。
3. 通常来说，在 Sample Entropy 的参数选择中，可以选择 $m = 2, r = 0.2 \cdot std$.

### 时间序列的分段特征

1. 分段线性逼近（Piecewise Linear Approximation）
2. 分段聚合逼近（Piecewise Aggregate Approximation）
3. 分段常数逼近（Piecewise Constant Approximation）

#### 分段聚合逼近（Piecewise Aggregate Approximation）— 类似 Riemann 积分

$\overline{x}_{i} = \frac{w}{N} \cdot \sum_{j=\frac{N}{w}(i-1)+1}^{\frac{N}{w}i} x_{j}$.

#### 符号逼近（Symbolic Approximation）— 类似 Riemann 积分

SAX 方法的流程如下：

Step 1. 正规化（normalization）：也就是该时间序列被映射到均值为零，方差为一的区间内。

Step 2. 分段表示（PAA）：$\{x_{1},\cdots, x_{N}\} \Rightarrow \{\overline{x}_{1},\cdots,\overline{x}_{w}\}$

Step 3. 符号表示（SAX）：如果 $\overline{x}_{i}，那么 $\hat{X}_{i}=l_{1}$；如果 $z_{(j-1)/\alpha}\leq \overline{x}_{i}，那么 $\hat{X}_{i} = l_{j}$；如果 $\overline{x}_{i}\geq z_{(\alpha-1)/\alpha}$，那么 $\hat{X}_{i} = l_{\alpha}$

# 时序数据与事件的关联分析

### 关联关系的挖掘分成三个部分：

（1）是否存在关联性（Existence of Dependency）：在事件（E）与时间序列（S）之间是否存在关联关系。

（2）关联关系的因果关系（Temporal Order of Dependency）：是事件（E）导致了时间序列（S）的变化还是时间序列（S）导致了事件（E）的发生。

（3）关联关系的单调性影响（Monotonic Effect of Dependency）：用于判断时间序列（S）是发生了突增或者是突降。

### 基本概念：

$e_{i}$来表示某个事件，$\ell_{k}^{rear}(S,e_{i})$表示序列S在事件$e_{i}$之后的长度为k的子序列，$\ell_{k}^{front}(S,e_{i})$表示序列S在事件$e_{i}$之前的长度为k的子序列。如果事件E与时间序列S之间存在关联关系，那么

$\Gamma^{front}=\{\ell_{k}^{front}(S,e_{i}), i=1,\cdots,n\}$

$\Gamma^{rear}=\{\ell_{k}^{rear}(S,e_{i}),i=1,\cdots,n\}$应该是不一样的。

### 方法论：

$\Gamma^{front}$来做例子，$\Gamma^{front}=\{\ell_{k}^{front}(S,e_{i}), i=1,\cdots,n\}$$\Theta =\{\theta_{1},\cdots,\theta_{\tilde{n}}\}$ 是随机选择的，$Z=\Gamma \cup \Theta$，可以标记为$Z_{1},\cdots,Z_{p}$，其中$p=n$+$\tilde{n}$$Z_{i}=\ell_{k}^{front}(S,e_{i})$ when $1\leq i\leq n$$Z_{i}=\theta_{i-n}$ when $n$+$1\leq i\leq p$。可以使用记号$A=A_{1}\cup A_{2}$，其中$A_{1}=\Gamma^{front}$$A_{2}=\Theta=\{\theta_{1},\cdots,\theta_{\tilde{n}}\}$是随机选择的。

$I_{r}(x,A_{1},A_{2})=1$ when $x\in A_{i} \&\& NN_{r}(x,A)\in A_{i}$,

$I_{r}(x,A_{1},A_{2})=0$ when otherwise.

$T_{r,p}=\frac{1}{pr}\sum_{i=1}^{p}\sum_{j=1}^{r}I_{j}(x_{i},A_{1},A_{2})$,

$\lambda_{1}=n/p=n/(n$+$\tilde{n})$, $\lambda_{2}=\tilde{n}/(n$+$\tilde{n})$

$\alpha = 1.96$ for $P=0.025$

$\alpha = 2.58$ for $P=0.001$

$t_{score}=\frac{\mu_{\Gamma^{front}} - \mu_{\Gamma^{rear}}}{\sqrt{\frac{\sigma_{\Gamma^{front}}^{2}+\sigma_{\Gamma^{rear}}^{2}}{n}}}$.

$\alpha = 1.96$ for $P=0.025$

$\alpha = 2.58$ for $P=0.001$

### 算法综述：

7-13行是 $E\rightarrow S$ 的情形，因为$\Gamma^{rear}$ 异常，同时 $\Gamma^{front}$ 正常，说明事件导致了时间序列的变化。7-13行是为了计算 $t_{score}$ 的范围，判断是显著的提升还是下降。

14-20行是 $S\rightarrow E$ 的情形，因为$\Gamma^{front}$ 异常，就导致了事件的发生。14-20行是为了计算 $t_{score}$ 的范围，判断是显著的提升还是下降。

（1）Pearson Correlation

（2）J-Measure Correlation

# Opprentice: Towards Practical and Automatic Anomaly Detection Through Machine Learning

### 系统遇到的挑战：

Definition Challenges: it is difficult to precisely define anomalies in reality.（在现实环境下很难精确的给出异常的定义）

Detector Challenges: In order to provide a reasonable detection accuracy, selecting the most suitable detector requires both the algorithm expertise and the domain knowledge about the given service KPI (Key Performance Indicators). To address the definition challenge and the detector challenge, we advocate for using supervised machine learning techniques. （使用有监督学习的方法来解决这个问题）

### 该系统的优势：

(i) Opprentice is the first detection framework to apply machine learning to acquiring realistic anomaly definitions and automatically combining and tuning diverse detectors to satisfy operators’ accuracy preference.

(ii) Opprentice addresses a few challenges in applying machine learning to such a problem: labeling overhead, infrequent anomalies, class imbalance, and irrelevant and redundant features.

(iii) Opprentice can automatically satisfy or approximate a reasonable accuracy preference (recall>=0.66 & precision>=0.66). （准确率和覆盖率的效果）

### 2. 背景描述：

KPIs and KPI Anomalies:

KPIs: The KPI data are the time series data with the format of (time stamp, value). In this paper, Opprentice pays attention to three kinds of KPIs: the search page view (PV), which is the number of successfully served queries; The number of slow responses of search data centers (#SR); The 80th percentile of search response time (SRT).

Anomalies: KPI time series data can also present several unexpected patterns (e.g. jitters, slow ramp ups, sudden spikes and dips) in different severity levels, such as a sudden drop by 20% or 50%.

### 问题和目标：

1-FDR（false discovery rate）：# of false anomalous points detected / # of anomalous points detected = 1 – precision

The quantitative goal of opprentice is precision>=0.66 and recall>=0.66.

The qualitative goal of opprentice is automatic enough so that the operators would not be involved in selecting and combining suitable detectors, or tuning them.

### 3. Opprentice Overview: （Opprentice系统的概况）

(i) Opprentice approaches the above problem through supervised machine learning.

(ii) Features of the data are the results of the detectors.（Basic Detectors 来计算出特征）

(iii) The labels of the data are from operators’ experience.（人工打标签）

(iv) Addressing Challenges in Machine Learning: （机器学习遇到的挑战）

(1) Label Overhead: Opprentice has a dedicated labeling tool with a simple and convenient interaction interface. （标签的获取）

(2) Incomplete Anomaly Cases:（异常情况的不完全信息）

(3) Class Imbalance Problem: （正负样本比例不均衡）

(4) Irrelevant and Redundant Features:（无关和多余的特征）

### 4. Opprentice’s Design:

Architecture: Operators label the data and numerous detectors functions are feature extractors for the data.

Label Tool:

Detectors:

(i) Detectors As Feature Extractors: （Detector用来提取特征）

Here for each parameter detector, we sample their parameters so that we can obtain several fixed detectors, and a detector with specific sampled parameters a (detector) configuration. Thus a configuration acts as a feature extractor:

data point + configuration (detector + sample parameters) -> feature,

(ii) Choosing Detectors: (Detector的选择，目前有14种较为常见的）

Opprentice can find suitable ones from broadly selected detectors, and achieve a relatively high accuracy. Here, we implement 14 widely-used detectors in Opprentice.

Opprentice has 14 widely-used detectors:

Diff“: it simply measures anomaly severity using the differences between the current point and the point of last slot, the point of last day, and the point of last week.

MA of diff“: it measures severity using the moving average of the difference between current point and the point of last slot.

The other 12 detectors come from previous literature. Among these detectors, there are two variants of detectors using MAD (Median Absolute Deviation) around the median, instead of the standard deviation around the mean, to measure anomaly severity.

(iii) Sampling Parameters: （Detector的参数选择方法，一种是扫描参数空间，另外一种是选择最佳的参数）

Two methods to sample the parameters of detectors.

(1) The first one is to sweep the parameter space. For example, in EWMA, we can choose $\alpha \in \{0.1,0.3,0.5,0.7,0.9\}$ to obtain 5 typical features from EWMA; Holt-Winters has three [0,1] valued parameters $\alpha,\beta,\gamma$. To choose $\alpha,\beta,\gamma \in \{0.2,0.4,0.6,0.8\}$, we have $4^3$ features; In ARIMA, we can estimate their “best” parameters from the data, and generate only one set of parameters, or one configuration for each detector.

Supervised Machine Learning Models:

Decision Trees, logistic regression, linear support vector machines (SVMs), and naive Bayes. 下面是决策树（Decision Tree）的一个简单例子。

Random Forest is an ensemble classifier using many decision trees. It main principle is that a group of weak learners (e.g. individual decision trees) can together form a strong learner. To grow different trees, a random forest adds some elements or randomness. First, each tree is trained on subsets sampled from the original training set. Second, instead of evaluating all the features at each level, the trees only consider a random subset of the features each time. The random forest combines those trees by majority vote. The above properties of randomness and ensemble make random forest more robust to noises and perform better when faced with irrelevant and redundant features than decisions trees.

Configuring cThlds: （阈值的计算和预估）

(i) methods to select proper cThlds: offline part

We need to figure cThlds rather than using the default one (e.g. 0.5) for two reasons.

(1) First, when faced with imbalanced data (anomalous data points are much less frequent than normal ones in data sets), machine learning algorithems typically fail to identify the anomalies (low recall) if using the default cThlds (e.g. 0.5).

(2) Second, operators have their own preference regarding the precision and recall of anomaly detection.

The metric to evaluate the precision and recall are:

(1) F-Score: F-Score = 2*precision*recall/(precision+recall).

(2) SD(1,1): it selects the point with the shortest Euclidean distance to the upper right corner where the precision and the recall are both perfect.

(3) PC-Score: （本文中采用这种评估指标来选择合适的阈值）

If r>=R and p>=P, then PC-Score(r,p)=2*r*p/(r+p) + 1; else PC-Score(r,p)=2*r*p/(r+p). Here, R and P are from the operators’ preference “recall>=R and precision>=P”. Since the F-Score is no more than 1, then we can choose the cThld corresponding to the point with the largest PC-Score.

(ii) EWMA Based cThld Prediction: （基于EWMA方法的阈值预估算法）

In online detection, we need to predict cThlds for detecting future data.

Use EWMA to predict the cThld of the i-th week ( or the i-th test set) based on the historical best cThlds. Specially, EWMA works as follows:

If $i=1$, then $cThld_{i}^{p}=cThld_{1}^{p}=$ 5-fold prediction

Else $i>1$, then $cThld_{i}^{p}=\alpha\cdot cThld_{i-1}^{b}$+$(1-\alpha)\cdot cThld_{i-1}^{p}$, where $cThld_{i-1}^{b}$ is the best cThld of the (i-1)-th week. $cThld_{i}^{p}$ is the predicted cThld of the i-th week, and also the one used for detecting the i-th week data. $\alpha\in [0,1]$ is the smoothing constant.

For the first week, we use 5-fold cross-validation to initialize $cThld_{1}^{p}$. As $\alpha$ increases, EWMA gives the recent best cThlds more influences in the prediction. We use $\alpha=0.8$ in this paper.

### 5. Evaluation（系统评估）

Opprentice has 14 detectors with about 9500 lines of Python, R and C++ code. The machine learning block is based on the scikit-learn library.

Random Forest is better than decision trees, logistic regression, linear support vector machines (SVMs), and naive Bayes.

# Focus: Shedding Light on the High Search Response Time in the Wild

### 问题描述：

To help search operators dubug HSRT (high search response time)，Focus is a search log analysis framework to answer the three questions:

(1) What is the HSRT condition?

(2) Which HSRT condition types are prevalent across days?

(3) How does each attribute affect SRT in those prevalent HSRT condition types?

### 解决方案：

Focus has one component for each of the above questions:

(1) A decision tree based classifier to identify HSRT conditions in search logs of each day;

(2) A clustering based condition type miner to combine similar HSRT conditions into one type, and find the prevalent condition types across days; following Occam’s razor principle.

(3) An attribute effect estimator to analyze the effect of each individual attribute of SRT within a prevalent condition type.

### 基础知识准备：

(A) Search Logs:

For each measured query, its search log records two types of data: SRT and SRT components, Query Attributes.

(1) SRT and SRT components:（特征层）

$t_{1}$ is when a query is submitted; $t_{2}$ is when the result HTML file has been downloaded; $t_{3}$ is when a brower finishes parsing the HTML; $t_{4}$ is when the page is completely rendered. SRT is measured by $t_{4}-t_{1}$, the user-received search response time.

$T_{server}$ is the server response time of the HTML file, which is recorded by servers; $T_{net}=t_{2}-t_{1}-T_{server}$ is the network transmission time of the HTML file; $T_{brower}=t_{3}-t_{2}$ is the browser parsing time of the HTML; $T_{other}=t_{4}-t_{3}$ is the remaining time spent before the page is rendered, e.g. download time of images from image servers.

(2) Query Attributes:（特征层）

The search logs record the following attributes for each measured query:

(i) Browser Engine: Webkit(e.g. Chrome, Safari and 360 Secure Browser), Gecko, Trident LEGC, Trident 4.0, Trident 5.0, and others.

(ii) ISP: China Telecom, China Unicom, China Mobile, China Netcom, CHina Tietong, others.

(iii) Localtion: Based on the client IP, convert IP to its geographic location. In total, there are 32 provinces.

(iv) #Image: the number of embedded images in the result page.

(vii) Background page views: On the service side, the search engine S also post-analyzes the logs and generates the background page views. The background PVs (page views) for a query q is measured by the number of queries served within 30 seconds before and after q is served.It reflects the average search request load where q is served. Due to confidentiality constraints, we normalize specific background PVs (page Views) by the maximum value.（事后分析，统计出一些必要的特征，输入 Focus 系统的机器学习模型中）

(B) HSRT and HSRT Conditions:（样本层）

Usually, we can use cumulative distribution fraction (CDF) of SRT in the search logs to determine the high search response time condition (HSRT condition). In this paper,  we define HSRT as the SRT longer than 1s.

Challenges of Identifying HSRT Conditions: In order to identify HSRT conditions in multi-dimensional search logs.（以下是这个系统的一些难点和挑战点）

(a) Naive Single Dimensional Based Methods: including pair-wise correlation analysis and so on, but is inefficient.

(b) Attributes can be potentially interdependent on each other: that means Naive Bayes Method may not applicable in this situation.

(c) Need to avoid output overlapping conditions: like {#image>30}, {ads=yes}, and {#image>20, ads=yes}.  （随着时间的推移，每天使用模型可能会推出类似或者重复的规则）。

### 关键思想和系统概况

Condition is a combination of attributes and specific values in search logs.

HSRT Condition is a condition that covers at least 1%\$ of total queries, and has the fraction of HSRT large than the global level:

(# of HSRT queries in a HSRT condition / #of queries in a HSRT condition) > (# of HSRT queries / # of queries). This is in order to assign to labels and we can change this definition in practice. （这只是用来打标签的定义，用于判断什么是HSRT，在实际的应用中，我们可以根据具体的场景采用不同的定义，例如返回码等指标）。

### ‘Focus’ System Overview:

Input: search logs（日志）

(i) Use a decision tree based classifier to identify HSRT conditions  in search logs every day; （每天可以使用决策树模型从日志中提取HSRT条件）。

(ii) Use a clustering based condition type miner to identify condition types of similar HSRT conditions, and fine prevalent condition types across days; （用于把类似的条件融合在一起）。

(iii) Use an attribute effect estimator to analyze how an attribute affects SRT and SRT components in each prevalent condition type. （用于判断哪些属性或者特征对这个标签影响更加深远）。

Output: prevalent condition types and their attributes effects on SRT.（第二步输出的条件以及第三步属性的重要性）。

Part (i): Decision Tree Based Classifier including ID3, C4.5, CART. It contains five important parts: (1) expressing attribute splits; (2) evaluate splits; (3) stopping tree growing; (4) assigning Labels: assign HSRT labels to the left nodes whose fraction of HSRT is larger than the global fraction of HSRT; (5) identify HSRT Branching Attribute Conditions. （这里是 Focus 系统所采用的机器学习算法）。

Part (ii): Condition Type Miner: group HSRT conditions according to (1) the same combination of attributes, (2) the same value from each category attribute, and (3) similar interval for each numeric attribute, using Jaccard Index to measure the similarity between intervals. （条件的融合）。

Part (iii): Attribute Effect Estimator: With each condition type

$C=\{c_{1}\wedge c_{2}\cdots \wedge c_{i} \wedge \cdots c_{n}\}$,

we design a method to understand how each attribute condition $c_{i}$ affects SRT.

For example, what is the HSRT fraction caused by $c_{i}$ in $C$? What SRT components (e.g. $T_{net}$ and $T_{server}$) are affected by $c_{i}$?

Main Idea: flip condition $c_{i}$ to the opposite $\overline{c}_{i}$ to get a variant condition type $C_{i}'=\{c_{1}\wedge c_{2}\cdots \wedge \overline{c}_{i} \wedge \cdots c_{n}\}$. In the past days, we have the number of HSRT events in total, the number of HSRT events in condition $C$ and the number of HSRT events in condition $C_{i}'$. As a result, we believe the historical data based comparison can provide a reasonable estimate of the attribute effects. The comparison between $C$ and $C_{i}'$ in these days is based on the specific HSRT conditions of these days. （用于判断哪些属性更能够引起 HSRT）。

In Table IV, the results are sorted by the variation of the fraction of HSRT in condition types (HSRT% column) caused by flipping an attribute condition.

(i) We highlight the variations greater than zero (getting worse after flipping an attribute condition).

(ii) We focus on that flipping the HSRT branching attribute conditions can yield improvements on HSRT%. For example, the condition #image>x are all ranked at the top. It means we need to reduce the impact of images on SRT and we can get the highest potential improvement of HSRT.

(iii) Table III and Table IV are the output of Focus to the operators for these months.

### Observations by Further Inverstigation

Table IV raises some interesting questions:（通过 Focus 输出的表格 Table IV 可以提出很多其余的问题，也许是人工经验不容易发现的问题）

(1) Why does reducing #images increase $T_{server}$, the time that servers prepare the result HTML (row 1, 2, 3, 4 of Table IV)?

(2) How do ads inflate SRT? Why do the pages with ads need more $T_{net}$ and $T_{brower}$ (row 7)?

(3) Why does Webkit engine perform better, especially greatly decreasing $T_{browser}$ (row 5, 10, 11, 12)?

(4) It is nature that switching ISPs can affect network transmission time $T_{net}$, but why does switching to China Telecom reduce $T_{server}$ by over 20% (row 6, 8, 9)?