【数模/预测】灰色预测

声明：文章参考数学建模清风的网课编写。

简介

 灰色预测是对既含有已知信息又含有不确定信息的系统进行预测，就是对在一定范围内变化的、与时间有关的灰色过程进行预测。

 灰色预测对原始数据进行生成处理来寻找系统变动的规律，并生成有较强规律性的数据序列，然后建立相应的微分方程模型，从而预测事物未来发展趋势的状况。

 GM(1, 1)模型的基本形式（也被称为灰色微分方程）：$$x^{(0)}(k)+a{z}^{(1)}(k)=b$$
 GM(1, 1)模型的白化方程：$$\frac{d{x}^{(1)}(t)}{dt}=-\hat{a}{x}^{(1)}(t)+\hat{b}$$

符号说明

符号	含义
$x^{(0)}$	原非负数据列
$x^{(1)}$	原数据的一次累加数据列
$z^{(1)}$	为$x^{(1)}$的紧邻均值数据列，下标从2开始
$a$	-$a$为发展系数
$b$	$b$为灰作用量

注意：变量上方加$\hat{}$表示变量的估计值。如：$\hat{b}$表示灰作用量的估计值。

适用条件

注意：使用GM(1, 1)模型。第一个‘1’表示微分方程是一阶的，后面的‘1’表示只有一个变量。

1.使用离散非负数据列（数据是以年份度量的非负数据）；
2.GM(1, 1)模型的本质是有条件的指数拟合，数据满足准指数规律；
3.数据的期数较短且和其他数据之间的关联性不强（小于等于10，也不能太短了，比如只有3期数据）
4.需检验GM(1, 1)模型对原数据的拟合程度。可以使用残差检验和级比偏差检验。

准指数规律检验

原序列的光滑比：$$\rho (k)=\frac{x^{(0)}(k)}{x^{(1)}(k-1)},k\ge 2$$原数据的光滑比小于0.5的个数占所有数据总个数的比例大于0.6表示可以通过准指数规律检验；或去掉前两个数据后，原数据的光滑比小于0.5的个数占所有数据总个数的比例大于0.9表示可以通过准指数规律检验。

注意：这里的0.6和0.9均为参考值，可以做出调整。必须注意，下标从2开始，因此数据总数为n-1。

残差检验

绝对残差：$$\epsilon (k)=x^{(0)}(k)-{\hat{x}}^{(0)}(k),k\ge 2$$
相对残差：$${\epsilon }_r(k)=\frac{|x^{(0)}(k)-{\hat{x}}^{(0)}(k)|}{x^{(0)}(k)}\times 100\%,k\ge 2$$
平均相对残差：$$\overline{{\epsilon }_r}=\frac{1}{n-1}\sum _{k=2}^n|{\epsilon }_r(k)|$$如果$\overline{{\epsilon }_r}<0.1$则认为GM(1, 1)模型对原数据的拟合效果非常不错，如果$\overline{{\epsilon }_r}<0.2$则认为GM(1, 1)模型对原数据的拟合效果达到了一般标准。

级比偏差检验

级比：$$\sigma (k)=\frac{x^{(0)}(k)}{x^{(0)}(k-1)},k\ge 2$$级比偏差：$$\eta (k)=\left|1-\frac{1-0.5\hat{a}}{1+0.5\hat{a}}\frac{1}{\sigma (k)}\right|,k\ge 2$$平均级比偏差：$$\bar{\eta }=\frac{1}{n-1}\sum _{k=2}^n\eta (k)$$$\eta (k)$越小说明$x^{(0)}(k)$和${\hat{x}}^{(0)}(k)$越接近。如果$\bar{\eta }<0.1$则认为GM(1, 1)模型对原数据的拟合效果非常不错，如果$\bar{\eta }<0.2$则认为GM(1, 1)模型对原数据的拟合效果达到了一般标准。

GM(1, 1)模型求解

 $\hat{b},\hat{a}$使用最小二乘法即可求出。

 求解GM(1, 1)模型的白化方程，即可用来预测：$$\frac{d{x}^{(1)}(t)}{dt}=-\hat{a}{x}^{(1)}(t)+\hat{b}$$
 求解微分方程得出：
