LSTM单变量时间序列预测及异常值检测

★★★ 本文源自AlStudio社区精品项目，【点击此处】查看更多精品内容 >>>

LSTM预测及异常值检测

本文主要参考了这篇文章：Anomaly Detection with LSTM in Keras，原文GitHub地址

对纽约市出租车需求量进行预测，并同时希望能检测出异常值。这篇文章运用了许多简单项目涉及不到的小技巧，因此我认为还是比较有价值的。

本文基本完整复现了原文的所有细节，还添加了个人的理解和一些新的细节（包括它文章提到的图但是代码没给出的图，顺便修复了一些小BUG）

原文是用keras训练神经网络的，这里用paddle2复现了一遍

对于原文中出现的难度重点提供了大量参考资料

欢迎评论讨论交流！~

如果paddle版本更新代码无法运行了可以私聊我更新代码，我也计划每年都复习一遍自己写的项目并更新。

但是有50%的概率会鸽

数据集介绍

纽约市出租车需求数据集

数据集共包含两列

• timestamp：时间戳，时间范围2014-07-01至2015-01-32。间隔30分钟
• value：这段时间内，纽约市出租车的需求量

任务：构建一个模型，实现出租车需求预测以及异常检测

数据探索

数据读取

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
1
2
3

df = pd.read_csv('data/data106368/nyc_taxi.csv')
df.head()
1
2

观察一天的波动情况

可以看到，每天的0点到6点载客需求量较小，晚上是需求高峰期

df.iloc[0:48].plot(y='value', x='timestamp', rot=45, figsize=(8, 3))
1

<matplotlib.axes._subplots.AxesSubplot at 0x7efe1a903e10>
1

观察一周的波动情况

可以看出，数据还是具有比较明显的日周期（白天高于晚上）的

df.iloc[0:48 * 7].plot(y='value', x='timestamp', rot=45, figsize=(8, 3))
1

<matplotlib.axes._subplots.AxesSubplot at 0x7efe17cea390>
1

全期数据观察

数据段内的时间段出现了5个异常。它们分别发生在纽约马拉松、感恩节、圣诞节、元旦和暴风雪

目的是提前发现这些异常的观测结果！

出租车的需求也受到每周趋势的驱动: 在一周的某几天，出租车的需求比其他几天更高。

df.plot(y='value', x='timestamp', rot=45, figsize=(16, 3))
1

<matplotlib.axes._subplots.AxesSubplot at 0x7efe1a3e77d0>
1

特征提取

提取一些月份、日期、小时、星期特征，方便后续可视化

# 提取一些日期特征
df['timestamp'] = pd.to_datetime(df['timestamp'])
df['year'] = df['timestamp'].dt.year
df['month'] = df['timestamp'].dt.month
df['day'] = df['timestamp'].dt.day
df['hour'] = df['timestamp'].dt.hour
df['weekday'] = df['timestamp'].dt.weekday

df.head()
1
2
3
4
5
6
7
8
9

# 自相关性
time_lags = np.arange(1,10*48*7)
auto_corr = [df.value.autocorr(lag=dt) for dt in time_lags]

plt.figure(figsize=(12, 4))
plt.plot(1.0 / (48 * 7) * time_lags, auto_corr)
plt.xlabel('time lag [weeks]')
plt.ylabel('correlation coeff', fontsize=12)
1
2
3
4
5
6
7
8

Text(0,0.5,'correlation coeff')
1

模型搭建

这里的模型作者使用的是分位数回归（quantile regression），我一开始看这个项目的时候也不太了解什么是分位数回归，所以感觉有必要简单介绍一下：

参考资料：分位数回归（Quantile Regression）

直观理解：例如0.9分位数回归，就是希望回归曲线之下能够包含90%的数据点（y）

作用：得到单条分位数回归曲线是不够的，进一步可以画出不同的分位数回归曲线，这样才能能更加明显地反映出随着x的增大，y的不同范围的数据是如何不同程度地变化的，而这个结论通过以前的回归分析是无法得到的，这就是分位数回归的作用。

