PHD的生涯, 互联网的生涯

2019年的博客

3 Comments

个人使用时间最长的博客平台就是 WordPress，从 2011 年开始算起，大约有 8 年时间了。

当年笔者还在南京大学本科读数学专业的时候，就见到大神 Terence Tao 建立了个人的网站：What’s new（terrytao.wordpress.com），并且后续数十年如一日的更新个人 Blog。在他的 Blog 上面，除了常见的数学基础知识和课程安排之外，更多的是当前数学界的一些新发展和新方向。由于数学系的博文或者论文除了文字之外，更多的是各种各样的数学公式，而这些数学公式使用公式编辑器一类的东西来处理是极其繁琐的，只能够使用 LaTex 等工具来写。恰好的是，WordPress 除了能够支持日常的文本编辑之外，还能够使用 LaTex 来对数学公式进行撰写，也就是说用户只需要在编辑框内写 LaTex，就能够编译成数学公式。因此，WordPress 对数学公式的支持是相对友好的，这对数学系爱好写数学博客的学生，工作人士提供了非常便利的条件。

除了整理数学公式比较容易之外，WordPress 上面还可以相对方便的选择各种各样的主题，这样的话刚注册的新用户也可以较为容易的上路，不用一开始就陷入编辑网站等一些繁琐的事情上面。同时，WordPress 上也有各种各样的小工具，包括日历（可以查看发表 Blog 的时间），文章的目录，比较受欢迎的文章或页面，文章的类别等内容。用户可以根据自己的爱好自行选择栏目，从而可以轻易地搭建出一个个人网站。另外，WordPress 可以统计出 Blog 的点击数量，包括每天，每周，每月，每年的具体点击数量。

第一次使用 WordPress（zr9558.com）正好是攻读博士学位一年，除了日常的科研工作之外，也打算写点东西来记录一下自己的成长。后续读博士期间也逐渐的写了一些数学方面的文章，不过后续回想起来其实应该在读本科的时候就开始写 Blog。如果这样的话，当年所学过的各种数学知识，整理过的各种资料都会更加清晰一些，也更加容易保存一些。毕竟写 Blog 的一个重要目的是给自己回顾用的，看看这段时间自己的积累是什么，自己学到了什么知识，相比去年成长了多少。其实，有的时候偶尔去浏览一下 Terence Tao 的博客，虽然也看不懂，但是可以明显地感受到“大神”的努力程度。牛人尚且如此努力，我等凡人有什么理由不努力呢。

当年在 NUS 读博士的时候也承担了一些教学的工作，并且也会给学生们讲解习题课。但是每年讲课的内容其实都差不多，于是在讲解习题课的时候，尤其是大学数学（MA1505），曾经做过一份 Tutorials，也就是把习题课的内容写在 WordPress 上面。不仅可以在课堂上方便地给学生们讲解，也可以让学生在课后方便的复习。时过境迁，当年听课的学生（还是大一）早已本科毕业，后面可能也会有学生看到这份资料，希望对大家有所帮助吧。链接在这里：MA 1505 Mathematics I，Tutorials。

后面进入公司之后，有的时候工作繁忙，整理 Blog 的时间就会减少许多。但是一份工作除了能够给人带来必要的薪资之外，更重要的是给自己不停地积攒经验。无论是现在还是将来都可以让自己在职场中更加值钱。因此，除了日常的搬砖工作之外，也要时刻注意自己的成长和经验的积累。而搬砖的经验会随着项目的结束和时间的迁移在记忆中逐渐淡忘而去，于是，适当的记录就成为了必要的工作。俗话说得好，“好记性不如烂笔头”。因此，隔一段时间（通常来说可以设定为一到两个月的时长）就整理一下经验就显得尤为重要，也是提升个人技术和经验的方法之一。但是在整理 Blog 的时候，一定要注意 Blog 的质量，只有不断地提炼自己 Blog 里面的内容才能够保证文章的质量。

不过现在的 Blog 远没有当年受欢迎，在各种各样 APP 横行的时代，已经比较少有人主动去看别人的 Blog 了。不过在使用搜索引擎来搜索某些关键字的时候，有的时候还是能够看到一些高质量的 Blog。在学习这些 Blog 的同时，其实也可以互相比较一下，取长补短才能够使自己的博客越做越好。其实，无论有没有人读自己所写的内容，都要坚持写下去，因为：

即使最后没有人为你鼓掌，也要优雅地谢幕，感谢自己的认真付出。

Advertisements
数据挖掘与机器学习

计算机视觉中的注意力机制

1 Comment

引言

在机器翻译（Machine Translation）或者自然语言处理（Natural Language Processing）领域，以前都是使用数理统计的方法来进行分析和处理。近些年来，随着 AlphaGo 的兴起，除了在游戏AI领域，深度学习在计算机视觉领域，机器翻译和自然语言处理领域也有着巨大的用武之地。在 2016 年，随着深度学习的进一步发展，seq2seq 的训练模式和翻译模式已经开始进入人们的视野。除此之外，在端到端的训练方法中，除了需要海量的业务数据之外，在网络结构中加入一些重要的模块也是非常必要的。在此情形下，基于循环神经网咯（Recurrent Neural Network）的注意力机制（Attention Mechanism）进入了人们的视野。除了之前提到的机器翻译和自然语言处理领域之外，计算机视觉中的注意力机制也是十分有趣的，本文将会简要介绍一下计算机视觉领域中的注意力方法。在此事先声明一下，笔者并不是从事这几个领域的，可能在撰写文章的过程中会有些理解不到位的地方，请各位读者指出其中的不足。

