灰色预测MATLAB程序资料讲解

灰色预测MATLAB程

序

灰色预测。设工◎=（工叫建工⑼工网⑺）为原始数列，其1次累分加生成数列为”=（/◎”©-/叫齐》,其中£'&）=工T⑨（立k=12…二退i-l。定义”）•的灰导数为成幻=工⑺的=户（后-戈田仇-1）一令产为数列/）的邻值生成数列，即工?左）=s叫此）十（1—ct）jra）（i-1）.于是定义GM（1.1）的灰微分方程模型为d（k）+口=⑴（左）=瓦即或 工⑼⑶+⑷⑴⑶=瓦 ⑴在式C1）中7⑺（衿称为灰导数，强称为发展系数，方®称为白化背景值，b称为灰作用量。将时刻表收=23,造代入⑴式有何Q）—里如（2）（■）也、_E8小引入矩阵向量记号于是GM(1-1)模型可表示为F=B/现在问题归结为求队b在值❷用一元线性回归，即最小二乘送求它们的拈计值为 &=[[=(£『团注二实际上回归分析中求去计值是用软件计算的，有标度程序求解，如出血疑等。GM(1.1)的白北型对于GM(1-1)的灰微分方程(D，如果将灰导救茎闹的时到上=工3</视为连续费里厂则工3视为时间t函敕产3-于是工I上)对应于导敕里级由“⑴，白化背景值冽体》对应于导颗一.于是Gy(1，1)的灰戳分方锭对应于的白制分方程为-加%)=2 (2)G31(1>1)G31(1>1)灰色独副的步霸L数据的检吩与处理为了保证GM(1,1)建攫方法的可行性，需要时已知数媚做必要的检哈处理<>被原始数提列为了丛0二世⑨①公力②,睛”3)；-计算孰列的组比

工⑸=工：潟]：比=，丸=如果斫有的策比都落在可喜覆盖区间 才-(:『斗」三)向津则数据列Hg可以建立GXICl.1)模型且可以遵行施色强娜•否则，对数据他适当的奏振处理，如平移变换;取。使潺救据列/叫左)二工期体)+Q#=12…；&的蛆比都蕉在可容覆盖内。