方式：本质上，这就是一个加权最小二乘法（虽然形式上有点不一样），给不同的y值（大于分位点和小于分位点的y）不同的权重，比如现在我们有一个数据集是1到10各整数，我们希望求0.7分位数回归曲线，则设q=0.7，因为希望模型希望回归曲线能包含70%的点，因此大于分为数的点会被赋予较高的损失权重，所有大于分位数q的数据点计算损失时赋上权重0.7，小于q的赋予权重0.3，然后最小化损失函数。但是实际上，可以定义不同的损失函数来满足上述要求。

本质：最根本的就是把损失函数从最小二乘法改成加权最小二乘法，通过不同的分位数得到不同的结果，再根据结果进行分析

参数设置和预处理

参数设置

参数主要有以下几个：

• seq_length预测时间窗（用多个个数据预测下一步数据）：这里设置为48，因为数据集时半小时采样一次数据点，48个数据点正好为一天
• train_size训练集长度：5000

num_samples = df.shape[0]  # 样本数
seq_length = 48  # 观测序列长度
train_size = 5000  # 训练集长度
1
2
3

预处理

首先构造了一个组合特征

weekday_hour是星期和小时的组合特征，例如下表第一行第一列的 1 0表示星期二的0点。（pandas提取的星期特征是从0开始的，0到6对应周一到周日）

m_weekday 是基于weekday_hour统计的租车载客流量均值特征，例如下表第一行第二列就表示：统计范围内所有星期二0点的平均载客流量是8774.433

这里的统计范围指数据集的前5000条数据，5000就是训练数据的大小train_size。

下面的特征构造涉及到的两个pandas函数对于不太了解熟悉的同学可能是个难点，下面附上了官网文档链接

pandas.Series.replace 官方文档

pandas.DataFrame.groupby 官方文档

df['weekday_hour'] = df['weekday'].astype(str) +' '+ df['hour'].astype(str)
df['m_weekday'] = df['weekday_hour'].replace(df[: train_size].groupby('weekday_hour')['value'].mean().to_dict())
df[['weekday_hour', 'm_weekday']].head(3)
1
2
3

然后对数据进行标准化，得到之后用于模型训练的标准化后的一维时间序列数据。

标准化：这里采用的标准化不是普通的标准化。刚刚得到一个特征m_weekday 表示某星期某小时的载客流量平均值

对原始数据载客流量 value 列 和 构造特征 m_weekday列先取对数然后做差，就得到了标准化后的数据。

为什么要这么做？

我的理解是，这样构造出来目标列既包含了星期和小时特征，又包含了均值统计特征，同时完成了归一化的任务

而且不同于普通均值归一化，对所有数据求均值，然后再做除法，这里做了更细粒度的均值归一化，一定程度上消除了每日波动的影响

这里的归一化计算方式是取对数做差，就相当于先做除法再取对数

# 取对数+做差
df['target'] = np.log(df['value']) - np.log(df['m_weekday'])
1
2

构建数据集

import paddle

# 数据集类
class MyDataset(paddle.io.Dataset):
    """
    步骤一：继承paddle.io.Dataset类
    """
    def __init__(self, data, num_features=48, num_labels=1):
        """
        步骤二：实现构造函数，定义数据集大小

        data: numpy.Array 1维数组
        """
        super(MyDataset, self).__init__()
        self.data = data
        self.num_features = num_features
        self.num_labels = num_labels

        x = []
        y = []
        for i in range(0, len(data) - num_features - num_labels + 1):
            x.append(data[i:i+num_features])
            y.append(data[i+num_features:i+num_features+num_labels])

        
        self.x = np.vstack(x).reshape(-1, self.num_features, 1)
        # self.y = np.vstack(y).reshape(-1, self.num_labels, 1)
        self.y = np.vstack(y)
        self.x = np.array(self.x, dtype="float32")
        self.y = np.array(self.y, dtype="float32")

        self.num_samples = len(x)

    def __getitem__(self, index):
        """
        步骤三：实现__getitem__方法，定义指定index时如何获取数据，并返回单条数据（训练数据，对应的标签）
        """
        data = self.x[index]
        label = self.y[index]

        return data, label

    def __len__(self):
        """
        步骤四：实现__len__方法，返回数据集总数目
        """
        return self.num_samples
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

