从R脚本复现GEFCOM14概率负荷预测：分位数回归与QRA集成-北京尧图网络

简介面向能源负荷预测学习者与R语言实践者的GEFCOM2014-EPFL资料包围绕多地区小时级电力负荷数据展开适合希望用统计模型和机器学习完成时间序列预测练习的中级读者。资料可用于理解负荷数据中的日期、时间、地理位置与能耗量结构并支撑季节性趋势分解、ARIMA与状态空间建模、随机森林等机器学习方法、温度节假日等特征工程以及MSE、RMSE、MAE和R²等指标评估。压缩包约111.53MB文件总数显示为0文件类型明细暂未提供具体构成需以实际包内为准已有1426人学习说明其在能源预测与R实践场景中具有一定关注度。读者可借助该资料熟悉从数据探索、模型比较到预测结果可视化与集成优化的分析链路积累负荷预测项目的建模与评估经验。1. 一份只有 R 脚本的仓库凭什么能复现 GEFCOM14 的概率负荷预测GEFCOM14 的能源负荷预测赛题把「明天 24 小时各用多少电」拆成了 99 个分位点来打分。这个仓库 gefcom2014-epfl 用 R 把整条链路串了起来逐时负荷与多站温度对齐、滞后与温度响应特征、按小时训练的分位数回归、QRA 集成最后落到 pinball 损失上。它真正的价值在于一套能拆开改的流程框架而非一份跑完就出成绩的脚本。数据读取、特征、模型、评分四层彼此独立你想换掉其中任何一层都不至于牵一发动全身。这也意味着它不会替你把参数调好拿到手该调的地方一个都少不了。适合两类人手上已经有几年逐时负荷数据、准备把点预测升级成分位数预测的同行以及刚接触概率预测、想找一个能在 R 里跑通、能对着改的例子的人。下面按数据、模型、集成、踩坑的顺序拆。2. 把逐时负荷和 11 站温度拼成训练表时间对齐与特征构造的三道关负荷预测的模型再花哨喂进去的表不对齐后面全是白搭。这个仓库里最值得抄的部分其实在最前面那两百行把三份原始文件拼成一张既能训练又能回测的宽表并且保证任何一个特征在预测时刻都是「当时真的拿得到」的。这一步做扎实后面换个模型分数就能往上走一截这一步糊弄回测能刷出很漂亮的数字换成真实数据立刻打回原形。2.1 三张表与 24 小时预测窗口原始数据很朴素就三份文件关键列粒度在模型里的角色负荷表时间戳、负荷值小时预测目标同时派生滞后特征温度表时间戳、站点编号、温度小时 × 站点采暖/制冷需求、平滑体感温度假期表日期、假期类型天区分法定假期与普通周末负荷跨度大约四年半逐小时中间有零星缺口。温度是预报值而不是实测值这一点很关键线上预测时你拿不到未来的实测温度只有数值天气预报所以训练阶段也必须用对应的预报序列否则模型学到的是「用未来真值解释现在」的假关系。气象站数量和编号在赛程后半段还调整过代码里千万别把站点列表写死按编号读取再对齐缺的补 NA。任务形态是每一轮只给到预测日零点为止的历史要求输出次日 0 到 23 时、τ 从 0.01 到 0.99 共 99 个分位点的预测值。一天就是 24 × 99 2376 个数而不是 24 个数。很多人第一次上手会本能地去做 24 个点预测然后发现评分函数完全不认——这是概率预测和点预测最直观的分水岭。2.2 时间对齐先把时区和缺口处理干净library(data.table) library(lubridate) # 1) 读入三份原始文件 load_dt - fread(data/load.csv) # 列timestamp, load temp_dt - fread(data/temperature.csv) # 列timestamp, station, temp holi_dt - fread(data/holiday.csv) # 列date, holiday # 2) 统一时间戳并强制 UTC避免本地时区把夏令时那一小时挪位 load_dt[, ts : as.POSIXct(timestamp, tz UTC)] temp_dt[, ts : as.POSIXct(timestamp, tz UTC)] # 3) 温度转宽表一行一小时每个站点一列 temp_dt[, st : paste0(t, station)] temp_wide - dcast(temp_dt, ts ~ st, value.var temp) # 4) 以负荷为主表左连接保证时间轴由负荷决定 df - merge(load_dt[, .(ts, load)], temp_wide, by ts, all.x TRUE) df[, date : as.Date(ts, tz UTC)] df - merge(df, holi_dt[, .(date as.Date(date), holiday)], by date, all.x TRUE) df[is.na(holiday), holiday : 0L] # 5) 补全缺失小时再做短缺口插值 full_idx - data.table(ts seq(min(df$ts), max(df$ts), by hour)) df - merge(full_idx, df, by ts, all.x TRUE) df[, load_fill : zoo::na.approx(load, x as.numeric(ts), maxgap 3, na.rm FALSE)]第 4 步必须用all.x TRUE以负荷为主表因为负荷才是时间轴的定义者如果反过来以温度为主表某些温度站晚开机的那几年会把负荷序列截短。第 5 步的maxgap 3是这套流程里最容易被忽略的参数只有连续缺失不超过 3 小时才做线性插值更长的缺口直接留 NA 让模型跳过。理由很直接一段二十小时的缺口被直线填平之后它下游所有lag24、lag168特征全都变成假数据而这些特征恰恰是模型最依赖的。na.approx的x参数要传数值时间戳而不是 POSIXct传日期对象在某些版本上会报错这是踩过一次就记住的小坑。2.3 特征构造滞后、滑动窗口与温度响应df[, :( hour hour(ts), dow wday(ts, week_start 1), # 1 周一 month month(ts), is_weekend as.integer(wday(ts, week_start 1) 6) )] # 滞后特征预测日零点能看到的最近值是 t-24所以 24/48/168 是安全边界 for (k in c(24, 48, 168)) { df[, paste0(lag, k) : shift(load_fill, k, type lag)] } # 七天滑动均值先整体后移 24 小时再滚动杜绝用到预测日的数据 df[, ma_7d : frollmean(shift(load_fill, 24), n 24 * 7, align right, na.rm TRUE)] # 多站温度先取均值再做 24 小时居中平滑近似体感温度 temp_cols - grep(^t[0-9]$, names(df), value TRUE) df[, temp_mean : rowMeans(.SD, na.rm TRUE), .SDcols temp_cols] df[, temp_smooth : frollmean(temp_mean, n 24, align center, fill NA)] # 分段温度响应低于 14 度开始采暖高于 22 度开始制冷 df[, :(hdd pmax(14 - temp_smooth, 0), cdd pmax(temp_smooth - 22, 0))]三个参数值得单独说。滞后阶选 24、48、168对应「昨天同一小时」「前天同一小时」「上周同一小时」这三条基本覆盖了日周期和周周期再往前加边际收益很低。ma_7d里那个shift(load_fill, 24)是防泄漏的关键动作直接对原始序列做frollmean会让预测日的负荷混进特征。temp_smooth用align center是想让温度曲线更平滑但居中窗口在序列末尾会产生 NA——如果预测日当天的温度预报拿不到完整窗口就把align改成right代价是响应略微滞后。温度响应做成折线hdd/cdd而不是直接丢原始温度是因为负荷对温度的关系是明显的 V 形太冷要供暖太热要制冷中间一段几乎平坦。直接把线性温度项塞进分位数回归模型会拟合出一条斜线两头都欠拟合。14 和 22 这两个拐点不是拍脑袋拿历史数据按温度分箱算平均负荷拐点位置一眼就能看出来不同地区要各自标定。2.4 训练集切分按预测日滚动不能随机抽概率预测的验证有个铁律不能随机抽 20% 的小时当测试集。相邻小时负荷的相关系数极高随机切分等于把答案拆成碎片撒进训练集回测分数会虚高到你自己都不信。正确做法是按预测日滚动用 D 日零点前全部历史训练预测 D 日 24 小时窗口往后推一天重训再预测。make_folds - function(dt, from, to) { days - seq(as.Date(from), as.Date(to), by day) lapply(days, function(d) { cut - as.POSIXct(d, tz UTC) list( train dt[ts cut], test dt[ts cut ts cut 24 * 3600] ) }) } folds - make_folds(df, 2010-01-01, 2010-12-31)from/to决定回测跨度一年大约 365 次重训用rq跑 24 个小时模型每个都要几秒整体下来是小时级。实际工程里常见的折中是每周重训一次、周内复用模型精度损失通常在小数点后第三位时间成本能砍掉八成。这个取舍在仓库里是可配置的别一上来就按最贵的配置跑。3. 99 个分位点怎么来的pinball 损失下的分位数回归与基线模型3.1 先写评分函数再谈模型评分函数不清楚模型就是玄学。GEFCOM14 用的是 pinball 损失也叫分位数损失对单个分位点 τL(y, q) τ·(y − q)当 y ≥ q否则 (1 − τ)·(q − y)写成 R 向量化形式特别短关键在把「τ 减指示函数」这个等价变形用上tau_grid - seq(0.01, 0.99, by 0.01) pinball_loss - function(y, Q, tau tau_grid) { # y长度 n 的实测值Qn × length(tau) 的预测矩阵第 j 列对应 tau[j] Q - as.matrix(Q) stopifnot(nrow(Q) length(y), ncol(Q) length(tau)) T_mat - matrix(tau, nrow nrow(Q), ncol length(tau)) # 按列填充每列一个 tau L - (y - Q) * (T_mat - (y Q)) mean(L) }matrix(tau, nrow(Q), length(tau))是按列填充的所以第 j 列整列都是tau[j]正好和Q的列对齐(y Q)里y是长度 n 的向量会按列自动循环成逻辑矩阵。最终mean(L)对所有小时、所有 τ 一视同仁地平均这就是赛题的官方口径。理解了这个函数模型形态就不用猜了——τ 取 0.9 时低估的惩罚是高估的 9 倍最优解自然就是条件 0.9 分位数而不是条件均值。这就是为什么不能拿一个点预测重复 99 遍那样每个分位点的损失都偏高平均下来分数很难看。3.2 按小时拆开的 24 组分位数回归一天 24 小时的负荷形态差异太大凌晨的分布窄、傍晚的分布宽塞进一个模型里互相打架。仓库的做法是按小时拆成 24 个模型每个模型内部再用 99 个 τ 一次性拟合。library(quantreg) fml - load_fill ~ lag24 lag48 lag168 ma_7d temp_smooth hdd cdd factor(dow) factor(month) holiday fit_one_hour - function(d, taus) { d - d[!is.na(load_fill)] # rq 支持一次拟合多个 tau返回对象可直接 predict 出矩阵 rq(fml, data d, tau taus, method fn) } models - lapply(1:24, function(h) { fit_one_hour(df[hour h], tau_grid) }) predict_one_hour - function(m, newd) { as.matrix(predict(m, newdata newd)) # n × 99 }几个参数要交代清楚。method fn是 Frisch-Newton 内点法样本量上万、自变量几十个时比默认单纯形法快一个量级。tau taus传向量是一次拟合全部 99 个分位点内部共享同一套设计矩阵比循环 99 次快得多。factor(dow)和factor(month)必须显式转因子否则会被当连续变量模型会学出「周五是 5比周二大 3」这种不存在的线性关系这个坑在后面的避坑章节还会展开。按小时拆分之后训练数据只剩原来的 1/24四年半也就 1600 来行。这个量级对rq刚好再少就该考虑把小时分组成早高峰、平峰、晚高峰几段用段内共享参数换样本量。3.3 两条必须跑的基线持续法与相似日任何模型出来之前先跑两条基线否则你不知道自己的提升是真的还是噪声。# 基线一持续法拿昨天同一小时当所有分位点的预测 naive_persist - function(newd) { matrix(newd$lag24, nrow nrow(newd), ncol length(tau_grid)) } # 基线二相似日按星期类型温度距离找 k 个历史日取同一小时的经验分位数 similar_day - function(hist, target, k 20) { d - hist[, .(date, dow, temp_smooth)] d - unique(d) setorder(d, temp_smooth) t0 - target$temp_smooth[1] d[, dist : abs(temp_smooth - t0) 5 * (dow ! target$dow[1])] pick - head(d[order(dist)], k) h - target$hour[1] vals - hist[date %in% pick$date hour h, load_fill] quantile(vals, probs tau_grid, na.rm TRUE, type 7) }持续法看着简陋但在负荷预测里极其难缠——它的 pinball 分数通常只比正经模型差 10% 到 20%任何新模型如果连这条线都甩不开说明特征工程或者数据对齐出了系统性错误。相似日那条里的5 * (dow ! ...)是给星期类型不匹配加的惩罚权重量纲上是「5 摄氏度」意思是宁可温度差 5 度也要找同为工作日的样本。这个权重按地区调商业负荷重的区域可以加到 8居民负荷为主的区域降到 3。quantile(..., type 7)是 R 的默认分位数算法样本量小的时候不同 type 出来的值能差几个点回测对比时最好三个基线统一用同一个 type。3.4 加一路梯度提升为后面的集成备料光靠线性分位数回归捕捉不了温度和负荷之间的交互——比如同样是 30 度工作日和周末的空调负荷曲线形状完全不同。补一路树模型是有必要的。library(xgboost) x_train - as.matrix(df[!is.na(load_fill), .(lag24, lag48, lag168, ma_7d, temp_smooth, hdd, cdd, hour, dow, month, holiday)]) y_train - df[!is.na(load_fill), load_fill] fit_gbm_quantile - function(X, y, tau) { dtrain - xgb.DMatrix(X, label y) xgb.train( params list(objective reg:quantileerror, quantile_alpha tau, eta 0.05, max_depth 6, subsample 0.8), data dtrain, nrounds 400, verbose 0 ) }reg:quantileerror加上quantile_alpha是指定分位点的官方写法eta 0.05配nrounds 400是我在几千行样本上比较稳的一组值样本更大就把轮数往上抬、eta 往下压。这路模型不需要多准它的价值在于和线性模型犯的错不一样——集成的前提是成员之间要有差异性两路高度相关的模型合起来毫无意义。4. QRA 集成与滚动回测让 24 小时预测向量不再各说各话4.1 QRA 在合成什么把三路预测当自变量再回归一次QRAQuantile Regression Averaging的思路朴素得有点反直觉把每个成员模型在同一个 τ 上的预测当作自变量真实负荷当作因变量再做一次 τ 分位数回归。得到的系数就是各成员的权重而且权重随 τ 变化——高峰时段的极端分位点可能更信树模型凌晨的窄分布更信线性模型。build_qra_matrix - function(pred_list) { # pred_list每个元素是 n × 99 的预测矩阵 do.call(cbind, pred_list) # 按列拼接成 n × (m*99) } qra_fit - function(X_list, y, taus) { # X 的每一列是一个模型在某一个 tau 上的预测逐 tau 单独回归 res - lapply(taus, function(t) { idx - seq_along(taus) (seq_along(X_list) - 1) * length(taus) # 上面这行只是示意列位实际按下面的方式取块更直观 NULL }) res }上面这段是骨架实际实现时按模型分块更清楚对每一个 τ从各成员矩阵里取同一列拼成一个 n × m 的矩阵命名成m1、m2、m3然后qra_one_tau - function(Xm, y, tau) { d - data.table(y y, Xm) rq(y ~ m1 m2 m3, data d, tau tau, method fn) }有三个细节决定了 QRA 是真提升还是假提升。第一训练集里的Xm必须是样本外预测也就是滚动回测产出的那一批预测值不能拿训练集拟合值凑——拟合值太准集成权重会严重偏向过拟合最厉害的那个模型上线就崩。第二系数允许为负这是 QRA 的正常行为不要强行lower 0负权重在成员高度相关时能起到去偏作用。第三y ~ m1 m2 m3里保留截距项允许集成有一个整体偏移去掉截距在成员有系统性偏差时会锁死误差。4.2 滚动原点回测把 15 轮任务还原成日常流程赛题只给了有限几轮任务拿来做模型选型样本太少。仓库里那套滚动回测才是日常用得上的东西每天推一次原点重训、预测、打分攒出几百天的样本外预测再在这个基础上做集成和参数对比。roll_backtest - function(dt, from, to, fit_fun, horizon 24) { origins - seq(as.POSIXct(from, tz UTC), as.POSIXct(to, tz UTC), by day) rbindlist(lapply(origins, function(o) { tr - dt[ts o] te - dt[ts o ts o horizon * 3600] if (nrow(te) horizon || nrow(tr) 24 * 60) return(NULL) P - fit_fun(tr, te) # horizon × 99 矩阵 data.table( ts rep(te$ts, times length(tau_grid)), y rep(te$load_fill, times length(tau_grid)), tau rep(tau_grid, each horizon), q as.vector(P) # 列优先与上面的 tau 展开顺序一致 ) }), fill TRUE) } res - roll_backtest(df, 2010-01-01, 2010-12-31, fit_fun my_pipeline) res[, mean((y - q) * (tau - (y q)))]as.vector(P)对矩阵是列优先展开所以 τ 的展开必须用rep(tau_grid, each horizon)两者顺序一旦错位分数会莫名其妙地变差而且看不出错在哪——这是这套代码里最容易翻车的一处。nrow(tr) 24 * 60那个守卫是给回测起点太靠前的情况兜底训练样本不足两个月的结果没有参考价值直接跳过比留着污染统计强。fit_fun这个参数是刻意留的接口传线性分位数回归就是一路模型传一个内部先跑三路再 QRA 合成的函数就是完整流程回测框架本身不用改。这种把「模型」和「评估」解耦的写法在这个仓库里贯彻得比较彻底也是它比一次性脚本耐改的原因。4.3 用 pinball 分解定位短板总分只有一个数看不出问题在哪。按小时、按月份拆开看才有信息量by_hour - res[, .(score mean((y - q) * (tau - (y q)))), by .(hour hour(ts))] setorder(by_hour, -score) by_tail - res[tau 0.05 | tau 0.95, .(score mean((y - q) * (tau - (y q)))), by .(tau)]by_hour里分数最高的那几个小时往往集中在傍晚负荷爬坡段和凌晨负荷低谷段。爬坡段的误差来自模型对「今天比昨天热多少」不敏感低谷段则是相对误差被放大——分母小同样的绝对误差在 pinball 上显得更刺眼。针对前者加温度差分特征今日均温减昨日均温针对后者考虑对负荷做对数变换再建模这是两条我验证过有效的改法。by_tail看的是尾部。如果 0.95 以上的分位点分数明显偏高说明模型对极端高温日的尖峰捕捉不足这时候该动的是训练样本加权给历史尖峰日更高权重而不是继续调树的深度。5. 避坑实录五个让 pinball 分数悄悄变差的细节5.1 时区解析错一小时滞后特征整体错位现象回测分数比预期差 15% 左右但模型残差图看起来很正常找不到明显异常点。原因时间戳按本地时区解析一年里有一次夏令时切换那天多出或少了一小时shift(load_fill, 24)之后所有滞后特征整体错位一小时。错位只在切换点附近出现摊到全年看不出来。解决所有as.POSIXct强制tz UTC并且在拼接完成后加一条断言检查每天恰好 24 行chk - df[, .N, by as.Date(ts, tz UTC)] stopifnot(all(chk$N 24))这条断言我现在固定在每个项目的第一段代码里五秒钟就能跑完能省掉半天排查。5.2 滑动窗口没后移信息泄漏让回测虚高现象回测 pinball 分数好得离谱比公开榜单上的水平还低一大截换一段数据立刻崩。原因frollmean(load_fill, n 168)直接对原始序列滚动预测日当天的负荷被算进了均值特征。模型在训练时看到了答案回测自然漂亮。解决先shift(load_fill, 24)再滚动并且写一个反向验证——把测试日的负荷列全部置为 NA重新生成特征对比两次结果是否一致。不一致就说明还有泄漏点没堵上。5.3 分位数交叉0.7 分位低于 0.6 分位现象分数比单模型还差检查发现某些小时里 τ0.7 的预测值小于 τ0.6。原因99 个 τ 是各自独立求解的线性规划的解不保证跨 τ 单调。交叉一旦发生pinball 损失会异常增大因为同一组预测里必然有分位点的方向是反的。解决对每一行预测做升序重排。这个操作叫重排rearrangement理论上不会增加 pinball 损失实现只有一行rearrange - function(Q) t(apply(as.matrix(Q), 1, sort))5.4 假期编码当连续变量丢进模型现象模型在春节、国庆前后预测偏差突然放大平时却很好。原因holiday列是整数编码直接进rq会被当连续变量模型学出「假期 3 比假期 1 更闲」这种凭空捏造的关系。解决factor(holiday)。更稳妥的做法是拆成两个特征is_holiday二值标记距离最近假期的天数前后各取 ±7 天后者能捕捉节前抢工和节后恢复的非对称效应。5.5 用拟合值训练 QRA权重全给了过拟合模型现象单模型回测 0.42集成之后 0.45越合越差。原因QRA 的训练矩阵用的是各成员在训练集上的拟合值。树模型在训练集上几乎完美权重被它吃光而这个权重在样本外完全不成立。解决QRA 的输入必须来自滚动回测产生的样本外预测哪怕这样只能拿到几百天、训练样本变少。宁可集成模型简单一点也不能用拟合值。这个坑我见过太多次本质上和 stacking 的做法是同一个道理。6. 提交前的自检习惯分位交叉修正与按小时拆解误差跑完一整轮回测别急着改模型先把下面这张表过一遍。这几项检查加起来不到两分钟但拦下的问题比我调参调一天发现的还多。检查项看什么通过标准时间轴完整性每天是否恰好 24 行全部为 24单调性每行 99 个分位是否升序交叉行数为 0与基线差距对比持续法的 pinball 分数优于基线 10% 以上分小时拆解哪几个小时分数最高无某小时异常突出尾部拆解τ ≥ 0.95 与 τ ≤ 0.05 的分数与整体分数同量级单调性那一项如果不通过就在输出前统一走一遍rearrange()这一步几乎是免费的午餐理论上有保证、实现上三行代码没有任何理由不做。分小时拆解那项如果发现某一个小时独占鳌头先怀疑特征里混进了时间相关的泄漏别急着给那个小时单独建模型。还有一个容易忽略的自检把预测矩阵按 τ 求积分得到隐含的预测均值和持续法的预测值画在同一张图上。如果两条线在某些日子分叉得离谱就去翻那几天的原始数据通常是传感器故障或者抄表异常留下的脏数据。这个方法不需要真值纯靠内部一致性就能提前发现大量问题。从那以后我每次跑完回测都强制先看按小时的 pinball 分解表确认没有单点异常之后才允许自己动模型参数——顺序反了很容易把数据问题当成模型问题越调越偏。这套流程现在已经是我的固定起手式。希望帮到你。本文还有配套的精品资源点击获取

