Rui - LSTM预测股价

回归问题实例：LSTM预测股价

RNN 在处理长期依赖（时间序列上距离较远的节点）时会遇到巨大的困难，因为计算距离较远的节点之间的联系时会涉及雅可比矩阵的多次相乘，会造成梯度消失或者梯度膨胀的现象。为了改善 RNN 的长期依赖问题，一个非常好的解决方式是引入门控机制来控制信息的累积速度，包括有选择地加入新的信息，并有选择地遗忘之前累积的信息。这一类网络可以称为基于门控的循环神经网络（Gate RNN），而 LSTM 就是 Gate RNN 中最著名的一种。

LSTM 使用说明

Keras 提供的 layer_lstm() 函数可以实现长短期记忆网络（LSTM），该函数表达形式为：

Code

layer_lstm(
  object, 
  units, 
  activation = "tanh",
  recurrent_activation = "hard_sigmoid", 
  use_bias = TRUE,
  return_sequences = FALSE, 
  return_state = FALSE,
  go_backwards = FALSE, 
  stateful = FALSE, 
  unroll = FALSE,
  kernel_initializer = "glorot_uniform",
  recurrent_initializer = "orthogonal", 
  bias_initializer = "zeros",
  unit_forget_bias = TRUE, 
  kernel_regularizer = NULL,
  recurrent_regularizer = NULL, 
  bias_regularizer = NULL,
  activity_regularizer = NULL, 
  kernel_constraint = NULL,
  recurrent_constraint = NULL, 
  bias_constraint = NULL, 
  dropout = 0,
  recurrent_dropout = 0, 
  input_shape = NULL,
  batch_input_shape = NULL, 
  batch_size = NULL, 
  dtype = NULL,
  name = NULL, 
  trainable = NULL, 
  weights = NULL
)

各参数描述如下：

object：模型或图层对象。
units：输出维度。
activation：激活函数，默认为双曲正切（tanh），如果传入 NULL，则不会使用任何激活函数（即线性激活a（x）=x）。
recurrent_activation：用于循环步骤的激活功能。
use_bias：布尔值（默认为TRUE），该层是否使用偏向量。
return_sequences：布尔值（默认为FALSE），若为FALSE，返回最后一层最后一个步长的隐藏状态；若为 TRUE，返回最后一层的所有隐藏状态。也就是说，当值为 TRUE 时，返回的是 c(samples, time_steps, output_dim) 的 3D 张量，当值为 FALSE 时，返回的是 (samples, output_dim)的 2D 张量。
return_state：布尔值（默认为 FALSE），若为 TRUE，除了返回输出值之外，还要返回最后一个单元状态。
go_backwards：布尔值（默认为 FALSE），如果为 TRUE，返回最后一层的最后一个步长的输出隐藏状态和输入单元状态。
stateful：布尔值（默认为 FALSE），如果为 TRUE，则批次中索引 i 处的每个样本的最后状态将用作后续批次中索引 i 的样本的初始状态。
unroll：布尔值（默认为 FALSE），如果为 TRUE，则将展开网络，否则将使用符号循环。网络展开可以加速 RNN，但往往会占用大量内存。展开仅适用于短序列。
kernel_initializer：kernel 权重矩阵的初始化器，用于输入的线性转换，默认为 glorot_uniform。
recurrent_initializer：recurrent_kernel 权重矩阵的初始化程序，用于循环状态的线性转换。
bias_initializer：偏置向量的初始化器。
kernel_regularizer：应用于 kernel 权重矩阵的正则化函数。
bias_regularizer：应用于偏置向量的正则化函数。
recurrent_regularizer：应用于 recurrent_kernel 权重矩阵的正则化函数。
activity_regularizer：应用于图层输出（它的激活值）的正则化函数。
kernel_constraints：应用于 kernel 权重矩阵的约束函数。
recurrent_constraints：应用于 recurrent_kernel 权重矩阵的约束函数。
bias_constraints：应用于偏置向量的约束函数。
dropout：0～1之间的浮点数，控制输入线性变换的神经元丢弃比例。

在进行多层 LSTM 时，需要注意以下两点：

需要对第一层的 LSTM 指定 input_shape 参数。
将前 N-1 层 LSTM 的 return_sequences 设置为 TRUE，保证每一层都会向下传播所有时间步长上的预测，同时保证最后一层的 return_sequences 设置为 FALSE。

股份数据处理

股价数据是非常常见的时序数据之一，本节将使用 LSTM 对股价进行预测。

导入数据

google_stock_prics.csv 文件记录了 Google 从 2010 年 1 月至 2018 年 5 月每日（非节假日）股价数据，共 2120 条记录。每条记录包含了当日的开盘价、最高价、最低价、收盘价、调整收盘价和成交量的信息。我们先将数据读入 R 中，并查看前 6 条记录。

