时间序列预测（二）基于LSTM的销售额预测

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import seaborn as sns
from fbprophet import Prophet
from sklearn.metrics import mean_squared_error
from math import sqrt
import datetime
from xgboost import XGBRegressor
from sklearn.metrics import explained_variance_score, mean_absolute_error, \
mean_squared_error, r2_score  # 批量导入指标算法

from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM
from tensorflow.keras.layers import Dense, Dropout
from sklearn.preprocessing import MinMaxScaler
from keras.wrappers.scikit_learn import KerasRegressor
from sklearn.model_selection import GridSearchCV

# 读取数据
raw_data = pd.read_csv('train.csv')
raw_data.set_index('datetime', inplace=True)
raw_data.head()

# 拆分数据集
num = 24*14 # 将最后2周划分为测试集
train, test = raw_data.iloc[:-num,:], raw_data.iloc[-num:,:]

# 数据缩放
scaler = MinMaxScaler(feature_range=(0,1))
train_scaled=scaler.fit_transform(train)
test_scaled=scaler.transform(test)

# 构造XY(样本数+时间步数+特征数)
def createXY(dataset, n_past, target_p=-1):
    '''
    将数据集构造成LSTM需要的格式
    dataset：数据集
    n_past：时间步数，利用过去n的时间作为特征,以下一个时间的目标值作为当前的y
    target_p：目标值在数据集的位置，默认为-1
    '''
    dataX = []
    dataY = []
    for i in range(n_past, len(dataset)):
        dataX.append(dataset[i - n_past:i, 0:dataset.shape[1]])
        dataY.append(dataset[i,target_p])
    return np.array(dataX),np.array(dataY)

X_train, Y_train=createXY(train_scaled,30)
X_test, Y_test=createXY(test_scaled,30)
X_train.shape

(10520, 30, 11)

# 定义LSTM
def build_model(optimizer):
    grid_model = Sequential()
    grid_model.add(LSTM(4, return_sequences=True, input_shape=(X_train.shape[1], X_train.shape[2])))
    grid_model.add(LSTM(4)) # 防止预测值为三维
    grid_model.add(Dropout(0.2))
    grid_model.add(Dense(1))
    grid_model.compile(loss = 'mse',optimizer = optimizer)
    return grid_model 

grid_model=KerasRegressor(build_fn=build_model,verbose=1,validation_data=(X_test,Y_test))
 
parameters = {'batch_size' : [16,20],
            'epochs' : [8,10],
            'optimizer' : ['adam','Adadelta'] }
 
grid_search = GridSearchCV(estimator = grid_model,
                            param_grid = parameters,
                            cv = 2)

