LSTM+ARIMA模型

使用LSTM+ARIMA模型，更好拟合数据

赵昀康

日期

长短期记忆网络 (Long Short-Term Memory, LSTM) 模型是一种循环神经网络的变种，它能够学习长期依赖信息。自回归积分滑动平均模型(Autoregressive Integrated Moving Average, ARIMA) 是一种经典的时间序列预测模型，用于分析和预测时间序列数据。将LSTM和ARIMA结合使用，可以利用LSTM处理非线性和复杂时间序列数据的能力，然后通过ARIMA模型对LSTM的残差进行校正，以提高整体模型的预测性能。

1. 数据和环境配置

Data and environment

请确保如下程序包都已安装：

  • 如果尚未安装，请使用以下命令进行安装：

pip install + 程序包名称

程序包作用解释：

pandas: 用于数据分析和操作的库
numpy: 用于数值计算的库
os: 提供了与操作系统交互的功能
time: 用于时间相关的操作
dotenv: 用于从.env文件加载环境变量
matplotlib.pyplot: 用于绘图和可视化
statsmodels.graphics.tsaplots: 用于绘制时间序列的自相关图和偏自相关图
sklearn.metrics: 用于评估模型性能的指标计算
tensorflow: 用于构建和训练深度学习模型的库
pmdarima: 用于自动ARIMA最优参数选择的库
statsmodels.tsa.arima.model: 用于构建ARIMA模型的库

import pandas as pd
import numpy as np
import os
import time
from dotenv import load_dotenv
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_pacf,plot_acf
from sklearn.metrics import mean_squared_error, mean_absolute_error
import tensorflow as tf
from pmdarima import auto_arima
from statsmodels.tsa.arima.model import ARIMA

使用数据库简介

全国新型冠状病毒感染疫情情况数据简介

数据概览

本数据集包含了2022年12月至2024年11月期间，全国哨点医院报告的流感样病例中新冠病毒的每日阳性率数据。这些数据由中国疾病预防控制中心提供，为研究和分析我国新冠疫情的变化趋势提供了重要参考。

数据内容

  • 时间范围：2022年12月12日至2024年11月3日
  • 数据类型：每日流感样病例中新冠病毒阳性率
  • 数据来源：中国疾病预防控制中心官方网站

数据用途

这些数据可以用于：

  • 监测和分析新冠疫情的发展趋势
  • 评估疫情防控措施的效果
  • 为公共卫生决策提供数据支持

数据访问

您可以通过以下链接访问原始数据： 中国疾病预防控制中心

注意事项

  • 数据使用时请遵守相关法律法规，尊重数据所有权和知识产权。
  • 数据分析结果仅供参考，具体情况请结合实际。

以上简介提供了数据的基本信息和用途，以及如何访问这些数据。如果您对其他类型的数据感兴趣，可以替换上述内容中的“数据内容”部分，以反映不同数据集的具体信息。

data = pd.read_excel("F:\\Others\\全国哨点医院流感样病例新冠和流感病毒阳性率.xlsx")

data.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 693 entries, 0 to 692
Data columns (total 2 columns):
 #   Column            Non-Null Count  Dtype         
---  ------            --------------  -----         
 0   date              693 non-null    datetime64[ns]
 1   Positive Rate(%)  693 non-null    float64       
dtypes: datetime64[ns](1), float64(1)
memory usage: 11.0 KB

2. 数据预处理

训练预测数据集的划分方法

数据集划分比例

在我们的模型训练过程中，我们采用了82开的方法来划分数据集，即80%的数据用于训练模型，剩余的20%用于模型的预测和验证。这种划分比例是一种常见的做法，旨在确保模型在未见过的数据上具有良好的泛化能力。

划分理由

  • 训练集（80%）：用于构建和优化模型，使模型学习到数据中的特征和模式。
  • 预测集（20%）：用于评估模型的性能，确保模型在新数据上的预测能力。

实际应用

您可以根据自己的实际需求调整训练预测集的比例。例如，如果您有更多的数据或希望模型在新数据上有更好的表现，您可能需要增加预测集的比例。相反，如果数据量较少，您可能需要增加训练集的比例以确保模型有足够的数据来学习。

调整方法

调整数据集划分比例的方法通常很简单，只需要更改数据划分的代码即可。
通过调整len(data)*0.2的比例和days相当即可。
比如37开，则调整为len(data)*0.3，修改days为208

print(len(data)*0.2)

days = 140


train_len_val = len(data) - days

train = data['Positive Rate(%)'].iloc[0:train_len_val]
test = data['Positive Rate(%)'].iloc[train_len_val:]

