简介一套基于Python和R的铁路货运量/客运量时序建模预测项目面向交通数据分析与工程应用。代码提供Python 3.7与R 4.0.2双入口覆盖简单线性、STL分解、Holt-Winters阻尼季节性、ARMA/ARIMA等经典时序模型并结合真实运输量数据完成建模与评估便于横向对比不同方法的预测效果。资源包共50个文件、约2.05MB其中包含17个Python脚本、4个R脚本、19张结果图以及csv/xlsx格式的运输量、汇率等数据集另附评估脚本和说明文档可支撑从数据读取到模型评估的完整流程。已有250人学习/浏览。适合需要掌握时序预测建模流程或对铁路货运量进行量化分析的研究者、数据分析师与开发人员借助该项目可直接运行模型代码并使用配套数据减少从零搭建与调参时间还可作为其他交通指标预测的基线参考。1. 把铁路货运量预测这件事做扎实从拿到数据集到跑通四个时序模型这篇笔记拆的是一个很务实的项目基于 Python 和 R 对铁路货运量、客运量做时序建模与预测压缩包里除了源代码还带着可直接用的数据集。源码结构很清晰Python 3.7 环境下run.py做入口model_evalute.py做模型评估算法目录里放了简单线性模型、STL、Holt-Winters 阻尼季节性模型和 ARMA/ARIMA 四个模型R 4.0.2 下也有对应的run.R入口。也就是说Python 和 R 各有一套完整流程拿来改改就能用到别的运输量、销量或流量预测场景。适合的对象包括做交通/物流方向课设的学生刚入手时序预测想找现成 baseline 的从业者以及想对比 R 与 Python 建模差异的工程师。这套代码最大的价值不是“跑通”而是四种模型放在同一份数据上做横向评估选型和取舍一目了然。2. 先摸数据再说模型运输量时序的预处理、差分与平稳性检验2.1 数据文件里到底有什么分清原始数据、规范化结果与额外数据压缩包解压后进入times_series-master/data目录能看到四类数据文件原始的铁路运输量 Excel、用于建模的规范化时序数据、评估指标用的结果数据还有一个附带的外部 CSV 数据。先用 pandas 快速读一遍弄清每个文件的行列规模和数据粒度这一步决定了后面所有模型的时间索引怎么设置。import pandas as pd # 读原始运输量数据 raw_df pd.read_excel(times_series-master/data/运输量.xlsx, sheet_nameNone) print(f工作表列表: {list(raw_df.keys())}) # 读规范化时序数据确认字段名 ts_df pd.read_excel(times_series-master/data/时序建模.xlsx) print(ts_df.head()) print(ts_df.columns.tolist()) print(ts_df.dtypes)这段代码里sheet_nameNone会把所有工作表一次性读进来能快速看出每个 sheet 里是什么指标时序建模.xlsx的列结构决定后面用什么字段做序列、什么字段做时间索引。我一般会在这一步顺手检查有没有缺失值和重复索引因为 STL 和 ARIMA 都要求连续且无重复的时间点缺失不处理后面全是坑。2.2 时间序列切分不能“乱劈柴”为什么不能直接用 train_test_split很多初学者习惯直接把 sklearn 的train_test_split用在时序数据上这是时序建模最典型的翻车姿势。随机切分会把未来的数据混进训练集模型学到的是“偷看答案”的模式评估结果虚高一上真实场景立刻原形毕露。train_size int(len(ts_df) * 0.8) train, test ts_df.iloc[:train_size], ts_df.iloc[train_size:] print(f训练集区间: {train.index[0]} ~ {train.index[-1]}, 样本数 {len(train)}) print(f测试集区间: {test.index[0]} ~ {test.index[-1]}, 样本数 {len(test)})常见的做法是按住时间顺序做 80/20 的前后切分严格保证训练集全部早于测试集。这份代码的思路同理用前 80% 的数据训练四类模型在最后 20% 的真实未来数据上做对比评估。这里为什么要多写一步输出索引区间因为确认切分点没有跨越季节周期很重要——如果数据是月度粒度80% 切分恰好把 12 月的冬季数据全部甩进测试集那模型在本该预测春运高峰的点上天然吃亏评估误差会虚高。2.3 差分与 ADF 检验判断序列是否满足平稳性要求ARMA 和 ARIMA 的适用前提是序列平稳或经过差分后平稳而 STL 和 Holt-Winters 主要靠趋势和季节项做分解对平稳性要求相对宽松。代码目录里的line_model.py、arma_arima_model.py等文件已经封装好了完整流程但自己手动做一遍差分和 ADF 检验能看懂模型选型背后的逻辑。from statsmodels.tsa.stattools import adfuller import numpy as np # 假设 ts_series 是货运量月度序列 d1 np.diff(ts_series, n1) adf_orig adfuller(ts_series, autolagAIC) adf_diff1 adfuller(d1, autolagAIC) print(f原始序列 ADF p值: {adf_orig[1]:.6f}) print(f一阶差分 ADF p值: {adf_diff1[1]:.6f})adfuller返回的第一个值是统计量第二个值是 p 值判断标准是 p 值小于 0.05 才拒绝“存在单位根”的原假设也就是序列平稳。这套代码里 ARIMA 的 d 参数就来自这一步的分析结果原始序列 ADF 检验不通过就做一阶差分再看差分后序列是否显著。跑完 ADF 后去翻arma_arima_model.py里的order设置你能看到模型选阶和这一步分析是对应的不是随手填的。2.4 季节性趋势分解从 STL 视角理解运输量的波动来源stl_model.py用 STL 做趋势和季节分解把原始序列拆成趋势项、季节项和残差项。铁路货运量和客运量的季节模式差异很大——货运受工业生产周期影响客运受节假日和学生寒暑假影响两者在同一模型里的季节分量明显不同。from statsmodels.tsa.seasonal import STL stl STL(ts_series, period12, robustTrue).fit() trend stl.trend seasonal stl.seasonal resid stl.resid # 看残差的标准差残差越小说明趋势季节分解越充分 print(f残差标准差: {resid.std():.4f})period12表示月度数据的年度季节性调整成 4 就是季度周期7 就是周周期。robustTrue会让分解对离群点不那么敏感运输量数据经常受疫情、政策调整等外生冲击用 robust 版本能减少个别极端值对趋势项的拉扯。跑完 STL 再看序列分解图能直接解释“为什么线性模型在这个数据集上表现有限”——运输量数据的趋势不是直线而是带明显季节摆动的曲线。3. 四个模型逐个拆解线性、STL、阻尼 Holt-Winters 与 ARIMA 的选型逻辑3.1 简单线性模型只能做趋势 baseline别指望它处理季节波动line_model.py实现的是对时间索引做线性回归本质上是 y a * t b 的结构只拟合递减或递增的整体趋势。这个模型在四个模型里担任 baseline 的角色用来衡量“引入季节项和自回归项到底提升了多少精度”。import numpy as np from sklearn.linear_model import LinearRegression t np.arange(len(ts_series)).reshape(-1, 1) model LinearRegression() model.fit(t, ts_series) trend_line model.predict(t) print(f斜率: {model.coef_[0]:.4f}, 截距: {model.intercept_:.4f})把t从 0 开始按顺序递增编码拟合出来的斜率代表每个时间步的平均增量。这个模型的残差里通常带着明显的季节摆幅这是正常的——它的作用就是给后面三个带季节能力的模型做参照。如果线性模型的误差和 STL/ARIMA 差不多说明这份数据的季节波动很小那选型思路就不一样了。3.2 STL 预测分解之后往前延拓趋势与季节分量STL 模型做预测的思路和直接对原始序列做拟合不同先把序列分解成趋势 季节 残差然后分别预测趋势部分和季节部分最后加总还原。stl_model.py里的做法是这个思路的标准封装。from statsmodels.tsa.holtwinters import ExponentialSmoothing # 对趋势项做二次指数平滑预测 trend_model ExponentialSmoothing( trend, trendadd, seasonalNone, initialization_methodestimated ).fit() # 季节项直接用最后一年同月均值外推或用 seasonal_naive forecast_trend trend_model.forecast(6) forecast_seasonal seasonal[-12:][:6] forecast forecast_trend forecast_seasonal.values print(f未来 6 期 STL 预测值: {forecast})常见做法是趋势项用简单指数平滑或线性外推季节项取近期季节分量继续沿用。这里我一般是把季节性当成确定的年度模式来外推不额外建模。对铁路客运量这种以年为周期、节假日导致波峰相对固定的数据STL 这种先分解再组合的预测方式对中期趋势变化的捕捉比单纯 ARIMA 更直观可控。3.3 阻尼 Holt-Winters给指数平滑加一道“刹车”Holt-Winters 在指数平滑里同时处理水平项、趋势项和季节项。damped 版本是在趋势项上乘一个阻尼系数 φ让趋势随预测步长增加逐渐衰减到水平值避免长期预测出现线性趋势无限外推的失控状态。holt_winters_damp_model.py这个名字里的 “damp” 就是这个作用。from statsmodels.tsa.holtwinters import ExponentialSmoothing damp_model ExponentialSmoothing( ts_series, trendadd, seasonaladd, seasonal_periods12, damped_trendTrue, initialization_methodestimated, ).fit() forecast_damp damp_model.forecast(6) print(forecast_damp)trendadd是加法趋势seasonaladd是加法季节项。货运量序列如果波动幅度随水平值增大而变大可以改成mul这需要看残差方差是否随趋势递增。damped_trendTrue是这份实现和传统 Holt-Winters 的核心差异——预测 12 个月时普通版本会按当前趋势一直冲上去阻尼版本到了第 6 个月后趋势增量明显变小两种预测曲线的分叉在长期预测中极其明显。3.4 ARIMA 与 ARMA从差分视角理解 d 参数以及节假日冲击的边界arma_arima_model.py同时封装了 ARMA 和 ARIMA 两个模型数据平稳时用 ARMA非平稳时先差分再套 ARMA就变成 ARIMA。这里的代码会基于 AIC/BIC 自动搜索 p、d、q 的范围或者直接读用户指定的 order。跑完第 2 章差分检验再看这个文件能明白 d 的取值不是拍脑袋的。from statsmodels.tsa.arima.model import ARIMA # 对一阶差分后序列拟合 ARMA(2,1) model ARIMA(ts_diff1, order(2, 0, 1)) result model.fit() print(result.summary()) print(fAIC: {result.aic:.2f}, BIC: {result.bic:.2f})order(2, 0, 1)中 p2 是自回归阶数d0 表示输入序列已经是平稳序列q1 是移动平均阶数。如果直接用原始序列建模d 要设成差分阶数。需要留意的是ARIMA 本质上靠序列自身的滞后相关性做预测它不感知“春节”“国庆”这类外部事件——如果某个月份的客运量因为假期调休出现异常尖峰而历史同期没有完全一致的模式ARIMA 的拟合残差会特别大。这种情况下引入外生变量的 SARIMAX 或干脆把节假日特征作为附加回归项是进阶方向。3.5 四个模型的适用边界与参数对比模型强项弱项关键参数在本项目中的预期定位简单线性整体趋势判断不含季节项残差大仅斜率截距baseline被其他模型吊打STL趋势/季节/残差可解释性强季节外推依赖历史同期模式period、robust可视化与解释优先阻尼 Holt-Winters加趋势衰减中短期预测稳季节模式固定变化剧烈时不敏感trend、seasonal、seasonal_periods、damped_trend月度运量的主力模型ARIMA/ARMA自回归结构捕捉惯性对季节性需手动差分/扩展order(p,d,q)平稳时序的标准做法这四组模型在evalute.xlsx里对应一批评估指标。跑完model_evalute.py后你会看到 MAE、RMSE 和 MAPE 三者一起比较不会只拿 RMSE 一个指标下结论——因为 RMSE 对个别异常月份极其敏感MAPE 对低基数月份的误差又会被放大。4. 把 run.py 到 model_evalute.py 整条链路跑通训练、预测、评估与结果解读4.1 从入口文件开始run.py 的流程和数据流向项目入口是run.py它做的事情是读入时序建模.xlsx→ 切分训练测试集 → 依次调用四个模型训练并预测 → 调用model_evalute.py计算评估指标 → 把预测值和真实值写入evalute.xlsx。这个结构很适合做两件事一是加自己写的模型比如 Prophet 或 LSTM进目录做横向比较二是换数据文件做别的行业时序预测。# 伪代码梳理 run.py 的调用链真实实现以仓库为准 from line_model import line_forecast from stl_model import stl_forecast from holt_winters_damp_model import hw_damp_forecast from arma_arima_model import arima_forecast models { linear: line_forecast, stl: stl_forecast, holtwinters_damp: hw_damp_forecast, arima: arima_forecast, } for name, func in models.items(): pred func(train_series, stepslen(test_series)) # 对比 pred 与 test_series 计算误差指标四个模型的函数签名统一成(train_series, steps)这种形式新增模型只要返回相同长度的 ndarray就能直接接入model_evalute.py的评估流程。另外注意一件事Python 3.7 在 long long 上处理 pandas 1.x 没有任何压力但 statsmodels 版本必须和 Python 版本匹配代码里requirements.txt已经把版本钉死了。4.2 model_evalute.py 的评估指标MAE、RMSE、MAPE 怎么算模型评估文件对每个模型的预测结果计算统一的误差指标。我一般不只是看数字大小还会把误差按月份拆开看模型到底在哪些时段集中翻车——这份代码的评估输出帮我把这个检查做得很快。import numpy as np def evaluate(y_true, y_pred): mae np.mean(np.abs(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)) * 100 return {MAE: mae, RMSE: rmse, MAPE: mape} # 逐月看误差定位集中翻车的时段 errors y_true - y_pred for idx, err in zip(test_index, errors): if abs(err) np.mean(abs(errors)) 2 * np.std(errors): print(f异常误差点位: {idx}, 误差 {err:.2f})单看 MAE 是拿“平均表现”说事RMSE 拉高了大幅偏差的权重MAPE 则对绝对值小的月份格外苛刻。用“平均值加两倍标准差”的方式筛异常误差点能找到具体月份比如每年 2 月春运峰值明显超出模型估计这代表模型的季节分量没有完全捕捉节假日加成效应。4.3 评估结果的合理解读方式别只看排名看差距与残差结构四个模型的评估结果写进evalute.xlsx后我建议做两个额外分析。一是看各模型评估指标的差距是否显著——如果 Holt-Winters 的 RMSE 只比线性模型低 5%那说明数据本身的季节规律很弱换模型解决的提升空间有限。二是看残差的滞后相关性对最优模型的残差再做一次 Ljung-Box 检验如果残差仍然存在显著自相关说明模型的时序信息还没提干净不是误差指标合格就能收工。from statsmodels.stats.diagnostic import acorr_ljungbox lb_test acorr_ljungbox(residuals, lags[6, 12], return_dfTrue) print(lb_test) # 如果 p 值 0.05残差还有自相关模型需要加阶数或换结构这里lags分别测 6 期和 12 期的滞后相关性p 值越小代表残差里残留的时间结构越强。如果最优模型的 Ljung-Box p 值低于 0.05我的习惯是回到arma_arima_model.py调高 p 或 q 的搜索范围重新跑一遍而不是直接接受现有模型。5. 这五个坑我踩过时序预测项目里最常见的翻车现场与排查路径5.1 一跑 model_evalute.py 就报维度不匹配现象测试集真实值长度是 24模型输出的预测值长度却是 23 或者 25直接报ValueError: operands could not be broadcast together。原因切分之后没有统一处理索引位移。forecast(steps)的长度选项不一致或者某个模型内部把最后一个观测值当作预测起点的基础上多预测了一期。解决在run.py的循环里加一个长度断言提前暴露问题检查 ARIMA 的forecast(steps)是否传入了和测试集完全一致的步数。assert len(pred) len(test_series), f长度不一致: {len(pred)} vs {len(test_series)}这个断言解决了一大部分排序混乱因为错误提示直接就把模型名和长度差打在控制台上了。5.2 原始数据里有缺失月份STL 直接抛错现象STL 在 statsmodels 里不接受缺失值运行stl_model.py时报ValueError: Naive seasonal decomp requires at least 2 complete cycles。原因运输量.xlsx里的某些月份没有记录。解决在 2.1 节的数据探索里顺手检查时间索引的连续性发现缺失就做前向填充或插值。ts_df ts_df.set_index(month).asfreq(MS) ts_df[value] ts_df[value].interpolate(methodlinear)asfreq(MS)把时间索引规范到每月开头插值则会对缺失位置用线性补全。千万别用dropna()删行时序模型对时间连续性极其敏感删掉一行就等于把相邻两个点的“间隔”从 1 个月拉伸成 2 个月季节周期全乱掉。5.3 ADF 检验 p 值怎么都过不了 0.05现象原始序列和一阶差分后 ADF 的 p 值都高于 0.05不管设order(2,1,1)还是(4,1,2)模型表现都稀烂。原因数据有异常跳变比如疫情初期的货运量探底这类结构性断点会让单位根检验的结论失真。解决先看序列折线图定位异常区间在建模窗口里去掉异常点或用干预变量标记时间充裕的话可以测试只用 2015-2019 的“正常周期”数据建模看模型在正常场景下的表现。5.4 用 sklearn 的 train_test_split 偷偷翻车现象模型评估指标异常地好好到每个模型的 RMSE 都比直觉低一个数量级预测曲线和真实曲线贴得严丝合缝。原因随机切分把未来数据混进了训练集模型“记住”了未来样本的位置和取值。这是时序预测里最害羞的错——指标越好越要警惕。解决一律前 80% 后 20% 切分并且在评估前打印训练集最后一个时间点和测试集第一个时间点确认没有交集。print(f训练集结束: {train_index[-1]}, 测试集开始: {test_index[0]})5.5 Holt-Winters 预测曲线最后变成一条直线现象阻尼 Holt-Winters 的 long-horizon 预测在末期完全平掉趋势项不再增长。原因阻尼系数 φ 的作用被拉满趋势项在 10 期之后被衰减到接近 0。damped_trendTrue的初始值由initialization_methodestimated自动估计如果数据整体趋势微弱估计出来的 φ 就会很小。解决手动给阻尼系数设一个下限或者跑两组对照——一组damped_trendFalse一组 True观察中短期和长期的交叉点在哪里。6. 用 R 再做一遍全流程验证run.R、tsibble 操作与双语言结果对照6.1 R 版本的结构差异与方法对应压缩包里同时有run.R和 R 4.0.2 的依赖文件这意味着同一个建模流程在 R 里也有一套完整实现。R 的时序生态里tsibble负责时间索引和切分fable包负责模型拟合与预测。和 Python 版本的四模型一一对应线性模型对应TSLM()STL 用decomposition_model()或model()加STL()Holt-Winters 对应ETS()里带 damped 的版本ARIMA 对应ARIMA()。library(fable) library(tsibble) ts_data - as_tsibble(data, index month) fit - ts_data %% model( linear TSLM(value ~ trend()), ets_damped ETS(value ~ error(A) trend(Ad) season(A)), arima ARIMA(value, stepwise FALSE) ) fit_forecast - fit %% forecast(h 12)ETS(value ~ error(A) trend(Ad) season(A))这一行准确对应 Python 里的trendadd, seasonaladd, damped_trendTrue三个成分分别用 A 表示加法、d 表示阻尼。R 的fable包最顺手的一点是forecast(h 12)自带置信区间和残差检验同一套流程可以快速验证 Python 侧结果的合理性。6.2 Python 与 R 结果对照的检查方法跑完 R 版本后把 Python 和 R 的预测结果按时间对齐画在一张图上。两个语言版本如果差异在 5% 以内基本可以确认代码链路没问题差异超过 10%优先查两边的缺失值处理方式是否一致其次查季节周期的设定是否都是 12。library(ggplot2) combined - bind_rows( python_pred %% mutate(source Python), r_pred %% mutate(source R) ) combined %% ggplot(aes(x month, y value, color source)) geom_line() labs(title Python vs R: 阻尼 Holt-Winters 预测对照)如果两版预测在小数点后一位才产生分支这趟双跑的最大意义已经完成了一半——说明代码迁移没有改变模型的数学本质。6.3 双语言交叉验证的最大收益从这个项目里拿到的最值钱经验是R 和 Python 各有各的优势但人为制造两边的“时差”没有意义。R 侧的重点是快速检验和可视化Python 侧则更适合接入生产环境做定时预测。从那以后我每次拿类似预测项目都强制走一遍Python 建模 → R 验证 → 对比差异 → 修正数据处理逻辑。这套流程看似多花了一个小时实际上把“数据预处理不一致”这类隐性 bug 掐死在源头。希望帮到你。本文还有配套的精品资源点击获取