ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

SEIR+LSTM残差学习:Python实现传染病预测时间序列模型

SEIR+LSTM残差学习:Python实现传染病预测时间序列模型 简介这份资源将传染病动力学模型SEIR与LSTM神经网络相结合用于2019新型冠状病毒肺炎COVID-19疫情预测适合计算机、数据科学、人工智能等相关专业的学生作为课程设计、期末大作业或初期立项演示使用。项目一方面通过SEIR模型控制接触率来模拟不同干预程度体现防控措施效果另一方面借助LSTM神经网络以前三天数据预测第四天感染趋势并预留干预值输入接口以优化长期预测。压缩包共31个文件其中12个Python脚本覆盖基础版、带干预的模型对比及累计预测等不同实验场景14张PNG图像展示模型拟合与预测结果另有2个Excel真实数据文件和2份Markdown项目说明整体仅1.79MB轻量易用。已有383人学习下载项目代码经验证可稳定运行注释详细、便于二次开发适合入门进阶或作为毕业设计、课程设计的参照方案。1. 传染病预测的难点为什么SEIR和LSTM要凑到一起2020年初我试着用公开数据预测疫情走势第一批模型全部翻车。纯SEIR模型对参数极其敏感β调一档预测峰值能差一倍纯LSTM则把早期平缓增长当成了常态预测峰值连真实的一半都不到。后来我改用串联结构SEIR先提供一个有物理约束的基线LSTM再去学习基线和真实数据之间的残差预测曲线才稳定下来。这就是标题里“SEIRLSTM”的核心思路也是这类python源码项目想让你复现的东西。适合正在做时间序列预测、传染病数据分析、公共卫生政策评估的人尤其是想搞懂“机制模型数据驱动”怎么真正结合的新手。2. 先把SEIR的数学账算清从仓室模型到模拟数据2.1 SEIR四个仓室与三个关键参数SEIR把人群分成四类易感者S、潜伏者E、感染者I、康复者R。模型假设康复者不再被感染对COVID-19这种短时间尺度来说基本成立。真正决定曲线形态的只有三个参数β是接触传染率σ是潜伏者转阳率等于1/平均潜伏期γ是感染者康复率等于1/平均感染周期。基本再生数R0 β / γ只看β和γσ不影响最终感染规模但决定峰值的早晚。微分方程组如下dS/dt -β S I / NdE/dt β S I / N - σ EdI/dt σ E - γ IdR/dt γ I这里省略了出生、死亡和自然免疫衰减也无症状感染者没有单独建仓。三个参数怎么理解β越大扩散越快σ越大潜伏期越短第一个峰值来得越早γ越大感染期越短整条曲线越扁平。这几个参数是所有预测的起点后面必须用真实数据去拟合不能拍脑袋。2.2 用Python把SEIR跑起来最小模拟代码import numpy as np from scipy.integrate import odeint def seir_simulate(params, days, N, E0, I0): params: [beta, sigma, gamma] days: 模拟天数 N: 总人口 E0: 初始潜伏者人数 I0: 初始感染者人数 beta, sigma, gamma params S0 N - E0 - I0 y0 [S0, E0, I0, 0] # S, E, I, R def deriv(y, t): S, E, I, R y dS -beta * S * I / N dE beta * S * I / N - sigma * E dI sigma * E - gamma * I dR gamma * I return [dS, dE, dI, dR] t np.arange(0, days, 1) result odeint(deriv, y0, t) return result # 示例总人口1000万初始潜伏者10人感染者5人 sim seir_simulate([0.8, 1/5.2, 1/14.0], 120, 10_000_000, 10, 5) daily_new_infected np.diff(sim[:, 2]) # 每日新增感染者约等于I的变化量这段代码把SEIR四个状态放进同一个向量用odeint解常微分方程组。注意每日新增没有直接出现在方程组里我在这里用np.diff(sim[:, 2])做近似更严谨的做法是用积分窗口内的净变化量但在日粒度数据下差异很小可以直接用。参数说明beta0.8 是早期疫情地区的估计值sigma1/5.2 对应约5.2天的平均潜伏期gamma1/14 对应平均14天从感染到康复。换毒株或改地区后这些值都要重新拟合。代码里odeint的返回值每一行是各仓室的人数列的顺序和y0一致所以取第3列是I。2.3 SEIR模型的三个硬伤为什么预测会失真第一β不是常数。封城、戴口罩、接种疫苗都会改变接触传染率经典SEIR却假设它全程不变。第二参数之间存在补偿效应不同的(β, γ)组合能拟合出几乎重合的确诊曲线但外推Future时差异巨大这是调参里最折磨人的黑匣子。第三报告数据有滞后确诊数不等于真实感染数直接用报告序列拟合会低估E和I峰值被拉平。这时LSTM的价值就出来了它不需要预设机制可以从残差里自动学到SEIR没建模的那部分动态。常见做法是用SEIR输出作为基线把真实新增数与SEIR预测值的差作为LSTM的训练目标。但要注意LSTM不是用来替代SEIR的。疫情数据通常只有几百个点用纯LSTM做多步外推几乎必然退化让LSTM只预测残差它的任务很轻模型容量不必很大反而更稳。这也是这个方案和坊间很多“直接拿LSTM预测确诊人数”的半吊子做法最大的区别。一句话总结SEIR负责把传染病常识写进预测LSTM负责修补SEIR没见过的东西。先把这个分工想清楚后面的代码才看得懂。3. 用LSTM修正SEIR残差时间序列预测的落地实现3.1 数据准备把每日新增确诊整理成监督学习样本真实数据通常是一列日期加一列每日新增数。第一步是把SEIR预测值和真实报告日对齐日期跨度不一致时要补零或插值。第二步是构造滑动窗口用过去14天的真实新增加上当前时刻的SEIR预测值去预测当天的残差。import numpy as np import torch def make_samples(real_new, seir_new, window14): X_real, X_seir, y [], [], [] for i in range(window, len(real_new)): X_real.append(real_new[i-window:i]) X_seir.append(seir_new[i]) # 当前时刻的SEIR预测 y.append(real_new[i] - seir_new[i]) # 残差 return (np.array(X_real, dtypenp.float32), np.array(X_seir, dtypenp.float32), np.array(y, dtypenp.float32)) # real_new来自公开数据seir_new来自seir_simulate()得到的每日新增 X_real, X_seir, y make_samples(real_new, seir_new, window14) print(X_real.shape, X_seir.shape, y.shape)逻辑说明X_real 是过去14天的真实新增序列X_seir 是当天SEIR模型的预测值不是序列是一个数y 是真实值减去SEIR预测值。窗口结束点就是预测点整个过程不混入未来信息这是时间序列预测不能破的底线。参数说明window 设为14一个潜伏期加一段传染期模型能看到一个相对完整的感染代际。数据只有200天时window14 会得到186个样本够训练一个轻量LSTM。若改用7样本数增加到193但短期波动更大模型更容易被周末效应干扰。可以先7后14对比。3.2 搭建LSTM模型两个输入分支的残差学习结构模型分两条路LSTM分支处理真实历史序列输出最后一个时刻的隐状态SEIR分支是一个很小的全连接网络把当前SEIR预测值映射成特征。两条路的输出拼起来再接全连接层输出残差。import torch.nn as nn class SEIRLSTM(nn.Module): def __init__(self, input_size1, hidden_size32, num_layers1): super(SEIRLSTM, self).__init__() self.lstm nn.LSTM(input_size, hidden_size, num_layers, batch_firstTrue) self.branch_seir nn.Linear(1, 8) self.fc_out nn.Linear(hidden_size 8, 1) def forward(self, x_real, x_seir): # x_real: (batch, window, 1) lstm_out, _ self.lstm(x_real) last lstm_out[:, -1, :] # 最后时刻隐状态 seir_feat torch.relu(self.branch_seir(x_seir)) # (batch, 8) merged torch.cat([last, seir_feat], dim1) return self.fc_out(merged)逻辑说明x_real的形状是(batch, window, 1)LSTM会输出每个时刻的隐状态我们只取最后一步last把它当作整段历史信息的压缩。SEIR分支用一层Linear加ReLU把单个数值转成8维向量这样SEIR判断可以作为外部条件参与预测。最后拼接后过一层Linear输出一个数就是残差的预测值。参数说明hidden_size32 对几百个样本的小数据集足够num_layers1 是默认首选加深LSTM几乎必过拟合。如果样本超过500个点可以把hidden_size加到64。branch_seir的输出维度我固定为8调它意义不大真正要调的是LSTM那部分。3.3 训练细节损失函数、优化器、早停import torch.optim as optim model SEIRLSTM(input_size1, hidden_size32, num_layers1) optimizer optim.Adam(model.parameters(), lr1e-3) loss_fn nn.MSELoss() # 按时间顺序切分前70%训练后30%测试 train_end int(len(X_real) * 0.7) X_real_t torch.tensor(X_real[:train_end]).unsqueeze(-1) X_seir_t torch.tensor(X_seir[:train_end]).unsqueeze(-1) y_t torch.tensor(y[:train_end]).unsqueeze(-1) # 归一化只用训练集的统计量 mu, std X_real_t.mean(), X_real_t.std() X_real_t (X_real_t - mu) / std y_t (y_t - mu) / std for epoch in range(50): model.train() optimizer.zero_grad() pred model(X_real_t, X_seir_t) loss loss_fn(pred, y_t) loss.backward() optimizer.step() if epoch % 10 0: print(fepoch {epoch}, loss {loss.item():.4f})逻辑说明两个关键点。一训练/测试必须按时间切不能随机打乱否则测试集的信息会通过乱序泄漏进训练。二归一化必须在切分之后只对训练集统计mu和std预测测试集时沿用训练集的统计量这是最容易出错的数据泄漏源。参数说明lr1e-3 是Adam的默认档数据量小时不用学习率调度器。50轮是起步值实际看loss曲线连续5轮不降就停。注意X_real_t.std()如果为0要防除零疫情数据一般不会但早期连续多天零新增时可能出现。损失函数用MSELoss因为残差是连续值且我们希望大误差被重点惩罚。4. 让SEIR-LSTM跑得更准参数调优与验证方法4.1 SEIR参数先行β、γ、σ怎么用最小二乘去拟合直接用默认参数预测很容易被真实数据甩开所以要先拿前N天数据拟合SEIR让基线贴近真实序列。from scipy.optimize import least_squares def seir_error(params, real_new, N, E0, I0): sim seir_simulate(params, len(real_new), N, E0, I0) seir_new np.diff(sim[:, 2]) seir_new seir_new[:len(real_new)] # 对齐长度 return seir_new - real_new # 用前30天拟合 best least_squares( seir_error, x0[0.5, 1/5.2, 1/14.0], args(real_new[:30], 10_000_000, 10, 5), bounds([0.01, 1/20, 1/30], [2.0, 1/1.0, 1/3.0]) ) beta, sigma, gamma best.x r0 beta / gamma print(f拟合结果: beta{beta:.3f}, sigma{sigma:.4f}, gamma{gamma:.4f}, R0{r0:.2f})逻辑说明least_squares 默认用Levenberg-Marquardt对3个参数收敛很快。误差函数直接返回“模拟值-真实值”让优化器在真值附近找参数。bounds 限制了物理合理范围比如β不能超过2γ对应的感染周期不能短于3天否则拟合出一个R020的模型曲线形状一定会翻车。参数说明x0 用常见文献值。建议多跑几组初值比如 [0.3, 1/7, 1/10] 和 [1.0, 1/3, 1/20]取残差平方和最小的一组。E0和I0也可以放进拟合变量里但那样自由度变高容易过拟合。我一般先固定E0和I0拟合三个核心参数如果曲线早期对不上再放开初始条件。4.2 LSTM超参数怎么配一张参数表和一组默认值我的经验是先调SEIR、后调LSTM不能反过来。SEIR基线不准LSTM会学出一套专门抵消SEIR错误的复杂逻辑一旦SEIR参数更新这套逻辑就作废了。常用参数范围如下参数默认值调试范围说明window147 ~ 21输入历史天数短数据用7hidden_size3216 ~ 64隐状态维度小数据用16num_layers11 ~ 22层以上极易过拟合lr1e-31e-4 ~ 3e-3数据少时1e-4更稳train_ratio0.70.6 ~ 0.8按时间切分不能用随机epochs5030 ~ 100必须配合早停batch_size3216 ~ 64样本少时可直接全批量表格里这些值是从几百个短时间序列案例里归纳的不是玄学。如果数据只有200天window14 会让有效样本不到200个此时hidden_size16、batch_size16 更稳。反过来数据超过500天window可以放大到21hidden_size上到64。另外一个实用流程先用默认参数跑一遍观察训练loss和测试loss的差距。如果训练loss很低、测试loss很高是过拟合优先降hidden_size或加早停如果两个loss都高先怀疑SEIR基线没拟合好再怀疑window太短。4.3 验证预测效果RMSE、MAPE和峰值误差import numpy as np def evaluate(y_true, y_pred): rmse np.sqrt(np.mean((y_true - y_pred) ** 2)) mape np.mean(np.abs((y_true - y_pred) / (y_true 1e-8))) * 100 peak_err (np.max(y_pred) - np.max(y_true)) / np.max(y_true) * 100 return rmse, mape, peak_err # y_pred_test来自测试集上的预测结果残差SEIR基线 y_pred_test residual_pred seir_baseline_test rmse, mape, peak_err evaluate(y_true_test, y_pred_test) print(fRMSE{rmse:.1f}, MAPE{mape:.2f}%, PeakErr{peak_err:.1f}%)逻辑说明RMSE量纲和真实值一样用来判断整体偏差MAPE是百分比但真实值接近0时会爆炸所以分母加1e-8保护。峰值误差是传染病预测里最该看的指标因为峰值决定了医疗资源峰值需求比平均误差更有业务意义。参数说明上面代码假设已经拿到了测试集的残差预测值。实际做多步预测时不能一次到位要用递归预测先把窗口往后滑每步生成一个残差再拼上SEIR基线得到最终预测。递归超过7天后误差会快速累积建议在验证时只报1、3、7天的结果不要无限外推。5. 避坑指南COVID-19预测里五个典型翻车现场这章写的是我复现类似项目时踩过的血泪经验按出现频率排序。前两条几乎每个做传染病预测的人都会遇到后三条属于数据质量问题和模型结构问题需要结合自己的数据去判断。5.1 数据泄漏归一化用了全局均值和方差现象测试集预测曲线离真实值非常近RMSE低到不可思议但把模型放到新数据上一预测就崩。原因在切分训练/测试之前就做了归一化mean和std混入了未来数据模型等于提前看到了测试集的尺度信息。解决把归一化放在切分之后只用训练集计算mu和std测试集和预测阶段沿用同一组统计量。另一个相关坑是特征工程里用了全序列的滑动平均同样会造成泄漏要警惕任何涉及未来时刻的计算。5.2 SEIR参数拟合出离谱R0现象least_squares 收敛到 R08 甚至 R020生成的曲线跟真实数据差异很大。原因早期报告病例数远小于实际感染人数数据不完备加上E0、I0的初值设得完全不靠谱拟合器只能靠极端参数去硬凑。解决把E0和I0也放进优化变量给初始条件一个合理范围或者固定σ和γ用已知潜伏期和感染期只拟合β和初始条件。更稳妥的是用多组初值多次拟合取拟合误差最小且R0落在1.5~6之间的一组。5.3 LSTM预测曲线在后段退化成一条直线现象多步递归预测时预测值逐渐趋于一个常数峰值完全消失曲线像一条被拉平的香肠。原因递归预测时每一步的误差都会作为下一步的输入误差不断累积模型输出往训练集均值回归这是LSTM做长期预测的通病。解决不要递归超过7天。每步都用最新的真实数据或SEIR基线做校准或者干脆只预测1~3天把更长期的输出当作趋势参考而不是正式预测值。另一种做法是在残差预测中引入不确定性区间至少让使用者知道长期预测不可靠。5.4 峰值低估训练集里根本没有这种突变现象模型预测峰值为真实值的60%误差远大于RMSE反映的水平。原因疫情峰值往往是政策干预、毒株变异共同作用的结果训练样本里没有出现过类似量级的突变LSTM学不到没见过的情况。解决把SEIR基线的峰值位置和峰值强度作为外部特征输入LSTM或者把峰值当作异常点处理在残差损失里对峰值区域单独加权。这里没有银弹更合理的心态是承认不确定性给出10%~90%的预测区间而不是追求一个确定值。5.5 日期口径混乱确诊日 vs 报告日现象预测曲线和真实数据总是错开一天某个波峰怎么调都对齐不上。原因很多公开数据的日期是“报告日期”模型输入却是“发病日期”两者存在一天或数天的报告延迟周末还会造成积压。解决统一口径后再用最省事的方法是用7日均线替代日更数据消除周末效应和报告延迟抖动。如果必须用日数据把滞后天数当成一个参数让SEIR去拟合不要手动猜。这个问题在数据清洗阶段解决比在模型阶段补救成本低得多。6. 进阶玩法把SEIR的约束塞进LSTM损失函数上面介绍的残差修正已经能跑但残差自由度太高模型可能学到违背传染病常识的曲线比如预测感染人数突然翻10倍。解决办法是把SEIR的动力学常识写成惩罚项加进损失函数这就是物理约束损失的简化做法。import torch def physics_penalty(seir_pred, lstm_residual): # seir_pred和lstm_residual都是张量 final_pred seir_pred lstm_residual # 感染人数不能为负 neg_penalty torch.relu(-final_pred).square().mean() # 相邻两天变化不能超过SEIR变化的3倍 diff final_pred[1:] - final_pred[:-1] seir_diff seir_pred[1:] - seir_pred[:-1] jump_penalty torch.relu( diff.abs() - 3 * seir_diff.abs().clamp(min1e-6) ).square().mean() return neg_penalty 0.1 * jump_penalty # 训练循环里 loss loss_fn(pred, y_t) 0.1 * physics_penalty(seir_baseline_t, pred)第一个惩罚让预测值保持非负第二个惩罚限制相邻日变化幅度不允模型拍脑袋跳出一个不合理的高峰。0.1是惩罚系数设置太大会让模型干脆不学残差设置太小约束形同虚设一般从0.1开始试观察训练曲线再调。这套做法可以看作“物理信息神经网络”在传染病预测里的轻量落地不需要懂PINN也能用。实际效果是预测曲线更平滑峰值位置更可信代价是对异常波动响应变钝。我个人后来的习惯是先把SEIR基线拟合到残差只有个位数再开LSTM最后加物理惩罚。换数据时先跑一遍完整流程再决定要不要动惩罚系数。整个方案跑下来你会发现最值得反复打磨的不是LSTM层而是数据对齐、切分、归一化和SEIR参数拟合那几十行代码90%的翻车都发生在那里。希望帮到你。本文还有配套的精品资源点击获取
返回列表