LSTM_1

注意力机制

顾名思义，注意力机制是本质上是为了模仿人类观察物品的方式。通常来说，人们在看一张图片的时候，除了从整体把握一幅图片之外，也会更加关注图片的某个局部信息，例如局部桌子的位置，商品的种类等等。在翻译领域，每次人们翻译一段话的时候，通常都是从句子入手，但是在阅读整个句子的时候，肯定就需要关注词语本身的信息，以及词语前后关系的信息和上下文的信息。在自然语言处理方向，如果要进行情感分类的话，在某个句子里面，肯定会涉及到表达情感的词语，包括但不限于“高兴”，“沮丧”，“开心”等关键词。而这些句子里面的其他词语，则是上下文的关系，并不是它们没有用，而是它们所起的作用没有那些表达情感的关键词大。

在以上描述下，注意力机制其实包含两个部分：

  1. 注意力机制需要决定整段输入的哪个部分需要更加关注；
  2. 从关键的部分进行特征提取，得到重要的信息。

通常来说，在机器翻译或者自然语言处理领域，人们阅读和理解一句话或者一段话其实是有着一定的先后顺序的，并且按照语言学的语法规则来进行阅读理解。在图片分类领域，人们看一幅图也是按照先整体再局部，或者先局部再整体来看的。再看局部的时候，尤其是手写的手机号，门牌号等信息，都是有先后顺序的。为了模拟人脑的思维方式和理解模式，循环神经网络（RNN）在处理这种具有明显先后顺序的问题上有着独特的优势，因此，Attention 机制通常都会应用在循环神经网络上面。

虽然，按照上面的描述，机器翻译，自然语言处理，计算机视觉领域的注意力机制差不多，但是其实仔细推敲起来，这三者的注意力机制是有明显区别的。

  1. 在机器翻译领域，翻译人员需要把已有的一句话翻译成另外一种语言的一句话。例如把一句话从英文翻译到中文，把中文翻译到法语。在这种情况下，输入语言和输出语言的词语之间的先后顺序其实是相对固定的，是具有一定的语法规则的；
  2. 在视频分类或者情感识别领域，视频的先后顺序是由时间戳和相应的片段组成的，输入的就是一段视频里面的关键片段，也就是一系列具有先后顺序的图片的组合。NLP 中的情感识别问题也是一样的，语言本身就具有先后顺序的特点；
  3. 图像识别，物体检测领域与前面两个有本质的不同。因为物体检测其实是在一幅图里面挖掘出必要的物体结构或者位置信息，在这种情况下，它的输入就是一幅图片，并没有非常明显的先后顺序，而且从人脑的角度来看，由于个体的差异性，很难找到一个通用的观察图片的方法。由于每个人都有着自己观察的先后顺序，因此很难统一成一个整体。

在这种情况下，机器翻译和自然语言处理领域使用基于 RNN 的 Attention 机制就变得相对自然，而计算机视觉领域领域则需要必要的改造才能够使用 Attention 机制。

LSTM_3

基于 RNN 的注意力机制

通常来说，RNN 等深度神经网络可以进行端到端的训练和预测，在机器翻译领域和或者文本识别领域有着独特的优势。对于端到端的 RNN 来说，有一个更简洁的名字叫做 sequence to sequence，简写就是 seq2seq。顾名思义，输入层是一句话，输出层是另外一句话，中间层包括编码和解码两个步骤。

而基于 RNN 的注意力机制指的是，对于 seq2seq 的诸多问题，在输入层和输出层之间，也就是词语（Items）与词语之间，存在着某种隐含的联系。例如：“中国” -> “China”，“Excellent” -> “优秀的”。在这种情况下，每次进行机器翻译的时候，模型需要了解当前更加关注某个词语或者某几个词语，只有这样才能够在整句话中进行必要的提炼。在这些初步的思考下，基于 RNN 的 Attention 机制就是：

  1. 建立一个编码（Encoder）和解码（Decoder）的非线性模型，神经网络的参数足够多，能够存储足够的信息；
  2. 除了关注句子的整体信息之外，每次翻译下一个词语的时候，需要对不同的词语赋予不同的权重，在这种情况下，再解码的时候，就可以同时考虑到整体的信息和局部的信息。

LSTM_4

注意力机制的种类

从初步的调研情况来看，注意力机制有两种方法，一种是基于强化学习（Reinforcement Learning）来做的，另外一种是基于梯度下降（Gradient Decent）来做的。强化学习的机制是通过收益函数（Reward）来激励，让模型更加关注到某个局部的细节。梯度下降法是通过目标函数以及相应的优化函数来做的。无论是 NLP 还是 CV 领域，都可以考虑这些方法来添加注意力机制。

LSTM_5

计算机视觉领域的 Attention 部分论文整理

下面将会简单的介绍几篇近期阅读的计算机视觉领域的关于注意力机制的文章。

Look Closer to See Better：Recurrent Attention Convolutional Neural Network for Fine-grained Image Recognition