train_data = df['target'][:train_size]
test_data = df['target'][train_size - seq_length:]

train_dataset = MyDataset(train_data)
test_dataset = MyDataset(test_data)

train_loader = paddle.io.DataLoader(train_dataset, batch_size=128, shuffle=False, drop_last=True)
1
2
3
4
5
6
7

模型搭建

这里使用的模型和原文一致，两层

import paddle

# 定义LSTM模型
class LSTM(paddle.nn.Layer):
    def __init__(self):
        super().__init__()
        self.rnn = paddle.nn.LSTM(input_size=1, hidden_size=64, direction='bidirectional')
        # input_size 为64 * 2 因为上一层的LSTM输出把输出变成64，而因为是双向LSTM，所以是64*2
        self.rnn2 = paddle.nn.LSTM(input_size=64*2, hidden_size=16, direction='bidirectional')

        self.dropout1 = paddle.nn.Dropout(p=0.5, mode='downscale_in_infer')
        self.dropout2 = paddle.nn.Dropout(p=0.5, mode='downscale_in_infer')

        # 同理，Linear层接收16*2的输入
        self.linear = paddle.nn.Linear(16*2, 50)

        self.out10 = paddle.nn.Linear(50, 1)  # 10%分位数回归输出
        self.out50 = paddle.nn.Linear(50, 1)  # 50%分位数回归输出
        self.out90 = paddle.nn.Linear(50, 1)  # 90%分位数回归输出

        # 也可以直接, 然后后续可以省略拼接的步骤
        # 写成上面那种方式是为了帮助理解
        # self.out = paddle.nn.Linear(50, 3)
        
    def forward(self, inputs):
        y, (hidden, cell) = self.rnn(inputs)
        y = self.dropout1(y)
        _, (out, _) = self.rnn2(y)  # 这一步输出的hidden作为下一步的输入
        
        out = self.dropout2(out)
        # 把前后向LSTM输出的 hidden 拼接起来
        out = paddle.concat([out[0], out[1]], axis=-1)
        # 以下两行代码效果同上一行代码
        # out = paddle.transpose(out, (1, 0, 2))
        # out = paddle.reshape(out, (out.shape[0], -1))
        out = self.linear(out)
        out10 = self.out10(out)
        out50 = self.out50(out)
        out90 = self.out90(out)

        out = paddle.concat([out10, out50, out90], axis=-1)
        return out
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

# 测试LSTM模型
input_size = 1  # 输入维度 单变量预测的话这个就是1 
seq_len = 48  # 用多少个数据预测后一个数据
batch_size = 4  # 一批4个样本长度为10

x = paddle.randn((batch_size, seq_len, input_size))
model = LSTM()

# x shape (batch_size, seq_len, input_size)
# 输出是三个分位数的回归预测
y_pred = model.forward(x)
y_pred
1
2
3
4
5
6
7
8
9
10
11
12

Tensor(shape=[4, 3], dtype=float32, place=CUDAPlace(0), stop_gradient=False,
       [[-0.07357128, -0.00638272, -0.00222786],
        [-0.00789048, -0.03552523,  0.04117572],
        [ 0.03548329,  0.08265428,  0.12370239],
        [-0.08020647,  0.02704403, -0.01321218]])
1
2
3
4
5

paddle.summary(model, (4, 64, 1))
1

-----------------------------------------------------------------------------------------------
 Layer (type)       Input Shape                    Output Shape                   Param #    
