MATLAB学习系列29. 灰色预测

即为灰因白果律事件预测。 （3） 具有可检验性， 包括： 建模可行性的级比检验 （事前检验） ， 建模精度检验（模型检验），预测的滚动检验（预测检验）。
二、GM(1,1)模型 1. 模型理论 GM(1,1)模型适合具有较强的指数规律的数列，只能描述单调的 变化过程。已知元素序列数据：
X (0) ( x (0) (1), x (0) (2), , x (0) (n ))
(2) 建立 GM(1,1)模型，计算出预测值列。 (3) 检验预测值： ① 相对残差检验，计算
(k )
ˆ (0) ( k ) x (0) ( k ) x , x (0) ( k ) k 1, , n
若 (k ) 0.2 ，则认为达到一般要求，若 (k ) 0.1 ，则认为达到较高要 求； ② 级比偏差值检验 根据前面计算出来的级比 (k ) , 和发展系数 a, 计算相应的级比 偏差：
j 1
由于 R 为非负对称矩阵，故存在最大特征值 及对应的（非负）
]。 特征向量 V [v1 , , vm 由特征向量与特征值的性质， 可取 wi vi / v j
j 1 m
融合后得到 a wi ai 。
i 1
பைடு நூலகம்
m
Matlab 实现：
function [a,w]=DataFusion(x) %x为长度≥2的行向量, 返回a为融合值, w为权重向量 Lx=length(x); if Lx==2 a=mean(x); w=[1/2,1/2]; else IndCom=nchoosek(1:Lx,2); 方式 d=abs(x(IndCom(:,1))-x(IndCom(:,2))); 间的距离 maxd=max(d); [Y,X]=meshgrid(1:Lx,1:Lx); R=cos(pi*(x(X)-x(Y))/(2*maxd)); %构造支持度矩阵 [V,D]=eig(R); w=V(:,Lx)/sum(V(:,Lx)); a=x*w; %计算任意两值之 %x中元素下标索引的所有两两组合
年份 销售额(亿元) 1999 2.67 2000 3.13 2001 3.25 2002 3.36 2003 3.56 2004 3.72
GM(1,1)函数：
function [pre,f]=GM11(x) %返回第n+1个预测值pre，和预测函数f %级比检验 n=length(x); lambda=x(1:n-1)./x(2:n); %计算级比
(k )
x (0) (k 1) , x (0) (k ) k 2, , n
若所有的级比 (k ) 都落在可容覆盖 ( e 2/( n1) , e 2/( n 2) ) 内， 则可进行灰色 预测；否则需要对 X (0) 做平移变换， Y (0) X (0) c ，使得 Y (0) 满足级比 要求。
建立支持度矩阵： R {rij }mm 。 为了从 a1 ,, am 融合得到最终的 a，需要确定每个 ai 的权重 wi 满足
w
i 1
m
i
1 ，注意到 wi 应综合包含 ri1 , , rim 的信息，从而要寻找一组非负
m
数 v1 ,, v m 使得 wi v j rij ，其矩阵形式为 W RV 。
即为预测方程。
注 1：原始序列数据不一定要全部使用，相应建立的模型也会不 同，即 a 和 b 不同； 注 2：原始序列数据必须要等时间间隔、不间断。
2. 算法步骤 (1) 数据的级比检验
为了保证灰色预测的可行性， 需要对原始序列数据进行级比检验。 对原始数据列 X (0) ( x (0) (1), x (0) (2), , x (0) (n )) ，计算序列的级比：
3. 数据融合算法 为了进一步提高灰色预测的精度，可以：先用原始数据列的 n 个 数据、后 n-1 个数据、……、后 n-k 个数据，进行多次 GM(1,1)预测， 得到 k+1 个预测值，再将它们进行数据融合得到最终的预测值。 数据融合算法原理： 设 a1 ,, am 为 m 个由 GM(1,1)模型预测得到的第 n+1 时刻的预测 值，定义任意两个值之间的距离：
Delta=abs(x_p-x); %绝对残差序列 phi=Delta./x; %相对残差序列 if max(phi) >= 0.2 disp('相对残差检验未通过'); else disp('相对残差检验通过'); end mean(phi); %平均相对残差 rho=1-(1-0.5*A(1))/(1+0.5*A(1))*lambda; %计算级比偏差值 if max(rho) >= 0.2 disp('级比偏差检验未通过'); else disp('级比偏差检验通过'); end %预测第n+1个数据 f=@(t) (x(1)-A(2)/A(1))*(1-exp(A(1)))*exp(-A(1)*t); % 预测公式 pre=f(n);
运行结果： pre = 3.9205 x_pre = 3.8948 3.8943 3.8756
4. 灰色预测的应用——灾变预测 灰色灾变预测的是给出下一个或几个异常值出现的时刻， 以便人 们提前防备，采取对策，减少损失。 基本思路步骤： 1. 对于给定的原始时间序列数据 X x (1), x (2), , x (n ) ， 根据灾变 界值选取出灾变子列：
range=[exp(-2/(n+1)),exp(2/(n+2))]; %可容覆盖的范围 if range(1) < min(lambda) <= max(lambda) < range(2) disp('级比检验通过'); else disp('级比检验未通过'); end %GM(1,1)建模 x1=cumsum(x); %一次累加 n=length(x1); z1=(x1(1:n-1)+x1(2:n))/2; Y=x(2:n)'; %构造矩阵Y B=[-z1',ones(n-1,1)]; %构造矩阵B A=(B'*B)\B'*Y; %计算模型的参数a,b k=1:n; x1=(x(1)-A(2)/A(1))*exp(-A(1)*(k-1))+A(2)/A(1); %利用模 型计算累加值的预测值 x_p=[x(1),diff(x1)]; %累减还原到预测值
end
例 2 对于例 1 的数据， 分别对近 6 年 （1999-2004） 、 近5年 （2000-2004） 、 近4年 （2001-2004） 的数据用 GM(1,1)预测出 2005 年的 3 个预测值， 再进行数据融合，得到最终的预测值。
代码：
x0=[2.67,3.13,3.25,3.36,3.56,3.72]; %原始数列 Lx=length(x0); pre=arrayfun(@(k) GM11(x0(Lx-k:end)),3:5) %最后一个预测值 是用全部6个数据 x_pre=DataFusion(pre)
29. 灰色预测
一、灰色预测概述 灰色预测是用灰色模型 GM(1,1)来进行定量分析的，通常分为以 下几类： (1) 灰色时间序列预测。用等时距观测到的反映预测对象特征的 一系列数量（如产量、销量、人口数量、存款数量、利率等）构造灰 色预测模型， 预测未来某一时刻的特征量， 或者达到某特征量的时间。 (2) 畸变预测（灾变预测）。通过模型预测异常值出现的时刻， 预测异常值什么时候出现在特定时区内。 (3) 波形预测，或称为拓扑预测，它是通过灰色模型预测事物未 来变动的轨迹。 (4) 系统预测，对系统行为特征指标建立一族相互关联的灰色预 测理论模型， 在预测系统整体变化的同时， 预测系统各个环节的变化。 上述灰色预测方法的共同特点是： （1）允许少数据预测； （2）允许对灰因果律事件进行预测，例如 灰因白果律事件： 在粮食生产预测中， 影响粮食生产的因子很多， 多到无法枚举，故为灰因，然而粮食产量却是具体的，故为白果。粮 食预测即为灰因白果律事件预测。 白因灰果律事件：在开发项目前景预测时，开发项目的投入是具 体的，为白因，而项目的效益暂时不很清楚，为灰果。项目前景预测
X x[ q(1)], x[ q(2)], , x[ q( m )]
dij | ai a j |, i, j 1,, m
构造两个数据间的支持度函数：
rij cos(
d ij ), 2 max{d ij }
i , j 1, , m
满足：(1) rij 与相对距离成反比，即两个值相差越大，彼此间的支持 程度越小； (2) rij [0,1] , 使数据的处理能够利用模糊集理论中隶属度函数的 优点，避免数据之间相互支持度的绝对化。
取 x (1) (0) x (0) (1) ，则
b b ˆ (1) (k +1) [ x (0) (1) ]e ak , x a a k 1,, n 1
再做累减还原可得
ˆ (0) (k 1) x ˆ (1) (k 1) x ˆ(1) (k ) x b [ x (0) (1) ](1 e a ) e ak , a k 1, , n 1
(k ) 1 (
1 0.5a ) ( k ) 1 0.5a
若 (k ) 0.2 , 则认为达到一般要求， 若 (k ) 0.1 , 则认为达到较高要求。 (4) 利用模型进行预测。
例 1 某大型企业 1999 年至 2004 年的产品销售额如下表，试建立
GM(1,1)预测模型，并预测 2005 年的产品销售额。
再建立灰色微分方程的白化方程（也叫影子方程）：
dx (1) ax (1) b dt
白化方程的解（也叫时间响应函数）为