使用LSTM进行多步预测

使用LSTM对PM2.5进行2步预测

1.项目简介

该项目的目标为：使用前1天24小时内影响PM2.5的相关因素数据，对第二天两小时内PM2.5进行预测分析
本项目数据和Jupyter代码：数据及代码
数据简介：北京市2010年1月1日~2014年12月31日的PM2.5相关数据
数据频率：每小时记录一次PM2.5相关数据
数据来源：UCI机器学习库中北京PM2.5数据集
数据指标简介：
　　No：序列号
　　year：年份
　　month：月份
　　day：日
　　hour：小时
　　pm2.5：PM2.5浓度
　　DEWP：露点温度
　　TEMP：温度
　　PRES：压力
　　cbwd：风向
　　Iws：风速
　　Is：积雪的时间
　　Ir：累积的下雨时数

2.导入所需库

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
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

3.导入数据

filepath = 'E:/Jupyter/Mojin/CSDN/LSTM_PM2.5多步预测/PRSA_data_2010.1.1-2014.12.31.csv'  # 文件路径
data = pd.read_csv(filepath, index_col=0)
data.shape
data.head(10)  # 输出前10行数据

运行结果：
从上面的输出数据看到：
1.数据存在缺失值
2.cbwd指标为分类数据，需要进行整数编码

4.数据清洗

4.1缺失值处理

# 查看缺失值信息
data.isnull().sum()

运行结果：1.看到只有pm2.5指标存在2067个缺失数据
2.打开数据excel文件，观看发现存在连续2天及以上的缺失情况，因此这里采用月平均填充缺失数据。

data["pm2.5"] = data.groupby(['year', 'month'])['pm2.5'].transform(lambda x: x.fillna(x.mean()))
# 查看填充后数据缺失情况
data.isnull().sum()  # 可以发现不存在缺失值了

4.2分类数据处理

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

运行结果：

4.3构造数据集

1.明确项目目标：使用前1天24小时内PM2.5的相关数据，对第二天两小时内PM2.5进行预测分析
2.确定影响因素，即X：pm2.5、DEWP、TEMP、PRES、cbwd、Iws、Is和Ir共8个指标；目标变量：第二天2小时内的pm2.5
3.明确处理后的数据shape，以一条样本为例：其特征变量shape为：(1,24,8)，目标变量shape为：(1,2)
4.原始数据有43824个样本，共43824/24=1826天，那么处理后的样本量为(43824/24)-1=1825个样本。故处理后的数据shape：X的shape为(1825,24,8)，Y的shape为(1825,2)
5.注意：特征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数据

X = np.zeros((X_data.shape[0] // 24 - 1, 
             24, 
              X_data.shape[-1]))
Y = np.zeros((Y_data.shape[0] // 24 - 1,
              2))
rows = range(0, X_data.shape[0] - 24, 24)
for i, row in enumerate(rows):
    X[i, :, :] = X_data.iloc[row : row + 24]
    Y[i, :] = [Y_data.iloc[row + 24],  Y_data.iloc[row + 24 + 1]]
X.shape
Y.shape

运行结果：
(1825,24,8)
(1825,2)

4.4拆分数据集

# 这里将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)::]

4.5数据标准化

注意：对于测试集的数据标准化只能使用训练集的均值和标准差，因为对于测试集的信息我们是一无所知的

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

5.建模

5.1构造模型

# 使用3层LSTM，输出层为2输出的Dense层
model = Sequential()
model.add(LSTM(32, 
               input_shape=(X_train_norm.shape[1], X_train_norm.shape[-1]), 
               return_sequences=True))
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.compile(loss='mse', 
              optimizer='rmsprop')
model.summary()

运行结果：

5.2训练模型

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();

运行结果：

5.3结果评估

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();

运行结果：

拟合R方分别为：92.52%和86.17%，感觉还可以！