 算法实现)
简介这份PDF文档围绕多变量灰色预测模型的建模方法与Matlab算法实现展开面向需处理小样本、多变量相互影响数据的科研人员、研究生与工程技术人员。内容从一次累加生成序列切入推导动态微分方程组给出离散化后基于最小二乘法的参数估计公式并说明时间响应函数求解、还原预测值以及残差均方差比值检验的完整流程。包内仅含1个PDF文件约122KB篇幅精炼核心是算法步骤与可直接运行的M程序片段涵盖数据累加、矩阵D构造、参数辨识与预测值计算等环节便于读者对照复现。目前已有154人学习。读者可借此掌握多变量灰色模型的编程实现思路并将其迁移至经济预测、工程数据分析等实际场景快速完成建模与精度评估。1. 多变量灰色预测模型适合什么样的数据场景做季度经营预算、区域用电量申报、连锁门店销量铺货这类任务的共同点是要预测的那个指标从来不孤立。气温拉动用电、客流量牵引销量、政策因素改变节奏单变量 GM(1,1) 只看一条曲线的自身惯性一旦外部变量出现拐点整条预测曲线就会整体跑偏。MGM(1,n) 把 n 个相关指标放进同一个一阶微分方程组用矩阵形式让变量之间互相提供约束小样本条件下比逐条建模更有抗扰动能力。在 MATLAB 里落地这套多变量灰色预测模型算法门槛不在数学推导而在四步矩阵维度对齐累加生成、最小二乘辨识、时间响应式求解、精度检验。维度错一位expm吐出来的就是一堆垃圾数而且不会报错。适合读下去的有两类人一类是已经用过 GM(1,1)想把单变量扩成多变量的人另一类手里压着几列强相关的历史数据正在犹豫先上 BP 神经网络拟合曲线还是先拿灰色模型把趋势踩准。2. MGM(1,n) 的建模链路从累加生成到白化方程2.1 为什么多变量灰色模型要先做一次累加生成灰色预测的底层假设是原始序列虽然看上去杂乱但它的累加序列近似服从指数规律。道理不玄乎——对任意弱增长序列做一次累加随机波动被平均掉曲线会变得单调平滑接近指数形状而指数形状正好是一阶线性微分方程的解这就给建模找到了数学落点。MGM(1,n) 沿用这个思路只是把单条序列的累加推广成矩阵运算。给定 m 个时刻、n 个变量的原始矩阵 X⁽⁰⁾尺寸 m×n沿时间维做 cumsum得到累加矩阵 X⁽¹⁾。方向必须写对MATLAB 里是cumsum(X0, 1)第二个参数 1 表示沿行方向累加。如果数据是按列排的每列一个时刻就得先转置再累加写成cumsum(X0)。这类错误不报错只会让后面的 A、B 全部失真。2.2 准光滑性与级比检验建模前的两道门槛累加能不能用得先让数据自己说话。四个指标最常看检验项计算式通过参考区间不通过时的处理准光滑比ρ(k)x⁽⁰⁾(k)/x⁽¹⁾(k-1)k3 时 ρ(k)0.5延长序列或加常数平移级比σ(k)x⁽¹⁾(k)/x⁽¹⁾(k-1)落在 (1, 1.5)加常数 c 使数据整体上移数据符号全序列非负x(k) ≥ 0平移到正数域再建模采样间隔等时间距Δt 恒定先插值或重采样成等距准光滑比直觉上回答的是“新增量相对历史累积量够不够小”——越小越说明序列趋近指数形态。级比则衡量累加序列的相邻比值是否落在指数函数的合理增长区间。下面这段代码把两项检验一次性算出来% 数据检验准光滑比与级比 X0 [ ...... ]; % m×n 原始矩阵行是时刻列是变量 X1 cumsum(X0, 1); % 一次累加生成沿时间维 rho X0(3:end,:) ./ X1(2:end-1,:); % 准光滑比长度 m-2 sigma X1(2:end,:) ./ X1(1:end-1,:); % 级比长度 m-1 fprintf(准光滑比最大值%.3f\n, max(rho(:))); fprintf(级比区间[%.3f, %.3f]\n, min(sigma(:)), max(sigma(:))); if max(rho(:)) 0.5 c 0.1 * max(abs(X0(:))); % 平移量一般取量级的一成 warning(准光滑性不满足建议 X0 X0 %.3f 后重试, c); end逻辑说明rho从第 3 个时刻开始算因为前两个时刻的累加值太小比值失真。平移量取0.1*max(abs(X0))是个保守经验值太小起不到作用太大会把原始差异压平。参数说明X0(3:end,:)与X1(2:end-1,:)长度都是 m-2维度刚好对齐这一步如果报“矩阵维度不一致”八成是原始数据第一列不是时刻标识、被误当成了变量。2.3 白化微分方程组与离散化的对应关系MGM(1,n) 的白化方程写成矩阵形式dX⁽¹⁾/dt A·X⁽¹⁾ BA 是 n×n 的发展系数矩阵B 是 n×1 的灰作用量向量。A 的对角元反映每个变量自身的增长惯性非对角元反映变量之间的耦合强度——这也是多变量灰色预测模型和逐条 GM(1,1) 最本质的区别后者的 A 只有对角线根本表达不了变量间的牵引。离散化时对导数项用前向差分对 X⁽¹⁾ 用紧邻均值X⁽¹⁾(k1) − X⁽¹⁾(k) A·Z⁽¹⁾(k) B其中 Z⁽¹⁾(k) (X⁽¹⁾(k1) X⁽¹⁾(k)) / 2这一步是整个算法的枢纽连续方程被换成一个线性方程组每个时刻贡献一行于是能一次性用最小二乘解出 A 和 B。紧邻均值这一步的直觉是——用区间中点的平均值来代表这段区间的“平均水平”比直接用左端点更贴合积分效果。3. 用 MATLAB 辨识参数矩阵 A 与 B3.1 构造紧邻均值矩阵与设计矩阵把离散方程对第 i 个变量展开x_i⁽¹⁾(k1) − x_i⁽¹⁾(k) a_i1·z_1(k) a_i2·z_2(k) … a_in·z_n(k) b_i对 k1..m-1 共 m-1 个时刻就得到 m-1 个方程、n1 个未知数。方程数必须不少于未知数也就是 m-1 ≥ n1即样本时刻数至少比变量数多两个。这是很多人建模卡住的第一道坎手上只有 6 期数据却想建 5 变量模型方程数根本不够\会给一个无意义的基础解。把 n 个变量拼起来设计矩阵 D 是所有时刻的紧邻均值再加一列 1D [Z1, ones(m-1,1)]尺寸 (m-1)×(n1)。待求矩阵 Θ 尺寸 (n1)×nY 的尺寸是 (m-1)×n每列对应一个变量的相邻差分。三者的维度必须死记行是时刻列是变量或系数。3.2 最小二乘辨识 A、B 与维度对齐function [A, B, X1] mgm_fit(X0) % MGM_FIT 拟合多变量灰色模型 MGM(1,n) 的参数 % 输入X0 m×n 原始数据m 为时刻数n 为变量数 % 输出A n×n 发展系数矩阵 % B n×1 灰作用量向量 % X1 m×n 一次累加生成序列 [m, n] size(X0); if m n 2 error(样本时刻数 m%d 不足至少需要 n2%d 个时刻, m, n2); end X1 cumsum(X0, 1); % 沿时间维累加 Z1 (X1(1:end-1,:) X1(2:end,:)) / 2; % 紧邻均值(m-1)×n Y X1(2:end,:) - X1(1:end-1,:); % 差分(m-1)×n D [Z1, ones(m-1, 1)]; % 设计矩阵(m-1)×(n1) Theta D \ Y; % 最小二乘(n1)×n A Theta(1:n, :); % 转置成 n×n B Theta(n1, :); % 1×n 转成 n×1 % 条件数检查D 接近奇异时参数不可信 c cond(D); if c 1e10 warning(设计矩阵条件数 %.2e 过大变量间可能存在强共线性, c); end end逻辑说明D \ Y走 QR 分解比手写inv(D*D)*D*Y数值上稳得多后者会把条件数平方掉。Theta(1:n,:)这一步转置最容易写错——最小二乘按列求解第 j 个变量的系数恰好摆在 Theta 的第 j 列上转置后才是按行摆放的 A。参数说明X0要求全为正数出现零或负数先整体平移。A的规模直接决定模型能容纳多少变量n越大、需要的样本时刻越多。cond(D)超过 1e10 时A 的非对角元会出现数量级上的抖动这时候应该先做标准化。3.3 时间响应式求解expm 与 A 的奇异性处理连续方程在初值 X⁽¹⁾(1) 下的解是矩阵指数形式X⁽¹⁾(k) e^{A(k−1)}·X⁽¹⁾(1) A⁻¹·(e^{A(k−1)} − I)·BMATLAB 里必须用expm而不是exp。后者是逐元素指数前者才是矩阵指数写错一个字母结果完全不同而且不报任何错。这是灰色模型算法在 MATLAB 实现里最常见的一个隐形坑。函数含义用错后的表现exp(A)逐元素 e^{a_ij}数值看似正常预测全线偏移expm(A)矩阵指数 Σ A^k/k!正确A\I求逆等价写法A 奇异时报警告pinv(A)伪逆A 奇异时的兜底方案function X1p mgm_predict(A, B, X1_1, steps) % MGM_PREDICT 由时间响应式外推累加序列 % X1_1 1×n 第一期累加值 % steps 需要外推的步数不含起始点 % X1p (steps1)×n 累加序列预测值 n size(A, 1); X1p zeros(steps 1, n); X1p(1, :) X1_1; I eye(n); % A 接近奇异时改用伪逆避免 A\... 报奇异矩阵警告 if rcond(A) 1e-12 Ainv pinv(A); % 伪逆兜底 else Ainv A \ I; % 等价于求逆数值更稳 end for k 1:steps E expm(A * k); % 直接算矩阵指数避免反复累乘的舍入漂移 X1p(k1, :) (E * X1_1 Ainv * (E - I) * B); end end逻辑说明每一步都重新计算expm(A*k)比反复左乘expm(A)累积的舍入误差更小尤其在 A 的特征值跨度大时。rcond(A) 1e-12判断的是 A 相对单位阵的接近奇异程度比直接比行列式可靠。参数说明steps是外推步数做 3 期预测就传 3。X1_1必须传入原始数据第一期的累加值也就是X1(1,:)不是第一期原始值。这是另一个高频错点。4. 预测还原、精度检验与残差修正4.1 累减还原与维度补齐模型输出的是累加序列 X⁽¹⁾必须累减还原到原始量纲x⁽⁰⁾(k) x⁽¹⁾(k) − x⁽¹⁾(k−1)X1p mgm_predict(A, B, X1(1,:), 5); % 向后外推 5 步 X0p [X1p(1,:); diff(X1p)]; % 累减还原(steps1)×n逻辑说明diff会把行数减 1所以必须手动补回第一行。X1p(1,:)对应的就是原始序列第一期的值丢弃它会让整条预测序列错位一步后面所有残差检验跟着全错。参数说明还原后的X0p与原始X0处于同一量纲可直接做残差运算如果之前在 3.2 里对数据做了标准化这里要按原均值和标准差反变换回去。4.2 后验差比值 C 与小误差概率 P两个核心指标是 C 和 PC 是残差标准差比原始序列标准差越大说明波动没被模型吃掉P 是残差落在 ±0.6745σ₁ 内的比例越大说明误差分布集中。等级小误差概率 P后验差比值 C模型状态一级P 0.95C 0.35可直接外推二级0.80 P ≤ 0.950.35 ≤ C 0.50适合短期预测三级0.70 P ≤ 0.800.50 ≤ C 0.65只做趋势参考四级P ≤ 0.70C ≥ 0.65参数不可用需重构function [C, P, level] mgm_check(X0, X0p) % 后验差检验 % X0 m×n 原始序列 % X0p m×n 模型拟合值与 X0 同长度 E X0 - X0p; % 残差矩阵 S1 std(X0, 0, 1); % 原始序列标准差1×n S2 std(E, 0, 1); % 残差标准差1×n C S2 ./ S1; % 后验差比值1×n thr 0.6745 * S1; % 小误差阈值 P mean(abs(E - mean(E,1)) thr, 1); % 小误差概率1×n level strings(1, size(X0,2)); for j 1:size(X0,2) if P(j) 0.95 C(j) 0.35; level(j) 一级; elseif P(j) 0.80 C(j) 0.50; level(j) 二级; elseif P(j) 0.70 C(j) 0.65; level(j) 三级; else; level(j) 四级; end end end逻辑说明std(X, 0, 1)的第二个参数 0 表示按样本标准差除以 N-1第三个参数 1 表示沿行方向统计也就是沿时间维。两个 1 都不写会默认沿第一维统计当 n1 时结果就是错的。mean(abs(...) thr, 1)用逻辑矩阵求均值得到的就是满足条件的时刻占比。参数说明X0和X0p必须包含相同的时刻范围。做样本内检验时用全部拟合值做样本外检验时X0 要裁到对应的外推区间否则 S1 会被训练期的波动稀释C 值虚高。注意C 和 P 是逐变量给出的只要有一个变量掉到四级整个模型的联合外推就不可信因为 A 的非对角元会让误差在各变量间传导。4.3 量纲差异、采样间隔与负值数据的排错实测中四个坑反复出现量纲差异过大。气温用摄氏度、GDP 用亿元紧邻均值里大数会把小数淹掉A 直接病态。做法是先zscore标准化再建模预测完反标准化。这在多变量场景下几乎是必做步骤。非等距采样。灰色模型隐含等时间距假设跨月缺测、跨年数据直接用会让 Z1 失真。先用interp1插成等距序列插值后要复核准光滑比重跑一遍。负值或零值。需求下滑、利润为负时累加序列不再单调级比检验直接失败。加常数 c 平移c 一般取1.1*max(abs(X0))预测完再减回去。变量数逼近样本数。m-1 与 n1 只差一两个时即便能解出来参数的方差也极大换一期数据 A 就翻脸。经验做法是变量数控制在(m-1)/3以内。对应的排查代码% 标准化后再建模避免量纲差异导致 A 病态 mu mean(X0, 1); sg std(X0, 0, 1); Xs (X0 - mu) ./ sg; % 逐列标准化 [As, Bs, X1s] mgm_fit(Xs 1); % 平移到正数域后再拟合 X1p mgm_predict(As, Bs, X1s(1,:), 5); X0p [X1p(1,:); diff(X1p)]; X0p (X0p - 1) .* sg mu; % 反标准化回原量纲逻辑说明加 1 平移是为了让标准化后的数据落进正数域减 1 时注意只在还原阶段减一次。sg和mu必须用训练期数据算出来不能混入外推期的统计量否则就是样本泄露。5. 背景值优化、滚动预测与 MATLAB 优化工具箱配合默认紧邻均值用 λ0.5也就是 Z1 0.5·X1(k1) 0.5·X1(k)。这个取值是几何中点但在序列增长较快时会引入系统性偏差。把 λ 当参数、以样本外 MAPE 为目标函数用fminsearch或粒子群算法搜索通常能把 C 值压下去一档。% 以样本外 MAPE 为目标搜索背景值权重 lambda obj (lam) mgm_mape(X0, lam); lam_opt fminsearch(obj, 0.5, ... optimset(TolX, 1e-4, MaxFunEvals, 200)); fprintf(最优背景值权重 lambda %.4f\n, lam_opt); function m mgm_mape(X0, lam) [~,~,~,A,B,X1] mgm_fit_lambda(X0, lam); % 内部使用 lam 构造 Z1 X1p mgm_predict(A, B, X1(1,:), size(X0,1)-1); X0p [X1p(1,:); diff(X1p)]; m mean(abs((X0(2:end,:) - X0p(2:end,:)) ./ X0(2:end,:)), all) * 100; endfminsearch属于 MATLAB 优化工具箱的基础函数无约束、无梯度对小参数问题够用变量多、目标函数非凸时换particleswarm但那需要全局优化工具箱。目标函数里算的是样本外一步误差注意别把训练区间一起丢进去算 MAPE否则 λ 会过拟合到训练期。另一个工程化习惯是滚动预测每来一期真实数据就整体重拟合一次 A、B而不是拿一个 A 硬撑十二期。多变量灰色预测模型的 A 是耦合矩阵任一变量出新拐点时非对角元会跟着变滚动重拟合是唯一稳健的做法。最后一段可复现的验证方式把最后两期留作样本外用前 m-2 期拟合、向后外推两步比对 MAPE 是否落在训练期 MAPE 的 1.5 倍以内。超出这个比例说明 A 的非对角元已经不稳定先回头检查标准化和变量数是否过密。本文还有配套的精品资源点击获取
MATLAB 多变量灰色预测 MGM(1,n) 算法实现
简介这份PDF文档围绕多变量灰色预测模型的建模方法与Matlab算法实现展开面向需处理小样本、多变量相互影响数据的科研人员、研究生与工程技术人员。内容从一次累加生成序列切入推导动态微分方程组给出离散化后基于最小二乘法的参数估计公式并说明时间响应函数求解、还原预测值以及残差均方差比值检验的完整流程。包内仅含1个PDF文件约122KB篇幅精炼核心是算法步骤与可直接运行的M程序片段涵盖数据累加、矩阵D构造、参数辨识与预测值计算等环节便于读者对照复现。目前已有154人学习。读者可借此掌握多变量灰色模型的编程实现思路并将其迁移至经济预测、工程数据分析等实际场景快速完成建模与精度评估。1. 多变量灰色预测模型适合什么样的数据场景做季度经营预算、区域用电量申报、连锁门店销量铺货这类任务的共同点是要预测的那个指标从来不孤立。气温拉动用电、客流量牵引销量、政策因素改变节奏单变量 GM(1,1) 只看一条曲线的自身惯性一旦外部变量出现拐点整条预测曲线就会整体跑偏。MGM(1,n) 把 n 个相关指标放进同一个一阶微分方程组用矩阵形式让变量之间互相提供约束小样本条件下比逐条建模更有抗扰动能力。在 MATLAB 里落地这套多变量灰色预测模型算法门槛不在数学推导而在四步矩阵维度对齐累加生成、最小二乘辨识、时间响应式求解、精度检验。维度错一位expm吐出来的就是一堆垃圾数而且不会报错。适合读下去的有两类人一类是已经用过 GM(1,1)想把单变量扩成多变量的人另一类手里压着几列强相关的历史数据正在犹豫先上 BP 神经网络拟合曲线还是先拿灰色模型把趋势踩准。2. MGM(1,n) 的建模链路从累加生成到白化方程2.1 为什么多变量灰色模型要先做一次累加生成灰色预测的底层假设是原始序列虽然看上去杂乱但它的累加序列近似服从指数规律。道理不玄乎——对任意弱增长序列做一次累加随机波动被平均掉曲线会变得单调平滑接近指数形状而指数形状正好是一阶线性微分方程的解这就给建模找到了数学落点。MGM(1,n) 沿用这个思路只是把单条序列的累加推广成矩阵运算。给定 m 个时刻、n 个变量的原始矩阵 X⁽⁰⁾尺寸 m×n沿时间维做 cumsum得到累加矩阵 X⁽¹⁾。方向必须写对MATLAB 里是cumsum(X0, 1)第二个参数 1 表示沿行方向累加。如果数据是按列排的每列一个时刻就得先转置再累加写成cumsum(X0)。这类错误不报错只会让后面的 A、B 全部失真。2.2 准光滑性与级比检验建模前的两道门槛累加能不能用得先让数据自己说话。四个指标最常看检验项计算式通过参考区间不通过时的处理准光滑比ρ(k)x⁽⁰⁾(k)/x⁽¹⁾(k-1)k3 时 ρ(k)0.5延长序列或加常数平移级比σ(k)x⁽¹⁾(k)/x⁽¹⁾(k-1)落在 (1, 1.5)加常数 c 使数据整体上移数据符号全序列非负x(k) ≥ 0平移到正数域再建模采样间隔等时间距Δt 恒定先插值或重采样成等距准光滑比直觉上回答的是“新增量相对历史累积量够不够小”——越小越说明序列趋近指数形态。级比则衡量累加序列的相邻比值是否落在指数函数的合理增长区间。下面这段代码把两项检验一次性算出来% 数据检验准光滑比与级比 X0 [ ...... ]; % m×n 原始矩阵行是时刻列是变量 X1 cumsum(X0, 1); % 一次累加生成沿时间维 rho X0(3:end,:) ./ X1(2:end-1,:); % 准光滑比长度 m-2 sigma X1(2:end,:) ./ X1(1:end-1,:); % 级比长度 m-1 fprintf(准光滑比最大值%.3f\n, max(rho(:))); fprintf(级比区间[%.3f, %.3f]\n, min(sigma(:)), max(sigma(:))); if max(rho(:)) 0.5 c 0.1 * max(abs(X0(:))); % 平移量一般取量级的一成 warning(准光滑性不满足建议 X0 X0 %.3f 后重试, c); end逻辑说明rho从第 3 个时刻开始算因为前两个时刻的累加值太小比值失真。平移量取0.1*max(abs(X0))是个保守经验值太小起不到作用太大会把原始差异压平。参数说明X0(3:end,:)与X1(2:end-1,:)长度都是 m-2维度刚好对齐这一步如果报“矩阵维度不一致”八成是原始数据第一列不是时刻标识、被误当成了变量。2.3 白化微分方程组与离散化的对应关系MGM(1,n) 的白化方程写成矩阵形式dX⁽¹⁾/dt A·X⁽¹⁾ BA 是 n×n 的发展系数矩阵B 是 n×1 的灰作用量向量。A 的对角元反映每个变量自身的增长惯性非对角元反映变量之间的耦合强度——这也是多变量灰色预测模型和逐条 GM(1,1) 最本质的区别后者的 A 只有对角线根本表达不了变量间的牵引。离散化时对导数项用前向差分对 X⁽¹⁾ 用紧邻均值X⁽¹⁾(k1) − X⁽¹⁾(k) A·Z⁽¹⁾(k) B其中 Z⁽¹⁾(k) (X⁽¹⁾(k1) X⁽¹⁾(k)) / 2这一步是整个算法的枢纽连续方程被换成一个线性方程组每个时刻贡献一行于是能一次性用最小二乘解出 A 和 B。紧邻均值这一步的直觉是——用区间中点的平均值来代表这段区间的“平均水平”比直接用左端点更贴合积分效果。3. 用 MATLAB 辨识参数矩阵 A 与 B3.1 构造紧邻均值矩阵与设计矩阵把离散方程对第 i 个变量展开x_i⁽¹⁾(k1) − x_i⁽¹⁾(k) a_i1·z_1(k) a_i2·z_2(k) … a_in·z_n(k) b_i对 k1..m-1 共 m-1 个时刻就得到 m-1 个方程、n1 个未知数。方程数必须不少于未知数也就是 m-1 ≥ n1即样本时刻数至少比变量数多两个。这是很多人建模卡住的第一道坎手上只有 6 期数据却想建 5 变量模型方程数根本不够\会给一个无意义的基础解。把 n 个变量拼起来设计矩阵 D 是所有时刻的紧邻均值再加一列 1D [Z1, ones(m-1,1)]尺寸 (m-1)×(n1)。待求矩阵 Θ 尺寸 (n1)×nY 的尺寸是 (m-1)×n每列对应一个变量的相邻差分。三者的维度必须死记行是时刻列是变量或系数。3.2 最小二乘辨识 A、B 与维度对齐function [A, B, X1] mgm_fit(X0) % MGM_FIT 拟合多变量灰色模型 MGM(1,n) 的参数 % 输入X0 m×n 原始数据m 为时刻数n 为变量数 % 输出A n×n 发展系数矩阵 % B n×1 灰作用量向量 % X1 m×n 一次累加生成序列 [m, n] size(X0); if m n 2 error(样本时刻数 m%d 不足至少需要 n2%d 个时刻, m, n2); end X1 cumsum(X0, 1); % 沿时间维累加 Z1 (X1(1:end-1,:) X1(2:end,:)) / 2; % 紧邻均值(m-1)×n Y X1(2:end,:) - X1(1:end-1,:); % 差分(m-1)×n D [Z1, ones(m-1, 1)]; % 设计矩阵(m-1)×(n1) Theta D \ Y; % 最小二乘(n1)×n A Theta(1:n, :); % 转置成 n×n B Theta(n1, :); % 1×n 转成 n×1 % 条件数检查D 接近奇异时参数不可信 c cond(D); if c 1e10 warning(设计矩阵条件数 %.2e 过大变量间可能存在强共线性, c); end end逻辑说明D \ Y走 QR 分解比手写inv(D*D)*D*Y数值上稳得多后者会把条件数平方掉。Theta(1:n,:)这一步转置最容易写错——最小二乘按列求解第 j 个变量的系数恰好摆在 Theta 的第 j 列上转置后才是按行摆放的 A。参数说明X0要求全为正数出现零或负数先整体平移。A的规模直接决定模型能容纳多少变量n越大、需要的样本时刻越多。cond(D)超过 1e10 时A 的非对角元会出现数量级上的抖动这时候应该先做标准化。3.3 时间响应式求解expm 与 A 的奇异性处理连续方程在初值 X⁽¹⁾(1) 下的解是矩阵指数形式X⁽¹⁾(k) e^{A(k−1)}·X⁽¹⁾(1) A⁻¹·(e^{A(k−1)} − I)·BMATLAB 里必须用expm而不是exp。后者是逐元素指数前者才是矩阵指数写错一个字母结果完全不同而且不报任何错。这是灰色模型算法在 MATLAB 实现里最常见的一个隐形坑。函数含义用错后的表现exp(A)逐元素 e^{a_ij}数值看似正常预测全线偏移expm(A)矩阵指数 Σ A^k/k!正确A\I求逆等价写法A 奇异时报警告pinv(A)伪逆A 奇异时的兜底方案function X1p mgm_predict(A, B, X1_1, steps) % MGM_PREDICT 由时间响应式外推累加序列 % X1_1 1×n 第一期累加值 % steps 需要外推的步数不含起始点 % X1p (steps1)×n 累加序列预测值 n size(A, 1); X1p zeros(steps 1, n); X1p(1, :) X1_1; I eye(n); % A 接近奇异时改用伪逆避免 A\... 报奇异矩阵警告 if rcond(A) 1e-12 Ainv pinv(A); % 伪逆兜底 else Ainv A \ I; % 等价于求逆数值更稳 end for k 1:steps E expm(A * k); % 直接算矩阵指数避免反复累乘的舍入漂移 X1p(k1, :) (E * X1_1 Ainv * (E - I) * B); end end逻辑说明每一步都重新计算expm(A*k)比反复左乘expm(A)累积的舍入误差更小尤其在 A 的特征值跨度大时。rcond(A) 1e-12判断的是 A 相对单位阵的接近奇异程度比直接比行列式可靠。参数说明steps是外推步数做 3 期预测就传 3。X1_1必须传入原始数据第一期的累加值也就是X1(1,:)不是第一期原始值。这是另一个高频错点。4. 预测还原、精度检验与残差修正4.1 累减还原与维度补齐模型输出的是累加序列 X⁽¹⁾必须累减还原到原始量纲x⁽⁰⁾(k) x⁽¹⁾(k) − x⁽¹⁾(k−1)X1p mgm_predict(A, B, X1(1,:), 5); % 向后外推 5 步 X0p [X1p(1,:); diff(X1p)]; % 累减还原(steps1)×n逻辑说明diff会把行数减 1所以必须手动补回第一行。X1p(1,:)对应的就是原始序列第一期的值丢弃它会让整条预测序列错位一步后面所有残差检验跟着全错。参数说明还原后的X0p与原始X0处于同一量纲可直接做残差运算如果之前在 3.2 里对数据做了标准化这里要按原均值和标准差反变换回去。4.2 后验差比值 C 与小误差概率 P两个核心指标是 C 和 PC 是残差标准差比原始序列标准差越大说明波动没被模型吃掉P 是残差落在 ±0.6745σ₁ 内的比例越大说明误差分布集中。等级小误差概率 P后验差比值 C模型状态一级P 0.95C 0.35可直接外推二级0.80 P ≤ 0.950.35 ≤ C 0.50适合短期预测三级0.70 P ≤ 0.800.50 ≤ C 0.65只做趋势参考四级P ≤ 0.70C ≥ 0.65参数不可用需重构function [C, P, level] mgm_check(X0, X0p) % 后验差检验 % X0 m×n 原始序列 % X0p m×n 模型拟合值与 X0 同长度 E X0 - X0p; % 残差矩阵 S1 std(X0, 0, 1); % 原始序列标准差1×n S2 std(E, 0, 1); % 残差标准差1×n C S2 ./ S1; % 后验差比值1×n thr 0.6745 * S1; % 小误差阈值 P mean(abs(E - mean(E,1)) thr, 1); % 小误差概率1×n level strings(1, size(X0,2)); for j 1:size(X0,2) if P(j) 0.95 C(j) 0.35; level(j) 一级; elseif P(j) 0.80 C(j) 0.50; level(j) 二级; elseif P(j) 0.70 C(j) 0.65; level(j) 三级; else; level(j) 四级; end end end逻辑说明std(X, 0, 1)的第二个参数 0 表示按样本标准差除以 N-1第三个参数 1 表示沿行方向统计也就是沿时间维。两个 1 都不写会默认沿第一维统计当 n1 时结果就是错的。mean(abs(...) thr, 1)用逻辑矩阵求均值得到的就是满足条件的时刻占比。参数说明X0和X0p必须包含相同的时刻范围。做样本内检验时用全部拟合值做样本外检验时X0 要裁到对应的外推区间否则 S1 会被训练期的波动稀释C 值虚高。注意C 和 P 是逐变量给出的只要有一个变量掉到四级整个模型的联合外推就不可信因为 A 的非对角元会让误差在各变量间传导。4.3 量纲差异、采样间隔与负值数据的排错实测中四个坑反复出现量纲差异过大。气温用摄氏度、GDP 用亿元紧邻均值里大数会把小数淹掉A 直接病态。做法是先zscore标准化再建模预测完反标准化。这在多变量场景下几乎是必做步骤。非等距采样。灰色模型隐含等时间距假设跨月缺测、跨年数据直接用会让 Z1 失真。先用interp1插成等距序列插值后要复核准光滑比重跑一遍。负值或零值。需求下滑、利润为负时累加序列不再单调级比检验直接失败。加常数 c 平移c 一般取1.1*max(abs(X0))预测完再减回去。变量数逼近样本数。m-1 与 n1 只差一两个时即便能解出来参数的方差也极大换一期数据 A 就翻脸。经验做法是变量数控制在(m-1)/3以内。对应的排查代码% 标准化后再建模避免量纲差异导致 A 病态 mu mean(X0, 1); sg std(X0, 0, 1); Xs (X0 - mu) ./ sg; % 逐列标准化 [As, Bs, X1s] mgm_fit(Xs 1); % 平移到正数域后再拟合 X1p mgm_predict(As, Bs, X1s(1,:), 5); X0p [X1p(1,:); diff(X1p)]; X0p (X0p - 1) .* sg mu; % 反标准化回原量纲逻辑说明加 1 平移是为了让标准化后的数据落进正数域减 1 时注意只在还原阶段减一次。sg和mu必须用训练期数据算出来不能混入外推期的统计量否则就是样本泄露。5. 背景值优化、滚动预测与 MATLAB 优化工具箱配合默认紧邻均值用 λ0.5也就是 Z1 0.5·X1(k1) 0.5·X1(k)。这个取值是几何中点但在序列增长较快时会引入系统性偏差。把 λ 当参数、以样本外 MAPE 为目标函数用fminsearch或粒子群算法搜索通常能把 C 值压下去一档。% 以样本外 MAPE 为目标搜索背景值权重 lambda obj (lam) mgm_mape(X0, lam); lam_opt fminsearch(obj, 0.5, ... optimset(TolX, 1e-4, MaxFunEvals, 200)); fprintf(最优背景值权重 lambda %.4f\n, lam_opt); function m mgm_mape(X0, lam) [~,~,~,A,B,X1] mgm_fit_lambda(X0, lam); % 内部使用 lam 构造 Z1 X1p mgm_predict(A, B, X1(1,:), size(X0,1)-1); X0p [X1p(1,:); diff(X1p)]; m mean(abs((X0(2:end,:) - X0p(2:end,:)) ./ X0(2:end,:)), all) * 100; endfminsearch属于 MATLAB 优化工具箱的基础函数无约束、无梯度对小参数问题够用变量多、目标函数非凸时换particleswarm但那需要全局优化工具箱。目标函数里算的是样本外一步误差注意别把训练区间一起丢进去算 MAPE否则 λ 会过拟合到训练期。另一个工程化习惯是滚动预测每来一期真实数据就整体重拟合一次 A、B而不是拿一个 A 硬撑十二期。多变量灰色预测模型的 A 是耦合矩阵任一变量出新拐点时非对角元会跟着变滚动重拟合是唯一稳健的做法。最后一段可复现的验证方式把最后两期留作样本外用前 m-2 期拟合、向后外推两步比对 MAPE 是否落在训练期 MAPE 的 1.5 倍以内。超出这个比例说明 A 的非对角元已经不稳定先回头检查标准化和变量数是否过密。本文还有配套的精品资源点击获取
