MATLAB实战WOA-SVM多变量时间序列预测鲸鱼优化SVR超参数、时序验证、误差诊断与完整代码多输入单输出 · Gaussian SVR · Whale Optimization Algorithm · 严格时序切分 · 可复现工程流程摘要多变量时间序列预测的难点并不只是“选一个回归器”而是如何在不泄露未来信息的前提下将多源变量的滞后关系转化为稳定的监督学习样本并为非线性模型找到合适的超参数。本文构建一套完整的 WOA-SVM严格来说是 WOA-SVR预测流程首先使用滑动窗口把多变量序列展开为固定维度特征随后按时间顺序划分训练集、验证集和测试集再由鲸鱼优化算法WOA在对数空间搜索支持向量回归的 BoxConstraint、KernelScale 与 Epsilon最后对预测曲线、RMSE、MAE、R² 及残差分布进行联合评估。文末给出可直接运行的 MATLAB 完整脚本与本地函数并补充参数边界、数据泄露、RBF 参数映射、MAPE 失效场景以及真实业务数据替换方法。关键词MATLAB鲸鱼优化算法WOA支持向量机SVR多变量时间序列超参数优化滑动窗口运行环境建议 MATLAB R2021a 及以上版本需要 Statistics and Machine Learning Toolbox。核心建模函数使用 fitrsvm。若 MATLAB 版本支持 tspartition也可将文中的手工前向验证替换为 expanding-window 或 sliding-window 交叉验证。1. 问题定义为什么“WOA SVM”适合中小样本多变量时序在电力负荷、设备状态、交通流量、环境监测、价格需求等场景中目标变量通常同时受多个外部因素影响而且这些因素往往存在滞后、耦合、周期与噪声。直接对原始序列做普通回归容易遗漏时间依赖直接使用深度网络则可能在样本量有限时带来更高的训练成本和调参复杂度。支持向量回归SVR对小样本、高维和非线性关系具有较好的适应性但它对惩罚系数 C、核尺度 σMATLAB 中为 KernelScale以及 ε 不敏感区间非常敏感。参数过强会把噪声也拟合进去参数过弱又会产生明显欠拟合。WOA 的作用不是“替代模型”而是把这组超参数选择转化为连续优化问题通过群体搜索降低人工试参的盲目性。SVM/SVR 负责建立“历史多变量特征 → 未来目标值”的非线性映射WOA 负责在连续超参数空间中搜索更合适的模型配置严格的时间顺序切分负责阻断未来信息泄露标准化、窗口构造和残差诊断负责保证工程链路可解释、可复现。图 1 WOA-SVR 多变量时间序列预测完整工程链路2. 模型原理SVR、RBF 核与 WOA 如何连接起来2.1 ε-SVR 的核心目标在回归任务中SVM 通常以 ε-SVR 的形式使用。它希望找到一个尽量“平坦”的函数同时允许预测值落在真实值上下 ε 范围内而不产生惩罚。经典原始问题可以写为min 1/2 ||w||² C Σ(ξᵢ ξᵢ*)s.t. yᵢ - (wᵀφ(xᵢ)b) ≤ ε ξᵢ(wᵀφ(xᵢ)b) - yᵢ ≤ ε ξᵢ*ξᵢ, ξᵢ* ≥ 0其中 C 决定“模型平滑性”和“训练误差惩罚”之间的权衡ε 决定多大的误差可以被视为可接受范围。C 太大时模型可能过度追随训练样本C 太小时则可能欠拟合ε 太大时模型会忽略过多细节太小时支持向量数量增加、模型更容易受噪声影响。2.2 MATLAB 的 Gaussian KernelScale 不要与 γ 混用RBFGaussian核可写成K(xᵢ, xⱼ) exp( -||xᵢ - xⱼ||² / (2σ²) )很多资料习惯用 γ 表示 exp(-γ||xᵢ-xⱼ||²)此时 γ 1/(2σ²)。但 MATLAB fitrsvm 直接暴露的是 KernelScale也就是 σ。工程实现中最稳妥的做法是让 WOA 直接搜索 KernelScale而不是先搜索 γ 再做二次换算。这样既避免参数定义歧义也让代码与 MATLAB 官方接口完全一致。2.3 WOA 的三类位置更新机制鲸鱼优化算法由 Mirjalili 与 Lewis 于 2016 年提出模拟座头鲸气泡网捕食过程。每一条“鲸鱼”代表一组候选超参数适应度函数则对应这组参数在验证集上的预测误差。经典 WOA 主要包含以下三种更新行为包围猎物当 |A| 1 时个体向当前最优解靠拢增强局部开发能力随机搜索当 |A| ≥ 1 时个体参考随机鲸鱼更新位置扩大搜索范围螺旋气泡网以螺旋轨迹围绕当前最优解移动在局部区域继续精细搜索。A 2ar - a, C 2rD |C·X* - X|X(t1) X* - A·D螺旋更新X(t1) D · exp(bl) · cos(2πl) X*这里 a 会从 2 逐步衰减到 0使搜索过程从更强调全局探索逐渐过渡到局部开发。需要强调的是WOA 是随机元启发式算法它能改善搜索覆盖并降低人工调参成本但并不提供“必然找到全局最优”的数学保证。2.4 为什么在对数空间搜索更合理C、KernelScale、Epsilon 往往跨越多个数量级。如果在 [0.01, 1000] 这样的线性区间直接均匀初始化大量搜索点会堆积在大数值区域小数值区间被严重稀释。因此本文让 WOA 搜索 log10 参数z₁ log10(C) ∈ [-2, 3]z₂ log10(KernelScale) ∈ [-2, 2]z₃ log10(Epsilon) ∈ [-3, -0.7]真实参数C 10^z₁, KernelScale 10^z₂, Epsilon 10^z₃这样每一个数量级获得近似均衡的搜索机会通常比在原始尺度上“让鲸鱼在超大海域随机游动”更有效。3. 数据设计把多变量时间序列转成监督学习样本为了让示例可以独立运行本文构造 4 个输入变量与 1 个目标变量。x1 表示主要周期驱动x2 表示较慢周期因子x3 表示较高频的辅助波动x4 表示缓慢趋势与随机扰动目标 y 同时包含线性项、二次非线性项、差分项与噪声。图 2 多变量合成时间序列及目标变量示例3.1 滑动窗口用过去 6 个时刻预测当前时刻假设共有 4 个输入变量历史窗口 lag6则每个监督学习样本包含 6×424 个特征。对第 i 个目标时刻输入为 i-1 到 i-6 的四个变量输出为当前时刻 y(i)。这种做法把“序列预测”转换成普通回归问题同时显式保留局部时间依赖。[x1(i-1), x2(i-1), x3(i-1), x4(i-1),x1(i-2), x2(i-2), x3(i-2), x4(i-2),...x1(i-6), x2(i-6), x3(i-6), x4(i-6)] → y(i)3.2 时间序列不能随机打乱普通机器学习常使用随机 K 折交叉验证但时间序列存在严格因果顺序。如果把未来样本随机混入训练集就可能让模型间接看到未来分布得到过于乐观的误差。本文采用两层时间切分外层前 70% 样本作为训练段后 30% 作为最终测试段内层在外层训练段中再取前 80% 作为 WOA 训练子集后 20% 作为验证子集WOA 的适应度只看内层验证误差最终测试集在参数搜索期间完全不可见。更严格的方案当样本量足够时建议用 expanding-window 或 sliding-window 的多折前向验证替代单一验证段。其代价是每一个候选参数都要训练多次 SVR计算量会显著增加但评价方差更小。3.3 标准化必须只使用训练段统计量标准化本身也会泄露信息。正确流程是先切分再用训练段计算均值和标准差然后用同一组统计量变换验证集与测试集。绝不能先对全量数据求均值/标准差再切分训练测试。本文同时标准化 X 和 y使不同尺度的特征更适合 RBF 核并让 ε 的搜索范围更稳定。4. WOA-SVR 优化流程与适应度设计每个鲸鱼位置包含 3 个维度[log10(C), log10(KernelScale), log10(Epsilon)]。对于一个候选位置先还原真实参数再训练 SVR最后以验证集标准化 RMSE 作为适应度。适应度越小参数组合越好。参数设置说明种群规模 pop24兼顾搜索多样性与计算量最大迭代 maxIter45足够观察收敛过程log10(C)[-2, 3]C ∈ [0.01, 1000]log10(KernelScale)[-2, 2]σ ∈ [0.01, 100]log10(Epsilon)[-3, -0.7]ε 约 ∈ [0.001, 0.20]适应度验证集 RMSE不使用测试集图 3 WOA 搜索过程中验证集最优 RMSE 的变化如果收敛曲线前期快速下降、后期基本稳定说明搜索已经找到较稳定区域如果直到最后仍持续大幅下降通常需要增加迭代次数如果从第一代开始几乎不变则应检查搜索边界、适应度函数、数据标准化和参数尺度是否合理。5. 等价复现实验结果不只看一条预测曲线结果说明下表与图 3图 6 使用与文末 MATLAB 脚本一致的数据构造、滑动窗口、时间切分和超参数定义在等价实现中固定随机种子生成。不同随机数生成器、MATLAB 版本与优化轨迹会使最优参数和小数点后的指标略有差异因此运行本地脚本时应以实际输出为准。等价复现实验中WOA 搜索得到 C≈2.0409、KernelScale≈6.0921、Epsilon≈0.1880。模型RMSEMAER²SMAPE(%)持续性基线0.24610.19440.742362.04线性 SVR0.16780.12930.880246.87固定参数 RBF-SVR0.40610.33790.2987117.46WOA-SVR0.15970.12280.891543.79可以看到WOA-SVR 相比固定参数 RBF-SVR 的 RMSE 降低约 60.7%说明“RBF 核本身并不自动等于高精度”核尺度与惩罚系数的配置会显著改变结果。与此同时线性 SVR 在这一合成数据上也是一个很强的基线WOA-SVR 的 RMSE 进一步降低约 4.8%。这说明对比实验不能只放一个弱基线更不能把“智能优化”写成必然碾压所有模型。图 4 WOA-SVR 测试集预测曲线图 5 WOA-SVR 测试集残差分布图 6 持续性基线、线性 SVR、固定参数 RBF-SVR 与 WOA-SVR 的误差对比5.1 为什么这里不把 MAPE 作为主指标本示例目标 y 可能接近 0甚至出现负值。MAPE 的分母是 |y|当真实值接近 0 时少量绝对误差会被放大成极高百分比当 y0 时更是直接失去定义。因此本文以 RMSE、MAE、R² 为主并补充对零值更稳定的 SMAPE。对于电力负荷、销量等严格为正且远离 0 的目标MAPE 才更有解释价值。6. MATLAB 完整代码一份脚本直接跑通下面代码包含数据生成、滑动窗口、时间切分、WOA 三参数优化、最终模型训练、四类基线对比、误差指标与图形输出。将整段保存为一个 .m 文件即可运行本地函数放在脚本末尾。clear; clc; close all;rng(2025,twister);%% 1. 构造多变量时间序列n 600;t (1:n);x1 sin(2*pi*t/24) 0.15*randn(n,1);x2 cos(2*pi*t/48) 0.12*randn(n,1);x3 0.5*sin(2*pi*t/12) 0.08*randn(n,1);x4 0.03*t/n 0.10*randn(n,1);dx1 [0; diff(x1)];y 0.6*x1 0.25*x2.^2 - 0.3*x3 ...0.4*x4 0.05*dx1 0.10*randn(n,1);data [x1,x2,x3,x4,y];%% 2. 滑动窗口过去 lag 个时刻预测当前 ylag 6;numSamples n - lag;numInput 4;X zeros(numSamples, lag*numInput);Y zeros(numSamples, 1);targetTime (lag1:n);for k 1:numSamplesi k lag;feat zeros(1, lag*numInput);c 1;for j 1:lagfeat(c:cnumInput-1) data(i-j,1:numInput);c c numInput;endX(k,:) feat;Y(k) data(i,5);end%% 3. 外层时间切分前 70% 训练后 30% 测试trainRatio 0.70;numTrain floor(numSamples*trainRatio);XTrain X(1:numTrain,:);YTrain Y(1:numTrain);XTest X(numTrain1:end,:);YTest Y(numTrain1:end);TTest targetTime(numTrain1:end);%% 4. WOA 在对数空间优化 C、KernelScale、Epsilonpop 24;maxIter 45;% z [log10(C), log10(KernelScale), log10(Epsilon)]lb [-2, -2, -3];ub [ 3, 2, -0.7];dim 3;fitnessFcn (z) svrFitnessTimeHoldout(z, XTrain, YTrain, 0.80);[bestZ, bestFit, curve] WOA(fitnessFcn, pop, maxIter, lb, ub, dim);bestC 10.^bestZ(1);bestKernelScale 10.^bestZ(2);bestEpsilon 10.^bestZ(3);fprintf(Best validation RMSE (standardized) %.6f\n, bestFit);fprintf(Best BoxConstraint %.6f\n, bestC);fprintf(Best KernelScale %.6f\n, bestKernelScale);fprintf(Best Epsilon %.6f\n, bestEpsilon);%% 5. 只用完整训练集统计量做标准化[muX, sigX] trainZScoreStats(XTrain);XTrainN applyZScore(XTrain, muX, sigX);XTestN applyZScore(XTest, muX, sigX);muY mean(YTrain);sigY std(YTrain,0,1);if sigY 0, sigY 1; endYTrainN (YTrain - muY)./sigY;%% 6. 训练最优 WOA-SVRmdlWOA fitrsvm(XTrainN, YTrainN, ...KernelFunction,gaussian, ...BoxConstraint,bestC, ...KernelScale,bestKernelScale, ...Epsilon,bestEpsilon, ...Standardize,false);YPredWOAN predict(mdlWOA, XTestN);YPredWOA YPredWOAN*sigY muY;%% 7. 三个对照模型% 7.1 固定参数 RBF-SVRmdlFixed fitrsvm(XTrainN, YTrainN, ...KernelFunction,gaussian, ...BoxConstraint,1, ...KernelScale,1, ...Epsilon,0.05, ...Standardize,false);YPredFixed predict(mdlFixed,XTestN)*sigY muY;% 7.2 线性 SVRmdlLinear fitrsvm(XTrainN, YTrainN, ...KernelFunction,linear, ...BoxConstraint,1, ...Epsilon,0.05, ...Standardize,false);YPredLinear predict(mdlLinear,XTestN)*sigY muY;% 7.3 持续性基线用上一时刻真实 y 预测当前 yYPredPersistence y(TTest-1);%% 8. 计算指标M1 regressionMetrics(YTest, YPredPersistence);M2 regressionMetrics(YTest, YPredLinear);M3 regressionMetrics(YTest, YPredFixed);M4 regressionMetrics(YTest, YPredWOA);Result table(...[M1.RMSE;M2.RMSE;M3.RMSE;M4.RMSE], ...[M1.MAE; M2.MAE; M3.MAE; M4.MAE], ...[M1.R2; M2.R2; M3.R2; M4.R2], ...[M1.SMAPE;M2.SMAPE;M3.SMAPE;M4.SMAPE], ...VariableNames,{RMSE,MAE,R2,SMAPE}, ...RowNames,{Persistence,Linear_SVR,Fixed_RBF_SVR,WOA_SVR});disp(Result);%% 9. 可视化figure(Color,w,Name,WOA 收敛曲线);plot(curve,-o,LineWidth,1.2,MarkerSize,3);xlabel(迭代次数); ylabel(验证集标准化 RMSE);title(WOA 超参数搜索收敛曲线); grid on;figure(Color,w,Name,预测结果);plot(YTest,k,LineWidth,1.3); hold on;plot(YPredWOA,r--,LineWidth,1.3);legend(真实值,WOA-SVR 预测值,Location,best);xlabel(测试样本序号); ylabel(目标值);title(WOA-SVR 多变量时间序列预测); grid on;figure(Color,w,Name,残差分布);residual YTest - YPredWOA;histogram(residual,30);xline(0,--);xlabel(残差); ylabel(频数);title(WOA-SVR 测试集残差分布); grid on;figure(Color,w,Name,模型误差对比);bar([M1.RMSE M1.MAE; M2.RMSE M2.MAE; ...M3.RMSE M3.MAE; M4.RMSE M4.MAE]);set(gca,XTickLabel,{Persistence,Linear SVR,Fixed RBF,WOA-SVR});ylabel(误差); legend(RMSE,MAE,Location,best);title(不同模型误差对比); grid on;%% -------------------- 本地函数 --------------------function f svrFitnessTimeHoldout(z, X, Y, innerTrainRatio)% 只在外层训练集内部做时间顺序 holdoutN size(X,1);nTr floor(N*innerTrainRatio);Xtr X(1:nTr,:);Ytr Y(1:nTr);Xval X(nTr1:end,:);Yval Y(nTr1:end);[muX, sigX] trainZScoreStats(Xtr);XtrN applyZScore(Xtr, muX, sigX);XvalN applyZScore(Xval, muX, sigX);muY mean(Ytr);sigY std(Ytr,0,1);if sigY 0, sigY 1; endYtrN (Ytr-muY)./sigY;YvalN (Yval-muY)./sigY;C 10.^z(1);kernelScale 10.^z(2);epsilonSVR 10.^z(3);trymdl fitrsvm(XtrN,YtrN, ...KernelFunction,gaussian, ...BoxConstraint,C, ...KernelScale,kernelScale, ...Epsilon,epsilonSVR, ...Standardize,false);pred predict(mdl,XvalN);f sqrt(mean((pred-YvalN).^2));if ~isfinite(f), f 1e6; endcatchf 1e6;endendfunction [bestPos, bestFit, curve] WOA(fitnessFcn, pop, maxIter, lb, ub, dim)Xpos zeros(pop,dim);fit zeros(pop,1);for i 1:popXpos(i,:) lb rand(1,dim).*(ub-lb);fit(i) fitnessFcn(Xpos(i,:));end[bestFit,idx] min(fit);bestPos Xpos(idx,:);curve zeros(maxIter,1);for iter 1:maxItera 2 - iter*(2/maxIter);for i 1:popr1 rand();r2 rand();A 2*a*r1 - a;Ccoef 2*r2;p rand();b 1;l -1 2*rand();if p 0.5if abs(A) 1D abs(Ccoef*bestPos - Xpos(i,:));Xnew bestPos - A.*D;elserandIdx randi(pop);Xrand Xpos(randIdx,:);D abs(Ccoef*Xrand - Xpos(i,:));Xnew Xrand - A.*D;endelseD abs(bestPos - Xpos(i,:));Xnew D.*exp(b*l).*cos(2*pi*l) bestPos;endXnew max(Xnew,lb);Xnew min(Xnew,ub);newFit fitnessFcn(Xnew);if newFit fit(i)Xpos(i,:) Xnew;fit(i) newFit;endif fit(i) bestFitbestFit fit(i);bestPos Xpos(i,:);endendcurve(iter) bestFit;endendfunction [muX, sigX] trainZScoreStats(X)muX mean(X,1);sigX std(X,0,1);sigX(sigX0) 1;endfunction XN applyZScore(X, muX, sigX)XN (X-muX)./sigX;endfunction M regressionMetrics(yTrue, yPred)err yTrue-yPred;M.RMSE sqrt(mean(err.^2));M.MAE mean(abs(err));denom sum((yTrue-mean(yTrue)).^2);if denom 0M.R2 NaN;elseM.R2 1-sum(err.^2)/denom;endM.SMAPE mean(2*abs(err)./(abs(yTrue)abs(yPred)1e-8))*100;end7. 真实业务数据替换从合成数据迁移到 CSV / Excel真正落地时最重要的不是“把 data 换成 readmatrix 一行就结束”而是明确时间列、输入列、目标列、采样频率与缺失处理。推荐先把数据整理成第一列时间戳后续若干列为外生输入最后一列为预测目标。T readtable(your_timeseries.csv);% 示例假设第 1 列为时间第 2:5 列为输入第 6 列为目标Xraw T{:,2:5};yraw T{:,6};% 缺失值示例线性插值真实项目应结合业务含义选择方法Xraw fillmissing(Xraw,linear);yraw fillmissing(yraw,linear);data [Xraw,yraw];n size(data,1);如果多个传感器采样频率不同应先按时间戳重采样和对齐如果目标存在明显日周期、周周期可加入 hour、weekday、holiday 等时间特征如果某些输入变量只在未来已知例如天气预报要明确预测时刻能否真实获得这些变量如果存在概念漂移应采用滚动训练或定期重估超参数而不是永久固定一次 WOA 结果。7.1 多步预测怎么改本文代码是一阶 ahead 的单步预测。如果要预测未来 12 个时刻可以选择两种常见策略Direct为每个 horizon 分别训练一个 SVR例如 h1、2、…、12 各一个模型。优点是不传播前一步预测误差缺点是模型数量增加Recursive只训练一个一步模型然后把前一步预测值回填继续预测。模型简单但误差会逐步累积。8. 常见问题与避坑这些细节比“多迭代 20 次”更重要问题 1先归一化全数据再划分训练测试。这是最隐蔽的数据泄露之一。均值和标准差必须只从训练段计算再应用到验证和测试。问题 2为了“样本均匀”把时间序列随机打乱。这会破坏因果顺序。最终测试必须是时间上晚于训练段的数据。问题 3把 γ 直接当成 MATLAB KernelScale。两者定义不同。本文直接优化 KernelScale避免不必要的换算。问题 4WOA 的适应度直接用测试集 RMSE。这样测试集就参与了模型选择最终误差失去独立性。适应度必须使用训练段内部验证集。问题 5搜索范围越大越好。范围过大反而浪费搜索预算。先根据标准化后的数据给出合理数量级再在对数空间搜索。问题 6只报告一次随机种子的结果。WOA 是随机算法。严谨实验应使用多个随机种子重复运行报告均值、标准差与最优值。问题 7只看 R²不看误差曲线与残差。R² 高并不代表每个关键时段都预测稳定。应同时检查峰值段、残差偏移、异常点与趋势滞后。问题 8认为 WOA 能弥补差的特征工程。优化器只能在给定特征与模型空间内搜索。窗口长度、变量质量、时序对齐和异常处理仍然决定上限。9. WOA-SVR 的适用边界什么时候值得用什么时候不必用方法更适合的场景主要优势主要代价WOA-SVR中小样本、非线性、超参数空间连续实现直观、可解释、无需梯度随机性、计算成本随 pop×iter 增长网格搜索 SVR参数范围很窄、维度少结果稳定、容易复现网格细时计算爆炸连续空间利用率低Bayesian Optimization单次训练较贵、超参数不多样本效率高、工具链成熟实现逻辑比 WOA 更复杂LSTM/GRU/Transformer大样本、长依赖、多步序列表达能力强数据量、训练成本和调参复杂度更高如果数据只有几百到几千个样本、特征维度中等、希望快速得到一个稳健的非线性回归基线同时又不希望手工尝试大量 C/KernelScale/Epsilon 组合那么 WOA-SVR 很有价值。反过来如果样本已经达到几十万甚至更大规模标准核 SVM 的训练复杂度可能成为瓶颈此时应考虑近似核方法、树模型或更适合大数据的深度学习方案。10. 进一步增强从“能跑”到“更稳健”的四个方向1.把窗口长度 lag 也加入优化变量但必须把整数变量做取整处理并注意搜索维度增加会显著放大计算量2.用 expanding-window 多折验证代替单一内层 holdout使适应度更能代表跨时间段稳定性3.对特征加入滞后差分、移动均值、周期编码或领域变量并用消融实验验证每类特征的真实贡献4.对 WOA 做多随机种子重复实验并把“平均 RMSE ± 标准差、最优 RMSE、运行时间”一起报告避免只展示一次幸运搜索。11. 总结WOA-SVM 多变量时间序列预测的关键不在于把两个算法名字拼在一起而在于把整个建模链路做对先用滑动窗口显式表达历史依赖再用严格时间切分隔离未来数据标准化只拟合训练段WOA 在对数空间优化 fitrsvm 真正使用的 BoxConstraint、KernelScale 与 Epsilon最终用独立测试集的 RMSE、MAE、R²、残差和预测曲线共同判断模型质量。对于中小样本、变量较多、关系非线性的预测任务这套方案结构清晰、实现成本可控也容易扩展到电力负荷、设备状态、气象、交通、需求预测等实际数据。更重要的是模型效果应由可复现的验证设计和对照实验决定而不是由“优化算法”这个标签决定。只有把数据、验证、超参数和误差诊断四个环节一起做好WOA-SVR 才真正具有工程价值。参考资料[1] Mirjalili, S.; Lewis, A. The Whale Optimization Algorithm. Advances in Engineering Software, 2016, 95: 51–67.[2] MathWorks. fitrsvm — Fit a support vector machine regression model. MATLAB Statistics and Machine Learning Toolbox Documentation.[3] MathWorks. tspartition — Partition time series data for cross-validation. MATLAB Statistics and Machine Learning Toolbox Documentation. MATLAB实战｜WOA-SVM多变量时间序列预测：鲸鱼优化SVR超参数、时序验证、误差诊断与完整代码 多输入单输出 严格时序切分 · 可复现工程流程
MATLAB实战WOA-SVM多变量时间序列预测鲸鱼优化SVR超参数、时序验证、误差诊断与完整代码多输入单输出 · Gaussian SVR · Whale Optimization Algorithm · 严格时序切分 · 可复现工程流程摘要多变量时间序列预测的难点并不只是“选一个回归器”而是如何在不泄露未来信息的前提下将多源变量的滞后关系转化为稳定的监督学习样本并为非线性模型找到合适的超参数。本文构建一套完整的 WOA-SVM严格来说是 WOA-SVR预测流程首先使用滑动窗口把多变量序列展开为固定维度特征随后按时间顺序划分训练集、验证集和测试集再由鲸鱼优化算法WOA在对数空间搜索支持向量回归的 BoxConstraint、KernelScale 与 Epsilon最后对预测曲线、RMSE、MAE、R² 及残差分布进行联合评估。文末给出可直接运行的 MATLAB 完整脚本与本地函数并补充参数边界、数据泄露、RBF 参数映射、MAPE 失效场景以及真实业务数据替换方法。关键词MATLAB鲸鱼优化算法WOA支持向量机SVR多变量时间序列超参数优化滑动窗口运行环境建议 MATLAB R2021a 及以上版本需要 Statistics and Machine Learning Toolbox。核心建模函数使用 fitrsvm。若 MATLAB 版本支持 tspartition也可将文中的手工前向验证替换为 expanding-window 或 sliding-window 交叉验证。1. 问题定义为什么“WOA SVM”适合中小样本多变量时序在电力负荷、设备状态、交通流量、环境监测、价格需求等场景中目标变量通常同时受多个外部因素影响而且这些因素往往存在滞后、耦合、周期与噪声。直接对原始序列做普通回归容易遗漏时间依赖直接使用深度网络则可能在样本量有限时带来更高的训练成本和调参复杂度。支持向量回归SVR对小样本、高维和非线性关系具有较好的适应性但它对惩罚系数 C、核尺度 σMATLAB 中为 KernelScale以及 ε 不敏感区间非常敏感。参数过强会把噪声也拟合进去参数过弱又会产生明显欠拟合。WOA 的作用不是“替代模型”而是把这组超参数选择转化为连续优化问题通过群体搜索降低人工试参的盲目性。SVM/SVR 负责建立“历史多变量特征 → 未来目标值”的非线性映射WOA 负责在连续超参数空间中搜索更合适的模型配置严格的时间顺序切分负责阻断未来信息泄露标准化、窗口构造和残差诊断负责保证工程链路可解释、可复现。图 1 WOA-SVR 多变量时间序列预测完整工程链路2. 模型原理SVR、RBF 核与 WOA 如何连接起来2.1 ε-SVR 的核心目标在回归任务中SVM 通常以 ε-SVR 的形式使用。它希望找到一个尽量“平坦”的函数同时允许预测值落在真实值上下 ε 范围内而不产生惩罚。经典原始问题可以写为min 1/2 ||w||² C Σ(ξᵢ ξᵢ*)s.t. yᵢ - (wᵀφ(xᵢ)b) ≤ ε ξᵢ(wᵀφ(xᵢ)b) - yᵢ ≤ ε ξᵢ*ξᵢ, ξᵢ* ≥ 0其中 C 决定“模型平滑性”和“训练误差惩罚”之间的权衡ε 决定多大的误差可以被视为可接受范围。C 太大时模型可能过度追随训练样本C 太小时则可能欠拟合ε 太大时模型会忽略过多细节太小时支持向量数量增加、模型更容易受噪声影响。2.2 MATLAB 的 Gaussian KernelScale 不要与 γ 混用RBFGaussian核可写成K(xᵢ, xⱼ) exp( -||xᵢ - xⱼ||² / (2σ²) )很多资料习惯用 γ 表示 exp(-γ||xᵢ-xⱼ||²)此时 γ 1/(2σ²)。但 MATLAB fitrsvm 直接暴露的是 KernelScale也就是 σ。工程实现中最稳妥的做法是让 WOA 直接搜索 KernelScale而不是先搜索 γ 再做二次换算。这样既避免参数定义歧义也让代码与 MATLAB 官方接口完全一致。2.3 WOA 的三类位置更新机制鲸鱼优化算法由 Mirjalili 与 Lewis 于 2016 年提出模拟座头鲸气泡网捕食过程。每一条“鲸鱼”代表一组候选超参数适应度函数则对应这组参数在验证集上的预测误差。经典 WOA 主要包含以下三种更新行为包围猎物当 |A| 1 时个体向当前最优解靠拢增强局部开发能力随机搜索当 |A| ≥ 1 时个体参考随机鲸鱼更新位置扩大搜索范围螺旋气泡网以螺旋轨迹围绕当前最优解移动在局部区域继续精细搜索。A 2ar - a, C 2rD |C·X* - X|X(t1) X* - A·D螺旋更新X(t1) D · exp(bl) · cos(2πl) X*这里 a 会从 2 逐步衰减到 0使搜索过程从更强调全局探索逐渐过渡到局部开发。需要强调的是WOA 是随机元启发式算法它能改善搜索覆盖并降低人工调参成本但并不提供“必然找到全局最优”的数学保证。2.4 为什么在对数空间搜索更合理C、KernelScale、Epsilon 往往跨越多个数量级。如果在 [0.01, 1000] 这样的线性区间直接均匀初始化大量搜索点会堆积在大数值区域小数值区间被严重稀释。因此本文让 WOA 搜索 log10 参数z₁ log10(C) ∈ [-2, 3]z₂ log10(KernelScale) ∈ [-2, 2]z₃ log10(Epsilon) ∈ [-3, -0.7]真实参数C 10^z₁, KernelScale 10^z₂, Epsilon 10^z₃这样每一个数量级获得近似均衡的搜索机会通常比在原始尺度上“让鲸鱼在超大海域随机游动”更有效。3. 数据设计把多变量时间序列转成监督学习样本为了让示例可以独立运行本文构造 4 个输入变量与 1 个目标变量。x1 表示主要周期驱动x2 表示较慢周期因子x3 表示较高频的辅助波动x4 表示缓慢趋势与随机扰动目标 y 同时包含线性项、二次非线性项、差分项与噪声。图 2 多变量合成时间序列及目标变量示例3.1 滑动窗口用过去 6 个时刻预测当前时刻假设共有 4 个输入变量历史窗口 lag6则每个监督学习样本包含 6×424 个特征。对第 i 个目标时刻输入为 i-1 到 i-6 的四个变量输出为当前时刻 y(i)。这种做法把“序列预测”转换成普通回归问题同时显式保留局部时间依赖。[x1(i-1), x2(i-1), x3(i-1), x4(i-1),x1(i-2), x2(i-2), x3(i-2), x4(i-2),...x1(i-6), x2(i-6), x3(i-6), x4(i-6)] → y(i)3.2 时间序列不能随机打乱普通机器学习常使用随机 K 折交叉验证但时间序列存在严格因果顺序。如果把未来样本随机混入训练集就可能让模型间接看到未来分布得到过于乐观的误差。本文采用两层时间切分外层前 70% 样本作为训练段后 30% 作为最终测试段内层在外层训练段中再取前 80% 作为 WOA 训练子集后 20% 作为验证子集WOA 的适应度只看内层验证误差最终测试集在参数搜索期间完全不可见。更严格的方案当样本量足够时建议用 expanding-window 或 sliding-window 的多折前向验证替代单一验证段。其代价是每一个候选参数都要训练多次 SVR计算量会显著增加但评价方差更小。3.3 标准化必须只使用训练段统计量标准化本身也会泄露信息。正确流程是先切分再用训练段计算均值和标准差然后用同一组统计量变换验证集与测试集。绝不能先对全量数据求均值/标准差再切分训练测试。本文同时标准化 X 和 y使不同尺度的特征更适合 RBF 核并让 ε 的搜索范围更稳定。4. WOA-SVR 优化流程与适应度设计每个鲸鱼位置包含 3 个维度[log10(C), log10(KernelScale), log10(Epsilon)]。对于一个候选位置先还原真实参数再训练 SVR最后以验证集标准化 RMSE 作为适应度。适应度越小参数组合越好。参数设置说明种群规模 pop24兼顾搜索多样性与计算量最大迭代 maxIter45足够观察收敛过程log10(C)[-2, 3]C ∈ [0.01, 1000]log10(KernelScale)[-2, 2]σ ∈ [0.01, 100]log10(Epsilon)[-3, -0.7]ε 约 ∈ [0.001, 0.20]适应度验证集 RMSE不使用测试集图 3 WOA 搜索过程中验证集最优 RMSE 的变化如果收敛曲线前期快速下降、后期基本稳定说明搜索已经找到较稳定区域如果直到最后仍持续大幅下降通常需要增加迭代次数如果从第一代开始几乎不变则应检查搜索边界、适应度函数、数据标准化和参数尺度是否合理。5. 等价复现实验结果不只看一条预测曲线结果说明下表与图 3图 6 使用与文末 MATLAB 脚本一致的数据构造、滑动窗口、时间切分和超参数定义在等价实现中固定随机种子生成。不同随机数生成器、MATLAB 版本与优化轨迹会使最优参数和小数点后的指标略有差异因此运行本地脚本时应以实际输出为准。等价复现实验中WOA 搜索得到 C≈2.0409、KernelScale≈6.0921、Epsilon≈0.1880。模型RMSEMAER²SMAPE(%)持续性基线0.24610.19440.742362.04线性 SVR0.16780.12930.880246.87固定参数 RBF-SVR0.40610.33790.2987117.46WOA-SVR0.15970.12280.891543.79可以看到WOA-SVR 相比固定参数 RBF-SVR 的 RMSE 降低约 60.7%说明“RBF 核本身并不自动等于高精度”核尺度与惩罚系数的配置会显著改变结果。与此同时线性 SVR 在这一合成数据上也是一个很强的基线WOA-SVR 的 RMSE 进一步降低约 4.8%。这说明对比实验不能只放一个弱基线更不能把“智能优化”写成必然碾压所有模型。图 4 WOA-SVR 测试集预测曲线图 5 WOA-SVR 测试集残差分布图 6 持续性基线、线性 SVR、固定参数 RBF-SVR 与 WOA-SVR 的误差对比5.1 为什么这里不把 MAPE 作为主指标本示例目标 y 可能接近 0甚至出现负值。MAPE 的分母是 |y|当真实值接近 0 时少量绝对误差会被放大成极高百分比当 y0 时更是直接失去定义。因此本文以 RMSE、MAE、R² 为主并补充对零值更稳定的 SMAPE。对于电力负荷、销量等严格为正且远离 0 的目标MAPE 才更有解释价值。6. MATLAB 完整代码一份脚本直接跑通下面代码包含数据生成、滑动窗口、时间切分、WOA 三参数优化、最终模型训练、四类基线对比、误差指标与图形输出。将整段保存为一个 .m 文件即可运行本地函数放在脚本末尾。clear; clc; close all;rng(2025,twister);%% 1. 构造多变量时间序列n 600;t (1:n);x1 sin(2*pi*t/24) 0.15*randn(n,1);x2 cos(2*pi*t/48) 0.12*randn(n,1);x3 0.5*sin(2*pi*t/12) 0.08*randn(n,1);x4 0.03*t/n 0.10*randn(n,1);dx1 [0; diff(x1)];y 0.6*x1 0.25*x2.^2 - 0.3*x3 ...0.4*x4 0.05*dx1 0.10*randn(n,1);data [x1,x2,x3,x4,y];%% 2. 滑动窗口过去 lag 个时刻预测当前 ylag 6;numSamples n - lag;numInput 4;X zeros(numSamples, lag*numInput);Y zeros(numSamples, 1);targetTime (lag1:n);for k 1:numSamplesi k lag;feat zeros(1, lag*numInput);c 1;for j 1:lagfeat(c:cnumInput-1) data(i-j,1:numInput);c c numInput;endX(k,:) feat;Y(k) data(i,5);end%% 3. 外层时间切分前 70% 训练后 30% 测试trainRatio 0.70;numTrain floor(numSamples*trainRatio);XTrain X(1:numTrain,:);YTrain Y(1:numTrain);XTest X(numTrain1:end,:);YTest Y(numTrain1:end);TTest targetTime(numTrain1:end);%% 4. WOA 在对数空间优化 C、KernelScale、Epsilonpop 24;maxIter 45;% z [log10(C), log10(KernelScale), log10(Epsilon)]lb [-2, -2, -3];ub [ 3, 2, -0.7];dim 3;fitnessFcn (z) svrFitnessTimeHoldout(z, XTrain, YTrain, 0.80);[bestZ, bestFit, curve] WOA(fitnessFcn, pop, maxIter, lb, ub, dim);bestC 10.^bestZ(1);bestKernelScale 10.^bestZ(2);bestEpsilon 10.^bestZ(3);fprintf(Best validation RMSE (standardized) %.6f\n, bestFit);fprintf(Best BoxConstraint %.6f\n, bestC);fprintf(Best KernelScale %.6f\n, bestKernelScale);fprintf(Best Epsilon %.6f\n, bestEpsilon);%% 5. 只用完整训练集统计量做标准化[muX, sigX] trainZScoreStats(XTrain);XTrainN applyZScore(XTrain, muX, sigX);XTestN applyZScore(XTest, muX, sigX);muY mean(YTrain);sigY std(YTrain,0,1);if sigY 0, sigY 1; endYTrainN (YTrain - muY)./sigY;%% 6. 训练最优 WOA-SVRmdlWOA fitrsvm(XTrainN, YTrainN, ...KernelFunction,gaussian, ...BoxConstraint,bestC, ...KernelScale,bestKernelScale, ...Epsilon,bestEpsilon, ...Standardize,false);YPredWOAN predict(mdlWOA, XTestN);YPredWOA YPredWOAN*sigY muY;%% 7. 三个对照模型% 7.1 固定参数 RBF-SVRmdlFixed fitrsvm(XTrainN, YTrainN, ...KernelFunction,gaussian, ...BoxConstraint,1, ...KernelScale,1, ...Epsilon,0.05, ...Standardize,false);YPredFixed predict(mdlFixed,XTestN)*sigY muY;% 7.2 线性 SVRmdlLinear fitrsvm(XTrainN, YTrainN, ...KernelFunction,linear, ...BoxConstraint,1, ...Epsilon,0.05, ...Standardize,false);YPredLinear predict(mdlLinear,XTestN)*sigY muY;% 7.3 持续性基线用上一时刻真实 y 预测当前 yYPredPersistence y(TTest-1);%% 8. 计算指标M1 regressionMetrics(YTest, YPredPersistence);M2 regressionMetrics(YTest, YPredLinear);M3 regressionMetrics(YTest, YPredFixed);M4 regressionMetrics(YTest, YPredWOA);Result table(...[M1.RMSE;M2.RMSE;M3.RMSE;M4.RMSE], ...[M1.MAE; M2.MAE; M3.MAE; M4.MAE], ...[M1.R2; M2.R2; M3.R2; M4.R2], ...[M1.SMAPE;M2.SMAPE;M3.SMAPE;M4.SMAPE], ...VariableNames,{RMSE,MAE,R2,SMAPE}, ...RowNames,{Persistence,Linear_SVR,Fixed_RBF_SVR,WOA_SVR});disp(Result);%% 9. 可视化figure(Color,w,Name,WOA 收敛曲线);plot(curve,-o,LineWidth,1.2,MarkerSize,3);xlabel(迭代次数); ylabel(验证集标准化 RMSE);title(WOA 超参数搜索收敛曲线); grid on;figure(Color,w,Name,预测结果);plot(YTest,k,LineWidth,1.3); hold on;plot(YPredWOA,r--,LineWidth,1.3);legend(真实值,WOA-SVR 预测值,Location,best);xlabel(测试样本序号); ylabel(目标值);title(WOA-SVR 多变量时间序列预测); grid on;figure(Color,w,Name,残差分布);residual YTest - YPredWOA;histogram(residual,30);xline(0,--);xlabel(残差); ylabel(频数);title(WOA-SVR 测试集残差分布); grid on;figure(Color,w,Name,模型误差对比);bar([M1.RMSE M1.MAE; M2.RMSE M2.MAE; ...M3.RMSE M3.MAE; M4.RMSE M4.MAE]);set(gca,XTickLabel,{Persistence,Linear SVR,Fixed RBF,WOA-SVR});ylabel(误差); legend(RMSE,MAE,Location,best);title(不同模型误差对比); grid on;%% -------------------- 本地函数 --------------------function f svrFitnessTimeHoldout(z, X, Y, innerTrainRatio)% 只在外层训练集内部做时间顺序 holdoutN size(X,1);nTr floor(N*innerTrainRatio);Xtr X(1:nTr,:);Ytr Y(1:nTr);Xval X(nTr1:end,:);Yval Y(nTr1:end);[muX, sigX] trainZScoreStats(Xtr);XtrN applyZScore(Xtr, muX, sigX);XvalN applyZScore(Xval, muX, sigX);muY mean(Ytr);sigY std(Ytr,0,1);if sigY 0, sigY 1; endYtrN (Ytr-muY)./sigY;YvalN (Yval-muY)./sigY;C 10.^z(1);kernelScale 10.^z(2);epsilonSVR 10.^z(3);trymdl fitrsvm(XtrN,YtrN, ...KernelFunction,gaussian, ...BoxConstraint,C, ...KernelScale,kernelScale, ...Epsilon,epsilonSVR, ...Standardize,false);pred predict(mdl,XvalN);f sqrt(mean((pred-YvalN).^2));if ~isfinite(f), f 1e6; endcatchf 1e6;endendfunction [bestPos, bestFit, curve] WOA(fitnessFcn, pop, maxIter, lb, ub, dim)Xpos zeros(pop,dim);fit zeros(pop,1);for i 1:popXpos(i,:) lb rand(1,dim).*(ub-lb);fit(i) fitnessFcn(Xpos(i,:));end[bestFit,idx] min(fit);bestPos Xpos(idx,:);curve zeros(maxIter,1);for iter 1:maxItera 2 - iter*(2/maxIter);for i 1:popr1 rand();r2 rand();A 2*a*r1 - a;Ccoef 2*r2;p rand();b 1;l -1 2*rand();if p 0.5if abs(A) 1D abs(Ccoef*bestPos - Xpos(i,:));Xnew bestPos - A.*D;elserandIdx randi(pop);Xrand Xpos(randIdx,:);D abs(Ccoef*Xrand - Xpos(i,:));Xnew Xrand - A.*D;endelseD abs(bestPos - Xpos(i,:));Xnew D.*exp(b*l).*cos(2*pi*l) bestPos;endXnew max(Xnew,lb);Xnew min(Xnew,ub);newFit fitnessFcn(Xnew);if newFit fit(i)Xpos(i,:) Xnew;fit(i) newFit;endif fit(i) bestFitbestFit fit(i);bestPos Xpos(i,:);endendcurve(iter) bestFit;endendfunction [muX, sigX] trainZScoreStats(X)muX mean(X,1);sigX std(X,0,1);sigX(sigX0) 1;endfunction XN applyZScore(X, muX, sigX)XN (X-muX)./sigX;endfunction M regressionMetrics(yTrue, yPred)err yTrue-yPred;M.RMSE sqrt(mean(err.^2));M.MAE mean(abs(err));denom sum((yTrue-mean(yTrue)).^2);if denom 0M.R2 NaN;elseM.R2 1-sum(err.^2)/denom;endM.SMAPE mean(2*abs(err)./(abs(yTrue)abs(yPred)1e-8))*100;end7. 真实业务数据替换从合成数据迁移到 CSV / Excel真正落地时最重要的不是“把 data 换成 readmatrix 一行就结束”而是明确时间列、输入列、目标列、采样频率与缺失处理。推荐先把数据整理成第一列时间戳后续若干列为外生输入最后一列为预测目标。T readtable(your_timeseries.csv);% 示例假设第 1 列为时间第 2:5 列为输入第 6 列为目标Xraw T{:,2:5};yraw T{:,6};% 缺失值示例线性插值真实项目应结合业务含义选择方法Xraw fillmissing(Xraw,linear);yraw fillmissing(yraw,linear);data [Xraw,yraw];n size(data,1);如果多个传感器采样频率不同应先按时间戳重采样和对齐如果目标存在明显日周期、周周期可加入 hour、weekday、holiday 等时间特征如果某些输入变量只在未来已知例如天气预报要明确预测时刻能否真实获得这些变量如果存在概念漂移应采用滚动训练或定期重估超参数而不是永久固定一次 WOA 结果。7.1 多步预测怎么改本文代码是一阶 ahead 的单步预测。如果要预测未来 12 个时刻可以选择两种常见策略Direct为每个 horizon 分别训练一个 SVR例如 h1、2、…、12 各一个模型。优点是不传播前一步预测误差缺点是模型数量增加Recursive只训练一个一步模型然后把前一步预测值回填继续预测。模型简单但误差会逐步累积。8. 常见问题与避坑这些细节比“多迭代 20 次”更重要问题 1先归一化全数据再划分训练测试。这是最隐蔽的数据泄露之一。均值和标准差必须只从训练段计算再应用到验证和测试。问题 2为了“样本均匀”把时间序列随机打乱。这会破坏因果顺序。最终测试必须是时间上晚于训练段的数据。问题 3把 γ 直接当成 MATLAB KernelScale。两者定义不同。本文直接优化 KernelScale避免不必要的换算。问题 4WOA 的适应度直接用测试集 RMSE。这样测试集就参与了模型选择最终误差失去独立性。适应度必须使用训练段内部验证集。问题 5搜索范围越大越好。范围过大反而浪费搜索预算。先根据标准化后的数据给出合理数量级再在对数空间搜索。问题 6只报告一次随机种子的结果。WOA 是随机算法。严谨实验应使用多个随机种子重复运行报告均值、标准差与最优值。问题 7只看 R²不看误差曲线与残差。R² 高并不代表每个关键时段都预测稳定。应同时检查峰值段、残差偏移、异常点与趋势滞后。问题 8认为 WOA 能弥补差的特征工程。优化器只能在给定特征与模型空间内搜索。窗口长度、变量质量、时序对齐和异常处理仍然决定上限。9. WOA-SVR 的适用边界什么时候值得用什么时候不必用方法更适合的场景主要优势主要代价WOA-SVR中小样本、非线性、超参数空间连续实现直观、可解释、无需梯度随机性、计算成本随 pop×iter 增长网格搜索 SVR参数范围很窄、维度少结果稳定、容易复现网格细时计算爆炸连续空间利用率低Bayesian Optimization单次训练较贵、超参数不多样本效率高、工具链成熟实现逻辑比 WOA 更复杂LSTM/GRU/Transformer大样本、长依赖、多步序列表达能力强数据量、训练成本和调参复杂度更高如果数据只有几百到几千个样本、特征维度中等、希望快速得到一个稳健的非线性回归基线同时又不希望手工尝试大量 C/KernelScale/Epsilon 组合那么 WOA-SVR 很有价值。反过来如果样本已经达到几十万甚至更大规模标准核 SVM 的训练复杂度可能成为瓶颈此时应考虑近似核方法、树模型或更适合大数据的深度学习方案。10. 进一步增强从“能跑”到“更稳健”的四个方向1.把窗口长度 lag 也加入优化变量但必须把整数变量做取整处理并注意搜索维度增加会显著放大计算量2.用 expanding-window 多折验证代替单一内层 holdout使适应度更能代表跨时间段稳定性3.对特征加入滞后差分、移动均值、周期编码或领域变量并用消融实验验证每类特征的真实贡献4.对 WOA 做多随机种子重复实验并把“平均 RMSE ± 标准差、最优 RMSE、运行时间”一起报告避免只展示一次幸运搜索。11. 总结WOA-SVM 多变量时间序列预测的关键不在于把两个算法名字拼在一起而在于把整个建模链路做对先用滑动窗口显式表达历史依赖再用严格时间切分隔离未来数据标准化只拟合训练段WOA 在对数空间优化 fitrsvm 真正使用的 BoxConstraint、KernelScale 与 Epsilon最终用独立测试集的 RMSE、MAE、R²、残差和预测曲线共同判断模型质量。对于中小样本、变量较多、关系非线性的预测任务这套方案结构清晰、实现成本可控也容易扩展到电力负荷、设备状态、气象、交通、需求预测等实际数据。更重要的是模型效果应由可复现的验证设计和对照实验决定而不是由“优化算法”这个标签决定。只有把数据、验证、超参数和误差诊断四个环节一起做好WOA-SVR 才真正具有工程价值。参考资料[1] Mirjalili, S.; Lewis, A. The Whale Optimization Algorithm. Advances in Engineering Software, 2016, 95: 51–67.[2] MathWorks. fitrsvm — Fit a support vector machine regression model. MATLAB Statistics and Machine Learning Toolbox Documentation.[3] MathWorks. tspartition — Partition time series data for cross-validation. MATLAB Statistics and Machine Learning Toolbox Documentation.