用于数据量较少的模型
模型引入什么是灰色一个系统通常分为黑，白，灰。白色系统，即数据完全清楚已知的系统。黑色系统完全未知，而灰色系统就是介于二者之间，仅有部分是已知的。
时间序列灰色模型通常用于时间序列的预测，但是灰色预测更加适合已知的历史数据较少的情况
核心过程灰色预测模型中，GM(1,1)（Gray Model (1,1)） 是最基础、最常用的单变量一阶灰色模型，适用于“单指标、短期、小样本”预测。括号中“1”分别表示“1个变量”和“1阶微分方程”。
模型的核心假设原始数据序列需满足“光滑性”和“指数增长趋势”（可通过级比检验验证），否则需进行数据预处理（如平移变换、对数变换）。
数学推导设原始非负数据序列为：$$ X^{(0)} = {x^{(0)}(1), x^{(0)}(2), …, x^{(0)}(n)}, \quad n \geq 4 \quad (n为样本数) $$
累加生成（AGO）对原始序列进行一次累加生成（1-AGO），得到生成序列 $ X^{(1)} $，目的是增强数据的规律性（让杂乱数据呈现指数增长趋势）：
$$ X^{(1)} = {x^{(1)}_1, x^{(1)}_2, …, x^{(1)}_n} $$
其中，$ x^{(1)}(k) = \sum_{i=1}^k x^{(0)}(i) \quad (k=1,2,…,n) $ (前缀和)
构建白化微分方程假设 $ X^{(1)} $ 满足一阶线性微分方程（白化方程）：$$ \frac{dx^{(1)}(t)}{dt} + a x^{(1)}(t) = b $$
$ a $：发展系数（反映序列增长趋势的强弱，$ |a|$ 越小，趋势越平稳）；
$ b $：灰作用量（反映系统的外部输入或固有特性）。
最小二乘估计参数将微分方程离散化，得到矩阵形式 $ Y = B\hat{\alpha} $，其中：
观测向量 $ Y = \begin{bmatrix} x^{(0)}(2) \ x^{(0)}(3) \ … \ x^{(0)}(n) \end{bmatrix} $
数据矩阵 $ B = \begin{bmatrix} -\frac{1}{2}[x^{(1)}(1)+x^{(1)}(2)] & 1 \ -\frac{1}{2}[x^{(1)}(2)+x^{(1)}(3)] & 1 \ … & … \ -\frac{1}{2}[x^{(1)}(n-1)+x^{(1)}(n)] & 1 \end{bmatrix} $
参数向量 $ \hat{\alpha} = \begin{bmatrix} a \ b \end{bmatrix} $
通过最小二乘法求解参数：$$ \hat{\alpha} = (B^T B)^{-1} B^T Y $$
求解白化方程与还原预测白化方程的解析解（时间响应函数）为：$$ x^{(1)}(t) = \left( x^{(0)}(1) - \frac{b}{a} \right) e^{-a(t-1)} + \frac{b}{a} $$
对 $ x^{(1)}(t) $ 进行累减生成（差分）,（IG），还原得到原始序列的预测值 $ \hat{x}^{(0)}(k) $：$$ \hat{x}^{(0)}(1) = x^{(0)}(1) \quad $$$$ \hat{x}^{(0)}(k) = \hat{x}^{(1)}(k) - \hat{x}^{(1)}(k-1) \quad (k \geq 2) $$
精度检验通过以下指标验证模型可靠性：
平均相对误差（MAPE）：$ \text{MAPE} = \frac{1}{n} \sum_{k=1}^n \left| \frac{x^{(0)}(k) - \hat{x}^{(0)}(k)}{x^{(0)}(k)} \right| \times 100% $
一级精度：$ \text{MAPE} < 10% $；二级精度：$ 10% \leq \text{MAPE} < 20% $；三级精度：$ 20% \leq \text{MAPE} < 30% $
后验差检验：
原始序列标准差 $ S_1 = \sqrt{\frac{1}{n} \sum_{k=1}^n (x^{(0)}(k) - \bar{x}^{(0)})^2} $
残差序列 $ \varepsilon(k) = x^{(0)}(k) - \hat{x}^{(0)}(k) $，残差标准差 $ S_2 = \sqrt{\frac{1}{n} \sum_{k=1}^n (\varepsilon(k) - \bar{\varepsilon})^2} $
后验差比 $ C = \frac{S_2}{S_1} $，小误差概率 $ P = P(|\varepsilon(k) - \bar{\varepsilon}| < 0.6745 S_1) $
一级精度：$ C < 0.35, P > 0.95 $；二级精度：$ C < 0.5, P > 0.8 $；三级精度：$ C < 0.65, P > 0.7 $
完整步骤数据预处理与检验
数据要求：原始序列 $ X^{(0)} $ 为非负序列（若有负数，进行平移变换 $ x’(0)(k) = x^{(0)}(k) + c $，$ c $ 为常数使所有数据非负）。
级比检验：验证数据是否适合GM(1,1)模型，计算级比 $ \lambda(k) = \frac{x^{(0)}(k-1)}{x^{(0)}(k)} \quad (k=2,…,n) $，若所有 $ \lambda(k) \in (e^{-2/(n+1)}, e^{2/(n+1)}) $，则适合建模；否则需进行数据变换（如对数变换、平移变换）。
累加生成（1-AGO）计算 $ X^{(1)} $，增强数据规律性。
构建矩阵 $ B $ 和 $ Y $按公式构造矩阵，注意维度匹配（$ B $ 为 $ (n-1) \times 2 $ 矩阵，$ Y $ 为 $ (n-1) \times 1 $ 向量）。
参数估计（求解 $ a, b $）通过最小二乘法计算 $ \hat{\alpha} = (B^T B)^{-1} B^T Y $。
预测与还原
短期预测：计算 $ \hat{x}^{(1)}(n+1), \hat{x}^{(1)}(n+2), … $（竞赛中预测步长建议≤3，步长越长误差越大）。
累减还原：得到原始序列的预测值 $ \hat{x}^{(0)}(n+1), … $。
精度检验计算MAPE、后验差比 $ C $、小误差概率 $ P $，若精度不达标，需优化模型（如残差修正、组合模型）。
代码实现：Python版本12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879import numpy as npimport matplotlib.pyplot as pltdef gm11(x0, predict_num=1): """ GM(1,1)灰色预测模型 :param x0: 原始序列（list或array，非负） :param predict_num: 预测步数（默认1步，竞赛建议≤3） :return: 原始序列预测值、未来预测值、模型参数(a,b)、平均相对误差(mape) """ x0 = np.array(x0).reshape(-1, 1) n = len(x0) # 1. 累加生成（1-AGO） x1 = np.cumsum(x0, axis=0) # 2. 构建矩阵B和Y B = np.zeros((n-1, 2)) Y = np.zeros((n-1, 1)) for i in range(n-1): B[i, 0] = -0.5 * (x1[i] + x1[i+1]) B[i, 1] = 1 Y[i] = x0[i+1] # 3. 最小二乘法估计参数a, b alpha = np.dot(np.dot(np.linalg.inv(np.dot(B.T, B)), B.T), Y) a, b = alpha[0, 0], alpha[1, 0] # 4. 时间响应函数与预测（x1_hat） def x1_hat(k): return (x0[0, 0] - b/a) * np.exp(-a*(k-1)) + b/a # 5. 还原预测值（x0_hat） x0_hat = np.zeros((n, 1)) x0_hat[0] = x0[0, 0] # 初始值不变 for k in range(1, n): x0_hat[k] = x1_hat(k) - x1_hat(k-1) # 6. 未来预测 future_x0 = [] for k in range(n+1, n+1+predict_num): future_x0.append(x1_hat(k) - x1_hat(k-1)) future_x0 = np.array(future_x0).reshape(-1, 1) # 7. 精度检验（平均相对误差MAPE） mape = np.mean(np.abs((x0 - x0_hat)/x0)) * 100 return x0_hat, future_x0, (a, b), mapeif __name__ == "__main__": # 原始数据：某地区2018-2022年GDP（单位：亿元） x0 = [500, 550, 610, 680, 760] predict_num = 2 # 预测2023、2024年GDP # 建模预测 x0_hat, future_x0, (a, b), mape = gm11(x0, predict_num) # 结果输出 print(f"模型参数：a={a:.4f}, b={b:.4f}") print(f"平均相对误差MAPE：{mape:.2f}%") print(f"原始序列预测值：{x0_hat.flatten().round(2)}") print(f"2023年预测GDP：{future_x0[0,0]:.2f}亿元") print(f"2024年预测GDP：{future_x0[1,0]:.2f}亿元") plt.rcParams['font.sans-serif'] = ['SimHei'] # 支持中文 plt.figure(figsize=(10, 6)) years = np.arange(2018, 2018+len(x0)) predict_years = np.arange(2018+len(x0), 2018+len(x0)+predict_num) plt.plot(years, x0, 'o-', label='原始GDP', linewidth=2, markersize=8) plt.plot(years, x0_hat.flatten(), 's-', label='拟合GDP', linewidth=2, markersize=8) plt.plot(predict_years, future_x0.flatten(), '^-', label='预测GDP', linewidth=2, markersize=8, color='red') plt.xlabel('年份', fontsize=12) plt.ylabel('GDP（亿元）', fontsize=12) plt.title(f'GM(1,1)灰色预测模型：某地区GDP预测（MAPE={mape:.2f}%）', fontsize=14) plt.legend(fontsize=10) plt.grid(True, alpha=0.3) plt.show()
Matlab版本123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172function [x0_hat, future_x0, a, b, mape] = gm11(x0, predict_num) % GM(1,1)灰色预测模型 % 输入：x0-原始序列（行向量），predict_num-预测步数 % 输出：x0_hat-原始序列预测值，future_x0-未来预测值，a-发展系数，b-灰作用量，mape-平均相对误差 x0 = x0(:); % 转为列向量 n = length(x0); % 1. 累加生成（1-AGO） x1 = cumsum(x0); % 2. 构建矩阵B和Y B = zeros(n-1, 2); Y = zeros(n-1, 1); for i = 1:n-1 B(i, 1) = -0.5 * (x1(i) + x1(i+1)); B(i, 2) = 1; Y(i) = x0(i+1); end % 3. 最小二乘法估计参数 alpha = (B' * B) \ B' * Y; a = alpha(1); b = alpha(2); % 4. 时间响应函数与还原预测 x0_hat = zeros(n, 1); x0_hat(1) = x0(1); for k = 2:n x1_hat_k = (x0(1) - b/a) * exp(-a*(k-1)) + b/a; x1_hat_k_1 = (x0(1) - b/a) * exp(-a*(k-2)) + b/a; x0_hat(k) = x1_hat_k - x1_hat_k_1; end % 5. 未来预测 future_x0 = zeros(predict_num, 1); for k = 1:predict_num t = n + k; x1_hat_t = (x0(1) - b/a) * exp(-a*(t-1)) + b/a; x1_hat_t_1 = (x0(1) - b/a) * exp(-a*(t-2)) + b/a; future_x0(k) = x1_hat_t - x1_hat_t_1; end % 6. 计算MAPE mape = mean(abs((x0 - x0_hat) ./ x0)) * 100;end% ---------------------- 测试代码 ----------------------x0 = [500, 550, 610, 680, 760];predict_num = 2;[x0_hat, future_x0, a, b, mape] = gm11(x0, predict_num);fprintf('模型参数：a=%.4f, b=%.4f\n', a, b);fprintf('平均相对误差MAPE：%.2f%%\n', mape);fprintf('原始序列预测值：');disp(x0_hat');fprintf('2023年预测GDP：%.2f亿元\n', future_x0(1));fprintf('2024年预测GDP：%.2f亿元\n', future_x0(2));% 可视化years = 2018:2022;predict_years = 2023:2024;plot(years, x0, 'o-', 'LineWidth', 2, 'MarkerSize', 8, 'DisplayName', '原始GDP');hold on;plot(years, x0_hat', 's-', 'LineWidth', 2, 'MarkerSize', 8, 'DisplayName', '拟合GDP');plot(predict_years, future_x0', '^-', 'LineWidth', 2, 'MarkerSize', 8, 'Color', 'red', 'DisplayName', '预测GDP');xlabel('年份', 'FontSize', 12);ylabel('GDP（亿元）', 'FontSize', 12);title(sprintf('GM(1,1)灰色预测模型：某地区GDP预测（MAPE=%.2f%%）', mape), 'FontSize', 14);legend('FontSize', 10);grid on;hold off;
总结与应用场景灰色预测的核心优势是“小样本、贫信息”，在数学建模中常用于以下场景：
短期经济预测（如GDP、居民收入、物价指数）；
公共卫生预测（如疫情初期确诊人数、疫苗接种需求）；
工程预测（如设备故障时间、原材料消耗、能耗预测）；
环境预测（如空气质量指数、污染物浓度短期预测）。
在竞赛中，若遇到“数据少、无明确分布、短期预测”需求，可优先考虑GM(1,1)模型，并通过“残差修正+组合模型”提升精度。