LSTM时间序列预测

1、tensorflow张量的形状

tensorflow中张量(tensor)的属性——维数(阶)、形状和数据类型
在形状的中括号中有多少个数字，就代表这个张量是多少维的张量。

形状的第一个元素要看张量最外边的中括号中有几个元素(被最外边的中括号里边的内中括号括起来的所有数据算作一个元素)被逗号隔开，有n1个则这个张量就是n1维的，形状的第一个元素就是n1;
形状的第二个元素要看张量中最左边的第二个中括号中有几个被逗号隔开的元素，有n2个则shape的第二个元素就是n2;形状的第二个元素之后的第3,4…n个元素依次类推，分别看第n个中括号中有几个元素即可：

1    # 形状为[]

[1,2,3]   # 形状为[3]

[[1,2],[3,4]]   # 形状为[2,2]

[[[1,2],[3,4]],[[1,2],[3,4]]]   # 形状为[2,2,2]

2、RNN

知乎大佬，写的很棒
LSTM物理架构

3、LSTM

长短期记忆（Long short-term memory, LSTM）是一种特殊的RNN，主要是为了解决长序列训练过程中的梯度消失和梯度爆炸问题。简单来说，就是相比普通的RNN，LSTM能够在更长的序列中有更好的表现。
LSTM结构（图右）和普通RNN的主要输入输出区别如下所示。

3.1 深入了解LSTM结构

4、keras的LSTM

LSTM一般堆叠2-3层就可以了

4.1 参数return_sequences含义

如果设置return_sequences = True，该LSTM层会返回每一个time step的h，那么该层返回的就是1个由多个h组成的2维数组了，如果下一层不是可以接收2维数组的层，就会报错。所以一般LSTM层后面接LSTM层的话，设置return_sequences = True，如果接全连接层的话，设置return_sequences = False。

5、LSTM实现PM2.5的预测

用24h前的特征数据，预测未来2h的pm2.5

# -*- coding: utf-8 -*- 
# @Time : 2021/9/9 14:20 
# @Author : Orange
# @File : mylstm.py

import warnings

warnings.filterwarnings("ignore")

import pandas as pd
import numpy as np
from sklearn import metrics
from keras.models import Sequential
from keras.layers import LSTM
from keras.layers import Dense, Activation
from keras.layers import Dropout
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

filepath = 'PRSA_data.csv'  # 文件路径
data = pd.read_csv(filepath, index_col=0)
print("data.shape", data.shape)
print(data.head(10))  # 输出前10行数据

# 查看每一列的缺失值的个数
print("每一列的缺失值的个数:", data.isnull().sum())

# 按月求平均并填充按月填充空值
data["pm2.5"] = data.groupby(['year', 'month'])['pm2.5'].transform(lambda x: x.fillna(x.mean()))
# 查看填充后数据缺失情况
print("每一列的缺失值的个数:", data.isnull().sum())

# 将cbwd指标数据类型转换成category
cbwd_category = data['cbwd'].astype('category')
# 使用标签的编码作为真正的数据
data['cbwd'] = cbwd_category.cat.codes

# 特征X中的pm2.5数据为前1天的24小时PM2.5浓度，而目标变量是第二天0时和1时的PM2.5浓度
X_data = data[['pm2.5', 'DEWP', 'TEMP', 'PRES', 'cbwd', 'Iws', 'Is', 'Ir']]  # 提取特征数据
Y_data = data[['pm2.5']]  # 提取pm2.5数据

name = ['pm2.5', 'DEWP', 'TEMP', 'PRES', 'cbwd', 'Iws', 'Is', 'Ir']

# 按天聚合，用前一天的24h数据预测今天的前两个小时数据
# X第1天到1825天的24h特征数据
X = np.array(
    data.groupby(by=['year', 'month', 'day']).apply(lambda x: np.array([x[i].values for i in name]).T).values.tolist())[
    :-1]

data_new = data[data['hour'].apply(lambda x: x == 0 or x == 1)]  # 只保留小时为0和1的数据
# Y第2天到第1826天的前2h的pm2.5数据
Y = np.array(
    data_new.groupby(by=['year', 'month', 'day']).apply(lambda x: x['pm2.5'].values).values.tolist())[1:]

# 这里将80%数据作为训练集，20%数据作为测试集，则训练集有1460个样本，测试集有365个样本
X_train = X[:int(X.shape[0] * 0.8)]
Y_train = Y[:int(X.shape[0] * 0.8)]

X_val = X[int(X.shape[0] * 0.8)::]
Y_val = Y[int(X.shape[0] * 0.8)::]

# 标准化
X_mean = X_train.mean(axis=0)
X_std = X_train.std(axis=0)
Y_mean = Y_train.mean(axis=0)
Y_std = Y_train.std(axis=0)

X_train_norm = (X_train - X_mean) / X_std
Y_train_norm = (Y_train - Y_mean) / Y_std
X_val_norm = (X_val - X_mean) / X_std
Y_val_norm = (Y_val - Y_mean) / Y_std

# 建模
# 使用3层LSTM，输出层为2输出的Dense层
model = Sequential()
model.add(LSTM(32,
               input_shape=(X_train_norm.shape[1], X_train_norm.shape[-1]),
               return_sequences=False))
# model.add(Dropout(0.2))  # 防止过拟合
# model.add(LSTM(32, return_sequences=True))
# model.add(Dropout(0.2))
# model.add(LSTM(32))
model.add(Dropout(0.2))

model.add(Dense(2))
model.add(Activation('linear'))
model.compile(loss='mse',
              optimizer='rmsprop')
model.summary()

## 训练模型
history = model.fit(X_train_norm, Y_train_norm,
                    epochs=60,
                    batch_size=128,
                    validation_data=(X_val_norm, Y_val_norm))

loss = history.history['loss']
val_loss = history.history['val_loss']

plt.plot(range(len(loss)), loss, 'b-', label='训练集损失')
plt.plot(range(len(loss)), val_loss, 'r-', label='测试集损失')
plt.legend(loc='best')
plt.show()

model_pred = model.predict(X_val_norm)
val_pred = model_pred * Y_std + Y_mean  # 别忘了，数据进行了标准化处理，因此预测值需要处理，再计算R方

# 计算R2
R_2_0 = metrics.r2_score(Y_val[:, 0], val_pred[:, 0])  # 计算0时预测的R方
R_2_1 = metrics.r2_score(Y_val[:, 1], val_pred[:, 1])  # 计算1时预测的R方

# 实际值与预测值对比图
fig = plt.subplots(figsize=(12, 6))
gs = gridspec.GridSpec(2, 1, height_ratios=[1, 1], hspace=0)

ax1 = plt.subplot(gs[0])
plt.plot(range(Y_val.shape[0]), Y_val[:, 0], 'b-', label='0时PM2.5实际图')
plt.plot(range(Y_val.shape[0]), val_pred[:, 0], 'r-', label='0时PM2.5预测图')
plt.legend(loc='best')
plt.text(150, 400, '拟合R2：{0}%'.format(round(R_2_0 * 100, 2)))

ax2 = plt.subplot(gs[1], sharex=ax1)
plt.plot(range(Y_val.shape[0]), Y_val[:, 1], 'b-', label='1时PM2.5实际图')
plt.plot(range(Y_val.shape[0]), val_pred[:, 1], 'r-', label='1时PM2.5预测图')
ax2.set_xlabel(' ')
plt.legend(loc='best')
plt.text(150, 400, '拟合R2：{0}%'.format(round(R_2_1 * 100, 2)))
plt.show()