Code

library(keras)
library(tictoc)

# 读入Google股价数据
googl <- read.csv("data/google_stock_prics.csv")
head(googl)
##        Date     Open     High      Low    Close Adj.Close   Volume
## 1  2010/1/4 313.7888 315.0701 312.4324 313.6887  313.6887  3908400
## 2  2010/1/5 313.9039 314.2342 311.0811 312.3073  312.3073  6003300
## 3  2010/1/6 313.2433 313.2433 303.4835 304.4344  304.4344  7949400
## 4  2010/1/7 305.0050 305.3053 296.6216 297.3474  297.3474 12815700
## 5  2010/1/8 296.2963 301.9269 294.8499 301.3113  301.3113  9439100
## 6 2010/1/11 302.5325 302.5325 297.3173 300.8559  300.8559 14411300

Code

plot(googl[ , "Close"], 
     col = "black", 
     type = "l", 
     lwd = 2,
     xlab = "Time", 
     ylab = "Google Stock Price",
     main = "Google Stock Price")

拆分训练集测试集

我们将数据集划分为两部分：前1760条记录作为训练集，用于模型训练，后360条记录作为测试集，用于模型验证。由于希望利用历史收盘价格（Close）的数据来预测未来收盘价格，所以在做数据分区时只提取收盘价格。

Code

train <- googl[1:1760, "Close"]
test <- googl[1761:2120, "Close"]

定义数据处理函数

Code

# 自定义归一化及反归一化函数
normalize <- function(vec, min, max) {
  (vec-min) / (max-min)
}
denormalize <- function(vec, min, max) {
  vec * (max-min) + min
}

Code

# 定义一个将数据转换为“时间步长形式”的单个矩阵
build_matrix <- function(tseries, overall_timesteps) {
  X <- t(sapply(1:(length(tseries) - overall_timesteps + 1), #!!!!!!!!!! +1
                function(x) tseries[x:(x + overall_timesteps - 1)]))
  
  cat("\nBuilt matrix with dimensions: ", dim(X))
  
  return(X)
}

Note

sapply() 函数的作用是：将列表、向量或数据帧作为输入，对其中每一个元素使用函数，并以向量或矩阵的形式输出结果。

Code

build_X <- function(tseries, num_timesteps) {
  X <- if (num_timesteps > 1) {
    t(sapply(1:(length(tseries) - num_timesteps),
             function(x) tseries[x:(x + num_timesteps - 1)]))
  } else {
    tseries[1:length(tseries) - num_timesteps]
  }
  
  if (num_timesteps == 1) dim(X) <- c(length(X), 1)
  
  cat("\nBuilt X matrix with dimensions: ", dim(X))
  
  return(X)
}

Code

build_y <- function(tseries, num_timesteps) {
  y <- sapply((num_timesteps + 1):(length(tseries)), function(x) tseries[x])
  
  cat("\nBuilt y vector with length: ", length(y))
  
  return(y)
}

Tip

可以将以上自定义函数单独放在一个 R 的脚本中，在必要时通过 source() 函数调用。例如 source("functions.R")。

处理数据

利用 normalize() 函数对训练集和测试集进行标准化处理，记住测试集需要用训练集的最小值和最大值进行转换。

Code

minval <- min(train)
maxval <- max(train)
train_normalize <- normalize(train, minval, maxval)
test_normalize <- normalize(test, minval, maxval)

我们想利用最近 60 天的数据来预测下一天的收盘股价，故将 build_X()、build_y() 函数 num_timesteps 参数设置为 60，构建 X、y。

Code

num_timesteps <- 60

X_train <- build_X(train_normalize, num_timesteps)
## 
## Built X matrix with dimensions:  1700 60
y_train <- build_y(train_normalize, num_timesteps)
## 
## Built y vector with length:  1700

X_test <- build_X(test_normalize, num_timesteps)
## 
## Built X matrix with dimensions:  300 60
y_test <- build_y(test_normalize, num_timesteps)
## 
## Built y vector with length:  300

因为LSTM要求输入矩阵的形状为 (sample, timesteps, input_dim)，所以需要对 X_train、X_test 进行形状处理。

Code

X_train <- array_reshape(X_train, c(dim(X_train), 1))
X_test <- array_reshape(X_test, c(dim(X_test), 1))
dim(X_train)
## [1] 1700   60    1
dim(X_test)
## [1] 300  60   1

至此，数据预处理工作已经完成。

无状态 LSTM

