2022年2022年灰色预测模型GM的matlab运行代码 .pdf

灰色预测模型GM(1,1)的 matlab 运行代码例 由 19902001 年中国蔬菜产量，建立模型预测2002年中国蔬菜产量，并对预测结果作检验。分析建模：给定原始时间19902001 年资料序列X)0(k),对 X)0(k)生成 1-AGO(累加)序列 X)1(k)及 Yn。见下表K 1 2 3 4 5 6 7 8 9 10 11 12X)0(19519,19578,19637,19695,16602,25723,30379,34473,38485,40514,42400,48337 X)1(19519,39097,58734,264605,307005,355342Yn-19578 19637 40514 42400 48337其中X)1(k)=)(1i)0(iXk;Yn=TXXX)12(,),3(),2()0()0()0(对上述 X)0(k)的 GM(1,1),得到15236.2-1331173.5-1244348.01204848.5-1168369.5-1135943.5-1107892.5-186730.0-168581.5-148915.5-129308.0-112115.0111105.011095.01985.01875.01765.01655.01545.01435.01325.01215.01)12(1)11(1)10(1)9(1)8(1)7(1)6(1)5(1)4(1)3(1)2()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1()1(）（）（）（）（）（）（）（）（）（）（）（）（）（）（）（）（）（）（）（）（）（）（XXXXXXXXXXXXXXXXXXXXXXzzzzzzzzzzzB将 B 和 Yn代入辨识算式，有：10.1062105()13999.9TTnaB BBYb得灰色 GM(1,1)模型为名师资料总结-精品资料欢迎下载-名师精心整理-第 1 页，共 7 页 -（1）灰微分方程X)0(k)-0.1062105 Z)1(k)=13999.9(2)白化方程9.139991062105.0)1()1(XdtdX(3)白化方程的时间响应式5.1318135.151332)1()1(?1062105.0)0()1(tateabeabXtX(4)还原为原始数据预测方程：)(?)1(?)1(?)1()1()0(tXtXtX,即tetX1062105.0)0(15248.968)1(?(5)残差检验：残差 error1=e1=)()?)0()0(iXiX（,这里残差有 12 个。相对残差 error2=e2=)0()0()0()0(1)(|)()(?|XeiXiXiX,这里相对残差有 12 个。(6)后验差检验：C=21SS,其 中1)(1212)0(01niSi）（,S1 为 绝 对 误 差 序 列 的 标 准 差。1)()(?)()0()0()0(eiXiXi,)(121121)0()0(iiS2为原始数据系列标准差，1)(1212)0(02nXiXSi）（,)(121121)0()0(iXXiC0.35 好；C0.6 不合格。利用 matlab 做求解 a,b,B,并作残差分析 x0=19519,19578,19637,19695,16602,25723,30379,34473,38485,40514,42400,48337;format long;(表示设计精度)n=length(x0);（输入数据长度）x1=;（表示 x1 是一矩阵）x1(1)=x0(1);for i=2:n;x1(i)=x1(i-1)+x0(i);end 名师资料总结-精品资料欢迎下载-名师精心整理-第 2 页，共 7 页 -for i=1:n-1;B(i,1)=-0.5*(x1(i)+x1(i+1);（矩阵 B 的第一列）B(i,2)=1;（矩阵 B 的第二列）Y(i)=x0(i+1);（表示 Yn 数据）end alpha=(B*B)(-1)*B*Y;a=alpha(1,1);b=alpha(2,1);d=b/a;（计算时间响应函数参数）c=x1(1)-d;x2(1)=x0(1);x(1)=x0(1);for i=1:n-1;x2(i+1)=c*exp(-a*i)+d;（这里 x2(i+1)相当上面所讲的)1(?)1(tX）x(i+1)=x2(i+1)-x2(i);（这里 x(i+1)相当原来输入数据的预测数据)1(?)0(tX）end for i=2:12;x2(i)=c*exp(-a*(i-1)+d;（对上面刚引出的x2(i)进行说明及计算）x(i)=x2(i)-x2(i-1);end for i=1:n;error(i)=x(i)-x0(i);（残差）error1(i)=abs(error(i);（计算残差，abs表示绝对值）error2(i)=error1(i)/x0(i);（计算相对误差）end C=std(error1)/std(x0);（计算后验差检验数,std 表示标准差）k=1;（k 表示预测长度，这里每次预测下一年）a=-0.106210475032772 b b=1.399996741173038e+04 B B=名师资料总结-精品资料欢迎下载-名师精心整理-第 3 页，共 7 页 -1.0e+05*-0.293080000000000 0.000010000000000-0.489155000000000 0.000010000000000-0.685815000000000 0.000010000000000-0.867300000000000 0.000010000000000-1.078925000000000 0.000010000000000-1.359435000000000 0.000010000000000-1.683695000000000 0.000010000000000-2.048485000000000 0.000010000000000-2.443480000000000 0.000010000000000-2.858050000000000 0.000010000000000-3.311735000000000 0.000010000000000 C（求后检验数）C=0.163969348419772 x（原始数据)()0(iX的对应的预测数据)(?)0(iX，这里也是12 个）x=1.0e+04*Columns 1 through 3 1.951900000000000 1.695769385830782 1.885790370699694 Columns 4 through 6 2.097104330304606 2.332097268346191 2.593422554345592 Columns 7 through 9 2.884030883565190 3.207203594115656 3.566589717441855 Columns 10 through 12 3.966247180534730 4.410688625094428 4.904932361001964 eroor1（求残差）eroor1=名师资料总结-精品资料欢迎下载-名师精心整理-第 4 页，共 7 页 -1.0e+03*Columns 1 through 4 0 2.620306141692185 0.779096293003065 1.276043303046055 Columns 5 through 8 6.718972683461907 0.211225543455919 1.538691164348100 2.400964058843441 Columns 9 through 12 2.819102825581445 0.851528194652696 1.706886250944284 0.712323610019637 eroor2（求相对误差）eroor2=Columns 1 through 4 0 0.133839316666267 0.039674914345525 0.064790215945471 Columns 5 through 8 0.404708630494031 0.008211543888968 0.050649829301429 0.069647667996503 Columns 9 through 12 0.073251989751369 0.021018121998635 0.040256751201516 0.014736611912606 a=-0.106210475032772 b=1.399996741173038e+04 B B=名师资料总结-精品资料欢迎下载-名师精心整理-第 5 页，共 7 页 -1.0e+05*-0.293080000000000 0.000010000000000-0.489155000000000 0.000010000000000-0.685815000000000 0.000010000000000-0.867300000000000 0.000010000000000-1.078925000000000 0.000010000000000-1.359435000000000 0.000010000000000-1.683695000000000 0.000010000000000-2.048485000000000 0.000010000000000-2.443480000000000 0.000010000000000-2.858050000000000 0.000010000000000-3.311735000000000 0.000010000000000 方法二程序（1）一次累加生成序列的matlab 命令 x0=19519,19578,19637,19695,16602,25723,30379,34473,38485,40514,42400,48337;x1(1)=x0(1);x1(1)x1(1)=19519 for t=2:12;x1(t)=x1(t-1)+x0(t);end x1 回车x1=Columns 1 through 8 19519 39097 58734 78429 95031 120754 151133 185606 Columns 9 through 12 224091 264605 307005 355342（2）由一次累加生成序列紧邻均值生成Z)1(的 matlab 命令：x0=19519,19578,19637,19695,16602,25723,30379,34473,38485,40514,42400,48337;x1(1)=x0(1);for t=2:12;x1(t)=x1(t-1)+x0(t);名师资料总结-精品资料欢迎下载-名师精心整理-第 6 页，共 7 页 -z1(t)=(1/2)*(x1(t)+x1(t-1);end z1=1.0e+05*Columns 1 through 7 0 0.2931 0.4892 0.6858 0.8673 1.0789 1.3594 Columns 8 through 12 1.6837 2.0485 2.4435 2.8581 3.3117(3)由于 GM(1,1)的灰微分方程为X)0(t)-+aZ)1(k)=b 设为待估参数，ba,利用最小二乘法得到YBBB1)（,matlab 程序为B=-z1(2:12),ones(11,1);Y=(x0(2:12);alpha=inv(B*B)*B*Y;alpha alpha=1.0e+04*-0.0000 1.4000 名师资料总结-精品资料欢迎下载-名师精心整理-第 7 页，共 7 页 -