MATLAB学习系列29.-灰色预测

29.灰色预测一、灰色预测概述灰色预测是用灰色模型GM(1,1)来进行定量分析的，通常分为以下几类：(1)灰色时间序列预测。用等时距观测到的反映预测对象特征的一系列数量（如产量、销量、人口数量、存款数量、利率等）构造灰色预测模型，预测未来某一时刻的特征量，或者达到某特征量的时间。(2)畸变预测（灾变预测）。通过模型预测异常值出现的时刻，预测异常值什么时候出现在特定时区内。(3)波形预测，或称为拓扑预测，它是通过灰色模型预测事物未来变动的轨迹。(4)系统预测，对系统行为特征指标建立一族相互关联的灰色预测理论模型，在预测系统整体变化的同时，预测系统各个环节的变化。上述灰色预测方法的共同特点是：（1）允许少数据预测；（2）允许对灰因果律事件进行预测，例如灰因白果律事件：在粮食生产预测中，影响粮食生产的因子很多，多到无法枚举，故为灰因，然而粮食产量却是具体的，故为白果。粮食预测即为灰因白果律事件预测。白因灰果律事件：在开发项目前景预测时，开发项目的投入是具体的，为白因，而项目的效益暂时不很清楚，为灰果。项目前景预测即为灰因白果律事件预测。（3）具有可检验性，包括：建模可行性的级比检验（事前检验），建模精度检验（模型检验），预测的滚动检验（预测检验）。二、GM(1,1)模型1.模型理论GM(1,1)模型适合具有较强的指数规律的数列，只能描述单调的变化过程。已知元素序列数据：(0)(0)(0)(0)((1),(2),,())Xxxxn做一次累加生成（1-AGO）序列：(1)(1)(1)(1)((1),(2),,())Xxxxn其中，(1)(0)1()(),1,,kixkxikn令(1)Z为(1)X的紧邻均值生成序列：(1)(1)(1)(1)((2),(3),,())Zzzzn其中，(1)(1)(1)()0.5()0.5(1)zkxkxk建立GM(1,1)的灰微分方程模型为：bkazkx)()()1()0(其中，a为发展系数，b为灰色作用量。设ˆ为待估参数向量，即ˆ(,)Tab，则灰微分方程的最小二乘估计参数列满足1ˆ()TTnBBBY其中B=(1)(1)(1)(2)1(3)1......()1zzzn，nY=(0)(0)(0)(2)(3)...()xxxn再建立灰色微分方程的白化方程（也叫影子方程）：(1)(1)dxaxbdt白化方程的解（也叫时间响应函数）为(1)(1)ˆ()((0))atbbxtxeaa那么相应的GM(1,1)灰色微分方程的时间响应序列为：(1)(1)ˆ(+1)[(0)],1,,akbbxkxeknaa取(1)(0)(0)(1)xx，则(1)(0)ˆ(+1)[(1)],1,,1akbbxkxeknaa再做累减还原可得(0)(1)(1)(0)ˆˆˆ(1)(1)()[(1)](1),1,,1aakxkxkxkbxeekna即为预测方程。注1：原始序列数据不一定要全部使用，相应建立的模型也会不同，即a和b不同；注2：原始序列数据必须要等时间间隔、不间断。2.算法步骤(1)数据的级比检验为了保证灰色预测的可行性，需要对原始序列数据进行级比检验。对原始数据列(0)(0)(0)(0)((1),(2),,())Xxxxn，计算序列的级比：(0)(0)(1)(),2,,()xkkknxk若所有的级比()k都落在可容覆盖2/(1)2/(2)(,)nnee内，则可进行灰色预测；否则需要对(0)X做平移变换，(0)(0)YXc，使得(0)Y满足级比要求。(2)建立GM(1,1)模型，计算出预测值列。(3)检验预测值：①相对残差检验，计算(0)(0)(0)ˆ()()(),1,,()xkxkkknxk若()0.2k，则认为达到一般要求，若()0.1k，则认为达到较高要求；②级比偏差值检验根据前面计算出来的级比()k,和发展系数a,计算相应的级比偏差：10.5()1()()10.5akka若()0.2k,则认为达到一般要求，若()0.1k,则认为达到较高要求。(4)利用模型进行预测。例1某大型企业1999年至2004年的产品销售额如下表，试建立GM(1,1)预测模型，并预测2005年的产品销售额。年份199920002001200220032004销售额(亿元)2.673.133.253.363.563.72GM(1,1)函数：function[pre,f]=GM11(x)%返回第n+1个预测值pre，和预测函数f%级比检验n=length(x);lambda=x(1:n-1)./x(2:n);%计算级比range=[exp(-2/(n+1)),exp(2/(n+2))];%可容覆盖的范围ifrange(1)min(lambda)=max(lambda)range(2)disp('级比检验通过');elsedisp('级比检验未通过');end%GM(1,1)建模x1=cumsum(x);%一次累加n=length(x1);z1=(x1(1:n-1)+x1(2:n))/2;Y=x(2:n)';%构造矩阵YB=[-z1',ones(n-1,1)];%构造矩阵BA=(B'*B)\B'*Y;%计算模型的参数a,bk=1:n;x1=(x(1)-A(2)/A(1))*exp(-A(1)*(k-1))+A(2)/A(1);%利用模型计算累加值的预测值x_p=[x(1),diff(x1)];%累减还原到预测值Delta=abs(x_p-x);%绝对残差序列phi=Delta./x;%相对残差序列ifmax(phi)=0.2disp('相对残差检验未通过');elsedisp('相对残差检验通过');endmean(phi);%平均相对残差rho=1-(1-0.5*A(1))/(1+0.5*A(1))*lambda;%计算级比偏差值ifmax(rho)=0.2disp('级比偏差检验未通过');elsedisp('级比偏差检验通过');end%预测第n+1个数据f=@(t)(x(1)-A(2)/A(1))*(1-exp(A(1)))*exp(-A(1)*t);%预测公式pre=f(n);主程序：x0=[2.67,3.13,3.25,3.36,3.56,3.72];%原始数列x_pre=GM11(x0)运行结果：级比检验通过相对残差检验通过级比偏差检验通过x_pre=3.8756注1：GM(1,1)模型的检验，还可以通过检验预测序列与原始序列的关联度R，分辨率为0.5时，R0.6结果可以接受；以及后验差检验，即对残差分布的统计特性进行检验（略）。注2：当原始数据序列建立的GM(1,1)模型检验不合格或结果不够精确时，可以用GM(1,1)残差模型（即对残差数列使用GM(1,1)模型）来修正或提高精度。3.数据融合算法为了进一步提高灰色预测的精度，可以：先用原始数据列的n个数据、后n-1个数据、……、后n-k个数据，进行多次GM(1,1)预测，得到k+1个预测值，再将它们进行数据融合得到最终的预测值。数据融合算法原理：设1,,maa为m个由GM(1,1)模型预测得到的第n+1时刻的预测值，定义任意两个值之间的距离：||,,1,,ijijdaaijm构造两个数据间的支持度函数：cos(),,1,,2max{}ijijijdrijmd满足：(1)ijr与相对距离成反比，即两个值相差越大，彼此间的支持程度越小；(2)[0,1]ijr,使数据的处理能够利用模糊集理论中隶属度函数的优点，避免数据之间相互支持度的绝对化。建立支持度矩阵：{}ijmmRr。为了从1,,maa融合得到最终的a，需要确定每个ia的权重iw满足11miiw，注意到iw应综合包含1,,iimrr的信息，从而要寻找一组非负数1,,mvv使得1mijijjwvr，其矩阵形式为WRV。由于R为非负对称矩阵，故存在最大特征值及对应的（非负）特征向量1[,,]mVvv。由特征向量与特征值的性质，可取1/miijjwvv融合后得到1miiiawa。Matlab实现：function[a,w]=DataFusion(x)%x为长度≥2的行向量,返回a为融合值,w为权重向量Lx=length(x);ifLx==2a=mean(x);w=[1/2,1/2];elseIndCom=nchoosek(1:Lx,2);%x中元素下标索引的所有两两组合方式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;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)运行结果：pre=3.92053.89433.8756x_pre=3.89484.灰色预测的应用——灾变预测灰色灾变预测的是给出下一个或几个异常值出现的时刻，以便人们提前防备，采取对策，减少损失。基本思路步骤：1.对于给定的原始时间序列数据(1),(2),,()Xxxxn，根据灾变界值选取出灾变子列：[(1)],[(2)],,[()]Xxqxqxqm则得到灾变日期序列：(0)(1),(2),,()Qqqqm2.对灾变日期序列(0)(1),(2),,()Qqqqm进行GM(1,1)预测，得到未来的灾变日期ˆˆ(1),(2),qmqm例3某地区平均降水量（毫米）的原始数据为：X={386.6,514.6,434.1,484.1,647.0,399.7,498.7,701.6,254.5,463.0,745.0,398.3,554.5,471.1,384.5,242.5,671.7,374.7,458.9,511.3,530.8,586.0,387.1,454.4},规定年降水量390(毫米)为旱灾年，试作旱灾预测。代码：X0=[386.6514.6434.1484.1647.0,399.7498.7701.6254.5463.0745.0398.3...554.5471.1384.5242.5671.7374.7458.9511.3530.8586.0387.1454.4];Xq=X0(X0=390)Q0=find(X0=390)%提取灾变日期序列[Q1,f]=GM11(Q0);%GM(1,1)预测Q1%预测未来第1个灾变日期m=length(Q0);f(m+1:m+3)%预测未来第1-3个灾变日期运行结果：Xq=386.6000254.5000384.5000242.5000374.7000387.1000Q0=1915161823Q1=27.4879ans=27.487933.187340.0684可见，GM(1,1)模型预测的未来3个灾变日期分别是第27年，第33年，第40年。三、GM(2,1)模型简述GM(2,1)模型为：(1)(0)(0)(1)12()()()xkaxkazkb其中，(1)(0)()xk为1次累减生成序列，(1)()zk为紧邻均值生成序列。其白化方程为：2(1)(1)