灰色預(yù)測模型matlab程序

1、灰色預(yù)測模型matlab程序%by alle n 紅嘴海鷗液色模型預(yù)測是在數(shù)據(jù)不呈現(xiàn)一定規(guī)律下可以采取的一種建模和預(yù)測方法，其預(yù)測數(shù)據(jù)與 原始數(shù)據(jù)存在一定的規(guī)律相似性蝦面程序是灰色模型GM(1,1)程序二次擬合和等維新陳代謝改進(jìn)預(yù)測程序，matlab6.5 ,用本程序請注明，程序存儲為gml.m%x = 5999,5903,5848,5700,7884;gm1(x);測試數(shù)據(jù)爐次擬合預(yù)測 GM(1,1)模型fun ctio n gmcal=gm1(x)sizexd2 = size(x,2);%求數(shù)組長度k=0;for y1=xk=k+1;if k>1x1(k)=x1(k_1)+x(k);

2、% 累加生成z1(k-1)=-0.5*(x1(k)+x1(k-1);%z1 維數(shù)減1,用于計(jì)算Byn 1(k-1)=x(k);elsex1(k)=x(k);endend%x1,z1,k, yn1sizez1=size(z1,2);%size(y n1);z2 = z1'z3 = on es(1,sizez1)'YN = yn 1' % 轉(zhuǎn)置%YNB=z2 z3;au0=i nv(B'*B)*B'*YN;au = au0'%B,au0,auafor = au(1);ufor = au(2);ua = au(2)./au(1);%afor,ufor,

3、ua%俞出預(yù)測的 a u 和u/a的值con sta nt1 = x(1)-ua;aforl = -afor;xltl = 'x1(t+1)'estr = 'exp'tstr = 't'leftbra = '('rightbra = ')' %constant1,afor1,x1t1,estr,tstr,leftbra,rightbra strcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,rightbra,'

4、+ ',leftbra,num2str(ua),rightbra)%輸出時(shí)間響應(yīng)方程%*%二次擬合k2 = 0;for y2 = x1k2 = k2 + 1;if k2 > kelseze1(k2) = exp(-(k2-1)*afor);endend%ze1sizeze1 = 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);%Ava

5、l,Bval%輸出預(yù)測的 A,B 的值 strcat(x1t1,'=',num2str(Aval),estr,leftbra,num2str(afor1),tstr,rightbra,'+',lef tbra,num2str(Bval),rightbra)%輸出時(shí)間響應(yīng)方程nfinal = sizexd2-1 + 1;%決定預(yù)測的步驟數(shù) 5 這個(gè)步驟可以通過函數(shù)傳入%nfinal = sizexd2 - 1 + 1;%預(yù)測的步驟數(shù) 1for k3=1:nfinalx3fcast(k3) = constant1*exp(afor1*k3)+ua;end%x3fca

6、st%一次擬合累加值for k31=nfinal:-1:0if k31>1x31fcast(k31+1) = x3fcast(k31)-x3fcast(k31-1); elseif k31>0x31fcast(k31+1) = x3fcast(k31)-x(1);elsex31fcast(k31+1) = x(1);endendend x31fcast %一次擬合預(yù)測值for k4=1:nfinalx4fcast(k4) = Aval*exp(afor1*k4)+Bval; end%x4fcastfor k41=nfinal:-1:0if k41>1x41fcast(k41+

7、1) = x4fcast(k41)-x4fcast(k41-1); elseif k41>0x41fcast(k41+1) = x4fcast(k41)-x(1);elsex41fcast(k41+1) = x(1);endend endx41fcast,x%二次擬合預(yù)測值%*精度檢驗(yàn) p C*/ k5 = 0;for y5 = xk5 = k5 + 1;if k5 > sizexd2elseerr1(k5) = x(k5) - x41fcast(k5);endend%err1%絕對誤差 xavg = mean(x);%xavg%x平均值err1avg = mean(err1);%

8、err1avg%err1 平均值k5 = 0; s1total = 0 ;for y5 = xk5 = k5 + 1;if k5 > sizexd2elsesltotal = sltotal + (x(k5) - xavg)A2;endends1suqare = s1total ./ sizexd2;s1sqrt = sqrt(s1suqare);%s1suqare,s1sqrt%s1suqare 殘差數(shù)列 x 的方差 s1sqrt 為 x 方差的平方根 S1 k5 = 0;s2total = 0 ;for y5 = xk5 = k5 + 1;if k5 > sizexd2elses2total = s2total + (err1(k5) - err1avg)A2;endends2suqare = s2total ./ sizexd2;%s2suqare 殘差數(shù)列 err1 的方差 S2Cval = sqrt(s2suqare ./ s1suqare);Cval%nnn = 0.6745 * s1sqrt%Cval C 檢驗(yàn)值k5 = 0;pnum = 0 ;for y5 = xk5 = k5 + 1;if abs( err1(k5) - err1avg ) < 0.6745 * s1sqrtp