ARIMA + LSTM模型

使用R实现ARIMA + LSTM结合模型，对疾病进行综合拟合与预测。

吴嘉栋

日期

基本介绍

ARIMA（AutoRegressive Integrated Moving Average，自动回归积分滑动平均）模型和LSTM（Long Short-Term Memory，长短时记忆网络）模型都是用于时间序列预测的强大工具。ARIMA模型主要用于捕捉时间序列中的线性依赖关系，适用于平稳或经过差分处理的非平稳数据。而LSTM则是一种神经网络，能够处理时间序列中的非线性模式，特别适用于长时间依赖的数据。通过将ARIMA与LSTM结合，可以同时利用它们各自的优势，进一步提高时间序列预测的准确性。

下载讲解视频

ARIMA

AR（自回归）：通过历史数据的线性组合来预测当前值，即当前时间点的数据依赖于其前几个时间点的值。

I（差分）：用于处理非平稳时间序列，使其变得平稳，从而便于建模。差分操作是通过计算序列当前值与前一个值的差异来实现的。

MA（滑动平均）：考虑到预测误差的影响，MA部分通过历史误差的线性组合来预测当前值。

ARIMA模型适用于平稳时间序列，若数据非平稳，通过差分（I部分）将其转换为平稳序列。该模型的核心思想是结合自回归（AR）和滑动平均（MA）过程来捕捉数据中的依赖性，并通过差分来解决非平稳问题。

更详细内容可见https://blog.csdn.net/weixin_62586597/article/details/120626460

LSTM

LSTM是一种深度学习模型，是改进版的RNN（递归神经网络），能够有效处理时间序列数据中的长期依赖问题。与传统的RNN不同，LSTM通过引入遗忘门、输入门和输出门的机制，有效避免了梯度消失问题，从而能够记住长期的依赖信息。LSTM尤其适用于具有非线性、复杂模式的数据，能够捕捉到数据中的长期趋势和周期性变化。与ARIMA不同，LSTM不需要假设数据的平稳性，因此对于非平稳且含有复杂非线性模式的时间序列表现优异。

更详细内容可见https://blog.csdn.net/mary19831/article/details/129570030

ARIMA + LSTM

ARIMA模型擅长处理时间序列中的线性关系，而LSTM则能很好地捕捉非线性关系。在时间序列数据中，许多现象具有复杂的非线性模式，单纯依靠ARIMA或LSTM往往难以取得最佳效果。因此，将ARIMA和LSTM结合起来，可以充分发挥两者的优势。

具体来说，结合方法通常为：首先，使用ARIMA对时间序列数据进行建模，得到拟合值和残差。ARIMA能够有效地拟合数据中的线性趋势和季节性变化。然后，将ARIMA的残差作为输入数据，使用LSTM进一步建模，捕捉数据中未被ARIMA捕捉到的非线性部分。最后，结合ARIMA的拟合结果和LSTM的预测值，得到最终的预测结果。通过这种结合方式，可以大幅度提高时间序列预测的准确性，特别是在非线性模式复杂的情况下，ARIMA+LSTM模型能够显著改进预测性能。

代码示例

1、数据导入

#arima以及数据处理相关包
library(tseries)
library(forecast)
library(ggplot2)
library(dplyr)

# LSTM需要python环境
#install.packages("keras")
#install_keras() 没装记得装
library(keras)
library(keras)
library(tensorflow)

