灰色模型算法文档.doc

灰色预测算法及相关程序1 引言32算法的基本原理32.1 GM（1，1）模型：32.2生成数42.2.1累加生成42.2.2累减生成53算法的具体实现流程63.1 算法流程图63.2 实现步骤83.3 数据准备与预处理104 算法程序实现104.1 程序使用说明104.2 程序源代码114.3 程序运行164.3.1程序运行及运行环境说明164.3.2 输入数据164.3.3 输出数据165 参考文献16灰色预测算法1 引言灰色预测（grey prediction）是利用灰色系统理论就灰色系统所作的预测灰色系统理论认为，尽管系统表象复杂，数据散乱，信息不充分，但作为系统，它必然有整体功能和内在规律，必然是有序的现有的分析方法大多依据过去的大量数据，按照统计方法分析其规律，这样不仅受数据量的限制，而且准确程度不高而灰色系统理论把随机量看作是在一定范围内变化的灰色量，对灰色量的处理不是寻求它的统计规律和概率分布，而是对原始数据加以处理，将杂乱无章的原始数据变为规律性较强的生成数据，通过对生成数据建立动态模型，来挖掘系统内部信息并充分利用信息进行分析预测目前，灰色系统理论用于预测主要通过GM（m，n）模型，该模型是灰色系统理论的量化体现，可用于以下几个方面的预测：（1）数列预测：对某个事物发展变化的大小与时间进行预测（2）灾变预测：预测灾变发生的时间或者说是异常值出现时区的分布如人体的血压过高或过低的时间预测（3）季节性灾变预测：对发生在每年特定时区的事件和命题作预测（4）拓扑预测：即事物整体的预测，亦称波形预测其特点是对于预先给定的多组数值建立GM（1，1）模型群，根据预测结果构造出整个波形（5）系统预测：对系统中众多变量间相互协调关系的发展变化所进行的预测2算法的基本原理2.1 GM（1，1）模型：灰色模型GM（1，1） GM（1，1）的含义为1阶，1个变量的灰色模型，它是在数据生成的基础上建立如下灰微分方程： 式中为原始序列，a称为发展系数，它反映和的发展态势；b称为灰作用量，它的大小反映数据变化的关系对序列，因为为与的平均值，故记为MEAN，即MEAN的白化型为：初始值用，则其解为：该式用于预测时称为时间响应函数，表示为 累减还原：其中（a,b）可通过最小二乘求解。2.2生成数灰色模型是将随机数经生成后变为有序的生成数据，然后建立微分方程，寻找生成数据的规律，再将运算结果还原的一种方法，其基础是数据的生成常用的生成方式有累加生成和累减生成2.2.1累加生成累加生成是将原始数据通过累加以生成新的数列记原始数列为:记的生成数列AGO为：其中称为的一次累加生成，记为1AGO定义的2AGO为：=AGO一般定义的r次AGO为：=AGO原始数据经累加生成后，其随机性明显减小，而规律性将增加对于非负的数据序列，累加生成的是一个递增的序列2.2.2累减生成累减生成是累加生成的逆运算，它是通过将原始序列前后两个数据相减生成新的数据序列累减生成可将累加生成还原为非生成数列，在建模中获得增量信息即：IAGO3算法的具体实现流程3.1 算法流程图3.2 实现步骤（1）数据检验与处理：具体算法：根据原始序列的级比的大小，判断GM（1，1）建模的可行性的定义为：可揭示数列的平滑情况或指数律的符合情况，若为常数，则具有白指数律（即确切的指数律），当然，此时无建立GM（1，1）模型的价值GM（1，1）建模希望的是的取值区间长度（覆盖的测度）接近于零而不等于零，其可容区为（01353，7389），即参数或级比的可容区是GM（1，1）建模的基本条件，但不是实用条件，要想建立满意有效的GM（1，1）模型，参数或级比应落于界区内设n是原始数据的个数，则我们有：a的界区：；的界区：，若，则认为是可作GM（1，1）建模的（2）建立模型GM(1,1)具体算法：利用算法思想的介绍建立模型，则可得累加预测值和预测值：模型的白化型为：初始值用，则其解为：该式用于预测时称为时间响应函数，表示为 累减还原；（3）检验预测值：残差检验：设为原始数列设为通过模型（11-52）和（11-53）得到的预测数据列则称为残差，称为GM（1，1）的残差相对值，若残差小于0.1，则称精度达到较高要求；若残差介于0.1与0.2之间，则称精度达到一般要求；若残差大于0.2，则模型不可靠。级比偏差值检验：与原始数列级比一样，我们定义模型级比推导可知，模型级比与k无关（因为符合指数律）其中a为发展系数定义级比偏差一般要求，最好；（4）由模型GM(1,1)作相应的预测。3.3 数据准备与预处理1）输入原始数据列。2）对原始数据做级比检验，若不满足要求，需要对原始数据作平移变换，直到满足要求为止。4 算法程序实现4.1 程序使用说明1程序采用Matlab语言2平移量由预测的人自行输入，输入值应逐渐增大3Matlab函数： GM(A)灰色预测函数；A初始输入数据向量nA的长度lamda级比向量P需要做平移的常数B一次累加向量alpha为分辨率系数，一般取0.5DB的灰导数a发展系统b灰作用量,且u=（a,b）Y提取A第2列后的所有值所构成的向量ZB的加权均值，也称为均值数列Z（白化背景值）T由Y扩展得到C原始累加数据预测值，且满足C(1)=A(1)C1原始数据(A)预测值,且满足C1(1)=A(1)C(i+1) 累加数据第i次预测值C1（j）原始数据第j次预测值eps(i) 预测值C1(i)与原始数据A（i）的相对误差，即残差ruo(i) 关于lamda的级比检验表达式m输入的值t预测天数，且t=m-nCt第t次累加数据预测值Ct_1第t-1次累加数据预测值pre_Ct第t次原始数据预测值，且pre_Ct=Ct-Ct_14.2 程序源代码function GM(A)%数据检验与处理n=length(A);%lamda为级比向量a=A(1);for j=1:n-1 lamda(j)=A(j)/A(j+1); %判断能否通过级比检验 if exp(-2/(n+1)lamda(j)&lamda(j)lamda(j) | lamda(j)exp(2/(n+1) P=input(平移系数c=); for i=1:n A(i)=A(i)+P; end for j=1:n-1 lamda(j)=A(j)/A(j+1); end end P=A(1)-a; end A;lamda;%对处理后的数据做一次累加，B为累加数列B=zeros(1,n);B(1)=A(1);for i=2:n j=i-1; B(i)=A(i)+B(j);endB;%求B的加权均值，即得均值数列Z（白化背景值）%alpha为分辨率系数，一般取0.5alpha=0.5;for i=1:n-1 Z(i)=alpha*B(i+1)+(1-alpha)*B(i);end%求B的灰导数D，由B数列做差可得D=A;%求发展系统a与灰作用量b,建立方程化简可得%提取A第2列后的所有值，称为数据向量YY=A(2:n);%将Y扩展为数据矩阵TT=(cat(1,-Z,ones(1,n-1);%由最小二乘法可得u=（a,b）u=(inv(T*T)*T*Y;%建立GM(1,1)模型求预测值%C是原始累加模型预测值，C1是原始数据预测值（均从原始数据第二项开始比较），可计算误差C(1)=A(1);for i=1:n-1 C(i+1)=(A(1)-u(2,1)/u(1,1)*exp(-u(1,1)*i)+u(2,1)/u(1,1);endC1(1)=A(1);for j=1:n-1 C1(j+1)=C(j+1)-C(j);endC;C1;%残差检验disp(以下为残差检验);for i=1:n eps(i)=(A(i)-C1(i)/A(i); if eps(i)0.1 disp(精度达到较高要求); elseif eps(i)0.1 disp(精度达到一般要求); elseif eps(i)0.2 disp(不能通过检验); endend%级比偏差检验disp(以下为级比偏差检验);for i=1:n-1 ruo(i)=1-(1-0.5*u(1,1)/(1+0.5*u(1,1)*lamda(i); if ruo(i)0.1 disp(精度达到较高要求); elseif ruo(i)0.1 disp(精度达到一般要求); elseif ruo(i)0.2 disp(不能通过检验); endend%数据检验到此结束*%*%* %改进预测算法等位新信息递补法，用于中长期动态预测%预测其后第t天%Ct为累加数列的预测值%pre_Ct为原始数据的预测值 m=input(请输入x=); while m0 disp(第,num2str(t),天的预测值为：); pre_Ct=pre_Ct-P elseif P=0 disp(第,num2str(t),天的预测值为：); pre_Ct end %等位新信息递补，即除去第一个元素，引入一个预测的新元素 for i=1:n-1 A(i)=A(i+1); end A(n)=pre_Ct+P; %判断是否结束预测 if t=m-n break; else%未结束，继续检验新得到的数据 %判断能否通过级比检验 a=A(1); for j=1:n-1 lamda(j)=A(j)/A(j+1); if exp(-2/(n+1)lamda(j)&lamda(j)lamda(j) | lamda(j)exp(2/(n+1) P1=input(平移系数c=);%整体平移数据 for i=1:n A(i)=A(i)+P1; end%计算级比 for j=1:n-1 lamda(j)=A(j)/A(j+1); end end%此循环所得到的平移总量 P1=A(1)-a; P1; end end %对处理后的数据做一次累加，B为累加数列 B(1)=A(1); for i=2:n j=i-1; B(i)=A(i)+B(j); end B; %求B的加权均值，即得均值数列Z（白化背景值）%alpha为分辨率系数，一般取0.5 alpha=0.5; for i=1:n-1 Z(i)=alpha*B(i+1)+(1-alpha)*B(i); end %求B的灰导数D，由B数列做差可得 D=A; %求发展系统a与灰作用量b,建立方程化简可得%提取A第2列后的所有值，称为数据向量Y Y=A(2:n); %将Y扩展为数据矩阵T T=(cat(1,-Z,ones(1,n-1); %由最小二乘法可得u=（a,b） u=(inv(T*T)*T*Y; end4.3 程序运行4.3.1程序运行及运行环境说明运行环境：windows操作系统运行软件：MATLAB适用范围：只能进行短期或中期的预测4.3.2 输入数据输入实际数据向量：a=2.28,2.98,3.39,4.24,6.86,8.64,11.85,12.15,12.71;输入平移量：c1=3,c2=4,c3=5;4.3.3 输出数据原始数据2.282.983.394.246.868.6411.8512.1512.71预测数据2.28002.89594.02135.27476.67068.22539.956811.885114.0328残差00.0084-0.0608-0.0920.0137 0.0265 0.10040.0138-0.0671级比偏差值检验值-0.0357-0.0699-0.02960.0967 0.01290.0758-0.0964 -0.0822由以上两表可见，残差值均小于0.1，模型精度较高；级比偏差值小于0.1，达到较高要求下面运用此模型预测并输出其后两天的值。第1天的预测值为：pre_Ct =16.4247，最终平移量：P =