时间序列分段算法 [Time series Breakout Detection]

在时间序列分析中，断点检测（breakout detection）是一个很基本的问题。

通过捕捉时序数据中的断点（breakout），来发现时序数据所表示的系统在过去是否发生了某种事件（event）。进而为系统诊断提供必要的数据支持。

为了实现对时序断点的检测，我们首先需要对时序的整体时序做拟合。

这里我们通过一条直线来拟合一段时序，如果时序的趋势发生了变化，则用多条直线来拟合整条时序数据。

如下是对一条波动规律明显的时序做拟合之后的结果。



每个红色线条的转折点，就是我们找到的断点。

以上数据是我们在实验环境下，为了检测算法效果而人工构造的一条时序。

那么，该算法在实际情况下表现如何？

一下是一条实际的股票价格时序数据。我们通过该算法进行断点检测，并将断点红红色线条连起来的效果：



更进一步，将拟合之后的线段图分段画出如下：



其中黑色表示原始时序，红色表示划分得到的断点，蓝线表示断点之间子序列的线性拟合结果。

 

算法介绍：

算法所使用的关键技术：

1. 单变量线性回归，用来拟合某一段时序

2. 动态规划算法，  用来全局最大化断点检测效果。

算法核心代码：

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "lsp.h"

static double loss(double * s, int n){
int i;
double t;
double g0   = 0.0, g1   = 0.0;
double h00  = 0.0, h01  = 0.0, h10  = 0.0, h11  = 0.0;
double hv00 = 0.0, hv01 = 0.0, hv10 = 0.0, hv11 = 0.0;
double l0, l1;

// grad and hessian matrix
for (i = 0; i < n; i++){
t    = s[i];
g0  += t;
g1  += t * (1.0 + i);
h00 += 1.0;
h01 += 1.0 + i;
h11 += (1.0 + i) * (1.0 + i);
}
h10 = h01;

// inverse of hessian
t = h00 * h11 - h01 * h10;
hv00 = h11 / t;
hv01 = hv10 = -h01 / t;
hv11 = h00 /t;

// the theta
l0 = hv00 * g0 + hv01 * g1;
l1 = hv10 * g0 + hv11 * g1;

// sqare loss
t = 0.0;
for (i = 0; i < n; i++){
t += (l0 + l1 * (i + 1) - s[i]) * (l0 + l1 * (i + 1) - s[i]);
}
return t;
}

int * lsp(double * ts, int n, int min_size, double beta, int *ol){

if (!ts || min_size < 2 || n < 2 * min_size || !ol){
return NULL;
}

// prev breakout point
int * prev = (int*)malloc(sizeof(int) * (n + 1));
memset(prev, 0, sizeof(int) * (n + 1));

// number of breakout point
int * num  = (int*)malloc(sizeof(int) * (n + 1));
memset(num, 0, sizeof(int) * (n + 1));

// F scores
double * F = (double*)malloc(sizeof(double) * (n + 1));
memset(F, 0, sizeof(double) * (n + 1));

// loss
double * lossv = (double*)malloc(sizeof(double) * (n + 1));
memset(lossv, 0, sizeof(double) * (n + 1));

for (int s = 2 * min_size; s < n + 1; ++s){
for (int t = min_size; t < s - min_size + 1; ++t){
//double ls = loss(ts + prev[t], t - prev[t]);
double ls = lossv[t];
double rs = loss(ts + t, s - t);
double as = loss(ts + prev[t], s - prev[t]);
double score = (as - ls - rs) * (t - prev[t]) * (s - t) /    \
((s - prev[t]) * (s - prev[t])) - num[t] * beta;
score += F[t];
if (score > F[s]){
num[s] = num[t] + 1;
F[s] = score;
prev[s] = t;
lossv[s] = rs;
}
}
}

int k = num
;
*ol = k;
int * re = (int*)malloc(sizeof(int) * k);
memset(re, 0, sizeof(int) * k);
int i = n;
while(i > 0){
if (prev[i])
re[--k] = prev[i];
i = prev[i];
}

free(prev);  prev  = NULL;
free(num);   num   = NULL;
free(F);     F     = NULL;
free(lossv); lossv = NULL;
return re;
}

算法复杂度上限为：O(n * n * n)