<ipython-input-41-930ab96e3919>:11: DeprecationWarning: KerasRegressor is deprecated, use Sci-Keras (https://github.com/adriangb/scikeras) instead. See https://www.adriangb.com/scikeras/stable/migration.html for help migrating.
  grid_model=KerasRegressor(build_fn=build_model,verbose=1,validation_data=(X_test,Y_test))

# 模型拟合
grid_search = grid_search.fit(X_train,Y_train)
grid_search.best_params_

# 输出最优模型参数
print(grid_search.best_params_)
# 获取最优模型
model_lstm=grid_search.best_estimator_.model
# 预测值
pre_y=model_lstm.predict(X_test)

{'batch_size': 16, 'epochs': 8, 'optimizer': 'adam'}

# 逆缩放
# 构造同等宽度的pre_y，即生成与缩放同等列数
pre_y_repeat = np.repeat(pre_y, X_train.shape[2], axis=-1)
# 逆缩放并获取预测值
pred=scaler.inverse_transform(np.reshape(pre_y_repeat,(len(pre_y), X_train.shape[2])))[:,-1] # 选取目标列所在位置

# 对Y_test进行逆缩放
Y_test_repeat = np.repeat(Y_test, X_train.shape[2], axis=-1)
Y_test_original=scaler.inverse_transform(np.reshape(Y_test_repeat,(len(Y_test),X_train.shape[2])))[:,-1] # 选取目标列所在位置

# 评估指标
model_metrics_functions = [explained_variance_score, mean_absolute_error, mean_squared_error,r2_score]  # 回归评估指标对象集
model_metrics_list = [[m(Y_test_original, pred) for m in model_metrics_functions]]  # 回归评估指标列表
regresstion_score = pd.DataFrame(model_metrics_list, index=['model_xgbr'],
                   columns=['explained_variance', 'mae', 'mse', 'r2'])  # 建立回归指标的数据框
regresstion_score  # 模型回归指标

# 模型效果可视化
fig = plt.figure(figsize=(16,6))
plt.title('True and LSTM result comparison', fontsize=20)
# 时间索引
ds_index = [datetime.datetime.strptime(x, "%Y-%m-%d %H:%M:%S") for x in test.index[30:]]
# 真实序列
true_s=pd.Series(Y_test_original, index=ds_index)
plt.plot(true_s, color='red')
# 预测序列
pre_s=pd.Series(pred, index=ds_index)
plt.plot(pre_s, color='green')
plt.xlabel('Hour', fontsize=16)
plt.ylabel('Number of Shared Bikes', fontsize=16)
plt.legend(labels=['True', 'Pre_y'], fontsize=16)
plt.grid()
plt.show()

# 预测未来值
# 历史30日数据作为构造第一条数据
df_30_days_past=raw_data.iloc[-30:,:]
# 读取未来数据
start_time = '2012-12-19 23:00:00' # 原始数据最后一条记录
end_time=[]
for i in range(30):
    
    time_temp = (datetime.datetime.strptime(start_time, "%Y-%m-%d %H:%M:%S") + datetime.timedelta(
        hours=i+1)).strftime("%Y-%m-%d %H:%M:%S")
    end_time.append(time_temp)

future=pd.read_csv("test.csv",parse_dates=["datetime"],index_col=[0])
df_30_days_future=future.loc[end_time,:]

# 补充缺失的列
df_30_days_future["registered"]=0
df_30_days_future["casual"]=0 # 测试样本缺失
df_30_days_future["count"]=0 # 测试样本缺失

# 构造旧新样本
old_scaled_array=scaler.transform(df_30_days_past)
new_scaled_array=scaler.transform(df_30_days_future)
new_scaled_df=pd.DataFrame(new_scaled_array)
new_scaled_df.iloc[:,-1]=np.nan # 对目标列的0修改为nan
full_df=pd.concat([pd.DataFrame(old_scaled_array),new_scaled_df]).reset_index().drop(["index"],axis=1)

full_df_scaled_array=full_df.values
all_data=[]
time_step=30
for i in range(time_step,len(full_df_scaled_array)):
    data_x=[]
    data_x.append(
        full_df_scaled_array[i-time_step :i , 0:full_df_scaled_array.shape[1]])
    data_x=np.array(data_x)
    prediction=model_lstm.predict(data_x)
    all_data.append(prediction)
    full_df.iloc[i,-1]=prediction # 替换目标列的nan

new_array=np.array(all_data)
new_array=new_array.reshape(-1,1)
prediction_copies_array = np.repeat(new_array,data_x.shape[2], axis=-1)
y_pred_future_30_days = scaler.inverse_transform(np.reshape(prediction_copies_array,(len(new_array),data_x.shape[2])))[:,-1]
print(y_pred_future_30_days)

# 因为测试样本casual和registered缺失，故预测不准确

[ 67.19963   28.726665  53.639095  77.41373   84.80948   97.550125
 105.34982  105.89546  105.45314  101.901405  98.7133    96.72031
 101.96429  107.16215  110.38648  111.0147   110.34529  109.8532
 107.00935  102.02435   91.798004  92.968605  93.628006  94.56274
 102.24291  113.515884  87.92808   78.71299   72.63906   70.64108 ]