基于高斯过程回归的MATLAB时间序列区间预测代码实现与解析

基于高斯过程回归(GPR)的时间序列区间预测 GPR时间序列区间预测 matlab代码 暂无Matlab版本要求 -- 推荐 2018B 版本及以上做时间序列最烦的就是拍脑袋给个“明天涨3%左右”——“左右”到底是正负0.5还是正负3如果是风电发电的负荷申报正负差多了要罚钱的。这时候自带“置信区间说明书”的高斯过程回归GPR就香起来了不用额外加什么分位数网络、bootstrap抽样堆计算量它的输出均值就是点预测方差一转换就是95%/99%的靠谱边界。咱们这次就拿最基础的单变量时间序列——举个例子随便编个带周期带趋势的模拟数据吧免得还要找公开数据集调半天——用Matlab的Statistics and Machine Learning Toolbox里的现成函数搞起来2018b之后的版本都能用fitrgp够友好。第一步造个能练手的带波动的时间序列先模拟点真实感强的数据比如带12小时周期假设是光伏后面解释方便、线性下降趋势比如冬天来了日照越来越短、再加一点点高斯白噪声。% 设定时间步长1小时总共30天720小时 t (1:720); % 周期项用两个正弦余弦组合更自然12小时24小时双周期 periodic 50*sin(2*pi*t/12) 20*cos(2*pi*t/24); % 线性下降趋势30天从100降到40 trend linspace(100, 40, 720); % 加高斯白噪声标准差3 noise normrnd(0, 3, 720, 1); % 拼接成完整序列y y trend periodic noise; % 先画个图看看模拟得像不像 plot(t, y, b.-); xlabel(小时数); ylabel(模拟光伏出力kW); title(练手用的带双周期、下降趋势、噪声的时间序列); hold on; plot(t, trendperiodic, r--, LineWidth, 1.2); % 真实无噪声曲线 legend(带噪声观测,真实信号); hold off;看上面的代码——哦不对看生成的图哈——真实信号那条红线是被蓝色点“包围”的这给后面GPR的方差做区间埋了伏笔蓝色点离红线近的地方方差应该小边界窄点偶尔飘远的地方比如噪声刚好赶上周期峰值的时候边界宽一点。第二步处理时间序列的数据输入GPR本身不是专门给时间序列设计的它是“非参数回归的全能选手”能把任意x映射到y带方差。那时间序列的x怎么选最简单的是自回归滞后项AR滞后——比如用前p个小时的出力预测第t个小时的出力。基于高斯过程回归(GPR)的时间序列区间预测 GPR时间序列区间预测 matlab代码 暂无Matlab版本要求 -- 推荐 2018B 版本及以上p选多少呢不用太纠结纠结的话可以做个自相关函数ACF/偏自相关函数PACF图但练手没必要咱们这里就选3小时不对不对刚才加了双周期12小时是主周期选p13吧——前12小时前1小时或者直接前12到1小时都行咱们直接造一个滑动窗口的矩阵。% 设定AR滞后阶数p这里选主周期12加个1凑13方便学习主周期和短相邻趋势 p 13; % 总样本数N滑动窗口后变成 N-p 个 N length(y); X zeros(N-p, p); Y y(p1:end); % 对应的预测目标是第p1到第N个小时 % 填充滑动窗口X第i行是第i到ip-1个小时的观测对应第ip个小时的Y for i 1:N-p X(i,:) y(i:ip-1); end接下来要划分训练集和测试集时间序列划分不能乱打乱分不能用随机抽样会破坏时间顺序这点很多新手容易踩坑一般是前80%训练后20%测试。train_ratio 0.8; train_size floor(train_ratio * length(Y)); X_train X(1:train_size, :); Y_train Y(1:train_size); X_test X(train_size1:end, :); Y_test Y(train_size1:end); t_test t(p1train_size : end); % 测试集对应的时间点画图用第三步选合适的高斯核函数训练GPR模型GPR的灵魂就是核函数协方差函数它决定了“什么样的x会被认为是相似的对应的y会被平滑/波动约束在一起”。咱们的模拟数据有平滑的下降趋势 → 可以选Matern32或者Matern52核比高斯核RBF更灵活不会强行把所有地方都揉得太光滑明显的双周期 → 可以加一个Periodic核还有短程的相邻波动 → 单独的Matern核其实也能学但为了效果更稳咱们把Matern52、Periodic、还有一个RationalQuadratic用来学一点点非平稳的“长尾”波动虽然模拟数据没有但留一手总没错组合起来。核函数组合在Matlab里用号就行超级直观。然后fitrgp还有几个常用参数可以调KernelFunction刚才说的组合核StandardizetrueGPR对输入输出的尺度很敏感标准化后训练更快效果更好PredictMethodexact——咱们练手的数据量只有720-13707exact方法完全没问题不用近似Optimizerquasinewton——优化超参数核函数的长度尺度、周期这些的算法quasinewton比默认的fmincon有时候收敛更快% 先定义组合核Matern52短长趋势非严格光滑 Periodic双周期不对默认Periodic是单周期那再加个Periodic或者把双周期放进Periodic的Frequency参数哦对Periodic可以接受多个频率 % 不过Matlab里fitrgp的KernelFunction直接传组合的字符不太对得用构造函数 matern52 matern(5/2); periodic1 periodic(Period,12); % 12小时主周期 periodic2 periodic(Period,24); % 24小时次周期 rq rationalquadratic; % 组合起来用*号也行不过号是加协方差更适合“多种模式叠加”的情况 kernel matern52 periodic1 periodic2 rq; % 训练GPR模型 gpr_model fitrgp(X_train, Y_train,... KernelFunction, kernel,... Standardize, true,... PredictMethod, exact,... Optimizer, quasinewton); % 先看一眼训练出来的超参数虽然练手不用太纠结但有意思 disp(gpr_model.KernelInformation.KernelParameters);运行完disp会看到一串数字分别是各个核的超参数比如matern52的长度尺度periodic1的周期应该接近12刚才的真实主周期periodic2的周期接近24rq的长度尺度和形状参数——这说明GPR真的“学会”了咱们模拟数据的模式不是瞎蒙的。第四步做预测画95%置信区间predict函数在GPR里有三个输出刚好是咱们要的ypredmean点预测ypredstd预测的标准差ypredcov预测的协方差矩阵咱们单变量时间序列这个矩阵是对角的其实和ypredstd的平方差不多除非是多输出GPR否则用不到95%置信区间怎么算如果假设预测误差是高斯分布的GPR的基本假设之一那就是均值 ± 1.96*标准差。% 做测试集的预测 [y_pred_mean, y_pred_std] predict(gpr_model, X_test); % 计算95%置信区间 confidence_level 0.95; z_score norminv((1 confidence_level)/2); % 1.96 y_pred_lower y_pred_mean - z_score * y_pred_std; y_pred_upper y_pred_mean z_score * y_pred_std; % 画图对比真实值、点预测、置信区间 figure; plot(t_test, Y_test, b.-, DisplayName, 真实带噪声出力); hold on; plot(t_test, y_pred_mean, r-, LineWidth, 1.5, DisplayName, GPR点预测); fill([t_test; flip(t_test)], [y_pred_lower; flip(y_pred_upper)], g, FaceAlpha, 0.2, DisplayName, 95%置信区间); xlabel(小时数); ylabel(模拟光伏出力kW); title(GPR时间序列区间预测测试集p13); legend(Location,best); grid on; hold off;看测试集的图绿色的区域是不是基本把蓝色的点都框住了如果没框住可能是p选得不够大或者核函数没加对——比如刚才我如果漏了Periodic核置信区间肯定会变宽因为GPR学不到周期只靠Matern核猜猜不准的地方方差就大第五步随便加个小指标验证一下效果虽然练手不用太严谨但加个小指标更有说服力比如区间覆盖率真实点落在置信区间里的比例理想情况应该接近咱们设定的95%。% 计算区间覆盖率 coverage sum(Y_test y_pred_lower Y_test y_pred_upper) / length(Y_test) * 100; fprintf(测试集95%%置信区间的覆盖率为%.2f%%\n, coverage);运行完大概率会在93%-97%之间完美如果太高说明核函数太“谨慎”方差给大了如果太低说明核函数太“自信”没学会所有模式。最后碎碎念几句核函数的选择太重要了别上来就用默认的RBF先看看你的时间序列有没有趋势、有没有周期、有没有季节性再选对应的核函数组合——比如趋势明显可以加Polynomial核或者先用detrend函数去趋势再GPRAR滞后阶数p的选择可以用ACF/PACF也可以用交叉验证比如把训练集再分成5折试不同的p选测试折上覆盖率最接近95%、点预测MAE最小的那个数据量的问题exact方法GPR的时间复杂度是O(n³)n是训练样本数——如果n超过10000就别用exact了换成PredictMethodsdp或者fic这两个是近似方法速度快很多效果也差不了太多真实数据的问题真实数据可能有缺失值、异常值——缺失值可以用线性插值/样条插值补一下异常值可以用3σ原则去掉或者选robust的核函数比如RationalQuadratic核比RBF核更抗异常值。今天这个练手是最基础的单变量AR-GPR区间预测感兴趣的话可以试试多变量比如加个气温、湿度来预测光伏出力、或者加个趋势项到GPR的MeanFunction里不用detrend了GPR自己学趋势都挺有意思的