===============================================================================================
    LSTM-2          [[4, 64, 1]]     [[4, 64, 128], [[2, 4, 64], [2, 4, 64]]]     34,304     
   Dropout-1       [[4, 64, 128]]                  [4, 64, 128]                      0       
    LSTM-3         [[4, 64, 128]]    [[4, 64, 32], [[2, 4, 16], [2, 4, 16]]]      18,688     
   Dropout-2        [[2, 4, 16]]                    [2, 4, 16]                       0       
   Linear-1          [[4, 32]]                       [4, 50]                       1,650     
   Linear-2          [[4, 50]]                        [4, 1]                        51       
   Linear-3          [[4, 50]]                        [4, 1]                        51       
   Linear-4          [[4, 50]]                        [4, 1]                        51       
===============================================================================================
Total params: 54,795
Trainable params: 54,795
Non-trainable params: 0
-----------------------------------------------------------------------------------------------
Input size (MB): 0.00
Forward/backward pass size (MB): 0.57
Params size (MB): 0.21
Estimated Total Size (MB): 0.78
-----------------------------------------------------------------------------------------------



/opt/conda/envs/python35-paddle120-env/lib/python3.7/site-packages/numpy/core/fromnumeric.py:87: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  return ufunc.reduce(obj, axis, dtype, out, **passkwargs)





{'total_params': 54795, 'trainable_params': 54795}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

自定义分位数损失函数

分位数拟合的损失函数定义得比较简单。

设真实值-预测值=error：

• 若真实值小于预测值, error < 0：误差就是(q-1) * error
• 若真实值大于预测值, error > 0：误差就是q * error

可以出，误差始终为正

为什么这么定义损失函数：

如果error < 0，则说明预测值比真实值大，假设 q=0.1, 表示希望只有10%的数据小于预测曲线，因此如果预测值比真实值大，就会给予较高的误差

为了给予loss就是0.9 * abs(error)

自定义loss函数的参数设置可以参考：

• 自定义Loss
• mse_loss

注意：自定义loss函数的参数固定为先预测值后真实值

# 自定义损失函数
def q_loss(y_pred, y_true):
    """ y_true shape (batch_size, 3) """
    q = paddle.to_tensor([0.1, 0.5, 0.9])  # 计算分位数
    # weight = 1 / y_pred.shape[-1]
    error = y_true - y_pred  # 计算误差
    # 原文给每个分位数都添加了一个0.3的权重，加在一起后权重就是0.9，可以不要这个0.9
    # 但为了准确还原我还是加上了
    return paddle.mean(paddle.maximum(q * error, (q-1) * error)) * 0.9
1
2
3
4
5
6
7
8
9

# 测试
y_true = paddle.rand((4, 3))
y_pred = paddle.rand((4, 3))
q_loss(y_true, y_pred)
1
2
3
4

/opt/conda/envs/python35-paddle120-env/lib/python3.7/site-packages/paddle/tensor/creation.py:125: DeprecationWarning: `np.object` is a deprecated alias for the builtin `object`. To silence this warning, use `object` by itself. Doing this will not modify any behavior and is safe. 
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  if data.dtype == np.object:





Tensor(shape=[1], dtype=float32, place=CUDAPlace(0), stop_gradient=True,
       [0.16217309])
1
2
3
4
5
6
7
8
9
10

配置优化器并训练

model = LSTM()
loss_fn = q_loss

optimizer = paddle.optimizer.Adam(learning_rate=0.001,
                                 parameters=model.parameters())

for epoch in range(50):
    for batch, (batch_x, batch_y) in enumerate(train_loader()):
        y_pred = model(batch_x)
        
        # loss = F.mse_loss(y_pred, batch_y, reduction='mean')
        loss = loss_fn(y_pred, batch_y)
        loss.backward()
        optimizer.step()
        optimizer.clear_grad()
    print("epoch {} loss {:.4f}".format(epoch+1, float(loss)))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

