其中df_train和df_test分別為模型的訓練集和測試集，數據不存在缺失值。
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時序圖中數據存在異常值。
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()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）的濾波器，旨在識別和去除時間序列數據中的異常值。
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.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模型。
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)
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()