在图像识别领域，通常都会遇到给图片中的鸟类进行分类，包括种类的识别，属性的识别等内容。为了区分不同的鸟，除了从整体来对图片把握之外，更加关注的是一个局部的信息，也就是鸟的样子，包括头部，身体，脚，颜色等内容。至于周边信息，例如花花草草之类的，则显得没有那么重要，它们只能作为一些参照物。因为不同的鸟类会停留在树木上，草地上，关注树木和草地的信息对鸟类的识别并不能够起到至关重要的作用。所以，在图像识别领域引入注意力机制就是一个非常关键的技术，让深度学习模型更加关注某个局部的信息。

RA_CNN_1

在这篇文章里面，作者们提出了一个基于 CNN 的注意力机制，叫做 recurrent attention convolutional neural network（RA-CNN），该模型递归地分析局部信息，从局部的信息中提取必要的特征。同时，在 RA-CNN 中的子网络（sub-network）中存在分类结构，也就是说从不同区域的图片里面，都能够得到一个对鸟类种类划分的概率。除此之外，还引入了 attention 机制，让整个网络结构不仅关注整体信息，还关注局部信息，也就是所谓的 Attention Proposal Sub-Network（APN）。这个 APN 结构是从整个图片（full-image）出发，迭代式地生成子区域，并且对这些子区域进行必要的预测，并将子区域所得到的预测结果进行必要的整合，从而得到整张图片的分类预测概率。

RA_CNN_2

RA-CNN 的特点是进行一个端到端的优化，并不需要提前标注 box，区域等信息就能够进行鸟类的识别和图像种类的划分。在数据集上面，该论文不仅在鸟类数据集（CUB Birds）上面进行了实验，也在狗类识别（Stanford Dogs）和车辆识别（Stanford Cars）上进行了实验，并且都取得了不错的效果。

RA_CNN_4

从深度学习的网络结构来看，RA-CNN 的输入时是整幅图片（Full Image），输出的时候就是分类的概率。而提取图片特征的方法通常来说都是使用卷积神经网络（CNN）的结构，然后把 Attention 机制加入到整个网络结构中。从下图来看，一开始，整幅图片从上方输入，然后判断出一个分类概率；然后中间层输出一个坐标值和尺寸大小，其中坐标值表示的是子图的中心点，尺寸大小表示子图的尺寸。在这种基础上，下一幅子图就是从坐标值和尺寸大小得到的图片，第二个网络就是在这种基础上构建的；再迭代持续放大图片，从而不停地聚焦在图片中的某些关键位置。不同尺寸的图片都能够输出不同的分类概率，再将其分类概率进行必要的融合，最终的到对整幅图片的鸟类识别概率。

因此，在整篇论文中，有几个关键点需要注意：

  1. 分类概率的计算，也就是最终的 loss 函数的设计；
  2. 从上一幅图片到下一幅图片的坐标值和尺寸大小。

只要获得了这些指标，就可以把整个 RA-CNN 网络搭建起来。

大体来说，第一步就是给定了一幅输入图片 X， 需要提取它的特征，可以记录为 W_{c}*X，这里的 * 指的是卷积等各种各样的操作。所以得到的概率分布情况其实就是 p(X) = f(W_{c}*X)，f 指的是从 CNN 的特征层到全连接层的函数，外层使用了 Softmax 激活函数来计算鸟类的概率。

第二步就是计算下一个 box 的坐标 (t_{x}, t_{y}) 和尺寸大小 t_{\ell}，其中 t_{x}, t_{y} 分别指的是横纵坐标，正方形的边长其实是 2*t_{\ell}。用数学公式来记录这个流程就是 [t_{x}, t_{y}, t_{\ell}] = g(W_{c}*X)。在坐标值的基础上，我们可以得到以下四个值，分别表示 x, y 两个坐标轴的上下界：

t_{x(t\ell)} = t_{x} - t_{\ell}, t_{x(br)} = t_{x} + t_{\ell},

t_{y(t\ell)} = t_{y} - t_{\ell}, t_{y(br)} = t_{y} + t_{\ell}.

局部注意力和放大策略（Attention Localization and Amplification）指的是：从上面的方法中拿到坐标值和尺寸，然后把图像进行必要的放大。为了提炼局部的信息，其实就需要在整张图片 X 的基础上加上一个面具（Mask）。所谓面具，指的是在原始图片的基础上进行点乘 0 或者 1 的操作，把一些数据丢失掉，把一些数据留下。在图片领域，就是把周边的信息丢掉，把鸟的信息留下。但是，有的时候，如果直接进行 0 或者 1 的硬编码，会显得网络结构不够连续或者光滑，因此就有其他的替代函数。

在激活函数里面，逻辑回归函数（Logistic Regression）是很常见的。其实通过逻辑回归函数，我们可以构造出近似的阶梯函数或者面具函数。

sigmoid_1

对于逻辑回归函数 \sigma(x) = 1/(1+e^{-kx}) 而言，当 k 足够大的时候，\sigma(x) \approx 1 当 x \geq 0；\sigma(x) \approx 0 当 x<0。此时的逻辑回归函数近似于一个阶梯函数。如果假设 x_{0}<x_{1}，那么 \sigma(x-x_{0}) - \sigma(x-x_{1}) 就是光滑一点的阶梯函数，\sigma(x-x_{0}) - \sigma(x-x_{1}) \approx 0 当 x < x_{0} \text{ or } x > x_{1}；\sigma(x-x_{0}) - \sigma(x-x_{1}) \approx 1 当 x_{0}\leq x\leq x_{1}。