接下来，我们构建一个两层无状态 LSTM 网络。两层 LSTM 的神经元数量均为 50，其中希望把第一层 LSTM 的隐藏层状态全部返回给下一层 LSTM，故将第一层的参数 return_sequences 设置为 TRUE。两层 LSTM 后均带有一个 Dropout 层，以防止过拟合。最后接一个密集层作为输出层，神经元数量为 1。因为是预测模型，所以不需要指定激活函数。编译模型时，采用 Adam 优化器，mean_squared_error 作为损失函数。

Code

tic("构建无状态 LSTM ")

model <- keras_model_sequential()
model %>% 
  layer_lstm(units = 50, return_sequences = TRUE, input_shape = dim(X_train)[-1]) %>%
  layer_dropout(rate = 0.2) %>%
  layer_lstm(units = 50) %>%
  layer_dropout(rate = 0.2) %>%
  layer_dense(units = 1) %>%
  compile(loss = "mean_squared_error", optimizer = "adam")

toc()
## 构建无状态 LSTM : 110.36 sec elapsed

无状态 LSTM 在输入样本后，默认会被打乱。也就是说，每个样本独立，样本之间无前后关系，适合输入一些没有关系的样本。在训练模型时，fit() 函数的参数 shuffle 默认为 TRUE，即 Keras 在训练时会默认打乱样本，导致序列之间的依赖性消失，样本和样本之间没有了时序关系，顺序被打乱。这时记忆参数在批量、小序列之间的传递就没有意义了，所以 Keras 要对记忆参数进行初始化。我们在训练过程中，设置一个回调函数，如果经过 5 次训练周期后验证集的损失函数没有明显改善，则自动停止训练。

Code

tic("训练模型 ")

history <- model %>%
  fit(X_train, 
      y_train, 
      batch_size = 50, 
      epochs = 50,
      validation_split = 0.2, 
      verbose = 2,
      callbacks = callback_early_stopping(patience = 5))

toc()
## 训练模型 : 39.17 sec elapsed

Code

plot(history)

经过 16 个周期后，模型停止了训练。利用 save_model_hdf5() 函数将训练好的模型保存到本地。

Code

model %>% save_model_hdf5("model/lstm_1v1_stateless.h5")

使用 predict() 函数对训练集和测试集进行预测。由于预测出来的结果是标准化后的数据，所以需要利用 denormalize() 函数进行转换。

Code

pred_train <- model %>% predict(X_train)
pred_test <- model %>% predict(X_test)

pred_train <- denormalize(pred_train, minval, maxval)
pred_test <- denormalize(pred_test, minval, maxval)

利用以下命令绘制训练集样本的实际值和预测值曲线。

Code

plot(train[61:length(train)], 
     col = "red", 
     type = "l", 
     lwd = 2,
     xlab = "Time", 
     ylab = "Google Stock Price",
     main = "Google Stock Price Prediction")
lines(pred_train, col = "blue", type = "l", lty = 3, lwd = 2)
legend("topleft", 
       legend = c("real_stock_price ","predicted_stock_price"),
       lty = c(1, 2),
       lwd = 2,
       col = c("red","blue"),
       bty = "n")

预测曲线与实际曲线基本保持一致，说明 LSTM 对训练集的拟合效果非常好。

利用以下命令绘制测试集样本的实际值和预测值曲线。

Code

plot(test[61:length(test)],
     col = "red",
     type = "l",
     lwd = 2,
     xlab = "Time",
     ylab = "Google Stock Price",
     main = "Google Stock Price Prediction")
lines(pred_test, col = "blue", type = "l", lty = 2, lwd = 2)
legend("topleft",
       legend = c("real_stock_price ","predicted_stock_price"),
       lty = c(1, 2),
       lwd = 2,
       col = c("red","blue"),
       bty = "n")

下面计算在测试集上的均方误差根。

Code

mean(test[61:length(test)])
## [1] 1012.324
(n <- length(pred_test))
## [1] 300
mse <- sum((test[61:length(test)]-pred_test)^2)/n
rmse <- sqrt(mse)
rmse
## [1] 60.26064

在平均股价为 1012 的 300 个测试样本里，无状态 LSTM 在测试集上的均方误差根为 60.2606367。

有状态 LSTM

有状态（stateful）LSTM 能够在训练中维护跨批次的有状态信息，即当前批次的训练数据计算的状态值，可以用作下一批次训练数据的初始隐藏状态。有状态 LSTM 能让模型学习到输入样本之间的时序特征，适合一些长序列的预测，此时样本的前后顺序对模型是有影响的。stateful 代表除了每个样本内的时间步内传递，每个样本之间还会有信息 (c, h) 传递。

当使用有状态 LSTM 时，需假定以下两种情况：