$${\hat{x}}^{(0)}(m+1)=(1-e^{\hat{a}})\left[x^{(0)}(1)-\frac{\hat{b}}{\hat{a}}\right]e^{-\hat{a}m},m\ge 1$$新陈代谢gm11模型即：调用普通gm11模型预测的同时，淘汰旧值并加入新值的预测模型。新陈代谢gm11预测代码（主函数main.m）：

clear;clc

% 输入应为两组向量 分别表示自变量与因变量且一一对应
% 注意：灰色预测的变量数不应小于3也不应大于10， 且不能包含负数

x = input('自变量：');
y = input('因变量：');

n = length(y);

% 画出原始数据的时间序列图
figure(1);
plot(x, y, 'o-');grid on;
set(gca, 'xtick', x(1:end)); % 设置x坐标间隔为1
title('标题描述');
xlabel('x坐标描述');
ylabel('y坐标描述');

% 检验是否存在负值
if sum(y < 0) > 0 
    disp('灰色预测不能对存在负值的数据进行预测');
end

% 准指数规律检验。 灰色预测原理为：指数拟合，需要数据具有指数规律。
% 通常要求：光滑比小于0.5的数据大于60%；去掉前两个数据后光滑比小于0.5的数据大于90%

cum_y = cumsum(y); % 一次累加和
% 光滑度
rho = y(2:end) ./ cum_y(1:end-1);
disp('光滑比小于0.5的数据占比为：' );
disp(sum(rho < 0.5)/(n-1));

disp('去掉前两个数据后，光滑比小于0.5的数据占比为：');
disp(sum(rho(3:end) < 0.5)/(n-3));

% 级比偏差和相对残差检验
check(y);

% 新尘代谢gm11模型

p = input('请输入预测的次数：');
res = zeros(p, 1);

yy = y;

for i = 1 : p
    res(i) = gm11(yy, 1);
    for j = 1 : n - 1
        yy(j) = yy(j + 1);
    end
    yy(n) = res(i);
end

res

级比偏差和相对残差检验（check.m）:

function [] = check(d)
    % 对待预测数据进行级比偏差和相对残差检验
    % d 要进行预测的数据

    n = length(d);
    cum_d = cumsum(d);
    % 紧邻均值序列
    z = (cum_d(1:end-1) + cum_d(2:end)) ./ 2;
    
    % 最小二乘
    y = d(2:end);
    x = z;
    m = n - 1;
    k = (m*sum(x.*y) - sum(x)*sum(y))/(m*sum(x.*x)-sum(x)*sum(x));
    % b为灰作用量， -a为发展系数
    b = (sum(x.*x)*sum(y)-sum(x)*sum(x.*y))/(m*sum(x.*x)-sum(x)*sum(x));
    a = -k;
    
    % 预测
    p_d = zeros(n, 1);
    for i = 1 : n - 1
        p_d(i + 1) = (1-exp(a))*(d(1)-b/a)*exp(-a*i);
    end
    
    % 级比偏差
    u = abs(d(2:end)-p_d(2:end))./d(2:end);
    disp('均级比偏差为(通常小于0.1则说明适合灰色预测)：');
    disp(mean(u));
    
    % 相对残差
    disp('均级相对残差为(通常小于0.1则说明适合灰色预测)：');
    e = abs(1-(1-0.5*a)/(1+0.5*a)*(d(1:end-1)./d(2:end)));
    disp(mean(e));
end

普通gm11代码：

function res = gm11(d, p)
    % d 要进行预测的数据
    % res 预测结果， 为大小为p的列向量
    % p 为预测次数

    n = length(d);
    res = zeros(p, 1);
    cum_d = cumsum(d);
    % 紧邻均值序列
    z = (cum_d(1:end-1) + cum_d(2:end)) ./ 2;
    
    % 最小二乘
    y = d(2:end);
    x = z;
    m = n - 1;
    k = (m*sum(x.*y) - sum(x)*sum(y))/(m*sum(x.*x)-sum(x)*sum(x));
    % b为灰作用量， -a为发展系数
    b = (sum(x.*x)*sum(y)-sum(x)*sum(x.*y))/(m*sum(x.*x)-sum(x)*sum(x));
    a = -k;
    
    disp(strcat('发展系数为', num2str(-a)));
    disp(strcat('灰作用量为', num2str(b)));
    
    % 预测
    for i = 1 : p
        res(i) = (1-exp(a))*(d(1)-b/a)*exp(-a*(n+i-1));
    end
  
end