138.6

数据转换

目的是将时间序列数据转换为适合输入到长短期记忆网络（LSTM）的格式。我们通过创建一个滑动窗口来实现的，该窗口用于生成输入序列和对应的目标值。

代码解释

  • apply_transform 函数接受两个参数：data（时间序列数据）和 n（滑动窗口的大小）。
  • middle_data 用于存储输入序列，而 target_data 用于存储目标值。
  • 通过循环，从索引 n 到数据长度，每次移动一步，创建长度为 n 的输入序列（input_sequence），并将当前索引处的值作为目标值。
  • middle_data 被转换成三维数组，其中第一维是样本数量，第二维是时间步长（n），第三维是特征数量（这里假设为1）。
  • target_data 被转换成一维数组，包含每个输入序列对应的目标值。
def apply_transform(data, n: int):
    middle_data = []
    target_data = []
    for i in range(n, len(data)):
        input_sequence = data[i-n:i]  
        middle_data.append(input_sequence) 
        target_data.append(data[i])
    middle_data = np.array(middle_data).reshape((len(middle_data), n, 1))
    target_data = np.array(target_data)
    return middle_data, target_data

# 假设 n 是滑动窗口大小，例如 n = 10
n = 30

# 将训练数据转换为适合LSTM的格式
middle_data, target_data = apply_transform(train.values, n)

3. 构建LSTM模型

模型参数说明

  • train: 训练集，因为是时间序列，所以是一维数组形式。
  • n: 滑动窗口大小，用于生成输入序列的长度。
  • number_nodes: LSTM模型中的节点数量。
  • learning_rate: 学习率，用于优化算法中调整权重的速率。
  • epochs: 训练周期数，即数据将被用于训练模型的次数。
  • batch_size: 每次梯度更新前用于训练模型的样本数量。

模型参数的重要性

在机器学习中，参数是学习过程中不可从数据中直接学习的参数。它们是影响算法性能的关键因素，通常需要通过经验或某种搜索策略来确定。超参数的设置直接影响模型的学习能力和泛化能力。

调整参数的意义

调整参数对于构建高效的机器学习模型至关重要，因为：

  1. 提高模型性能：合适的参数可以显著提高模型的准确度和泛化能力。
  2. 防止过拟合和欠拟合：通过调整参数，可以平衡模型的复杂度，避免过拟合（模型在训练数据上表现太好，在新数据上表现差）和欠拟合（模型在训练数据上表现就很差）。
  3. 优化训练过程：合理的参数设置可以加快模型的训练速度，减少计算资源的消耗。

参数调优步骤

1. 定义参数网格

我们首先定义一个字典 param_grid，其中包含了我们想要尝试的不同参数值。这些参数包括LSTM层的节点数、学习率、批次大小和训练周期数。

2. 防过拟合

防止模型过拟合，我们使用早停策略。这是一个回调函数，当模型在一定数量的周期内没有改善时，它会停止训练。

3. 选择评估指标

我们选择平均绝对误差（mean_absolute_error，mse）评估模型的拟合的好坏，mse越低，说明模型拟合的效果越好

4.网格搜索

我们将遍历超参数网格中的所有可能组合，并使用 LSTM 函数训练模型。对于每一组参数，我们都会评估模型的性能，并更新最佳分数和最佳参数。

5. 输出最佳参数

from tensorflow.keras.callbacks import EarlyStopping
param_grid = {
    'number_nodes': [20, 50, 80],
    'learning_rate': [0.001, 0.0005, 0.0001],
    'batch_size': [16, 32, 64],
    'epochs': [50, 100]
}

# 初始化最佳分数和最佳参数
best_score = np.inf
best_params = None

# 定义早停策略以避免过拟合
early_stopping = EarlyStopping(monitor='loss', patience=10, min_delta=0.001)

def LSTM(train, n: int, number_nodes, learning_rate, epochs, batch_size, early_stopping):
    middle_data, target_data = apply_transform(train, n)
    model = tf.keras.Sequential([
        tf.keras.layers.Input((n, 1)),
        tf.keras.layers.LSTM(number_nodes, input_shape=(n, 1)),
        tf.keras.layers.Dense(units=number_nodes, activation="relu"),
        tf.keras.layers.Dense(units=number_nodes, activation="relu"),
        tf.keras.layers.Dense(1)
    ])
    model.compile(loss='mse', optimizer=tf.keras.optimizers.Adam(learning_rate), metrics=["mean_absolute_error"])
    
    history = model.fit(middle_data, target_data, epochs=epochs, batch_size=batch_size, verbose=0, callbacks=[early_stopping])
    full_predictions = model.predict(middle_data).flatten()
    return model, history, full_predictions



