使用LSTM模型预测多特征变量的时间序列

数据集为印度气候天气预报数据，该数据集提供2013年1月1日至2017年4月24日在印度德里市的数据，存在 4 个参数分别是*meantemp, humidity, wind_speed, meanpressure，*其中meantemp为待预测变量，其余为特征变量

1. 实现流程

代码实现流程基于 LSTM 模型的时间序列预测任务，以下是其简要描述：

• 数据准备：

分别读取训练集和测试集的数据，进行异常值处理、归一化处理，最后确保它们在0-1的数 值范围内

将归一化后的数据转换成模型可用的格式，包括将时间序列数据转换成带有滑动窗口的输入特征矩阵和目标值数组

• 模型构建：

使用 Keras 库构建 LSTM 模型，该模型包含了一个 LSTM 层和若干个全连接层

输入形状是 (窗口大小, 特征数量)，其中窗口大小定义了时间序列的历史信息，特征数量表示每个时间步的特征数

• 模型训练：

编译模型时指定损失函数和优化器

最后拟合模型，将训练数据和标签传递给模型进行训练，并指定了训练的迭代次数和批次大小

• 模型评价：

计算

LSTM 模型的评价指标，包括均方误差（MSE）、均方根误差（RMSE）、平均绝对误差（MAE）和拟合优度（R²）

• 未来预测：

使用训练好的模型对测试集 中的最后一个时间窗口进行预测，并根据模型输出的结果更新序列

重复以上步骤，预测未来若干步的数值

通过这样的流程建立一个LSTM多特征变量时间序列预测模型，在划分时间窗口的步骤将给出 两种模式一种是把待预测的特征也纳入输入特征进行时间窗口划分，另外一种模式是不把待预测的特征纳入输入特征进行时间窗口划分

2. 代码实现

2.1 数据展示

          
import pandas as pd
          
import matplotlib.pyplot as plt
          
plt.rcParams['font.sans-serif'] = 'SimHei' # 设置中文显示
          
plt.rcParams['axes.unicode_minus'] = False
          

          
df_train = pd.read_csv('DailyDelhiClimateTrain.csv', index_col=0, parse_dates=['date'])
          
df_test = pd.read_csv('DailyDelhiClimateTest.csv', index_col=0, parse_dates=['date'])
          
df_train
      

其中df_train和df_test分别为模型的训练集和测试集，数据不存在缺失值

2.2 训练集时序图

          
plt.figure(figsize=(15, 7))
          
plt.subplot(2, 2, 1)
          
plt.plot(df_train['meantemp'], color='g',  alpha=0.3)
          
plt.title('meantemp时序图')
          
plt.grid(True)
          
plt.subplot(2, 2, 2)
          
plt.plot(df_train['humidity'], color='g',  alpha=0.3)
          
plt.title('humidity时序图')
          
plt.grid(True)
          
plt.subplot(2, 2, 3)
          
plt.plot(df_train['wind_speed'], color='g',  alpha=0.3)
          
plt.title('wind_speed时序图')
          
plt.grid(True)
          
plt.subplot(2, 2, 4)
          
plt.plot(df_train['meanpressure'], color='g',  alpha=0.3)
          
plt.title('meanpressure时序图')
          
plt.grid(True)
          
plt.show()
      

很明显在meanpressure时序图中数据存在异常值

2.3 测试集时序图

          
plt.figure(figsize=(20, 7))
          
plt.subplot(2, 2, 1)
          
plt.plot(df_test['meantemp'], color='g',  alpha=0.3)
          
plt.title('meantemp时序图')
          
plt.grid(True)
          
plt.subplot(2, 2, 2)
          
plt.plot(df_test['humidity'], color='g',  alpha=0.3)
          
plt.title('humidity时序图')
          
plt.grid(True)
          
plt.subplot(2, 2, 3)
          
plt.plot(df_test['wind_speed'], color='g',  alpha=0.3)
          
plt.title('wind_speed时序图')
          
plt.grid(True)
          
plt.subplot(2, 2, 4)
          
plt.plot(df_test['meanpressure'], color='g',  alpha=0.3)
          
plt.title('meanpressure时序图')
          
plt.grid(True)
          
plt.show()
      

2.4 Hampel滤波器异常值处理

          
def hampel(vals_orig, k=7, t0=3):
          
    vals_filt = np.copy(vals_orig)
          
    outliers_indices = []
          
    n = len(vals_orig)
          

          
    for i in range(k, n - k):
          
        window = vals_orig[i - k:i + k + 1]
          
        median = np.median(window)
          
        mad = np.median(np.abs(window - median))
          
        if np.abs(vals_orig[i] - median) > t0 * mad:
          
            vals_filt[i] = median
          
            outliers_indices.append(i)
          

          
    return vals_filt, outliers_indices
          
filtered_data, outliers_indices = hampel(df_train['meanpressure'])
          

          
go_over = df_train['meanpressure']
          