所有的批次都有相同数量的样本，即训练集和测试集的样本数量均能被 batch_size 整除。
如果 x1 和 x2 是连续批次的样本，则对于每一个 $i\in timesteps$，x2[i] 是 x1[i] 的后续序列。

有状态LSTM的实现步骤如下：

必须将参数 batch_size 显式地传递给模型的第一层。
在 LSTM 层中将 stateful 设置为 TRUE。
在调用 fit() 函数时需指定 shuffle 为 FALSE，因为打乱样本后，序列之间就没有依赖性了。
训练完一个周期后，需要使用 reset_state() 函数来重置状态。

接下来，将上一小节的无状态 LSTM 调整为有状态 LSTM，并使用训练集数据进行训练。 在构建和拟合模型时的关键代码如下：

将 batch_size 设置为 50，则 1700（训练样本量）/50=34，300（测试样本量）/50=6。
构建模型时，将第一层 LSTM 网络参数 return_sequences 设置为 TRUE，stateful 设置为 TRUE，batch_input_shape 设置为 c(batch_size, 60, 1))；第二层 LSTM 网络的 stateful 设置为 TRUE。
使用 for 循环进行模型训练，需指定 fit() 函数的 batch_size，且将循环周期次数设置为 1，shuffle 设置为 FALSE；在 fit() 函数后使用 reset_states() 函数重置模型状态。

以下是构建和训练有状态 LSTM 网络的程序代码。

Code

tic("构建模型 ")

batch_size <- 50
num_epochs <- 50

# 构建模型
# 期望输入数据尺寸: (batch_size, timesteps, input_dim)
stateful_model <- keras_model_sequential()
stateful_model %>%
  layer_lstm(units = 50, return_sequences = TRUE, stateful = TRUE,
             batch_input_shape = c(batch_size, 60, 1)) %>%
  layer_dropout(rate = 0.2) %>%
  layer_lstm(units = 50, stateful = TRUE) %>%
  layer_dropout(rate = 0.2) %>%
  layer_dense(units = 1) %>%
  compile(loss = "mean_squared_error", optimizer = "adam")

toc()
## 构建模型 : 0.78 sec elapsed

Code

tic("训练模型 ")

# 训练模型
for (i in seq(num_epochs)){
  stateful_model %>% 
    fit(
      X_train, 
      y_train, 
      batch_size = batch_size, 
      epochs = 1,
      shuffle = FALSE
    )
  
  stateful_model %>% reset_states()
}

toc()
## 训练模型 : 56.4 sec elapsed

模型训练好后，使用 save_model_hdf5() 函数将模型保存到本地。

Code

stateful_model %>% save_model_hdf5("model/lstm_stateful_model.h5")

接下来，我们对测试集和训练集进行预测，进行反归一化处理，并分别绘制训练集和测试集的实际值和预测值曲线。

Code

# 模型预测
pred_train_stateful <- stateful_model %>% predict(X_train, batch_size = batch_size)
pred_test_stateful <- stateful_model %>% predict(X_test, batch_size = batch_size)

pred_train_stateful <- denormalize(pred_train_stateful, minval, maxval)
pred_test_stateful <- denormalize(pred_test_stateful, minval, maxval)

Code

par(mfrow = c(1, 2))
# 绘制训练集的实际值与预测值曲线
plot(train[61:length(train)], 
     col = "red", 
     type = "l", lwd = 2,
     xlab = "Time", ylab = "Google Stock Price", 
     ylim = c(200, 1200),
     main = "Train Data Prediction")
lines(pred_train_stateful, col = "blue", type = "l", lty = 3, lwd = 2)
legend("topleft", 
       legend = c("real ", "predicted"),
       lty = c(1, 2), lwd = 2, 
       horiz = FALSE,
       col = c("red", "blue"), 
       bty = "n")

# 绘制测试集的实际值与预测值曲线
plot(test[61:length(test)], 
     col = "red", 
     type = "l", lwd = 2,
     xlab = "Time", ylab = "Google Stock Price", 
     ylim = c(200, 1200),
     main = "Test Data Prediction") 
lines(pred_test_stateful, col = "blue", type = "l", lty = 2, lwd = 2)
legend("topleft", 
       legend = c("real ", "predicted"),
       lty = c(1, 2), 
       lwd = 2, 
       horiz = FALSE,
       col = c("red", "blue"), 
       bty = "n") # 去除图例外部方框

Code

par(mfrow = c(1, 1))

下面计算有状态 LSTM 在测试集上的均方误差根。

Code

mean(test[61:length(test)])
## [1] 1012.324
(n <- length(pred_test_stateful))
## [1] 300
mse <- sum((test[61:length(test)]-pred_test_stateful)^2)/n
rmse <- sqrt(mse)
rmse
## [1] 93.55118