for nn in param_grid['number_nodes']:
    for lr in param_grid['learning_rate']:
        for bs in param_grid['batch_size']:
            for ep in param_grid['epochs']:
                #print(f"Training model with params: {nn} nodes, {lr} learning rate, {bs} batch size, {ep} epochs")
                model, history, predictions = LSTM(train, n, nn, lr, ep, bs, early_stopping)
                
                # 计算验证集上的损失，这里假设 target_data 是验证集
                val_loss = model.evaluate(middle_data, target_data, verbose=0)
                
                # 检查是否有更好的分数
                if val_loss[0] < best_score:  # 选择返回列表的第一个元素
                    best_score = val_loss[0]
                    best_params = (nn, lr, bs, ep)
                    best_model = model  # 保存最佳模型

# 输出最佳参数和模型
print(f"Best params: {best_params}, Best score: {best_score}")

Best params: (80, 0.001, 16, 50), Best score: 7.317050767596811e-05

使用最优参数的模型拟合

number_nodes = 80
learning_rate = 0.001
epochs = 16
batch_size = 50


def LSTM(train, n: int, number_nodes, learning_rate, epochs, batch_size):
    middle_data, target_data = apply_transform(train, n)
    model = tf.keras.Sequential([
        tf.keras.layers.Input((n, 1)),
        tf.keras.layers.LSTM(number_nodes, input_shape=(n, 1)),
        tf.keras.layers.Dense(units=number_nodes, activation="relu"),
        tf.keras.layers.Dense(units=number_nodes, activation="relu"),
        tf.keras.layers.Dense(1)
    ])
    model.compile(loss='mse', optimizer=tf.keras.optimizers.Adam(learning_rate), metrics=["mean_absolute_error"])
    
    history = model.fit(middle_data, target_data, epochs=epochs, batch_size=batch_size, verbose=0)
    full_predictions = model.predict(middle_data).flatten()
    return model, history, full_predictions

# 拟合模型
model, history, full_predictions = LSTM(train, n, number_nodes, learning_rate, epochs, batch_size)

模型效果评估

模型性能评估指标

在机器学习和统计建模中，评估模型性能是一个关键步骤。我们使用如下几种常用的性能评估指标，它们可以帮助我们了解模型预测的准确性和可靠性。

1. 均方误差（Mean Squared Error, MSE）

均方误差是预测值与实际值之差的平方的平均值。它能够量化模型预测误差的大小，并且对较大的误差给予更大的惩罚。

2. 均方根误差（Root Mean Squared Error, RMSE）

均方根误差是MSE的平方根，与原始数据在同一量纲上，更容易解释。它衡量了预测值与实际值之间的标准偏差。

3. 平均绝对误差（Mean Absolute Error, MAE）

平均绝对误差是预测值与实际值之差的绝对值的平均值。与MSE相比，MAE对异常值的敏感度较低。

4. 决定系数（R² Score, R²）

决定系数，也称为R²分数，衡量了模型预测值与实际值之间的相关程度。R²值越接近1，表示模型的解释能力越强。

from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score

def evaluate_performance(actual, predicted):
    # 计算 MSE
    mse = mean_squared_error(actual, predicted)
    # 计算 RMSE
    rmse = np.sqrt(mse)
    # 计算 MAE
    mae = mean_absolute_error(actual, predicted)
    # 计算 R²
    r2 = r2_score(actual, predicted)
    
    return mse, rmse, mae, r2

def Error_Evaluation(train_data, predict_train_data, n: int):
    errors = []
    for i in range(len(predict_train_data)):
        err = train_data[n + i] - predict_train_data[i]
        errors.append(err)
    return errors

errors = Error_Evaluation(train[n:], full_predictions, n)

last_sequence = train[-n:].values.reshape((1, n, 1))
predictions = []
for i in range(days):
    next_prediction = model.predict(last_sequence).flatten()[0]
    predictions.append(next_prediction)
    if i < len(test):
        actual_value = test.iloc[i]
        new_row = np.append(last_sequence[:, 1:, :], np.array([[[actual_value]]]), axis=1)
    else:
        new_row = np.append(last_sequence[:, 1:, :], np.array([[[next_prediction]]]), axis=1)        
    last_sequence = new_row.reshape((1, n, 1))

evaluate_performance(test,predictions )

(0.0006135117138462903,
 0.02476916861435382,
 0.02390339080614434,
 0.9185233816505528)