df_train['meanpressure'] = filtered_data
          

          
plt.figure(figsize=(20, 7))
          
plt.subplot(2, 1, 1)
          
plt.plot(go_over, color='g',  alpha=0.3)
          
plt.title('meanpressure异常时序图')
          
plt.grid(True)
          
plt.subplot(2, 1, 2)
          
plt.plot(df_train['meanpressure'], color='g',  alpha=0.3)
          
plt.title('meanpressure Hampel时序图')
          
plt.grid(True)
          
plt.show()
      

在这里利用Hampel滤波器进行异常值处理， Hampel滤波器是一种基于中值和中值绝对偏差（MAD）的滤波器，旨在识别和去除时间序列数据中的异常值，该滤波器详细讲解参考往期文章时间序列异常值检验神器——Hampel滤波器

2.5 数据0-1标准化

          
from sklearn.preprocessing import MinMaxScaler
          

          
def normalize_dataframe(train_df, test_df):
          
    scaler = MinMaxScaler()
          
    scaler.fit(train_df)  # 在训练集上拟合归一化模型
          
    train_data = pd.DataFrame(scaler.transform(train_df), columns=train_df.columns, index = df_train.index)
          
    test_data = pd.DataFrame(scaler.transform(test_df), columns=test_df.columns, index = df_test.index)
          
    return train_data, test_data
          

          
data_train, data_test = normalize_dataframe(df_train, df_test)
          
data_train
      

在这里归一化时只使用训练集的统计量，并将归一化后的转换应用于训练集和测试集，避免直接对所有数据集进行归一化处理从而产生信息泄露

2.6 滑动窗口定义

2.6.1 滑动窗口不包含待预测特征

          
import numpy as np
          
def prepare_data(data, win_size, target_feature_idx, exclude_features=[]):
          
    num_features = data.shape[1] - len(exclude_features)  # 更新特征数量
          
    X = []  # 用于存储输入特征的列表
          
    y = []  # 用于存储目标值的列表
          

          
    # 遍历数据，形成样本窗口
          
    for i in range(len(data) - win_size):
          
        temp_x = []
          
        for j in range(data.shape[1]):
          
            if j != target_feature_idx and j not in exclude_features:
          
                temp_x.append(data[i:i + win_size, j])  # 提取窗口大小范围内的输入特征
          
        temp_y = data[i + win_size, target_feature_idx]  # 提取对应的目标值
          
        X.append(temp_x)
          
        y.append(temp_y)
          

          
    # 转换列表为 NumPy 数组，并调整维度顺序
          
    X = np.asarray(X).transpose(0, 2, 1)  
          
    y = np.asarray(y)
          

          
    return X, y
          

          
win_size = 12
          
target_feature_idx = 0 # 指定待预测特征
          
exclude_features = [0] #  需要删除的自变量特征也就是不把待预测的特征纳入输入特征进行时间窗口划分
          
train_x, train_y = prepare_data(data_train.values, win_size, target_feature_idx)
          
test_x, test_y = prepare_data(data_test.values, win_size, target_feature_idx)
          
print("训练集形状:", train_x.shape, train_y.shape)
          
print("测试集形状:", test_x.shape, test_y.shape)
      

在这里模型的训练集的输入特征为3个、时间窗口为12、存在1449条数据，输入特征分别为humidity、wind_speed、meanpressure, 输出特征meantemp

2.6.2 滑动窗口包含待预测特征

          
def prepare_data(data, win_size, target_feature_idx):
          
    num_features = data.shape[1]
          
    X = []  
          
    y = []  
          
    for i in range(len(data) - win_size):
          
        temp_x = data[i:i + win_size, :]  
          
        temp_y = data[i + win_size, target_feature_idx]  
          
        X.append(temp_x)
          
        y.append(temp_y)
          
    X = np.asarray(X)
          
    y = np.asarray(y)
          

          
    return X, y
          

          
win_size = 12 # 时间窗口
          
target_feature_idx = 0 # 指定待预测特征列
          
train_x, train_y = prepare_data(data_train.values, win_size, target_feature_idx)
          
test_x, test_y = prepare_data(data_test.values, win_size, target_feature_idx)
          
print("训练集形状:", train_x.shape, train_y.shape)
          
print("测试集形状:", test_x.shape, test_y.shape)
      

在这里模型的训练集的输入特征为4个、时间窗口为12、存在1449条数据，输入特征分别为*meantemp、humidity、wind_speed、meanpressure, 输出特征meantemp，*接下来将采用这种形式的滑动窗口构建LSTM模型

2.7 模型建立及评价

          
from keras.layers import LSTM, Dense
          
from keras.models import Sequential
          

          
# 模型构建
          
model = Sequential()
          
model.add(LSTM(128, activation='relu', input_shape=(train_x.shape[1], train_x.shape[2])))
          