epoch 1 loss 0.0280
epoch 2 loss 0.0283
epoch 3 loss 0.0294
epoch 4 loss 0.0278
epoch 5 loss 0.0264
epoch 6 loss 0.0244
epoch 7 loss 0.0261
epoch 8 loss 0.0242
epoch 9 loss 0.0244
epoch 10 loss 0.0242
epoch 11 loss 0.0232
epoch 12 loss 0.0231
epoch 13 loss 0.0221
epoch 14 loss 0.0205
epoch 15 loss 0.0213
epoch 16 loss 0.0205
epoch 17 loss 0.0192
epoch 18 loss 0.0190
epoch 19 loss 0.0180
epoch 20 loss 0.0172
epoch 21 loss 0.0183
epoch 22 loss 0.0173
epoch 23 loss 0.0168
epoch 24 loss 0.0175
epoch 25 loss 0.0147
epoch 26 loss 0.0144
epoch 27 loss 0.0159
epoch 28 loss 0.0144
epoch 29 loss 0.0161
epoch 30 loss 0.0155
epoch 31 loss 0.0150
epoch 32 loss 0.0149
epoch 33 loss 0.0141
epoch 34 loss 0.0155
epoch 35 loss 0.0147
epoch 36 loss 0.0146
epoch 37 loss 0.0151
epoch 38 loss 0.0150
epoch 39 loss 0.0137
epoch 40 loss 0.0124
epoch 41 loss 0.0148
epoch 42 loss 0.0136
epoch 43 loss 0.0134
epoch 44 loss 0.0146
epoch 45 loss 0.0157
epoch 46 loss 0.0138
epoch 47 loss 0.0156
epoch 48 loss 0.0142
epoch 49 loss 0.0155
epoch 50 loss 0.0147
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

结果展示

提取测试集特征和标签

test_x = paddle.to_tensor(test_dataset.x, dtype="float32")
y_true = test_dataset.y
# 反归一化
y_true = np.exp(y_true + np.log(df['m_weekday'][5000:].values.reshape(-1, 1)))

model.train()  # 确保模型还处于训练阶段
# y_pred = model(test_x).numpy()

1
2
3
4
5
6
7
8

源代码采用的是重复预测100次结果

然后对100次结果按照不同的分位数等级，利用np.quantile获得最终的预测结果

为什么要这么做？

在应用神经网络训练模型时，经常出现由于内部权重初始化导致结果的不确定的问题。

分位数回归预测似乎特别受此类问题的困扰； 即计算分位数预测时分位数不能重叠，这没有意义！

为了避免这个问题，作者在预测阶段使用了自助采样的技巧?（bootstrapping）：

作者没有直接采用预测结果，而是重新激活了网络的 dropout层（可训练：模型中为 true），迭代预测 100 次

存储它们并最终计算所需的分位数（作者认为这是一个trick，有助于改善预测结果）。

import tqdm

pred_10 = []
pred_50 = []
pred_90 = []

# 重复预测100次结果
for i in tqdm.tqdm(range(0,100)):
    
    predd = model(test_x).numpy()
    pred_10.append(predd[:, 0])
    pred_50.append(predd[:, 1])
    pred_90.append(predd[:, 2])

pred_10 = np.stack(pred_10)
pred_50 = np.stack(pred_50)
pred_90 = np.stack(pred_90)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

100%|██████████| 100/100 [00:02<00:00, 37.92it/s]
1

下面以图形方式解释了这个过程，三条不同颜色的区域代表不同的预测结果。

设定好不同的分位数，给予90%分位数回归结果90%分位点，计算了它们的汇总度量（红线），避免分位数之间发生交叉。

# 还原回真实值 shape (100, test_size)
pred_10 = np.exp(pred_10.T + np.log(df['m_weekday'][5000:].values.reshape(-1, 1)))
pred_50 = np.exp(pred_50.T + np.log(df['m_weekday'][5000:].values.reshape(-1, 1)))
pred_90 = np.exp(pred_90.T + np.log(df['m_weekday'][5000:].values.reshape(-1, 1)))
1
2
3
4

