小样本预测利器：灰色预测GM(1,1)模型详解与Python实现

简介这是一个面向数据分析与机器学习学习者的Python灰色预测模型示例代码包围绕GM(1,1)一阶单变量模型演示完整流程涵盖数据预处理、一次累加生成、参数求解、模型预测和误差分析。适合需要快速上手灰色预测、处理小样本或不完全数据的读者。压缩包共六个文件包含五个Python脚本和一个txt测试数据文件整体约4KB脚本对应不同实验用例数据文件用于模型验证便于对照运行和二次修改。已有3907人学习下载关注度较高。通过研读这套代码可以掌握使用numpy、pandas等库实现灰色模型的具体写法理解最小二乘参数估计与累加逆变换等关键环节并可将该模型迁移到自身的时间序列预测场景。代码按功能拆分结构清晰特别适合数据量少、信息不完全的场景能够弥补传统统计方法对数据长度要求较高的不足可作为课程设计、期末作业或毕业设计前的入门参考。1. 小样本预测的灰色预测为什么我劝你先别看神经网络当你手里只有十几条历史数据却要预测未来三四期的走势时机器学习模型常常陷入两难训练集太小模型不是欠拟合就是过拟合。我做过不少次这种场景比如水位测站的月度观测、设备告警次数统计最后发现最管用的反而是“土办法”——灰色预测模型。它用一条原始数据序列就能建模GM(1,1) 靠一次累加生成加最小二乘求解不依赖大样本也不需要调超参数。这套 Python 源码包里共有五份脚本文件Pex15_1.py、Pex15_2.py、Pex15_3.py、Pex15_4_1.py、Pex15_4_2.py和一个测试数据文件 Pdata15_3.txt核心实现全部基于 numpy 与 pandas把 GM(1,1) 从建模到检验的完整链路都演示了一遍。适合数学建模参赛者、在小样本数据上做短期趋势判断的数据从业者以及想真正看懂灰色预测内部计算过程而非只会调库的 Python 学习者。2. GM(1,1) 建模与代码解读从累加生成到最小二乘参数求解2.1 为什么累加生成这一步是灰色预测的“魔法”灰色预测模型建立在对观测序列的再加工上。原始数据往往上下波动直接拿去做拟合容易被随机噪声带偏。GM(1,1) 的做法是先做一次累加生成把原始序列 x(0) 变成 x(1)公式是 x(1)(k) x(0)(1) x(0)(2) … x(0)(k)。累加的过程是一个逐步累积的过程会把原本起伏不定的数据压成一条基本单调的曲线。这条平滑曲线更接近指数型增长曲线也是后续建立一阶微分方程的数学前提。这一步不是可有可无的前处理而是灰色预测算法的本质特征用累加弱化随机性用累减还原预测值整个模型都围绕这一正一反两个变换构造。累加完成后要估计的是白化微分方程 dx(1)/dt a x(1) b 中的参数 a 和 b。这里的 a 叫发展系数代表序列的内在增长或衰减速度b 叫灰作用量相当于外部驱动力的拟合结果。离散化后方程变成与紧邻均值序列 z(1) 相关的线性关系x(0)(k) -a z(1)(k) b。当 n 个观测值组成 n-1 个方程时方程组是超定的所以要用最小二乘法求解。这里容易混的是 numpy 里的函数核心不是np.polyfit而是np.linalg.lstsq它会根据最小化误差平方和的原则一次性给出 a 和 b这也是 Pex15_1.py 里最核心的计算环节。2.2 Pex15_1.py 与 Pdata15_3.txt一个能跑通的灰色预测最小实现Pdata15_3.txt 里存的是示例数据格式通常是每行一个数值或者数字之间用空白分隔。读取方法我一般直接用 numpy 的 loadtxt因为这类数据文件很少带表头列数也统一。下面这段代码基本就是 Pex15_1.py 的核心逻辑我把它整理成一个可直接套用的函数。import numpy as np import pandas as pd # 读取数据Pdata15_3.txt 为空白分隔的文本数据无表头 x0 np.loadtxt(Pdata15_3.txt).reshape(-1) print(原始序列长度 n , len(x0)) print(前 5 条观测, x0[:5]) def gm11(x0, n_pred5): # x0: 一维原始数据序列, n_pred: 向后预测的期数 x0 np.asarray(x0, dtypenp.float64) n len(x0) # 一次累加生成序列 x1 x1 np.cumsum(x0) # 紧邻均值生成序列 z1长度比 x1 少 1 z1 (x1[:-1] x1[1:]) / 2.0 # B 阵和 Y 阵是 GM(1,1) 标准形式x0[k] -a*z1[k] b B np.column_stack((-z1, np.ones(n - 1))) Y x0[1:].reshape(-1, 1) # 最小二乘求解参数 a、b u np.linalg.lstsq(B, Y, rcondNone)[0].ravel() a, b u[0], u[1] # 时间响应式x1_hat(k1) (x0[0] - b/a) * exp(-a*k) b/a k np.arange(n n_pred) x1_hat (x0[0] - b / a) * np.exp(-a * k) b / a # 累减还原得到预测值 x0_hat np.r_[x1_hat[0], np.diff(x1_hat)] return x0_hat, a, b x0_hat, a, b gm11(x0, n_pred5) print(预测值含历史回代, x0_hat) print(发展系数 a , a, 灰作用量 b , b)这段代码有四个值得记住的细节。第一np.cumsum完成累加生成千万别自己写循环numpy 版本跑几十万数据也不卡。第二z1用的是“紧邻均值生成”也就是相邻两个累加值的平均这是灰色预测与普通指数拟合的关键区别直接拿 x1 本身去构建 B 阵参数估计会产生系统性偏差。第三x1_hat表达式里(x0[0] - b / a)的符号直接决定预测曲线方向如果原始序列首值相对 b/a 的差距过大外推值可能出现发散。第四np.diff(x1_hat)是累减还原它把累加空间的预测值逐项取差得到原始量纲下的预测结果长度刚好和 x1_hat 对齐。n_pred 是向后预测期数灰色预测适合短中期预测我自己在实际项目里一般控制在 3 到 5 期设太长相当于把指数推到看不见边界的方向。2.3 拆开压缩包看构成五个脚本分别示范了什么这个压缩包里脚本命名带明显的教材痕迹Pex 前缀加上章节号大概率是某本数学建模教材第 15 章的配套 Python 例子。按文件名称和灰色预测的知识结构对应关系如下表所示。文件大概率示范内容学习价值Pex15_1.pyGM(1,1) 基础建模与预测最小实现适合入门Pex15_2.py级比检验与模型适用性判断解决“能不能用灰色预测”的选题问题Pex15_3.py残差修正模型提高预测精度的常用扩展Pex15_4_1.py灰色灾变预测对超阈值时点建模类似异常预警Pex15_4_2.py后验差检验或滚动预测评估模型质量处理数据更新场景我建议的阅读顺序是从 Pex15_1.py 起步跑通后再把 Pex15_2.py 和 Pex15_3.py 对照着看。Pex15_2.py 里的级比检验是灰色预测里最容易被跳过的环节直接影响模型能不能用后面第 3 章会专门展开。Pex15_3.py 走的是残差修正路线实际场景里命中率很高先用基础 GM(1,1) 预测一次算出真实值与预测值的残差再对残差序列单独建一个 GM(1,1)两部分预测累加后作为最终输出。这种方法本质上是在误差序列里再挖一次灰色规律实战中常常能挽回几个百分点的精度。3. 把 GM(1,1) 代码落到自己的数据Pdata15_3.txt 读取、级比检验与步长设置3.1 数据预处理拿到 Pdata15_3.txt 之后第一步先做体检许多初学灰色预测的人拿到数据文件后就急着跑脚本结果输出全是 nan。Pdata15_3.txt 这类文件本身内容简单但换成自己的数据时格式问题就来了。我的固定做法是先用 pandas 读一遍、.describe()看分布再决定下一步。import pandas as pd import numpy as np # 先读成 DataFramesep 用正则匹配任意空白header 用 None df pd.read_csv(Pdata15_3.txt, sepr\s, headerNone) desc df.describe().T print(desc) # 检查缺失与无穷值 print(NaN 数量, df.isna().sum().sum()) print(无穷值数量, np.isinf(df.values).sum())这里有两个操作值得解释。sepr\s是把 Tab、空格、多个空白全部当作列分隔符比固定sep 健壮得多headerNone避免把第一行数据当成列名。灰色预测对数据有隐性要求原始序列最好为正数。数据里有零或负数时累加生成序列的均值 z1 可能跨过零点参数估计结果变得不稳定预测值甚至出现负值。常见做法是先做整体平移给整个序列加一个常数让最小值大于等于 1预测完再平移回来。平移后序列的水平和波动幅度都变了级比检验区间也会跟着变这一步在 Pex15_3.py 这样的残差修正例子里尤其重要。3.2 预测步长和参数输出a、b 到底在说什么gm11 函数跑完除了预测值还会返回 a 和 b。a 的绝对值越大说明序列自身的增长或衰减趋势越强a 为负时累减还原后的预测值通常在增长a 为正时预测趋势朝水平线收敛。b 是灰作用量相当于模型中的驱动项它的数值大小依赖原始数据的绝对量级单独看没有意义要和 a 一起配合使用。记住一个边界|a| 若大于 0.3模型外推的稳定性会明显下降。此时与其硬跑 GM(1,1)不如先对数据做对数变换或平移把趋势放缓再建模。预测期数 n_pred 是最容易拍脑袋决定的一个参数。灰色预测的本质是指数外推预测期拉长误差会成指数级放大。我的一般建议是训练数据有 n 个点最多向后预测 n/3 期如果原始序列呈现明显非线性波动最好只预测 1 到 2 期。下面的表格是我在项目里的参考值。原始序列长度建议预测期数理由5 ~ 81 ~ 2数据太少模型信息量有限9 ~ 152 ~ 4级比检验通过时可适当外推16 ~ 303 ~ 5样本充足滚动预测更合适需要强调表格里的“期数”不是季度或年而是与数据频率对应的自然单位。预测完以后至少留下最后几个观测点作为验证集别把全部数据都喂进模型否则后验差检验再好看也说明不了泛化能力。3.3 级比检验跑 Pex15_2.py 前先学会自己算级比检验是判断原始数据适不适合用 GM(1,1) 的“体检报告”。级比定义为当前观测除以前一观测公式是 λ(k) x(0)(k-1) / x(0)(k)。注意分子是更早的那一项很多资料这里符号写法不统一会导致检验结论恰好相反。合格区间是 (e^(-2/(n1)), e^(2/(n1)))n 是序列长度。所有级比都落在区间内说明数据指数趋势明显适合灰色预测。核心判断代码很短。def level_ratio_check(x0): n len(x0) # 级比前一项除以当前项 lam x0[:-1] / x0[1:] lb np.exp(-2 / (n 1)) ub np.exp(2 / (n 1)) all_pass np.all((lam lb) (lam ub)) return all_pass, lam, lb, ub ok, lam, lb, ub level_ratio_check(np.loadtxt(Pdata15_3.txt)) print(级比是否全部落在区间内, ok) print(区间边界, lb, ub)如果你的数据没有通过级比检验不用急着放弃灰色预测。先检查是不是存在极端异常值把异常点去除后再算一次如果还是不通过就对原始序列做一次方根变换或对数变换。做变换后记住最终预测结果要还原到原始量纲——对数变换就取指数平移变换就减回去。很多人在 Pex15_2.py 里看到“未通过”提示就直接换模型其实忽略了变换这一步错过了一次修正的机会。提示级比检验区间的上下界都依赖样本量 n改动数据长度后区间会变别拿固定常数当标准。4. 避坑/常见问题/排查GM(1,1) 复现时最容易翻车的 5 个细节4.1 复现数值对不上先查累加、均值与符号方向问题 1预测结果出现负值而且趋势明显发散。现象使用 Pdata15_3.txt 或其他自己采集的数据跑 Pex15_1.py预测的第三第四期直接变负输出看起来完全不可用。原因多数情况下原始序列里包含非正数或者数据量级跨度超过 10 倍。累加序列在这种情况下会出现 z1 穿越零点的情形最小二乘求出的 a 失去稳定性时间响应式的外推值被指数放大于是预测值翻车。少数情况是读取数据时把文件里的非数值字符读进了数组导致后续计算直接异常。解决先看原始序列最小值如果小于等于 0就给整个序列加一个正常数平移平移量取绝对值最小值的 1.5 倍如果数据跨度太大先取对数。预测完成后记得反向还原否则结果完全不可用。问题 2RuntimeWarning 不断刷屏预测值全是 nan。现象numpy 在运行 gm11 时不断弹出一串 “invalid value encountered in exp”print 出来的预测结果全是 nan 或 inf。原因np.linalg.lstsq求出的 a 为 0 或极其接近 0时间响应式里 b/a 直接除零。这种事最容易发生在原始序列为常数或近常数序列时比如整列数据都是同一个值。另外Pdata15_3.txt 里如果有空行loadtxt 可能把后一行的数值错位读进来同样造成运算异常。解决读入数值后立刻做一次合法性检查np.isfinite(x0).all()不通过就停止建模。对近常数序列先做差分变换或加微小噪声别拿原始序列硬刚。a 接近 0 时可以给 b/a 的除法加一个极小 epsilon 保护但更稳妥的思路是承认该序列不适用灰色预测。4.2 输出结果不合理级比方向、后验差指标和还原公式是重灾区问题 3同样的数据代码包里明明通过了级比检验自己一算却显示不通过。现象直接用列表推导计算级比结果与 Pex15_2.py 打印的结论相反模型适用性判断完全反过来。原因级比有“前除以后”和“后除以前”两种写法区间边界也对应不同方向。写成 x0[1:] / x0[:-1] 的话若数据呈整体下降趋势结论会从“通过”变成“不通过”。解决固定使用x0[:-1] / x0[1:]这种前一项除以后一项的写法并把区间边界与代码里的公式一一对照。同时确保读入序列顺序与时间方向一致不要因为排序操作把时间序列倒置。问题 4后验差比值算出来小得离谱或者根本算不出来。现象调用后验差检验时C 值接近 0或者代码报分母为零的错。原因后验差比 C S2 / S1S1 是原始序列标准差S2 是残差序列标准差。当数据本身几乎不波动S1 趋近 0比值就会变成畸形数值当模型回代误差极小S2 趋近 0则 C 显得“异常优秀”其实是用同一份数据既拟合又验证造成的假象。解决用独立的验证集计算误差不要把回代残差当作泛化误差。S1 过小时改用平均绝对百分比误差 MAPE 作为评价指标或者对数据做差分后再建模。后验差检验只是参考不能作为一个独立指标决定模型去留。问题 5预测曲线整体平移和真实数据差一个固定常数。现象回代预测的走势和真实值一致但整体错开一个固定量首期就不对。原因累减还原公式写错。常见错误是用循环从 k0 开始计算 x0_hat(k) x1_hat(k) - x1_hat(k-1)第一个还原值被算成了 x1_hat(0) 减一个不存在的前项而不是直接取 x1_hat(0)。解决用np.r_[x1_hat[0], np.diff(x1_hat)]完成还原或者循环从 i1 开始做差分。这个现象常发生在自己重写代码时也是“精度高但结果错”的典型警醒。以上五条是我在复现灰色预测代码时踩过的典型坑基本都对应代码包里至少一个脚本。你会发现排错的核心往往不是改模型而是确认数据读写、级比方向、参数计算和评估口径这四件事。5. 实际场景调参白化响应、滚动预测与 ARIMA 残差组合5.1 白化响应表达式里的决策点a 值大小决定预测能走多远灰色预测的参数估计只是开始真正会用起来的关键在白化响应表达式也就是 x1_hat 的通项公式。公式里 (x0[0] - b/a) 是初始偏移量exp(-a·k) 控制整个预测曲线的形状。a 对预测的影响可以分几个区间看。a 的取值范围数据形态实操建议a -0.5强增长且增速加快不宜做长期外推-0.5 a 0稳定增长型可做中短期预测a ≈ 0数据接近随机游走或平稳灰色预测意义有限0 a 0.5衰减趋稳型可预测收敛水平a 0.5快速衰减型先查是否过度处理数据实际建模时我习惯先把参数打出来看一眼再决定要不要加修正。a 的绝对值太大说明数据本身波动剧烈预测步长建议压缩到最小。碰到这种情况Pex15_3.py 的残差修正思路就派上用场第一次建模得到的误差序列里往往还有灰色规律对残差再建一个 GM(1,1)把两部分预测叠加。如果残差序列的趋势不再明显就不需要修正硬加反而会引入第二重过拟合。5.2 滚动 GM(1,1)Pex15_4_2.py 风格的新陈代谢模型很多实际业务预测进入运营阶段后数据会不断更新。这时最忌讳用固定模型一路外推因为老数据的信息会慢慢失效。Pex15_4_2.py 这个文件名后缀暗示它比 Pex15_4_1.py 多了一层扩展处理实际中常见形态就是“新陈代谢 GM(1,1)”。思路很简洁每预测完一期把真实观测并入窗口丢掉最早的一个观测保持窗口长度不变再重新拟合 GM(1,1)。这样能限制模型对陈旧信息的依赖预测精度通常比一次性外推要好。def rolling_gm11(x0, window6, n_pred3): # window: 每次建模使用的观测窗口长度, n_pred: 滚动预测期数 x0 np.asarray(x0, dtypenp.float64) result [] for i in range(n_pred): # 截取最近 window 个观测重新调用基础 GM(1,1) xw x0[-window:] xhat, a, b gm11(xw, n_pred1) x_ xhat[-1] result.append(x_) # 把预测值加入序列尾部下一次循环再滚动一次 x0 np.r_[x0, x_] return np.array(result)这段代码的逻辑核心是加一截窗口再滚动每次循环都重新调用 gm11而不是沿用第一次算出的 a、b。灰色的“适应性”正是体现在这里参数会随窗口移动而变化。window 取值过小模型容易抖动window 取值过大滚动就失去了意义。经验参考值是取原始序列长度的一半最低别少于 4 个观测。滚动预测在数据更新过程中需要人工重复触发更好的做法是把它包成一个定时任务预测新值后立刻追加到数据库形成自动更新闭环。5.3 与 ARIMA 组合灰色抓趋势、残差抓结构灰色预测单独用在强波动数据上会吃亏所以很多人想到了组合预测。灰色预测和 ARIMA 结合是常见扩展方向。我一般这样搭先用 GM(1,1) 预测趋势部分计算出历史回代误差序列再用误差序列拟合 ARIMA 模型。ARIMA 的阶数交给 AIC 自动选择但样本量小时别用高阶参数。组合结果等于灰色预测值加上 ARIMA 残差预测值整体公式就是最终预测 G_t R_t。from statsmodels.tsa.arima.model import ARIMA def hybrid_forecast(x0, n_pred3, arima_order(1, 0, 1)): # 第一步基础灰色预测 trend, a, b gm11(x0, n_predn_pred) n len(x0) # 第二步回代残差 resid x0[:n] - trend[:n] # 第三步对残差序列拟合 ARIMA model ARIMA(resid, orderarima_order).fit() resid_pred model.forecast(n_pred) # 第四步组合输出 return trend[n:n n_pred] resid_pred需要提醒的是残差序列本身可能已经是白噪声此时 ARIMA 拟合出来的外推值基本是零组合模型精度不会比纯灰色预测好太多。判断依据很简单对 resid 做 Ljung-Box 检验p 值大于 0.05 说明无明显自相关直接用灰色预测即可不用叠第二层。组合预测在提高稳定性的同时也把误差来源从一处变成两处任何一环节参数设错整体效果都会打折。所以我在项目里先跑单模型确认趋势部分可靠后再决定要不要加这一层。6. 预测完不检验等于白做后验差检验与最小验证三连灰色预测代码跑起来很容易难的是判断结果能不能交付。我从第一次复现 Pex15_4_2.py 时就养成了一个习惯任何灰色预测模型必须过三道验证才敢上报数据。第一道是级比检验确认数据具备灰色预测的前提条件第二道是后验差检验用 C 值衡量预测误差相对数据波动的大小第三道是曲线端点检查观察预测序列首值和历史数据末值是否平滑衔接。三道都过模型才算真正可用。后验差检验的计算极短核心就三行。resid x0[:n] - x0_hat[:n] S1 np.std(x0, ddof1) S2 np.std(resid, ddof1) C S2 / S1 print(后验差比值 C , C)C 越小说明预测误差相对原始波动越小。教材里的分级标准是 C 0.35 为优小于 0.5 算合格大于 0.65 不合格。当 C 值落在不合格区间时我通常回到 Pex15_2.py 的级比检验和 Pex15_3.py 的残差修正两条路去排查而不是直接改代码。端点检查则是个土办法把预测曲线和历史曲线画在同一坐标系里如果预测第一点和历史最后一点之间出现断崖式跳动基本可以断定累减还原或者平移还原出了问题检查一下 np.diff 用没用对、有没有忘减平移量往往一抓一个准。从那以后我每次拿到灰色预测相关代码都会强制把级比检验、后验差检验和端点连续性这三件事走一遍。数据不干净就先平移或变换预测步长压到样本量的三分之一以内实在不行再上残差修正和滚动重建。这套流程看着笨却让我躲过了不少预测翻车的尴尬。希望上面的拆解对你有用这套灰色预测 Python 源码包值得每个做小样本预测的人手跑一遍。本文还有配套的精品资源点击获取