#导入数据
data <- data.frame(
  t = c(1:120), 
  cases = c(157.366035757709, 171.130138146691, 161.693157234169, 135.374180545006, 114.007159079689, 111.082056488276, 119.201409558046, 95.6803503331459, 71.1906322070109, 35.598717573228, 60.2034334670303, 112.048883891603, 144.845902708285, 165.889144760438, 152.0816505714, 122.766314554685, 106.380414739982, 97.8680987913461, 108.02581570496, 111.821570958674, 75.0370942955287, 50.6659003170483, 82.0655650775372, 140.953367194482, 195.892170500805, 205.459477893884, 182.034976442226, 142.0568961126, 130.694989846144, 127.059814511306, 132.95483025752, 113.760595547534, 84.9911804004458, 63.0694971323979, 90.0000850405908, 131.101883588582, 195.680742185678, 220.807156372956, 199.504676188912, 182.840636455921, 169.473941911273, 169.899653426993, 180.722696441136, 157.949920670289, 121.773973179475, 115.681882005576, 121.686950594746, 184.678962922351, 227.012229122205, 251.875799188187, 232.125155418183, 199.032258191127, 198.285702378847, 191.950711312486, 201.131080443658, 190.094942226544, 164.339401112219, 150.475743218537, 159.624770150229, 234.153867274929, 287.156980668974, 310.589169876908, 296.510037547029, 269.726787903443, 243.033455321117, 251.481290333471, 253.896512319433, 238.564959508394, 211.098463303815, 191.882538635884, 206.554307921509, 266.889571574292, 328.465174202973, 350.754993853567, 338.149838426153, 322.733486292356, 311.394593710661, 317.045436273286, 314.319160862374, 306.373790888873, 286.2824039625, 271.294420906394, 289.564214535972, 338.61876054452, 392.831965356688, 419.760580817102, 401.763647709361, 376.57394082233, 369.104109829162, 374.128549311613, 382.159382916695, 370.228358626996, 348.937021133561, 343.86690133925, 359.684003622442, 422.110655747725, 480.976729832121, 504.073238092168, 503.568155353901, 479.991825540222, 471.450621177379, 475.095540047885, 472.060657681026, 474.68655194079, 431.146806261249, 420.567658383368, 443.848371373664, 503.134929966743, 565.632875634455, 591.392284612901, 581.346867955645, 563.328484217726, 554.556184160117, 536.846908386078, 522.688272498988, 511.096942609685, 501.894252847389, 516.030657951073, 521.970923505595, 548.937412697586 )) #生成的，120步

2、数据处理

#将数据被分为两部分：训练和验证
training_data<-subset(data,t >= 1 & t <= 108) 
validation_data<-subset(data,t >= 109 & t <= 120) 

# 将数据转化为时间序列,使用ts()函数
ts_data <- ts(training_data$cases, frequency =12, start = c(2000, 1))
test_data <- ts(validation_data$cases, frequency = 12, start = c(2009, 1))

# 拆分季节性、趋势、随机性(时间序列分解，看看各种特性)
decomposition_data <- decompose(ts_data)
plot(decomposition_data)

3、数据检查

#稳定性检验，P值小于0.01说明数据稳定，可以建模
adf.test(ts_data)
#白噪声检验，P值小于0.05，序列不是白噪声，可以建模
Box.test(ts_data,type = "Ljung-Box")

4、模型构建

#使用auto.arima，自动选取最优模型
auto.arima(ts_data,stepwise = F,trace = T,stationary = T,ic=c("aic"))
#找到上方代码给出的best model
fit <- arima(ts_data,order=c(1,0,2),seasonal = list(order=c(0,0,2),period=12))

5、模型检查

#正态性检验
shapiro.test(fit$residuals)#残差是否服从正态(0假设残差正态)
qqnorm(fit$residuals, main = "Q-Q Plot of Residuals")#QQ图
#自相关性检验(0假设残差是白噪声)
Box.test(fit$residuals,lag=24,type="Ljung-Box") #非季节性=10,季节性=period*2
tsdisplay(residuals(fit),lag.max=20,main='残差')#残差相关性

残差正态但有明显相关性，说明arima模型没有很好的将数据中存在的某些规律抓取出来，此时，可以考虑使用lstm模型对残差进行二次建模。

6、效果评估

#预测测试集
prediction<-forecast(fit,12)#向后预测12个月，也就是测试集的年份
# 提取预测值和实际值
predicted_values <- prediction$mean  # 对验证集的预测值
actual_values <- test_data  # 验证集的实际值

# 计算残差（实际-预测）
errors <- actual_values - predicted_values

# 计算指标
mean(abs(errors))  # 平均绝对误差（通常情况下越小越好）
mean(abs(errors / actual_values)) * 100

#预测2008年
plot(prediction)

mae 42.45393

mape 7.867773%（效果不错，但残差提示拟合有进一步优化的空间，LSTM启动）