因此，基于以上的分析和假设，我们可以构造如下的函数：X^{attr} = X \odot M(t_{x}, t_{y}, t_{\ell}), 其中，X^{attr} 表示图片需要关注的区域，M(\cdot) 函数就是 M(t_{x}, t_{y}, t_{\ell}) = [\sigma(x-t_{x(t\ell)}) - \sigma(x-t_{x(br)})]\cdot[\sigma(y-t_{y(t\ell)}) - \sigma(y-t_{y(br)})], 这里的 \sigma 函数对应了一个足够大的 k 值。

当然，从一张完整的图片到小图片，在实际操作的时候，需要把小图片继续放大，在放大的过程中，可以考虑使用双线性插值算法来扩大。也就是说：

X_{(i,j)}^{amp} = \sum_{\alpha,\beta=0}^{1}|1-\alpha-\{i/\lambda\}|\cdot|1-\beta-\{j/\lambda\}|\cdot X_{(m,n)}^{att},

其中 m = [i/\lambda] + \alpha, n = [j/\lambda] + \beta，\lambda 表示上采样因子，[\cdot], \{\cdot\} 分别表示一个实数的正数部分和小数部分。

在分类（Classification）和排序（Ranking）部分，RA-CNN 也有着自己的方法论。在损失函数（Loss Function）里面有两个重要的部分，第一个部分就是三幅图片的 LOSS 函数相加，也就是所谓的 classification loss，Y^{(s)} 表示预测类别的概率，Y 表示真实的类别。除此之外，另外一个部分就是排序的部分，L_{rank}(p_{t}^{(s)}, p_{t}^{(s+1)}) = \max\{0,p_{t}^{(s)}-p_{t+1}^{(s+1)}+margin\}, 其中 p^{(s)} 表示在第 s 个尺寸下所得到的类别 t 的预测概率，并且最大值函数强制了该深度学习模型在训练中可以保证 p_{t}^{(s+1)} > p_{t}^{(s)} + margin，也就是说，局部预测的概率值应该高于整体的概率值。

L(X) = \sum_{s=1}^{3}\{L_{cls}(Y^{(s)},Y^{*})\} + \sum_{s=1}^{2}\{L_{rank}(p_{t}^{(s)},p_{t}^{(s+1)})\}.

RA_CNN_3

在这种 Attention 机制下，可以使用训练好的 conv5_4 或者 VGG-19 来进行特征的提取。在图像领域，location 的位置是需要通过训练而得到的，因为每张图片的鸟的位置都有所不同。进一步通过数学计算可以得到，t_{\ell} 会随着网络而变得越来越小，也就是一个层次递进的关系，越来越关注到局部信息的提取。简单来看，

\frac{\partial L_{rank}}{\partial t_{x}} \propto D_{top} \odot \frac{\partial M(t_{x},t_{y},t_{\ell})}{\partial t_{x}},

这里的 \odot 表示元素的点乘，D_{top} 表示之前的网络所得到的导数。

当 x\rightarrow t_{x(t\ell)}，\frac{\partial M}{\partial t_{x}}<0;

当 x \rightarrow t_{x(br)}，\frac{\partial M}{\partial t_{x}}>0;

其余情况，\frac{\partial M}{\partial t_{x}}=0.

当 y\rightarrow t_{y(t\ell)}，\frac{\partial M}{\partial t_{y}}<0;

当 y \rightarrow t_{y(br)}，\frac{\partial M}{\partial t_{y}}>0;

其余情况，\frac{\partial M}{\partial t_{y}}=0.

当 x \rightarrow t_{x(t\ell)}\text{ or } x \rightarrow t_{x(br)}\text{ or } y \rightarrow t_{y(t\ell)}\text{ or } y \rightarrow t_{y(br)}, \frac{\partial M}{\partial t_{\ell}}>0;

其余情况，\frac{\partial M}{\partial t_{\ell}}<0.

因此，t_{\ell} 在迭代的过程中会越来越小，也就是说关注的区域会越来越集中。

RA-CNN 的实验效果如下：

RA_CNN_5.png

 

Multiple Granularity Descriptors for Fine-grained Categorization

这篇文中同样做了鸟类的分类工作，与 RA-CNN 不同之处在于它使用了层次的结构，因为鸟类的区分是按照一定的层次关系来进行的，粗糙来看，有科 -> 属 -> 种三个层次结构。

MC_CNN_1

因此，在设计网络结构的过程中，需要有并行的网络结构，分别对应科，属，种三个层次。从前往后的顺序是检测网络（Detection Network），区域发现（Region Discovery），描述网络（Description Network）。并行的结构是 Family-grained CNN + Family-grained Descriptor，Genus-grained CNN + Genus-grained Descriptor，Species-grained CNN + Species-grained Descriptor。而在区域发现的地方，作者使用了 energy 的思想，让神经网络分别聚焦在图片中的不同部分，最终的到鸟类的预测结果。

MC_CNN_2MC_CNN_3

Recurrent Models of Visual Attention

在计算机视觉中引入注意力机制，DeepMind 的这篇文章 recurrent models of visual attention 发表于 2014 年。在这篇文章中，作者使用了基于强化学习方法的注意力机制，并且使用收益函数来进行模型的训练。从网络结构来看，不仅从整体来观察图片，也从局部来提取必要的信息。

DeepMind_1

DeepMind_2DeepMind_3