在平均股价为 1012 的 300 个测试样本里，有状态 LSTM 在测试集上的均方误差根为 93.5511821。

--- title: "LSTM预测股价" author: "Rui" date: "2023-05-14" categories: [R, Keras] image: "FinancialTowers.jpg" format: html: code-fold: true code-tools: true --- ```{r setup, include = FALSE} # 设置默认参数 knitr::opts_chunk$set( echo = TRUE, fig.align = "center", message = FALSE, warning = FALSE, collapse = TRUE ) ``` ```{=html} <iframe frameborder="no" border="0" marginwidth="0" marginheight="0" width=330 height=86 src="//music.163.com/outchain/player?type=2&id=1931193123&auto=1&height=66"></iframe> ``` ## 回归问题实例：LSTM预测股价 RNN 在处理长期依赖（时间序列上距离较远的节点）时会遇到巨大的困难，因为计算距离较远的节点之间的联系时会涉及雅可比矩阵的多次相乘，会造成梯度消失或者梯度膨胀的现象。为了改善 RNN 的长期依赖问题，一个非常好的解决方式是引入门控机制来控制信息的累积速度，包括有选择地加入新的信息，并有选择地遗忘之前累积的信息。这一类网络可以称为基于门控的循环神经网络（Gate RNN），而 LSTM 就是 Gate RNN 中最著名的一种。 ### LSTM 使用说明 Keras 提供的 `layer_lstm()` 函数可以实现长短期记忆网络（LSTM），该函数表达形式为： ```{r eval=FALSE} layer_lstm( object, units, activation = "tanh", recurrent_activation = "hard_sigmoid", use_bias = TRUE, return_sequences = FALSE, return_state = FALSE, go_backwards = FALSE, stateful = FALSE, unroll = FALSE, kernel_initializer = "glorot_uniform", recurrent_initializer = "orthogonal", bias_initializer = "zeros", unit_forget_bias = TRUE, kernel_regularizer = NULL, recurrent_regularizer = NULL, bias_regularizer = NULL, activity_regularizer = NULL, kernel_constraint = NULL, recurrent_constraint = NULL, bias_constraint = NULL, dropout = 0, recurrent_dropout = 0, input_shape = NULL, batch_input_shape = NULL, batch_size = NULL, dtype = NULL, name = NULL, trainable = NULL, weights = NULL ) ``` **各参数描述如下：** - `object`：模型或图层对象。 - `units`：输出维度。 - `activation`：激活函数，默认为双曲正切（tanh），如果传入 NULL，则不会使用任何激活函数（即线性激活a（x）=x）。 - `recurrent_activation`：用于循环步骤的激活功能。 - `use_bias`：布尔值（默认为TRUE），该层是否使用偏向量。 - `return_sequences`：布尔值（默认为FALSE），若为FALSE，返回最后一层最后一个步长的隐藏状态；若为 TRUE，返回最后一层的所有隐藏状态。也就是说，当值为 TRUE 时，返回的是 `c(samples, time_steps, output_dim)` 的 3D 张量，当值为 FALSE 时，返回的是 `(samples, output_dim)`的 2D 张量。 - `return_state`：布尔值（默认为 FALSE），若为 TRUE，除了返回输出值之外，还要返回最后一个单元状态。 - `go_backwards`：布尔值（默认为 FALSE），如果为 TRUE，返回最后一层的最后一个步长的输出隐藏状态和输入单元状态。 - `stateful`：布尔值（默认为 FALSE），如果为 TRUE，则批次中索引 `i` 处的每个样本的最后状态将用作后续批次中索引 `i` 的样本的初始状态。 - `unroll`：布尔值（默认为 FALSE），如果为 TRUE，则将展开网络，否则将使用符号循环。网络展开可以加速 RNN，但往往会占用大量内存。展开仅适用于短序列。 - `kernel_initializer`：`kernel` 权重矩阵的初始化器，用于输入的线性转换，默认为 `glorot_uniform`。 - `recurrent_initializer`：`recurrent_kernel` 权重矩阵的初始化程序，用于循环状态的线性转换。 - `bias_initializer`：偏置向量的初始化器。 - `kernel_regularizer`：应用于 `kernel` 权重矩阵的正则化函数。 - `bias_regularizer`：应用于偏置向量的正则化函数。 - `recurrent_regularizer`：应用于 `recurrent_kernel` 权重矩阵的正则化函数。 - `activity_regularizer`：应用于图层输出（它的激活值）的正则化函数。 - `kernel_constraints`：应用于 `kernel` 权重矩阵的约束函数。 - `recurrent_constraints`：应用于 `recurrent_kernel` 权重矩阵的约束函数。 - `bias_constraints`：应用于偏置向量的约束函数。 - `dropout`：0～1之间的浮点数，控制输入线性变换的神经元丢弃比例。 **在进行多层 LSTM 时，需要注意以下两点：** - 需要对第一层的 LSTM 指定 `input_shape` 参数。 - 将前 `N-1` 层 LSTM 的 `return_sequences` 设置为 `TRUE`，保证每一层都会向下传播所有时间步长上的预测，同时保证最后一层的 `return_sequences` 设置为 `FALSE`。 ### 股份数据处理股价数据是非常常见的时序数据之一，本节将使用 LSTM 对股价进行预测。 #### 导入数据 google_stock_prics.csv 文件记录了 Google 从 2010 年 1 月至 2018 年 5 月每日（非节假日）股价数据，共 2120 条记录。每条记录包含了当日的开盘价、最高价、最低价、收盘价、调整收盘价和成交量的信息。我们先将数据读入 R 中，并查看前 6 条记录。 ```{r} library(keras) library(tictoc) # 读入Google股价数据 googl <- read.csv("data/google_stock_prics.csv") head(googl) ``` ```{r} plot(googl[ , "Close"], col = "black", type = "l", lwd = 2, xlab = "Time", ylab = "Google Stock Price", main = "Google Stock Price") ``` #### 拆分训练集测试集我们将数据集划分为两部分：前1760条记录作为训练集，用于模型训练，后360条记录作为测试集，用于模型验证。由于希望利用历史收盘价格（Close）的数据来预测未来收盘价格，所以在做数据分区时只提取收盘价格。 ```{r} train <- googl[1:1760, "Close"] test <- googl[1761:2120, "Close"] ``` #### 定义数据处理函数 ```{r} # 自定义归一化及反归一化函数 normalize <- function(vec, min, max) { (vec-min) / (max-min) } denormalize <- function(vec, min, max) { vec * (max-min) + min } ``` ```{r} # 定义一个将数据转换为“时间步长形式”的单个矩阵 build_matrix <- function(tseries, overall_timesteps) { X <- t(sapply(1:(length(tseries) - overall_timesteps + 1), #!!!!!!!!!! +1 function(x) tseries[x:(x + overall_timesteps - 1)])) cat("\nBuilt matrix with dimensions: ", dim(X)) return(X) } ``` ![](函数示意图.svg) ::: callout-note `sapply()` 函数的作用是：将列表、向量或数据帧作为输入，对其中每一个元素使用函数，并以向量或矩阵的形式输出结果。 ::: ```{r} build_X <- function(tseries, num_timesteps) { X <- if (num_timesteps > 1) { t(sapply(1:(length(tseries) - num_timesteps), function(x) tseries[x:(x + num_timesteps - 1)])) } else { tseries[1:length(tseries) - num_timesteps] } if (num_timesteps == 1) dim(X) <- c(length(X), 1) cat("\nBuilt X matrix with dimensions: ", dim(X)) return(X) } ``` ```{r} build_y <- function(tseries, num_timesteps) { y <- sapply((num_timesteps + 1):(length(tseries)), function(x) tseries[x]) cat("\nBuilt y vector with length: ", length(y)) return(y) } ``` ::: callout-tip 可以将以上自定义函数单独放在一个 R 的脚本中，在必要时通过 `source()` 函数调用。例如 `source("functions.R")`。 ::: #### 处理数据利用 `normalize()` 函数对训练集和测试集进行标准化处理，记住测试集需要用训练集的最小值和最大值进行转换。 ```{r} minval <- min(train) maxval <- max(train) train_normalize <- normalize(train, minval, maxval) test_normalize <- normalize(test, minval, maxval) ``` 我们想利用最近 60 天的数据来预测下一天的收盘股价，故将 `build_X()`、`build_y()` 函数 `num_timesteps` 参数设置为 60，构建 X、y。 ```{r} num_timesteps <- 60 X_train <- build_X(train_normalize, num_timesteps) y_train <- build_y(train_normalize, num_timesteps) X_test <- build_X(test_normalize, num_timesteps) y_test <- build_y(test_normalize, num_timesteps) ``` 因为LSTM要求输入矩阵的形状为 `(sample, timesteps, input_dim)`，所以需要对 X_train、X_test 进行形状处理。 ```{r} X_train <- array_reshape(X_train, c(dim(X_train), 1)) X_test <- array_reshape(X_test, c(dim(X_test), 1)) dim(X_train) dim(X_test) ``` 至此，数据预处理工作已经完成。 ### 无状态 LSTM 接下来，我们构建一个两层无状态 LSTM 网络。两层 LSTM 的神经元数量均为 50，其中希望把第一层 LSTM 的隐藏层状态全部返回给下一层 LSTM，故将第一层的参数 `return_sequences` 设置为 TRUE。两层 LSTM 后均带有一个 Dropout 层，以防止过拟合。最后接一个密集层作为输出层，神经元数量为 1。因为是预测模型，所以不需要指定激活函数。编译模型时，采用 `Adam` 优化器，`mean_squared_error` 作为损失函数。 ```{r} tic("构建无状态 LSTM ") model <- keras_model_sequential() model %>% layer_lstm(units = 50, return_sequences = TRUE, input_shape = dim(X_train)[-1]) %>% layer_dropout(rate = 0.2) %>% layer_lstm(units = 50) %>% layer_dropout(rate = 0.2) %>% layer_dense(units = 1) %>% compile(loss = "mean_squared_error", optimizer = "adam") toc() ``` ```{r eval=FALSE, include=FALSE} deepviz::plot_model(model) # 网络拓扑可视化 ``` 无状态 LSTM 在输入样本后，默认会被打乱。也就是说，每个样本独立，样本之间无前后关系，适合输入一些没有关系的样本。在训练模型时，`fit()` 函数的参数 `shuffle` 默认为 TRUE，即 Keras 在训练时会默认打乱样本，导致序列之间的依赖性消失，样本和样本之间没有了时序关系，顺序被打乱。这时记忆参数在批量、小序列之间的传递就没有意义了，所以 Keras 要对记忆参数进行初始化。我们在训练过程中，设置一个回调函数，如果经过 5 次训练周期后验证集的损失函数没有明显改善，则自动停止训练。 ```{r} tic("训练模型 ") history <- model %>% fit(X_train, y_train, batch_size = 50, epochs = 50, validation_split = 0.2, verbose = 2, callbacks = callback_early_stopping(patience = 5)) toc() ``` ```{r} plot(history) ``` 经过 `r length(history[["metrics"]][["val_loss"]])` 个周期后，模型停止了训练。利用 `save_model_hdf5()` 函数将训练好的模型保存到本地。 ```{r} model %>% save_model_hdf5("model/lstm_1v1_stateless.h5") ``` 使用 `predict()` 函数对训练集和测试集进行预测。由于预测出来的结果是标准化后的数据，所以需要利用 `denormalize()` 函数进行转换。 ```{r} pred_train <- model %>% predict(X_train) pred_test <- model %>% predict(X_test) pred_train <- denormalize(pred_train, minval, maxval) pred_test <- denormalize(pred_test, minval, maxval) ``` 利用以下命令绘制训练集样本的实际值和预测值曲线。 ```{r} plot(train[61:length(train)], col = "red", type = "l", lwd = 2, xlab = "Time", ylab = "Google Stock Price", main = "Google Stock Price Prediction") lines(pred_train, col = "blue", type = "l", lty = 3, lwd = 2) legend("topleft", legend = c("real_stock_price ","predicted_stock_price"), lty = c(1, 2), lwd = 2, col = c("red","blue"), bty = "n") ``` 预测曲线与实际曲线基本保持一致，说明 LSTM 对训练集的拟合效果非常好。利用以下命令绘制测试集样本的实际值和预测值曲线。 ```{r} plot(test[61:length(test)], col = "red", type = "l", lwd = 2, xlab = "Time", ylab = "Google Stock Price", main = "Google Stock Price Prediction") lines(pred_test, col = "blue", type = "l", lty = 2, lwd = 2) legend("topleft", legend = c("real_stock_price ","predicted_stock_price"), lty = c(1, 2), lwd = 2, col = c("red","blue"), bty = "n") ``` 下面计算在测试集上的均方误差根。 ```{r} mean(test[61:length(test)]) (n <- length(pred_test)) mse <- sum((test[61:length(test)]-pred_test)^2)/n rmse <- sqrt(mse) rmse ``` 在平均股价为 1012 的 300 个测试样本里，无状态 LSTM 在测试集上的均方误差根为 `r rmse`。 ### 有状态 LSTM 有状态（`stateful`）LSTM 能够在训练中维护跨批次的有状态信息，即当前批次的训练数据计算的状态值，可以用作下一批次训练数据的初始隐藏状态。有状态 LSTM 能让模型学习到输入样本之间的时序特征，适合一些长序列的预测，此时样本的前后顺序对模型是有影响的。`stateful` 代表除了每个样本内的时间步内传递，每个样本之间还会有信息 `(c, h)` 传递。 **当使用有状态 LSTM 时，需假定以下两种情况：** - 所有的批次都有相同数量的样本，即训练集和测试集的样本数量均能被 `batch_size` 整除。 - 如果 `x1` 和 `x2` 是连续批次的样本，则对于每一个 $i\in timesteps$，`x2[i]` 是 `x1[i]` 的后续序列。 **有状态LSTM的实现步骤如下：** - 必须将参数 `batch_size` 显式地传递给模型的第一层。 - 在 LSTM 层中将 `stateful` 设置为 `TRUE`。 - 在调用 `fit()` 函数时需指定 `shuffle` 为 `FALSE`，因为打乱样本后，序列之间就没有依赖性了。 - 训练完一个周期后，需要使用 `reset_state()` 函数来重置状态。接下来，将上一小节的无状态 LSTM 调整为有状态 LSTM，并使用训练集数据进行训练。 **在构建和拟合模型时的关键代码如下：** - 将 `batch_size` 设置为 50，则 1700（训练样本量）/50=34，300（测试样本量）/50=6。 - 构建模型时，将第一层 LSTM 网络参数 `return_sequences` 设置为 `TRUE`，`stateful` 设置为 `TRUE`，`batch_input_shape` 设置为 `c(batch_size, 60, 1))`；第二层 LSTM 网络的 `stateful` 设置为 `TRUE`。 - 使用 `for` 循环进行模型训练，需指定 `fit()` 函数的 `batch_size`，且将循环周期次数设置为 1，`shuffle` 设置为 `FALSE`；在 `fit()` 函数后使用 `reset_states()` 函数重置模型状态。以下是构建和训练有状态 LSTM 网络的程序代码。 ```{r} tic("构建模型 ") batch_size <- 50 num_epochs <- 50 # 构建模型 # 期望输入数据尺寸: (batch_size, timesteps, input_dim) stateful_model <- keras_model_sequential() stateful_model %>% layer_lstm(units = 50, return_sequences = TRUE, stateful = TRUE, batch_input_shape = c(batch_size, 60, 1)) %>% layer_dropout(rate = 0.2) %>% layer_lstm(units = 50, stateful = TRUE) %>% layer_dropout(rate = 0.2) %>% layer_dense(units = 1) %>% compile(loss = "mean_squared_error", optimizer = "adam") toc() ``` ```{r} tic("训练模型 ") # 训练模型 for (i in seq(num_epochs)){ stateful_model %>% fit( X_train, y_train, batch_size = batch_size, epochs = 1, shuffle = FALSE ) stateful_model %>% reset_states() } toc() ``` 模型训练好后，使用 `save_model_hdf5()` 函数将模型保存到本地。 ```{r} stateful_model %>% save_model_hdf5("model/lstm_stateful_model.h5") ``` 接下来，我们对测试集和训练集进行预测，进行反归一化处理，并分别绘制训练集和测试集的实际值和预测值曲线。 ```{r} # 模型预测 pred_train_stateful <- stateful_model %>% predict(X_train, batch_size = batch_size) pred_test_stateful <- stateful_model %>% predict(X_test, batch_size = batch_size) pred_train_stateful <- denormalize(pred_train_stateful, minval, maxval) pred_test_stateful <- denormalize(pred_test_stateful, minval, maxval) ``` ```{r} par(mfrow = c(1, 2)) # 绘制训练集的实际值与预测值曲线 plot(train[61:length(train)], col = "red", type = "l", lwd = 2, xlab = "Time", ylab = "Google Stock Price", ylim = c(200, 1200), main = "Train Data Prediction") lines(pred_train_stateful, col = "blue", type = "l", lty = 3, lwd = 2) legend("topleft", legend = c("real ", "predicted"), lty = c(1, 2), lwd = 2, horiz = FALSE, col = c("red", "blue"), bty = "n") # 绘制测试集的实际值与预测值曲线 plot(test[61:length(test)], col = "red", type = "l", lwd = 2, xlab = "Time", ylab = "Google Stock Price", ylim = c(200, 1200), main = "Test Data Prediction") lines(pred_test_stateful, col = "blue", type = "l", lty = 2, lwd = 2) legend("topleft", legend = c("real ", "predicted"), lty = c(1, 2), lwd = 2, horiz = FALSE, col = c("red", "blue"), bty = "n") # 去除图例外部方框 par(mfrow = c(1, 1)) ``` 下面计算有状态 LSTM 在测试集上的均方误差根。 ```{r} mean(test[61:length(test)]) (n <- length(pred_test_stateful)) mse <- sum((test[61:length(test)]-pred_test_stateful)^2)/n rmse <- sqrt(mse) rmse ``` 在平均股价为 1012 的 300 个测试样本里，有状态 LSTM 在测试集上的均方误差根为 `r rmse`。