资讯动态

时间序列回归预测实战:用高斯过程回归建模疫情新增病例

发布时间:2026/9/26 4:47:57 来源:尧图企业网站定制
前阵子整理旧项目翻到年初写的新冠回归预测代码顺手把整个流程重新跑了一遍。这个项目的内容很朴素拿公开的疫情历史数据用回归模型预测未来新增确诊数顺便把预测区间画出来。代码不算长但里面有大量和时间序列有关的问题比如数据顺序、滞后特征、时间序列交叉验证随便踩一个就能让结果看起来离谱。这篇就按我实际动手的顺序把项目从数据准备到模型训练、再到结果解读完整讲一遍适合刚接触回归预测、想用真实数据练手的朋友。看完你至少能照着自己写一套可运行的版本而不是停留在把sklearn的demo跑通。1. 回归预测项目的整体设计思路1.1 预测目标怎么定不要一上来就预测原始每日新增最开始我图省事直接把“每日新增确诊”当作回归目标结果模型训练得再好预测曲线也是锯齿状和真实数据对的乱七八糟。问题不在算法而在目标变量本身太脏。新冠疫情数据有很明显的报告周期效应周末和节假日通报量会明显偏低周一和节后第一天会堆积补报导致原始序列出现周期性尖刺。模型如果直接学这个序列它会花大量容量去拟合“周一反弹”这种统计噪声而不是真正影响疫情扩散的长期趋势。我最终把预测目标改成了“7日中心移动平均后的每日新增确诊数”。这一步相当于给信号加了一个低通滤波器把短周期的报告抖动抹平让模型集中去学趋势变化。原始序列仍然保留用于画图对比但喂给回归模型的一律是平滑值。数据源我用的是公开的OWID数据集里面按国家/地区提供了每日的 new_cases、new_deaths 字段。为了避免引入过多周期性扰动我只取了 location World 的全球汇总行这样得到的是一条连续时间线不需要再考虑各地区合并的问题。如果你拿到的是按地区分开的文件注意先用 groupby(date).sum() 汇总再去做移动平均否则 rolling 运算会跨地区串数据。还有一点我一开始也试过对目标做 log1p 对数变换让那些爆发期的大数字不要主导损失函数。但高斯过程回归预测的是对数空间里的均值和方差还原到原始尺度之后置信区间会变成不对称形状处理起来比较麻烦。如果你是第一次做这个项目建议先不取对数靠7日均值配合 GPR 的噪声项来吸收波动后期再慢慢优化。1.2 模型选型逻辑为什么是线性回归配合高斯过程回归回归模型的选择很多线性回归、多项式回归、岭回归、随机森林、XGBoost 都能做但我在这个项目里最后保留了两类线性回归作为基线高斯过程回归Gaussian Process RegressionGPR作为进阶方案。线性回归是必须有的。它的作用不是直接产出最佳结果而是当做一个“下限”参考。如果 GPR 跑到最后连线性回归都不如那大概率是数据预处理出了差错比如滞后特征没有对齐、训练集和测试集混在一起、或者特征尺度差太大。先把线性回归结果打出来心里就有底了。GPR 的引入理由很直接疫情历史数据本质上是小样本非线性轨迹整条时间线可能有上千个点但有效的“拐点”和“阶段性特征”很少。GPR 在小样本场景下表现非常稳定而且直接输出预测均值和方差不需要像随机森林那样用袋外偏差去近似不确定性也不像 XGBoost 那样需要调一堆树参数。它特别适合“样本量不大但希望预测结果带置信区间”的任务。可能有人会想到在 MATLAB 里用相关向量机RVM做多输出回归。RVM 的优势是可以在小样本下得到稀疏解概率输出和 GPR 共享部分贝叶斯思想。但 scikit-learn 没有原生 RVM 实现自己从零撸一个成本不低。所以我这个项目的正式代码里没有用 RVM只在对比讨论了它的适用性。如果你恰好是 MATLAB 用户又需要多输出回归RVM 完全可以尝试否则直接用 Python 的 GPR 会让整个流程更顺畅。1.3 开发环境和项目结构整个项目用 Python 3.9 常见数据科学库就能跑通不涉及复杂的工程框架。我实测过的版本是pandas 2.1、numpy 1.26、scikit-learn 1.3、matplotlib 3.8。安装直接用一行命令pip install pandas numpy matplotlib scikit-learn代码组织上我在 Jupyter Notebook 里按流程分成四块数据加载与预处理特征构造模型训练与评估预测未来与可视化如果你更习惯脚本工程也可以拆成 load_data.py、features.py、train.py 三个文件。对于这种规模的项目我觉得 Notebook 更方便因为每一步中间结果都能及时看到排查特征错位时尤其好用。2. 核心细节解析特征工程比模型更影响结果2.1 时间特征、滞后特征与移动平均回归模型不认识日期所以第一步必须把日期转换成机器能用的数字特征。最基础的是“当前日期是时间序列的第几天”直接用日期减去最小值得到天数差。这个特征让模型意识到时间在前进也是预测未来趋势最主要的线性信号。但光有天数不够。疫情的传播有明显自相关性今天的新增数字很大程度上取决于一周前和两周前的情况因为感染、发病、报告之间存在时间滞后。我在项目里构造了滞后特征world[day_idx] (world[date] - world[date].min()).dt.days world[target] world[new_cases].rolling(7, centerTrue).mean() world[lag_7] world[target].shift(7) world[lag_14] world[target].shift(14)这里需要注意target 用的是带 centerTrue 的7日均值也就是说第 t 天的 target 实际上是 t-3 到 t3 天的均值。它会引入很小的“未来信息”但对于训练和测试评估来说只要训练集和测试集按照相同方式切分这个问题影响不大。真正要留意的是预测未来时滞后特征 lag_7 和 lag_14 需要用到 target 的历史值而在未来段我们是没有真实数据的。我的处理方式是用最后一个已知的 target 值临时填充如果只预测未来14天内误差尚可接受更严格的滚动预测需要把 GPR 预测出来的均值回填到滞后特征里再继续预测。我还试过加入“星期几”的 one-hot 特征但数据已经做过7日均值平滑周期性报告信号基本被抹平了加入之后没有带来明显提升反而增加了特征维度后来就去掉了。2.2 用 TimeSeriesSplit 做时间序列交叉验证这是很多初学者踩得最惨的坑。普通的 KFold 交叉验证会随机打乱样本顺序把早先的数据混到验证集里让模型提前“看到”未来结果验证分数特别漂亮但一旦真正预测未来就崩掉。时间序列预测必须保住时间顺序。sklearn 提供了现成的 TimeSeriesSplit它保证训练集永远是验证集早期的时间片段。我用它来做模型评估from sklearn.model_selection import TimeSeriesSplit tscv TimeSeriesSplit(n_splits5, gap7, test_size14) for train_index, test_index in tscv.split(X): X_train, X_test X[train_index], X[test_index] y_train, y_test y[train_index], y[test_index]gap 参数表示训练集最后一天和验证集第一天之间要留出多少天间隔我设置为7天尽量模拟真实的预测场景避免相邻日期的强相关性让验证结果虚高。2.3 数据清洗中容易忽略的问题疫情数据的质量问题比想象中严重。常见的包括某一天数据缺失被补成0某天突然出现负值某些国家在改口径后一次性补报大量历史病例。这些都是原始序列中的“离群点”。我的处理顺序是这样先 clip 负值归零再计算7日均值。对于个别历史极端补报我选择保留不强行剔除。GPR 本身有一个 WhiteKernel 噪声项会把一部分离群值当作噪声处理效果比线性回归稳健得多。如果你发现模型被某个尖峰带偏可以检查该时间点前后14天的数据决定是否要把 target 直接置为前后两天的均值。另外训练集的尾端数据影响最大。因为近期样本直接决定了近期预测的走向如果你从完整数据里随机丢掉尾部几个点GPR 预测未来7天的效果会明显变差。切训练集时我最看重最后30天的样本是否干净。3. 实操过程从线性回归到高斯过程回归的完整实现3.1 数据加载与预处理把CSV变成训练矩阵这个项目完整流程我整理成了可以直接跑的代码。第一步是读取数据保留全球汇总行并构造特征和目标列import pandas as pd import numpy as np import matplotlib.pyplot as plt from sklearn.linear_model import LinearRegression from sklearn.gaussian_process import GaussianProcessRegressor from sklearn.gaussian_process.kernels import RBF, WhiteKernel, ConstantKernel from sklearn.model_selection import TimeSeriesSplit from sklearn.preprocessing import StandardScaler from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score df pd.read_csv(owid-covid-data.csv, parse_dates[date]) world df[df[location] World].sort_values(date).reset_index(dropTrue) world[new_cases] world[new_cases].clip(lower0) world[target] world[new_cases].rolling(7, centerTrue).mean() world[day_idx] (world[date] - world[date].min()).dt.days world[lag_7] world[target].shift(7) world[lag_14] world[target].shift(14) train world.dropna(subset[target, lag_7, lag_14]).copy() X train[[day_idx, lag_7, lag_14]].values y train[target].values注意 lag_7 和 lag_14 构造完以后数据前14行会因为 shift 产生 NaN必须先 dropna否则 sklearn 会直接报错。这一行代码虽然不起眼但能省掉后面大量排查时间。3.2 第一步线性回归基线模型先不急着上复杂模型用线性回归打个底。我这里把最后30天作为验证集前边全部当作训练集X_train, X_test X[:-30], X[-30:] y_train, y_test y[:-30], y[-30:] lr LinearRegression() lr.fit(X_train, y_train) lr_pred lr.predict(X_test) lr_mae mean_absolute_error(y_test, lr_pred) lr_rmse np.sqrt(mean_squared_error(y_test, lr_pred)) lr_r2 r2_score(y_test, lr_pred) print(fLinearRegression MAE{lr_mae:.2f} RMSE{lr_rmse:.2f} R2{lr_r2:.4f})在我手里的这份世界数据上线性回归的 R² 能做到 0.8 左右看起来不错但画出曲线就会发现它的问题整体趋势跟住了但遇到波动拐点时反应慢预测值经常比真实值矮一截或者多冲出一段。线性模型没法应对疫情传播的阶段性突变这就是需要非线性模型的原因。3.3 第二步高斯过程回归带置信区间GPR 的核心是核函数。我用了“常数核 RBF 核 白噪声核”的组合常数核 ConstantKernel 用来控制整体方差。RBF 核用来衡量相邻日期的相似度length_scale 控制了“多远算相关”。WhiteKernel 用于捕获观测噪声对应数据本身的随机波动。kernel 1.0 * RBF(length_scale30.0, length_scale_bounds(1e-2, 1e4)) \ WhiteKernel(noise_level1.0, noise_level_bounds(1e-3, 1e3)) scaler StandardScaler() X_train_scaled scaler.fit_transform(X_train) X_test_scaled scaler.transform(X_test) gpr GaussianProcessRegressor( kernelkernel, alpha1e-4, normalize_yTrue, n_restarts_optimizer5, random_state42 ) gpr.fit(X_train_scaled, y_train) gpr_pred, gpr_std gpr.predict(X_test_scaled, return_stdTrue)为什么先做 StandardScaler因为 day_idx 的数值是几百甚至上千而 lag_7 是几万级别的病例数两个特征尺度差距太大。GPR 对特征尺度很敏感如果不缩放length_scale 会被某个大尺度特征带偏。LinearRegression 不做缩放也能跑但 GPR 建议统一缩放。预测时我把 gpr_std 作为置信区间边界plt.figure(figsize(12, 6)) plt.plot(X_test[:, 0], y_test, labelactual, colorblack, linewidth2) plt.plot(X_test[:, 0], gpr_pred, labelGPR pred, colortab:blue) plt.fill_between( X_test[:, 0], gpr_pred - 1.96 * gpr_std, gpr_pred 1.96 * gpr_std, colortab:blue, alpha0.2, label95% interval ) plt.legend() plt.title(Gaussian Process Regression on New Cases) plt.show()我这里用了 1.96 倍标准差近似95%置信区间。真实情况是 GPR 假设高斯后验分布所以这个近似是合理的。画出来的区间宽度很有意思在数据密集的地方区间窄在预测远期或者数据稀疏的地方区间明显变宽这是一个非常直观的“模型有多少信心”的信号。3.4 模型对比与结果解读在我跑的这条全球汇总数据上最后30天验证结果大致是线性回归 MAE 在2万左右RMSE 在3万左右R² 只有0.8出头GPR 的 MAE 降到约1.2万RMSE约1.8万R² 接近0.94。GPR 的提升主要来自两个地方一是 RBF 核能捕捉非线性拐点二是预测置信区间让我能看出哪些位置模型不太确定。这里要提醒一下不同时间段跑出来的数字差异非常大。瘟疫数据不是平稳时间序列后面的数据取决于前面怎么传播。如果你用我这段代码下载的数据日期不同结果完全不同。对比模型的效果时不要只看均值指标还要看预测曲线在拐点处的拟合能力和置信区间是否合理。各类模型在这个场景下的特点我整理成一个表模型适用场景输出形式主要风险线性回归快速基线、单调趋势点预测无法拟合突变拐点多项式回归趋势简单的曲线点预测高次项容易过拟合振荡高斯过程回归小样本、非线性、需要不确定性均值方差计算量随样本量增长RVMMATLAB高维稀疏小样本均值概率Python生态不成熟4. 常见问题与排查技巧实录4.1 为什么预测结果看着像一条直线这个问题我遇到过当时差点把 GPR 丢掉。后来发现最核心的原因是 RBF 核的初始 length_scale 设得太大。length_scale 太大时核函数觉得所有样本之间都差不多最终学出来的函数非常平滑几乎退化成直线。解决办法是把初始长度设成和特征量级匹配的值。比如时间跨度是1500天设成30到60比较合理如果你已经标准化特征length_scale 的初始值可以设为1左右让优化器自己调。另外特征没有标准化也会导致类似问题。GPR 的结果依赖特征距离计算一组特征动辄上千、另一组特征只有几十距离基本被大尺度特征统治核函数的 length_scale 很难同时适应所有维度。先用 StandardScaler 把特征压到均值0、方差1再训练 GPR是成本最低的优化手段。4.2 置信区间过窄或过宽如何调整置信区间太窄通常是因为 WhiteKernel 的噪声水平没有被优化起来或者目标变量几乎没有随机波动。你可以调大初始化 noise_level比如从1改到10然后观察训练后核函数的 noise_level 最终值。另外一个常见原因是没有设置 normalize_yTrue。目标变量如果不在0均值附近GPR 先验假设会出问题导致方差被严重低估。这个参数只需要一行代码却能带来肉眼可见的改善。区间太宽多半是训练数据里有几个特别极端的离群点。GPR 会把它们解释成高噪声导致 WhiteKernel 的噪声水平调得很大整体预测区间被撑开。我的处理方式是画训练数据的分布图找到极端尖峰把它们手动替换成附近7天的均值再重新训练。注意不要随便删除因为时间序列的连续性要求数据点位置不能空缺用平滑值替换比删除更安全。4.3 滞后特征错位的隐蔽问题滞后特征 out-by-one 错误非常隐蔽。我一开始在构造 shift 之前忘了 sort_values结果数据顺序还是原始的按照地区汇总排列的顺序shift(7) 取到的根本不是7天前而是7行以前整个特征全是乱的。排查方式很简单print(train[[date, target, lag_7]].tail(20))人工看几行确认 lag_7 是否等于7天前的 target。这个问题很基础但一旦发生模型效果会变得非常奇怪而你不会想到是 shift 用错了地方。所以我的建议是每次构造完滞后特征先打印检查再往下走。另外TimeSeriesSplit 需要注意 test_size 和 gap 的单位是“样本数”不是“天数”。如果你的数据做了移平均清洗、中间又有缺失行那么样本间隔和时间间隔可能不一致。稳妥的做法是先连续按天补全索引再做特征和切分。4.4 实操问题排查速查表最后把项目里容易踩的坑整理成一张速查表方便你回头排查现象可能原因解决办法RMSE高得离谱目标没平滑、极端尖峰未处理先做7日均值极端点用局部均值替换验证集结果虚高普通KFold随机打乱时间改用TimeSeriesSplit并设置gaplag列出现大量NaNshift后没有dropna删除前14行或对NaN行过滤预测未来单调下坠用last value填滞后特征改用滚动填充或只预测短期步长GPR预测像直线length_scale初始值过大或特征未标准化调小初始值并先做standard scaling置信区间窄得离谱忘记normalize_yTrue加上normalize_yTrue检查噪声核做这个项目最大的体会是传染病数据回归预测难点多半不在算法而在目标定义和数据时间对齐。你只要把“预测的东西到底是什么、特征来自哪一天”写清楚模型用常见库就能跑出像样的效果。我自己现在做这类小项目已经养成一个习惯先花半小时把时间线画出来确认移动平均和滞后特征没问题再开始调模型。顺序对了后面就顺了。

读完文章,也想定制专属网站?

尧图设计师 24 小时内与您沟通定制方案

免费获取报价 →
↑