整体来看，其网络结构是 RNN，上一个阶段得到的信息和坐标会被传递到下一个阶段。这个网络只在最后一步进行分类的概率判断，这是与 RA-CNN 不同之处。这是为了模拟人类看物品的方式，人类并非会一直把注意力放在整张图片上，而是按照某种潜在的顺序对图像进行扫描。Recurrent Models of Visual Attention 本质上是把图片按照某种时间序列的形式进行输入，一次处理原始图片的一部分信息，并且在处理信息的过程中，需要根据过去的信息和任务选择下一个合适的位置进行处理。这样就可以不需要进行事先的位置标记和物品定位了。

DeepMind_4

正如上图所示，enc 指的是对图片进行编码，r_{i}^{(1)} 表示解码的过程，x_{i} 表示图片的一个子区域。而 y_{s} 表示对图片的预测概率或者预测标签。

 

Multiple Object Recognition with Visual Attention

这篇文章同样是 DeepMind 的论文，与 Recurrent Models of Visual Attention 不同之处在于，它是一个两层的 RNN 结构，并且在最上层把原始图片进行输入。其中 enc 是编码网络，r^{(1)}_{i} 是解码网络，r_{i}^{(2)} 是注意力网络，输出概率在解码网络的最后一个单元输出。

deep_recurrent_attention_model_1

在门牌识别里面，该网络是按照从左到右的顺序来进行图片扫描的，这与人类识别物品的方式极其相似。除了门牌识别之外，该论文也对手写字体进行了识别，同样取得了不错的效果。

deep_recurrent_attention_model_3deep_recurrent_attention_model_2

实验效果如下：

deep_recurrent_attention_model_4.png

总结

本篇 Blog 初步介绍了计算机视觉中的 Attention 机制，除了这些方法之外，应该还有一些更巧妙的方法，希望各位读者多多指教。

参考文献

  1. Look Closer to See Better：Recurrent Attention Convolutional Neural Network for Fine-grained Image Recognition，CVPR，2017.
  2. Recurrent Models of Visual Attention，NIPS，2014
  3. GitHub 代码：Recurrent-Attention-CNN，https://github.com/Jianlong-Fu/Recurrent-Attention-CNN
  4. Multiple Granularity Descriptors for Fine-grained Categorization，ICCV，2015
  5. Multiple Object Recognition with Visual Attention，ICRL，2015
  6. Understanding LSTM Networks，Colah’s Blog，2015，http://colah.github.io/posts/2015-08-Understanding-LSTMs/
  7. Survey on the attention based RNN model and its applications in computer vision，2016

 

时间序列

时间序列的聚类

Leave a comment

在机器学习领域，聚类问题一直是一个非常常见的问题。无论是在传统的机器学习（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
时间序列

时间序列的单调性

Leave a comment

在时间序列的众多研究方向上，除了时间序列异常检测，时间序列的相似性，时间序列的趋势预测之外，无论是在量化交易领域还是其余领域，时间序列的单调性都是一个重要课题。本文将会对时间序列的单调性作简单的介绍。

连续函数的单调性