7、LSTM模型构建

# 提取ARIMA模型的残差
arima_residuals <- fit$residuals

# 标准化残差
scaled_residuals <- scale(arima_residuals)

# 将数据分为训练集和测试集
train_size <- 80
train_data <- scaled_residuals[1:train_size]
test_data <- scaled_residuals[(train_size + 1):length(scaled_residuals)]

# 将数据重塑为LSTM输入格式 [samples, time steps, features]
create_dataset <- function(data, time_step = 12) {  # 时间步为12
  x <- list()
  y <- list()
  for(i in 1:(length(data) - time_step)) {
    x[[i]] <- data[i:(i + time_step - 1)]
    y[[i]] <- data[i + time_step]
  }
  return(list(array(unlist(x), dim = c(length(x), time_step, 1)), unlist(y)))
}

time_step <- 12  # 设置时间步为12
train_data_lstm <- create_dataset(train_data, time_step)
test_data_lstm <- create_dataset(test_data, time_step)

x_train <- train_data_lstm[[1]]
y_train <- train_data_lstm[[2]]
x_test <- test_data_lstm[[1]]
y_test <- test_data_lstm[[2]]


# 3.构建LSTM模型
model <- keras_model_sequential() %>%
  layer_lstm(units = 64, input_shape = c(time_step, 1)) %>% 
  # 单层LSTM，50个神经元
  layer_dense(units = 1)  # 直接输出结果


# 4. 编译模型，调整学习率
model %>% compile(
  loss = 'mean_squared_error',
  optimizer = optimizer_adam(lr = 0.001)  # 设置学习率为0.001
)

# 5. 设置早停机制（如果需要）
early_stop <- callback_early_stopping(monitor = "val_loss", patience = 50, restore_best_weights = TRUE)

# 6. 训练模型
history <- model %>% fit(
  x_train, y_train,
  epochs = 200,
  batch_size = 16,
  validation_data = list(x_test, y_test),
  callbacks = list(early_stop),  # 如果使用早停
  verbose = 1
)
save_model_hdf5(model, "arima_lstm_model.h5")

8、预测残差

#导入模型
model <- load_model_hdf5("arima_lstm_model.h5")

# 预测未来第109-120到步的残差
future_steps <- 12  # 预测12步

# 获取最后一段训练数据作为输入
last_train_data <- tail(train_data, time_step)

# 构建未来的预测输入
future_predictions <- numeric(future_steps)
current_input <- last_train_data

for (i in 1:future_steps) {
  # 预测当前步
  next_prediction <- model %>% predict(array(current_input, dim = c(1, time_step, 1)))
  future_predictions[i] <- next_prediction
  
  # 将预测结果更新到输入数据中（窗口滑动）
  current_input <- c(current_input[-1], next_prediction)
}

# 反归一
mu <- mean(arima_residuals)  # 原始残差数据的均值
sigma <- sd(arima_residuals)  # 原始残差数据的标准差
future_predictions_original <- future_predictions * sigma + mu

# 输出反归一化后的预测结果
future_predictions_df <- data.frame(
  step = 109:120,
  predicted_residuals = future_predictions_original
)

pre_al <- prediction$mean+future_predictions_df$predicted_residuals

pre_al <- data.frame(
  t = 109:120,  # 时间序列
  cases = as.numeric(pre_al))

9、效果评估

# 计算残差（实际-预测）
errors2 <- validation_data$cases - pre_al$cases
# 计算指标
mean(abs(errors2))  # 平均绝对误差（通常情况下越小越好）
mean(abs(errors2 / actual_values)) * 100

mae 41.818

mape 7.754341%（变好了！）

10、结合模型

##将arima和lstm的拟合值和预测值结合
# 获取ARIMA模型的拟合值和预测值
arima_fitted_values <- prediction$fit  # ARIMA模型在历史数据上的拟合值（1到108步）
arima_forecast_values <- prediction$mean  # ARIMA模型对未来步骤的预测值（109到120步）

# 获取LSTM模型的未来预测残差（第13步到第108步的残差）
future_steps_initial <- 96  # 从第13步到第108步，共96步
current_input <- train_data  # 使用训练集的全部数据（1到108步）作为输入