.建立GMCl»1)槿型不妨设m=0用炉冉…产叫项莉足上面的要求，以它为数据列建立GM。，1，模型W(k)+az.建立GMCl»1)槿型不妨设m=0用炉冉…产叫项莉足上面的要求，以它为数据列建立GM。，1，模型W(k)+azU>(A)=b,I用回归分析求得出b的估计值・于是相应的白化模型为成;十“Q"。一瓦解为婢的=(工冏Q)—士州7-9一匕 (3)<7 G于是得到颈皿值却\而+1)=。⑼(1)一日州-龙+白承=12…：网一1;aaM而相应地得到预则值二皿)(t+1)=xcl)(t+l)-ia)(i)3i=1,2,,n-L.检脸籁则值(1)残差桧缝;计算相对残差g留徜-货般,17一= 收团==12一一如果对斫有的<0.1"则认为达到较高的要求；否则，若对所有的|式©K口工，则认为达到一艇要求口(2)艇比偏差值横嗓二计算雄)=1-詈辔靛尬1-0.如果对斫有的|p(A)|<0.1 ，则认为达到较高的要求m否则若对所有的I四©%02，则认为达到一般要求。小例北方某城市1986〜1992年道路交通噪声平均声.据见知 序号年附219S771A319S872.44219S771A319S872.441OT372A3阚了1.4619917207199771.6第一步:纵比检骗建立胃通噪薛均声级数据时间序列如下：=(.叫1),钟⑵:钟⑺)<711,72,4,72.4,711n71.4,720,71.6)◎")求级比,出)M通留位-。叁，=(忒2),a⑶「1(7》"总=(0.982JJ.0042J,0098.0.9917J,0056)⑵级比判断由于所有的入®曰［0.982」.009间,k=23 7故可以用式0)作满意的GM(L1)建模二第二步:GM(1,1)建模(1)对原始数据祚一次累加，即王⑴=(71.1,143.5,215,9,288,359,4,431.4,503)(2)构造数据矩阵B及数据向量丫年求解得工切3=（工⑼0）--）e+£仪 仪mg+D= 当L+♦=—-T["31t+3IODO⑸求生成效列值 城喙县卷 ：产（f令k=1口340怎由上面的新瑜曲却可算得，其中期r由寸⑴（1）= 眈!3与乳5"冷春丁11：'涉®=迎®-通色-6（71.U7Z.417L2,72.L71J,71.7,71.6）逊=（£—，迎团…,＜第三步;模型检验＞模型的各种检验指标值的计算结果见表工：•表7 检验表19S671.171.1 0 01987 72.4 72.4 -0.0057 19S671.171.1 0 01987 72.4 72.4 -0.0057 0.01%19SS 72.4 72.2 0.163区 0.23%19S9 72.1 72.1 0.0329 0.05%1990 71.4 71.9 -0.4州40.7%1991 72.0 71.7 0.26990.37%1992 71.6 71.6 0.037S0.05%0.00230.0203-0.00180.00740.0107-0.00323经验证.该模型的精度较高.可进行f页测和预报.

计算的Matl疝程序如木厂…也仃顶测和预报。clc.clearx0-；711172,472,472,171,472n 丈i,.,n-lengthtxO)； , '「*注意这里为列向置lamda=xD(l：n-I),/<0(2：n"计算级比range-minmaxtlajnda*)%计算级比的范围xl=cumsum(x0)；%累加运算xl-0.5(xl(1 -1)+xl(2:n))tones(n^1,1)]Tu= %拟合参数=si,u(2)=bx=dsolve(,Dx+a'x=b\*x(0) %求微分方程的符号薜x=subs(xj Ju(1),u(2),xO(J)I)；%代入牯计参数值和初始位yuceL=subs(x/t\lOin-lJ);%求巳知数据的掂测便yvpa(xt6)%其中的6表示显不6位数字yuce= ,diff(yucel)]*差分运算，近原数据iIon~ —yuce%计算残差―彳普詈没：慧m’,计…….-ttlQ=1—(I=0.5U(I))/1+0"-a .k■■"inh0■STI^Zjprhulst作用：求累加数列、求ab的值、求预测方程、求残差clc%clc%清屏，以使结果独立显示x=[71.172.472.472.171.472.071.6];formatlong;formatlong;%设置计算精度iflength(x(:,1))==1 %对输入矩阵进行判断，如不是一维列矩阵，进行转置变换x=x';endn=length(x);z=0;%取输入数据的样本量fori=1:nz=z+x(i,:)be(i,:)=z;iflength(x(:,1))==1 %对输入矩阵进行判断，如不是一维列矩阵，进行转置变换x=x';endn=length(x);z=0;%取输入数据的样本量fori=1:nz=z+x(i,:)be(i,:)=z;%计算累加值，并将值赋予矩阵beendfori=2:n%对原始数列平行移位endfori=1:n-1%计算数据矩阵B的第一列数据c(i,:)=-0.5*(be(i,:)+be(i+1,:));endforj=1:n-1 %计算数据矩阵B的第二列数据e(j,:)=1;endfori=1:n-1 %构造数据矩阵BB(i,1)=c(i,:);B(i,2)=e(i,:);endalpha=inv(B'*B)*B'*y; %计算参数矩阵即ab的值fori=1:n+1 %计算数据估计值的累加数列，如改为n+1为n+m可预测后m-1个值ago(i,:)=(x(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1,:)*(i-1))+alpha(2,:)/alpha(1,:);%显示输出预测值的累加数列endvar(1,:)=ago(1,:) %显示输出预测值fori=1:n %如改n为n+m-1，可预测后m-1个值var(i+1,:)=ago(i+1,:)-ago(i,:); %估计值的累加数列的还原，并计算出下一预测值endfori=1:nerror(i,:)=x(i,:)-var(i,:); %计算残差endc=std(error)/std(x); %调用统计工具箱的标准差函数计算后验差的比值cago %显示输出预测值的累加数列alpha %显示输出参数数列var %显示输出预测值error %显示输出误差c %显示后验差的比值作用：数据处理判断是否可以用灰色预测、求级比、求累加数列、求ab的值、求预测方程clc,clearx0=[71.172.472.472.171.472.071.6],;%注意这里为列向量n=length(x0);lamda=x0(1:n-1)./x0(2:n) %计算级比range=minmax(lamda，)%计算级比的范围x1=cumsum(x0)%累加运算B=[-0.5*(x1(1:n-1)+x1(2:n)),ones(n-1,1)];Y=x0(2:n);u=B\Y%拟合参