4. 构建ARIMA模型

什么是ARIMA模型？

ARIMA模型，即自回归积分滑动平均模型（AutoRegressive Integrated Moving Average Model），是一种广泛应用于时间序列预测的统计模型。它结合了自回归（AR）、差分（I）和移动平均（MA）三个部分来建模时间序列数据。

ARIMA模型的参数：p、d、q

ARIMA模型的参数p、d、q分别代表：

  • p（自回归项）：模型中自回归项的阶数，即模型将使用多少个过去的值来预测未来的值。
  • d（差分次数）：使时间序列平稳所需的差分次数。差分是将时间序列转换为平稳序列的过程。
  • q（移动平均项）：模型中移动平均项的阶数，即模型将使用多少个过去的预测误差来预测未来的值。

参数调整的重要性

  • 提高预测准确性：正确的参数可以使模型更好地捕捉时间序列的数据结构，从而提高预测的准确性。
  • 模型稳定性：适当的差分次数（d）可以使非平稳时间序列变得平稳，这是ARIMA模型有效性的基础。
  • 避免过拟合和欠拟合：选择合适的p和q可以避免模型过于复杂（过拟合）或过于简单（欠拟合）。

使用auto_arima自动寻找最优参数

auto_arima函数来自pmdarima库，它自动寻找最佳的p、d、q参数组合，以拟合给定的时间序列数据。

本次我们使用ARIMA模型校准LSTM模型的残差，提高模型的性能 

def Parameter_calculation(data):
   finding = auto_arima(data,trace = True)
   ord = finding.order
   return ord

ord = Parameter_calculation(errors)

Performing stepwise search to minimize aic

Best model:  ARIMA(4,0,2)(0,0,0)[0] intercept
Total fit time: 136.189 seconds

使用最优的参数拟合ARIMA模型

def ARIMA_Model(train,len_test,ord):
   model = ARIMA(train, order = ord)
   model = model.fit()
   predictions = model.predict(start = len(train),end = len(train) + len_test ,type='levels')
   full_predictions = model.predict(start = 0,end = len(train)-1,type='levels')
   return model,predictions,full_predictions

Arima_Model,predictions_errors,full_predictions_errors = ARIMA_Model(errors,len(test),ord)

Arima_Model.summary()
SARIMAX Results
Dep. Variable: y No. Observations: 523
Model: ARIMA(4, 0, 2) Log Likelihood 2248.723
Date: Tue, 14 Jan 2025 AIC -4481.447
Time: 16:38:15 BIC -4447.370
Sample: 0 HQIC -4468.101
- 523
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
const -0.0046 0.006 -0.803 0.422 -0.016 0.007
ar.L1 1.5830 0.109 14.544 0.000 1.370 1.796
ar.L2 -0.0771 0.055 -1.410 0.159 -0.184 0.030
ar.L3 -0.9744 0.108 -9.063 0.000 -1.185 -0.764
ar.L4 0.4533 0.072 6.285 0.000 0.312 0.595
ma.L1 -0.0085 0.113 -0.075 0.940 -0.230 0.213
ma.L2 -0.6293 0.120 -5.260 0.000 -0.864 -0.395
sigma2 1.07e-05 1.11e-07 95.952 0.000 1.05e-05 1.09e-05

结合LSTM的结果和校准残差值

def Final_Predictions(predictions_errors,predictions):
   final_values = []
   for i in range(days):
      final_values.append(predictions_errors[i] + predictions[i])
   return final_values

final_predictions = Final_Predictions(predictions_errors,predictions)

5.前后效果评估对比

校准后，效果拔群 

evaluate_performance(test,predictions )

(0.0006135117138462903,
 0.02476916861435382,
 0.02390339080614434,
 0.9185233816505528)

evaluate_performance(test[:days], final_predictions[:days])

(0.0003946475114226618,
 0.019865737122560085,
 0.018279961643969753,
 0.9475893549462051)

6. 应用拓展

a. 趋势预测
b. 其他实践项目实践: 换用其他数据源分析
c. 其他模型组合项目实践: 探索将与其他时间序列模型（如Prophet、SARIMA）或机器学习模型（如XGBoost、Random Forest）结合，发挥各自的优势，提升预测效果

赵昀康
艾泽拉斯的的守护者、德玛西亚的护卫者
Ljung-Box (L1) (Q): 0.00 Jarque-Bera (JB): 2602302.20
Prob(Q): 0.95 Prob(JB): 0.00
Heteroskedasticity (H): 0.03 Skew: 16.42
Prob(H) (two-sided): 0.00 Kurtosis: 347.00