future_predictions_initial <- numeric(future_steps_initial)

for (i in 1:future_steps_initial) {
  # 预测当前步
  next_prediction <- model %>% predict(array(current_input, dim = c(1, time_step, 1)))
  future_predictions_initial[i] <- next_prediction
  
  # 更新输入数据（窗口滑动）
  current_input <- c(current_input[-1], next_prediction)
}

# 反归一化处理（这里依然使用原始残差的均值和标准差）
future_predictions_initial_original <- future_predictions_initial * sigma + mu

# 合并ARIMA拟合值和LSTM预测残差
# 对于第1到12步，直接使用ARIMA拟合值
final_predictions_1_to_12 <- arima_fitted_values[1:12]

# 对于第13到108步，将ARIMA的拟合值和LSTM的预测残差加和
final_predictions_13_to_108 <- arima_fitted_values[13:108] + future_predictions_initial_original

# 对于第109到120步，将ARIMA的预测值和LSTM的预测残差加和
final_predictions_109_to_120 <- arima_forecast_values + future_predictions_initial_original[(future_steps_initial - 11):future_steps_initial]

# 合并所有的预测值（1到120步）
final_predictions <- c(final_predictions_1_to_12, final_predictions_13_to_108, final_predictions_109_to_120)

# 创建最终的预测数据框
pre_al_df <- data.frame(
  t = 1:120,  # 时间序列
  cases = final_predictions  # 合并后的预测值作为 cases 数据
)


# 合并真实数据和预测数据（前12步用实际数据，后续步骤用预测值）
combined_data <- data.frame(
  t = 1:120,
  actual_cases = data$cases,  # 真实数据
  predicted_cases = pre_al_df$cases  # 合并后的预测数据
)

11、来个简图


# 绘图
ggplot(combined_data, aes(x = t)) +
   # 添加浅蓝色背景，放置在底层
  geom_rect(aes(xmin = 109, xmax = 120, ymin = -Inf, ymax = Inf), fill = "lightblue", alpha = 0.01) + 
  geom_line(aes(y = actual_cases), color = "black", size = 1) +  # 绘制真实数据线
  geom_line(aes(y = predicted_cases), color = "red", size = 1, linetype = "dashed") +  # 绘制预测数据线
  labs(
    title = "原始数据与拟合数据的比较",
    x = "时间",
    y = "病例数"
  ) +
  theme_minimal() +
  scale_x_continuous(breaks = seq(1, 120, by = 12))  # 设置x轴的刻度

分析讨论

  1. ARIMA模型擅长捕捉时间序列中的线性趋势和短期依赖性，但对于复杂的非线性关系和长期依赖性，它的表现较弱。LSTM通过其独特的门控机制，能够从过去的时序数据中学习到复杂的动态变化，因此结合ARIMA和LSTM的方法，能够同时处理线性和非线性的部分，提升模型的整体预测精度。

  2. ARIMA+LSTM方法同样存在局限性。首先，LSTM模型的训练需要大量数据和计算资源计算成本和时间开销较大；其次，LSTM的超参数调优过程比较复杂，需要通过网格搜索、随机搜索等方法反复调整参数；再次，LSTM不是万能的，数据不适合例如残差是白噪声（很多时候就会这样），LSTM表现会很差。再加上这个方法本身最大头的部分是依赖于arima，所以预测的好不好还要看arima表现。

  3. 为了改进ARIMA+LSTM组合模型的性能，首先，可以优化ARIMA部分，选择更好的模型，甚至可以直接使用单一的lstm看看效果；其次，针对LSTM的部分，可以通过超参数调优（如学习率、层数、单位数等）来提升预测精度，同时要注意避免过拟合；再次，可以通过采用GPU加速技术、分布式计算等手段提高训练效率。此外，还可以探索其他深度学习架构（如GRU、Transformer）来补充和替代LSTM，以期获得更好的处理能力和预测精度。

  4. 个人拙见：对于数据拟合预测，究竟使用什么模型，要不要用混合模型，都是要多试，找到最合适的才对，而不是复杂至上。

吴嘉栋
Still KD