导数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 1，M_{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}]，为了计算某个时间戳 n 下 x_{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
时间序列

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

Leave a comment

引言

在时间序列异常检测中，特征工程往往是非常繁琐而复杂的，怎样才能够减少时间序列的特征工程工作量一直是一个关键问题。在本文中，作者们提出了一个新的思路，使用深度学习的办法来进行端到端的训练，从而减少时间序列的特征工程。

提到深度学习，大家都能够想到卷积神经网络（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
Computer Science

非计算机专业学生如何转行 AI

Leave a comment

个人背景

笔者本科和博士期间都在数学系攻读基础数学，也就是那种跟工业界基本上挨不上边的东西。后来博士毕业之后进入互联网公司搬砖，于是就开始做机器学习方向。之前也写过关于转行的文章，不过近期看到知乎上有类似的问题，于是整理一下之前所写分享给大家。

转专业的困难

虽然现在很多人都会说数学学完之后转计算机有优势，学了数学之后学金融如鱼得水。但是这些人很可能既没学过数学，也没学过计算机和金融，只是看了网络或者报纸上的宣传就开始四处说这些观点。其实，作为一个数学系的学生，如果要想转专业的话，其实是需要付出很多时间和精力的。因为数学系所上的课程和计算机所上的课程是不一样的。通常来说：数学系和计算机系的不完全课表大致如下：
IMG_4671.PNG

数学系的课程：

数学分析，高等代数，解析几何，C++，离散数学，常微分方程，偏微分方程，抽象代数，复变函数，实变函数，泛函分析，数值计算，偏微分方程数值解，拓扑学，微分几何，概率论与数理统计，随机过程等。

计算机系的课程：

微积分，线性代数，离散数学，数据结构与算法，数字电路，计算机组成原理，操作系统，编译原理，计算机网络，数据库原理，软件工程，汇编语言等。
从这两个表格的对比情况来看，如果要想从数学系转行到计算机系，那么基本上要把计算机的一些基础知识课程都大致过一遍才行，否则企业为什么不直接招聘一个计算机系的，而需要一个跨专业的人呢？在这种情况下，对数学系的人其实提出了很高的挑战，因为在数学系繁重的课程下，想要同时兼顾数学系和计算机系两个专业的课程是比较困难的，需要同学耗费巨大的时间和精力才能够做好。

人工智能行业所带来的机遇和挑战

就这几年的人工智能发展情况和笔者的个人经验而言，人工智能可以大致分成以下几个方向：
  1. 计算机视觉方向；
  2. 自然语言处理方向；
  3. 语音识别方向；
  4. 机器学习方向。

IMG_4673

作为一个转专业的学生，如果要在一些比较成熟的行业里面去和科班的人去竞争，那么自然就会出现劣势，因为企业是非常看个人产出的。既然能够招聘到一个有经验的人，其实没有必要培养一个无经验者。这种时候，转专业的同学一定要找好自己的定位，也就是传说中的“打法”，需要凸显自己的优势，然后尽量避开劣势。举个例子，假设让一个数学系的人去做美术设计，在大多数情况下就是把自身的缺点暴露给别人，当然在数学系也有美术不错的同学。但是在大多数情况下，数学系的美术功底比艺术院校的美术功底其实是差很远的。在人工智能领域也是这样的，在一些成熟的领域，其实计算机系就能给工业界源源不断地提供人才，对于转行的人来说其实是不算特别友好的。计算机视觉方向（Computer Vision）无论是在学校还是在公司，都有着大量的从业者，并且 ImageNet 项目可以提供上千万的标注图片供大家使用。既然 ImageNet 是开源的数据集，那么无论是学校的教授还是学生，不管是大型互联网公司还是初创企业，都可以轻易地获取到这些数据集，不仅可以进行 CV 算法的研究工作，还可以进行相关的工程实践。由于计算机视觉方向的历史悠久，不管是计算机系，工程系，甚至数学系，都有着大量的老师和相应的学生从事该方向的研究工作，因此，学校或者研究所能够对工业界输出的计算机视觉人才数量也是可观的。其他两个，自然语言处理和语音识别的老师其实也是有的，可能相对于图像来说是少了一些。

IMG_4674
如果针对机器学习领域的话，就公司或者学术界的一些情况来看，其实机器学习领域的应用范围十分广泛。最经典的当然属于广告推荐和个性化推荐这一块，无论是今日头条，抖音，还是各个 APP，其实都包含了推荐系统，无论这个推荐系统是通过规则的形式做出来的，通过逻辑回归的方法做出来的，还是通过深度学习做出来的，都是可以在点击率和利润等方向上获得收益的。除了推荐系统之外，游戏 AI 也是一个不错的研究方向，几年前强化学习这个方向也是不温不火，但是在 AlphaGo 崛起之后，深度学习和强化学习就已经开始进入了大多数人的视野。随着围棋被攻克之后，德州扑克AI，或者其他的游戏 AI 也被很多学者和大型游戏公司所关注。DeepMind 也在 2017 年开放了星际争霸的研究平台，今年无论是在 Dota2 还是星际争霸上，游戏 AI 相比之前都有了巨大的突破。因此，如何在诸多业务线中，选择一个适合自己的研究方向，才是比较关键的问题。是选择一个成熟的领域努力奋斗，还是选择一个新兴领域开疆拓土，都是需要自己去考虑的。

机器学习如何入门

上一部分介绍了人工智能方向的一些情况，下面可能各位同学比较关心机器学习领域该如何入门。由于笔者是做机器学习方向的，对计算机视觉，自然语言处理，语音识别等方向不太了解，所以这次着重讲一下如何转行到机器学习。
IMG_4676
在公司里面工作通常都需要有提取数据的工具，在大多数情况下就是写 SQL。SQL 是为了从数据库中提取数据，然后进行必要的数据过滤，数据分析，数据提取。对于 SQL，需要掌握的内容有以下几点：聚合函数，数学函数，字符串函数，表格的连接函数，条件语句等。SQL 的经典教材有两本，分别是：《HIVE编程指南》，作者 Edward Capriolo；《SQL基础教程》，作者 Mick。个人特别喜欢《SQL基础教程》，极易上手，易学易通。
目前工业界的机器学习编程语言很多，就个人浅显的经验来看，现在比较常见的编程语言还是 Python。Python 的话包括各种各样的工具包，例如 Numpy，Scipy，Scikit-Learn，Tensorflow 等等。其中，Scikitlearn 的文档是非常详细的，特别适合初学者入门学习。至于 Python 教材的话，其实有很多，例如：《Python基础教程》，作者是 Magnus Lie Hetland，这本书特别适合初学者看。如果是网络教材的话，推荐参考廖雪峰的官方网站，地址是：http://www.liaoxuefeng.com/。开发环境的话，公司一般都会使用 Linux，而不是 Windows 系统。在这里，特别推荐转行的同学掌握 Python 和 Linux。
其实，除了 SQL 和 Python 之外，C++ 或者 Java 也要选择一个方向来学。因为做机器学习的时候，除了离线使用 SQL 提取各种数据，用 Python 或者各种大数据工具来进行必要的模型训练，另外一个需要做的事情就是模型的上线工作。而模型的上线的时候是需要根据实际的数据来进行模型的预测，在模型预测的时候，通常来说用 Python 的话效率会有瓶颈，这种时候都会换成 C++ 或者 Java，因此，如果大家有时间准备的话，可以考虑把 SQL，Python，C++/Java 一起准备了。
除了基础的工具之外，其实做机器学习必不可少的那就是大学数学。在一般情况下，如果要读机器学习的普通书籍，数学系二年级左右的数学课程基本上就够用了。在一些特殊的时候需要读数学系更高年级的课程，不过在工作中还是相对偏少一些。如果是数学系的优秀学生的话，这一块基本上不需要特别担心，因为机器学习的数学没有数学系所教的那么难，机器学习更强调的是应用。基本上，数学分析，高等代数，概率论与数理统计，离散数学等课程其实就够用了。除了数学之外，计算机系的数据结构和算法的课程也是需要学习的，同时也应该多做一些算法题目，毕竟面试的时候是很有可能考这些算法题目的。
除了这些计算机与数学的基础知识之外，下面就要开始讲机器学习的入门了。通常来说，随便翻开一本机器学习的教材，都能够看到以下内容。
  1. 线性回归
  2. 决策树
  3. 朴素 Bayes
  4. 神经网络
  5. 集成学习
  6. 强化学习
在这里，推荐给大家看的教材有三本：
  1. 《机器学习实战》，作者是 Peter Harrington，
  2. 《机器学习》，作者是周志华
  3. 《Scikit-learn 与 Tensorflow 机器学习实用指南》，作者是 Adrelien Geron。
通过这几本书的学习，其实按照书上的目录和内容把机器学习过一遍，基本上就能够掌握机器学习的绝大部分知识点了。如果想学计算机视觉方向或者自然语言处理方向的话，建议再去看相关的书籍和教材，这一方面的知识点的话，Stanford 有不少优秀的课程。

转专业的人如何求职机器学习

其实，绝大部分的人都是需要求职的，在这种情况下，如何在转专业的时候脱颖而出就是一个关键的问题。其实，上面所说的知识点，无论是自学机器学习，还是努力刷题目，其实都是可以通过自我学习得到的，下面来说一些不太可能通过自己就能够得到的技能。
IMG_4677
1. 机器学习的竞赛。通常来说，无论是计算机视觉领域，NLP 领域，还是机器学习领域，都会有着各种各样的比赛。也就是公司或者竞赛的主办方会在网上公开竞赛的题目，并且提供必要的数据，让大家在该数据集上进行比赛，并且得分高的队伍获胜。而在比赛的过程中，通常都会用到各种各样的机器学习知识，这个是练习机器学习能力的第一步，也就是通过比赛来检验之前的学习是否达到了一个不错的效果。而在比赛的途中，建议还是通过组队的方式，几个人共同完成一项比赛。而且优秀的比赛结果其实最终也是可以写到简历里面的。
2. 实习经历。其实打比赛这件事情，自己做也是可以的。但是实习这个经历通常来说还是比较重要的，在一般情况下，公司选人都会从实习生当中来选，如果有合适的实习生，通常就不需要继续在校园招聘中寻找人才了。所以，找实习对应届生来说应该是比较重要的事情，而且这件事情建议早点做，而不要等到最后找工作的那一年才开始。一般来说，第一份实习可能比较难找一些，但是如果在低年级的暑假或者寒假就有实习机会的话，其实对于后续的实习或者就业是有很大的帮助的。通过实习，可以了解公司的一些业务和数据的情况，在这种情况下，才能够逐步理解工业界的数据和学术界的数据的差异性。
其实在选择做机器学习方向的时候，可以考虑得更加清楚一些，因为无论是做 CV，NLP 还是其他，都面临着和计算机系的人竞争的场面。在这里提个醒，其实机器学习并不是这两年才发展起来的，很多年前机器学习技术就已经存在了。最理想的状况就是，把机器学习的技术运用到本领域中，因为本专业的领域知识是计算机系的人并不具备的，例如，做气象研究，做金融分析。这类研究方向其实对本专业的人更加友好，如果能够把机器学习的知识应用到这些方向，那么求职的时候在本方向就更有竞争力，也会更有优势。而计算机系的人如果要转行做金融之类的，其实也要学习金融方面的知识，不过计算机系的人应该还是会倾向于去更熟悉的环境或者企业找工作。
对于转专业的人来说，如果去互联网的企业寻找工作，其实就已经算跨行了。这种时候其实是占有一定劣势的。所以，需要做的事情就是尽快补齐一些知识上的不足，并且尽快找一份实习，通过实习来提升自己的实战经验，最差也要多做比赛项目。除此之外，如果要寻找工作的话，能够找到熟悉的人帮忙内部推荐一下，拿到一些面试的机会或者名额其实也是很关键的。内推并不是保送的意思，而是免除了筛选简历的过程，直接就能够拿到面试的机会。

转专业的工作感受

IMG_4679

1. 给自己压力。一般来说，转专业求职是一个艰苦的过程，但是入职之后的生活则更加辛苦。因为公司的考核是每半年甚至两个月就一次，所以，在这种情况下，任何人都需要有一个上手的速度。有的人因为在学校学过相关的内容，或者之前实习过，因此上手的时候比较快；但是有的人转专业就面临上手慢的情况。其实这些对于应届生来说都可以理解，毕竟所有的人都需要有一个适应的过程。在这种情况下，在工作的初期一定要给自己一定的压力。意思就是说：在刚工作的第一年，每三个月就要让自己有一个飞速的提升；在工作的第二年，每半年就要让自己有一个提升；后续的话，每一年都要让自己有提升才是关键。因此，无论是本专业还是转专业的同学，都建议在前两年工作的时候，多给自己一些压力，只有这样，才能够让自己有更好的进步空间。
2. 业务经验。公司里面有很多东西并不是直接使用开源代码就能够发挥作用的，在公司里面无论做什么事情，最重要的一点就是对业务的理解。在对业务的理解方面，老员工相对于新人来说确实有着不少的优势。其次，在做业务的过程中，通常都会经历很多的坑，无论是别人主动挖的，还是自己踩坑踩出来的，都是自身宝贵的财富和经验。而这些经验只能够通过靠做大量的业务来获得。如果要想长期保持自身的优势，通过长期的训练和学习确实是一个有效的办法。无论是天才还是普通人，要想提升自身的技术，不花一定的时间去学习是不可行的。因此，无论在任何时候都不能够放弃让自己学习和充电的机会。
3. 勇于接受新的挑战。公司里面除了已有的项目之外，通常来说都会开启各种各样的新项目，在这种情况下，如果有机会做新的项目，也就是别人没有做过的项目。这种机会已经要把握住，因为对于新人来说，能够接触全新的项目肯定是好过维护已有的项目的。但是几乎所有的人都是从维护旧的项目开始的，只有旧的项目做好了，才有机会拿到新的项目。
4. 不要永远抱着已有的方向不放手。在公司里面，业务方向总会或多或少的发生变化，随着部门的调整，方向的变化，所做的内容总会发生一些变化。在工作的时候，最好不要抱着我就是来做这个方向的，除了这个方向之外其他的内容我一概不想做。因为当时的工作岗位未必能够提供你想做的方向，但是说不定能够提供其他的研究方向。有的时候，在公司里面，根据方向的变化来调整自己的工作内容也是一个必要的技能。而且，在公司的时候，一定要多做一些有挑战的项目，只有通过这些项目，才能够让自己的技术壁垒更加深厚。当然，在求职的时候，每个人都有着自己的想法和选择，所以，在求职的时候，是可以选择一个自己喜欢的方向来做的。
 

机器学习还能持续多久

IMG_4681
如果是在数学界，要想成为一个数学工作者，一般来说都要经过以下的学习路程：
数学分析/高等代数/概率论/复变函数/实变函数/泛函分析/微分几何/抽象代数
等一系列越来越难的课程的学习，而以上的这些课程只是基础课，连数学科研的边还没碰到。所有想要从事数学工作的人，都必须一步一步地，从头建立自己的数学知识体系，完善自己的数学工具库。而且数学的学习路径没有捷径，除非人绝顶聪明之外，都是需要一步一步，一年一年的来花费时间学习数学，才能够逐步体会数学所带来的奥妙。
与之截然不同的是 AI 领域（Machine Learning, CV, NLP, 语音等），一般来说只需要学习微积分/线性代数/概率论就可以基本上看懂机器学习的相关课程，当然要想深入学习 AI 的话还是需要很多数学基础的。随着科技的发展，各种开源工具的层出不穷，很多学校的学生甚至工业界的人士都已经不需要从底层从头开始，一步一步地建立自己的工具库。根据各种丰富的文档和 Blog，不少人都可以快速上手各种 AI 的工作内容，无论是使用 Tensorflow 建立图像分类器，还是使用 XGBoost 刷竞赛的成绩。这种时候，从事 AI 相关工作的门槛将会比之前变得越来越低，毕竟从头开始手动写一个 BP 算法或者说 GBDT 算法还是有一定门槛的。
根据经济学的基础知识，供需关系与价格有着一定的关系。一旦人数过多，而市场上的蛋糕并没有那么大的时候，很多人就要降低自己之前的期望，甚至转行做其他的事情。
之前在学校读书的时候，就听一些老师说过，最近放出来一个助理教授的职位，但是收到了200-300封简历，全部都是北美欧洲或者国内名校的PHD。之前听说在1980年的时候，数学PHD还不需要做posdoc就可以找到工作；到了1990年，基本上都要做一两年的posdoc才可以找到下家；到了00年以后，回国的话马上还能找到一个不错的职位，虽然工作不高，但是对论文的要求也没那么高；等到了2010年以后，国家千人计划等项目的开启，没有在国外混到一个好职位的，没有几篇好文章的，基本上在国内就没法找到教职了。等到了2020年以后，还真不知道是什么样的行情了，进入好学校的要求肯定是越来越高，要求的论文数量也是越来越多，质量也是越来越高了。
如果现在有十个岗位，但是只有五个 AI 专业的人来应聘，当然这些人都能够找到工作；但是随着人工智能专业的开设，相关的本科生和研究生开始培养，AI 从业者将会变得越来越多，但是岗位是否能够得到相应的增加就不是特别清楚了。就之前的经验而言，数学系的学生之间在毕业的时候差距还是挺大的，有的很强，有的很差。相信在人工智能专业也会有类似的情况，优秀的学生总是占少数。
就笔者的经验实在是无法确定这一波 AI 浪潮能够持续多久，如果五六年之后这波浪潮还在，蛋糕越来越大，那么现在选择攻读 AI 相关专业的人将会是受益者。但是如果这波浪潮不在了，蛋糕保持稳定甚至缩小的时候，AI 相关专业的人的竞争将会变得更加激烈。无论是工业界还是学术界的竞争，将会比现在的情况变大很多倍。而这波浪潮退去之后，能够留在沙滩上继续前进的永远都是少数人。
时间序列

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

3 Comments

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/

 

zr9558's Blog