一.灰色预测
1.灰色系统下的灰色预测
<1>什么是灰色系统?
所谓的灰色系统其实就是夹杂在白色系统和黑色系统之中的一种系统,而白色系统就是全部信息已知的系统,黑色系统就是全部信息未知的系统。所以,夹在这两种系统中间的灰色系统就是部分信息已知,而部分信息也是未知的系统。
<2>什么是灰色预测?
灰色预测就是在灰色系统中所作的预测。灰色预测就是在部分信息已知,而部分信息也是未知的前提下做的一种预测分析。灰色预测通过鉴别个因素之间的差异程度,进行关联分析,对原始数据处理后生成一定规律性的序列,然后建立相应的微分方程模型,从而预测事物未来的发展趋势,最后得到其发展的模型。
2.如何实现灰色预测GM(1,1)
<1>了解灰色生成
①累加生成: 通过数列间各时刻数据的依个累加以得到新的数据与数列。
(最开始看的时候我还不信累加生成能增强数据规律性所以我用MATLAB画了个图,对于[6 3 8 10 7]这组数据好像确实累加后规律更明显,累加后的数列呈指数增长趋势)
②累减生成: 将原始序列前后两个数据相减得到累减生成序列。(将累加的数据累减就得到原始序列)
③均值生成: 将前后两个数据相加在除以2得到均值生成序列。
<2>实现GM(1,1)——一阶微分方程,且只含一个变量的灰色模型
原始序列
x(0)=(x(0)(1),x(0)(2),x(0)(3),...,x(0)(n))x^{(0)}=(x^{(0)}(1),x^{(0)}(2),x^{(0)}(3),...,x^{(0)}(n))x(0)=(x(0)(1),x(0)(2),x(0)(3),...,x(0)(n))
先累加操作得到序列
x(1)=(x(1)(1),x(1)(2),x(1)(3),...,x(1)(n))x^{(1)}=(x^{(1)}(1),x^{(1)}(2),x^{(1)}(3),...,x^{(1)}(n))x(1)=(x(1)(1),x(1)(2),x(1)(3),...,x(1)(n))x(1)(t)=∑k=1tx(1)(k)x^{(1)}(t)=\sum_{k=1}^{t}x^{(1)}(k)x(1)(t)=k=1∑tx(1)(k)
然后均值操作累加生成的数列
z(1)=(z(1)(2),z(1)(3),...,z(1)(n))z^{(1)}=(z^{(1)}(2),z^{(1)}(3),...,z^{(1)}(n))z(1)=(z(1)(2),z(1)(3),...,z(1)(n))z(1)(k)=0.5x(1)(k)+0.5x(1)(k)z^{(1)}(k)=0.5x^{(1)}(k)+0.5x^{(1)}(k)z(1)(k)=0.5x(1)(k)+0.5x(1)(k)
对累加生成序列建立微分方程——白化微分方程(这里要说一下因为刚刚提过累加后的数列呈指数增长,所以可以用指数函数进行拟合所以建立如下微分方程)如下微分方程用数学知识求得明显为指数关系
dx(1)dt+ax(1)=b\frac{dx^{(1)}}{dt}+ax^{(1)}=bdtdx(1)+ax(1)=b
然后建立灰微分方程这里是知乎网友解释的建立过程
x(0)(k)+az(1)(k)=bx^{(0)}(k)+az^{(1)}(k)=bx(0)(k)+az(1)(k)=b
aaa为发展系数,bbb为灰色作用量
根据最小二乘法求aaa,bbb其中a^\hat{a}a^=[aaa;bbb]
求解刚刚建立的白化微分方程
x(1)(t)=(x(1)(0)−ba)e−at+bax^{(1)}(t)=(x^{(1)}(0)-\frac{b}{a})e^{-at}+\frac{b}{a}x(1)(t)=(x(1)(0)−ab)e−at+ab
取x(1)(0)=x(0)(1)x^{(1)}(0)=x^{(0)}(1)x(1)(0)=x(0)(1),则
x(1)(k+1)=(x(0)(1)−ba)e−ak+bax^{(1)}(k+1)=(x^{(0)}(1)-\frac{b}{a})e^{-ak}+\frac{b}{a}x(1)(k+1)=(x(0)(1)−ab)e−ak+abk=1,2,3,⋯,n−1k=1,2,3,\cdots,n-1k=1,2,3,⋯,n−1
最后累减还原即为预测方程x(0)(k+1)=x(1)(k+1)−x(1)(k)=(x(0)(1)−ba)(1−ea)e−akx^{(0)}(k+1)=x^{(1)}(k+1)-x^{(1)}(k)=(x^{(0)}(1)-\frac{b}{a})(1-e^a)e^{-ak}x(0)(k+1)=x(1)(k+1)−x(1)(k)=(x(0)(1)−ab)(1−ea)e−akk=1,2,3,⋯,n−1k=1,2,3,\cdots,n-1k=1,2,3,⋯,n−1
注意:①原始序列数据不一定要全部使用,相应建立的模型也会不同,即aaa和bbb不同。②原始序列数据必须要等时间间隔、不间断。
<3>灰色模型检验
① 计算x(0)(t)x^{(0)}(t)x(0)(t)与x^(0)(t)\hat{x}^{(0)}(t)x^(0)(t)之间的残差e(0)(t)e^{(0)}(t)e(0)(t)和相对误差q(x)q(x)q(x):
e(0)(t)=x(0)(t)−x^(0)(t)e^{(0)}(t)=x^{(0)}(t)-\hat{x}^{(0)}(t)e(0)(t)=x(0)(t)−x^(0)(t)q(x)=e(0)(t)x(0)(t)q(x)=\frac{e^{(0)}(t)}{x^{(0)}(t)}q(x)=x(0)(t)e(0)(t)
② 求原始数据的均值x0x^0x0以及方差s1s_1s1。
③ 求e(0)(t)e^{(0)}(t)e(0)(t)的平均值qˉ\bar{q}qˉ以及残差的方差s2s_2s2。
④ 计算方差比C=s2s1C=\frac{s_2}{s_1}C=s1s2。
⑤ 求小误差概率P=P(∣e(t)∣<0.6745s1)P=P({|e(t)|<0.6745s_1})P=P(∣e(t)∣<0.6745s1)。
⑥ 灰色精度检验如下表
等级 | 相对误差qqq | 方差比CCC | 小误差概率PPP |
---|---|---|---|
ⅠⅠⅠ级 | <0.01<0.01<0.01 | <0.35<0.35<0.35 | >0.95>0.95>0.95 |
ⅡⅡⅡ级 | <0.05<0.05<0.05 | <0.50<0.50<0.50 | <0.80<0.80<0.80 |
ⅢⅢⅢ级 | <0.10<0.10<0.10 | <0.65<0.65<0.65 | <0.70<0.70<0.70 |
ⅣⅣⅣ级 | >0.20>0.20>0.20 | >0.80>0.80>0.80 | <0.60<0.60<0.60 |
在实际应用的过程中,检验模型精度的方法并不唯一。可以利用上述方法进行模型的检验,也可以根据q(x)q(x)q(x)的误差百分比进行结合预测数据与实际数据之间的测试结果酌情认定模型是否合理。
<4>利用灰色模型预测
x^(0)\hat{x}^{(0)}x^(0)= [x^(0)(1),x^(0)(2),⋯,x^(0)(n)⏟原数列的模拟,x^(0)(n+1),x^(0)(n+2),⋯,x^(0)(n+m)⏟未来数列的预测][ \begin{matrix} \underbrace{\hat{x}^{(0)}(1),\hat{x}^{(0)}(2),\cdots,\hat{x}^{(0)}(n)} \\ 原数列的模拟 \end{matrix},\begin{matrix} \underbrace{\hat{x}^{(0)}(n+1),\hat{x}^{(0)}(n+2),\cdots,\hat{x}^{(0)}(n+m)} \\ 未来数列的预测\end{matrix} ][x^(0)(1),x^(0)(2),⋯,x^(0)(n)原数列的模拟,x^(0)(n+1),x^(0)(n+2),⋯,x^(0)(n+m)未来数列的预测]
二.灰色预测例子及代码
1.长江污水排放预测
请以上表的数据为依据,预测2005-2014年长江的污水排放量(单位:亿吨)。
2.长江污水排放预测MATLAB代码
%2020/2/14,灰色模型学习进行预测2005-2014年长江的污水排放量(单位:亿吨)。
clc;clear;close all;
%% 数据准备
syms a b; %建立符号变量a(发展系数)和b(灰作用量)
c=[a,b]';
A = [174, 179, 183, 189, 207, 234, 220.5, 256, 270, 285];%原始数列 A
n = length(A); %求出A中元素个数
years=10; %预测未来年数
B = cumsum(A); %对原始数列 A 做累加得到数列 B
for i = 2:n %对数列 B 做紧邻均值生成C(i) = (B(i) + B(i - 1))/2;
end
C(1) = []; %这里注意C中元素个数为n-1所以删去第一列
B=[-C;ones(1,n-1)]'; %构造数据矩阵
Y=A(2:n)';
%% 最小二乘法计算参数并求出预测数据
c=inv(B'*B)*B'*Y; %使用最小二乘法计算参数 a(发展系数)和b(灰作用量)
a=c(1);b=c(2);
%预测未来数据
F = []; F(1) = A(1);
for i = 2:(n+years)F(i) = (A(1)-b/a)/exp(a*(i-1))+ b/a;
end
G = []; G(1) = A(1); %对数列 F 累减还原,得到预测出的数据
for i = 2:(n+years)G(i) = F(i) - F(i-1); %得到预测出来的数据
end
disp(['预测数据为:',num2str(G)]);
%% 模型检验
%计算残差
H = G(1:10);
epsilon = A - H;
%计算相对误差q
q= abs(epsilon./A);
Q=mean(q);
disp(['相对残差Q检验:',num2str(Q)]);
%方差比C检验
C = std(epsilon, 1)/std(A, 1);
disp(['方差比C检验:',num2str(C)]);
%小误差概率P检验
S1 = std(A, 1);
temp = find(abs(epsilon - mean(epsilon))< 0.6745 * S1); %找到数据中满足条件数据位置
P = length(temp)/n; %求出满足条件个数然后除以总个数求出概率
disp(['小误差概率P检验:',num2str(P)]);
%% 绘图
t1 = 1995:2004;
t2 = 1995:2014;
plot(t1, A,'r*'); hold on;
plot(t2, G, 'b-');
xlabel('年份'); ylabel('污水量/亿吨');
legend('实际污水排放量','预测污水排放量');
title('长江污水排放量增长曲线');
grid on; %显示网格
3.运行结果
预测数据为:174 172.809 183.9355 195.7785 208.3839 221.801 236.082 251.2825 267.4616 284.6825 303.0122 322.5221 343.2881 365.3912 388.9175 413.9585 440.6118 468.9812 499.1772 531.3174
相对残差Q检验:0.023399
方差比C检验:0.18697
小误差概率P检验:1
2020/2/14灰色预测学习笔记,代码参考《数学建模算法与应用》