灰色模型Matlab程序

1、%by alien 红嘴海鸥%灰色模型预测是在数据不呈现一定规律下可以采取的一种建模和预测方法，其预测数据与原始数据存在一定的规律相似性%下面程序是灰色模型GM(1,1)程序二次拟合和等维新陈代谢改进预测程序，matlab6.5 ,使用本程序请注测试数据明，程序存储为gml.m%x =5999,5903,5848,5700,7884;gm1(x);%二次拟合预测GM(1,1)模型 function gmcal=gm1(x) sizexd2 = size(x,2);%求数组长度 k=0;for y1=x k=k+1;ifk1 x1(k)=x1(k-1)+x(k);%累加生成 z1(k-1)=-0

2、.5*(x1(k)+x1(k-1);%z1维数减1,用于计算B yn1(k-1)=x(k);else x1(k)=x(k);end end %x1,z1,k,yn1 sizez1=size(z1,2);%size(yn1);z2 = z1;z3 = ones(1,sizez1);YN = yn1; % 转置%YNB=z2 z3;auO=inv(B*B)*B*YN;au = auO;%B,auO,au afor = au(1);ufor = au(2);ua = au(2)./au(1);%afor,ufor,ua%输出预测的a u和u/a的值 constant1 = x(1)-ua;afor1

3、 = -afor;x1t1 = x1(t+1);estr = ex p;tstr = t;leftbra =(;rightbra =);%constant1,afor1,x1t1,estr,tstr,leftbra,nghtbra strcat(x1t1,=,num2str(constant1),estr,leftbra,num2str(afor1),tstr,nghtbra,+,leftbra,num2str(ua).rightbra)%输出时间响应方程%*%二次拟合k2 = 0;for y2 = x1 k2 = k2 +1;if k2 k else ze1(k2) =ex p(-(k2-1

4、)*afor);end end %ze1 sizezel = size(ze1,2);z4 = ones(1,sizeze1);G=ze1 z4;X1 = x1;au20=inv(G*G)*G*X1;au2 = au20;%z4,X1,G,au20Aval = au2(1);Bval = au2(2);%Aval,Bval%输出预测的A,B的值 strcat(x1t1,=,num2str(Aval),estr,leftbra,num2str(afor1),tstr,nghtbra,+,leftbra,num2str(Bval),righ tbra)%输出时间响应方程 nfinal = size

5、xd2-1 + 1;%决定预测的步骤数5这个步骤可以通过函数传入%nfinal = sizexd2 - 1 + 1;%预测的步骤数1for k3=1:nfinalx3fcast(k3)= constant1*ex p( afor1*k3)+ua;end%x3fcast%一次拟合累加值for k31=nfinal:-1:0ifk311x31fcast(k31+1) = x3fcast(k31)-x3fcast(k31-1);elseif k310x31fcast(k31+1) = x3fcast(k31)-x(1);elsex31fcast(k31+1) = x(1);endendendx31f

6、cast%一次拟合预测值for k4=1:nfinalx4fcast(k4)= Aval*ex p( afor1*k4)+Bval;end%x4fcastfor k41=nfinal:-1:0ifk411x41fcast(k41+1) = x4fcast(k41)-x4fcast(k41-1);elseif k410x41fcast(k41+1) = x4fcast(k41)-x(1);elsex41fcast(k41+1) = x(1);endendendx41fcast,x%二次拟合预测值%* 精度检验 pC*/k5 = 0;for y5 = xk5 = k5 +1;if k5 sizex

7、d2elseerr1(k5) = x(k5) -x41fcast(k5);endend%err1%绝对误差xavg = mean(x);%xavg%x平均值err1avg = mean(err1);%err1avg%err1平均值k5 = 0;sitotal = 0 ;for y5 = x k5 = k5 +1;if k5 sizexd2 else sitotal = sitotal + (x(k5) -xavg)2;end end s1suqare = s1total ./ sizexd2;s1sqrt = sqrt(s1suqare);%s1suqare,s1sqrt %s1suqare

8、残差数列x的方差s1sqrt为x方差的平方根 S1 k5 = 0;s2total = 0 ;for y5 = x k5 = k5 +1;if k5 sizexd2 else s2total = s2total + (err1(k5) -err1avg)2;end end s2suqare = s2total ./ sizexd2;%s2suqare 残差数列erri的方差S2Cval = sqrt(s2suqare ./ sisuqare);Cval %nnn = 0.6745 * sisqrt %Cval C检验值k5 = 0;pnum = 0 ;for y5 = x k5 = k5 +1;

9、if abs(err1(k5) - errlavg ) 0.6745 * sisqrt pnum = pnum + 1;%ppp = abs( err1(k5) - err1avg) else end end p val = pnum ./ sizexd2;pval%p检验值%arr1 = x41fcast(1:6)%预测结果为区间范围预测步长和数据长度可调整程序参数进行改进程序为原创，引用请注明? gg=x1,x2,x3,x4,x5,zeros(1,m);%矩阵 gg,其中Xl,x2,x3,x4,x5是已知数据,前面说到灰色理论只要四个数据以上就可以建模 ，所以这里用五个比较合理,zeros

10、(1,m)中的m是你所想得到的预测数据的个数，由前面用于预测的已知数据是五个，所以m的值最好也在五个左右.在F面的语句中，实际操作的时候不要忘了把 m替换成实际的数值.? result=zeros(1,m);? for j=1:m? go=gg(:,j：j+4)；? n=5;? xO=gO;? x1=zeros(n-1,1);? x1(1)=x0(1);? for i=2:n,x1(i)=x1(i-1)+x0(i);? end? for k=1:n-1z1(k)=(x1(k)+x1(k+1)/2;? end? z=zeros(1,n-1);? for i=1: n-1z(i)=-z1(i);? end? B=z;ones(1,n-1);? B=B;? y=zeros(n-1,1);? for i=2:ny(i-i)=xO(i)；? end? u=inv(B*B)*B*y;? a=u(1);? b=u (2);? for k=0:5,? x2(k+