# 汇总：按照分位数等级获取最终的预测结果
pred_90_m = np.quantile(pred_90, 0.9,axis=-1)
pred_50_m = pred_50.mean(axis=-1)
pred_10_m = np.quantile(pred_10,0.1,axis=-1)
1
2
3
4

# 这一步绘图时间较久
start, end = 30, 35
plt.figure(figsize=(12,8))
l1 = plt.plot(pred_90[start:end], c='cyan', alpha=.5)
l2 = plt.plot(pred_50[start:end], c='blue', alpha=.5)
l3 = plt.plot(pred_10[start:end], c='green', alpha=.5)
l4 = plt.plot(pred_90_m[start:end], c='r', lw=3)
plt.plot(pred_50_m[start:end], c='r', lw=3)
plt.plot(pred_10_m[start:end], c='r', lw=3)
plt.legend([l1[0], l2[0], l3[0]], ['pred_90', 'pred_50', 'pred_10'])
plt.show()
1
2
3
4
5
6
7
8
9
10
11

评估指标和预测结果

模型在第 50 分位数预测出租车需求方面达到了很好的性能。

大约 0.055 均方对数误差是一个很好的结果！

这意味着 LSTM 网络能够理解驱动出租车需求的潜在规则。

所以我们的异常检测方法听起来很棒……我们计算了第 90 个分位数预测和第 10 个分位数预测之间的差异，看看发生了什么。

关于均方对数误差，可以看做是能够反映比例的误差
$$L(y,\hat{y})=\frac{1}{N}\sum _{i=0}^N(\log (y_i+1)-\log ({\hat{y}}_i+1){)}^2$$

参考资料: Mean squared logarithmic error (MSLE)

一开始，我训练出的模型计算得出的msle与原文差异较大(msle = 1.1108，原文msle = 0.055)

后来经排查发现，可能是dropout mode的设置问题

默认设置是：upscale_in_train, 在训练时增大输出结果。

我改成了downscale_in_infer, 在预测时减小输出结果

原本的msle: 1.1108 改之后 1.0958，并没有多大变化

经过再次排查后发现，是y_ture的维度有问题（之前并没有y_true.ravel()）

y_true 的维度: (test_size, 1), pred_50_m 的维度 (test_size, )

造成二者相减后的维度为: (test_size, test_size)

下面代码中改成y_true.ravel()后，问题解决，msle的计算恢复正常

msle = np.mean((np.log(y_true.ravel() + 1) - np.log(pred_50_m + 1))**2)
print('均方对数误差: ', msle)
1
2

均方对数误差:  0.053641254799958224
1

# 观察预测值和真实值之间的偏差确实也不大
plt.plot(pred_50_m[100:150], label='pred')
plt.plot(y_true[100:150], '.r', label='true')
plt.show()
1
2
3
4

# 预测结果绘图
test_date = df.timestamp.values[5000:]

plt.figure(figsize=(16,8))
plt.plot(test_date, pred_90_m, color='cyan')
plt.plot(test_date, pred_50_m, color='blue')
plt.plot(test_date, pred_10_m, color='green')
plt.show()
1
2
3
4
5
6
7
8

90分位数与10分位数差异绘图

对第 90 分位数预测值和第 10 分位数预测值做差，并绘制散点（蓝色），看看发生了什么

粉色线为真实值

在异常区间，10和90分位数差异更大（蓝点更高）。

5个异常时间点分别是：纽约马拉松、感恩节、圣诞节、元旦和暴风雪。

我就不标出这5个时间点了，具体可以看原文

plt.figure(figsize=(16,8))
plt.plot(test_date, y_true, color='red', alpha=0.4)
plt.scatter(test_date, pred_90_m - pred_10_m)
大（蓝点更高）。

5个异常时间点分别是：纽约马拉松、感恩节、圣诞节、元旦和暴风雪。

> 我就不标出这5个时间点了，具体可以看原文


```python
plt.figure(figsize=(16,8))
plt.plot(test_date, y_true, color='red', alpha=0.4)
plt.scatter(test_date, pred_90_m - pred_10_m)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15