model.add(Dense(64, activation='relu'))
          
model.add(Dense(32, activation='relu'))
          
model.add(Dense(1))
          

          
# 编译模型
          
model.compile(loss='mse', optimizer = 'adam')
          

          
# 模型拟合
          
history = model.fit(train_x, train_y, epochs=100, batch_size=32, validation_data=(test_x, test_y))
          

          
plt.figure()
          
plt.plot(history.history['loss'], c='b', label = 'loss')
          
plt.plot(history.history['val_loss'], c='g', label = 'val_loss')
          
plt.legend()
          
plt.show()
      

        
            

          model.summary()
        
      

          
from sklearn import metrics
          
y_pred = model.predict(test_x)
          
# 计算均方误差（MSE）
          
mse = metrics.mean_squared_error(test_y, np.array([i for arr in y_pred for i in arr]))
          
# 计算均方根误差（RMSE）
          
rmse = np.sqrt(mse)
          
# 计算平均绝对误差（MAE）
          
mae = metrics.mean_absolute_error(test_y, np.array([i for arr in y_pred for i in arr]))
          
from sklearn.metrics import r2_score # 拟合优度
          
r2 = r2_score(test_y, np.array([i for arr in y_pred for i in arr]))
          
print("均方误差 (MSE):", mse)
          
print("均方根误差 (RMSE):", rmse)
          
print("平均绝对误差 (MAE):", mae)
          
print("拟合优度:", r2)
      

2.8 向后预测并可视化

          
def predict_future(model, initial_sequence, steps):
          
    predicted_values = []  # 存储预测结果
          
    current_sequence = initial_sequence.copy()  # 初始序列
          
    for i in range(steps):
          
        # 使用模型进行单步预测
          
        predicted_value = model.predict(current_sequence.reshape(1, initial_sequence.shape[0], initial_sequence.shape[1]))
          
        # 将预测结果添加到列表中
          
        predicted_values.append(predicted_value[0][0])
          
        # 更新当前序列，删除第一个时间步并将预测值添加到最后一个时间步
          
        current_sequence[:-1] = current_sequence[1:]
          
        current_sequence[-1] = predicted_value
          
    return predicted_values
          

          
# 使用该函数进行预测
          
steps_to_predict =  11 # 要预测的步数
          
predicted_values = predict_future(model, test_x[-1], steps_to_predict)
          

          
train_max = np.max(df_train['meantemp'])
          
train_min = np.min(df_train['meantemp'])
          

          
series = np.array(predicted_values)*(train_max-train_min)+train_min
          
plt.figure(figsize=(15,4), dpi =300)
          
plt.subplot(2, 1, 1)
          
plt.plot(pd.date_range(start='2013-01-13', end='2016-12-31', freq='D')
          
         ,model.predict(train_x)*(train_max-train_min)+train_min, color = 'c', label = '训练集meantemp')
          
plt.plot(pd.date_range(start='2017-01-13', end='2017-04-24', freq='D')
          
         ,test_y*(train_max-train_min)+train_min, color = 'y', label = '测试集meantemp')
          
plt.plot(pd.date_range(start='2017-01-13', end='2017-04-24', freq='D')
          
         ,y_pred*(train_max-train_min)+train_min, color = 'r', label = '测试集预测meantemp')
          
plt.plot(pd.date_range(start='2017-04-24', end='2017-05-4', freq='D'), series, color = 'b', label = '向后预测10天')
          
plt.title('meantemp')
          
plt.grid(True)
          
plt.xlabel('time')
          
plt.ylabel('°C')
          
plt.legend()
          

          
plt.subplot(2, 1, 2)
          
plt.plot(pd.date_range(start='2017-01-13', end='2017-04-24', freq='D')
          
         ,test_y*(train_max-train_min)+train_min, color = 'y', label = '测试集meantemp')
          
plt.plot(pd.date_range(start='2017-01-13', end='2017-04-24', freq='D')
          
         ,y_pred*(train_max-train_min)+train_min, color = 'r', label = '测试集预测meantemp')
          
plt.plot(pd.date_range(start='2017-04-24', end='2017-05-4', freq='D'), series, color = 'b', label = '向后预测10天')
          

          
plt.grid(True)
          
plt.xlabel('time')
          
plt.ylabel('°C')
          
plt.legend()
          
plt.show()
      

3.**** 往期推荐

TCN时间序列卷积神经网络

基于VMD分解的VMD-CNN-LSTM时间序列预测模型实现

基于VMD分解的VMD-LSTM时间序列预测模型实现，大力提升预测精度！

ARIMA模型进阶——具有季节性的SARIMAX模型实现

时间序列——移动平均（SMA、EMA、WMA）

时间序列——平稳性检验

时间序列——季节性检验

长短期记忆网络LSTM在时序数据预测演示

如果你对类似于这样的文章感兴趣。

欢迎关注、点赞、转发~