从R脚本复现GEFCOM14概率负荷预测：分位数回归与QRA集成

相关推荐

基于傅里叶描述子与SVM的Python手势识别源码实践

HyperDbg 单机内核调试：VMM 加载、EPT Hook 与事件排错

傅里叶描述子与SVM驱动的低算力实时手势识别实践

更多资讯

完全二叉树、非完全二叉树与哈希表的实现与应用

Python中文词云7案例：从能出图到能控图

Conda环境清理全攻略：释放磁盘空间与提升性能

阶梯碳交易与电制氢综合能源系统24小时调度建模

冒泡排序与稀疏数组：算法优化与高效存储技术

拆开AI黑箱：从梯度下降到SHAP的模型可解释性实践

今日更新

PLFM_RADAR：从被动告警到主动感知的微服务健康监控系统

SPI隔离设计实战：MAX14431BAEE+与STM32F407工业级应用

数字隔离器如何根治SPI通信地环路干扰：MAX14431与TM4C129实战

本周精选

面向复杂光照与角度的宠物面部识别鲁棒性增强策略

基于DRV8818与STM32的双极步进电机驱动方案实战解析

DRV8818PWPR+PIC18F45K80双极步进电机控制方案实战

本月热门

面向复杂光照与角度的宠物面部识别鲁棒性增强策略

基于DRV8818与STM32的双极步进电机驱动方案实战解析

DRV8818PWPR+PIC18F45K80双极步进电机控制方案实战

读完文章，想聊聊您的网站？