灰色預測模型matlab程序精確版

1、%x=1019,1088,1324,1408,1601;gm1(x); 測試數(shù)據(jù)%1次擬合預測 GM(1,1)模型function gmcal=gm1(x)if nargin=0x=1019,1088,1324,1408,1601endformat long gsizex=length(x);%求數(shù)組長度k=0;for y1=xk=k+1;if k>1x1(k)=x1(k-1)+x(k);% 累加生成z1(k-1)=-0.5*(x1(k)+x1(k-1);%z1 維數(shù)減 1，用于計算Byn1(k-1)=x(k);elsex1(k)=x(k);endend%x1,z1,k,yn1sizez

2、1=length(z1);%size(yn1);z2 = z1'z3 = ones(1,sizez1)'YN = yn1' % 轉置%YNB=z2 z3;au0=inv(B'*B)*B'*YN;au = au0'%B,au0,auafor = au(1);ufor = au(2);ua = au(2)./au(1);%afor,ufor,ua%輸出預測的a u 和 u/a 的值constant1 = x(1)-ua;afor1 = -afor;x1t1 = 'x1(t+1)'estr = 'exp'tstr = &

3、#39;t'leftbra = '('rightbra =')'%constant1,afo1x1t1,estr,tstr,leftbra,rightbrastrcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,rightbra,'+',leftbra,num2str(ua),rightbra)%俞出時間響應方程%*%=次擬合k2 = 0;for y2 = x1k2 = k2 + 1;if k2 > kelseze1(k2) = exp(

4、-(k2-1)*afor);endend%ze1sizeze1=length(ze1);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,rightbra,'+',leftbr

5、a,num2str(Bval),rightbra)喻出時間響應方程nfinal = sizex-1 + 1; (其中+1可改為+5等其他數(shù)字，即可預測更多的數(shù)字)燭定預測的步驟數(shù) 5這個步驟可以通過函數(shù)傳入%nfinal = sizexd2 - 1 + 1;函測的步驟數(shù)1for k3=1:nfinalx3fcast(k3) = constant1*exp(afor1*k3)+ua;end%x3fcast獷次擬合累加值for k31=nfinal:-1:0if k31>1x31fcast(k31+1) = x3fcast(k31)-x3fcast(k31-1);elseif k31>

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

7、) = x(1);endendendx41fcast,x%二次擬合預測值%*精度檢驗 p C*/k5 = 0;for y5 = xk5 = k5 + 1;if k5 > sizexelseerr1(k5) = x(k5) - x41fcast(k5);endend%err1%絕對誤差xavg = mean(x);%xavg%x平均值err1avg = mean(err1);%err1avg%err1 平均值k5 = 0;s1total = 0 ;for y5 = xk5 = k5 + 1;if k5 > sizexelsesltotal = sltotal + (x(k5) - x

8、avg42;endends1suqare = s1total ./ sizex;s1sqrt = sqrt(s1suqare);%s1suqare,s1sqrt%s1suqare 殘差數(shù)列 x 的方差 s1sqrt 為 x 方差的平方根S1k5 = 0;s2total = 0 ;for y5 = xk5 = k5 + 1;if k5 > sizexelses2total = s2total + (err1(k5) - err1avg)A2;endends2suqare = s2total ./ sizex;%s2suqare 殘差數(shù)列 err1 的方差 S2Cval = sqrt(s2s

9、uqare ./ s1suqare);Cval%nnn = 0.6745 * s1sqrt%Cval C 檢驗值k5 = 0;pnum = 0 ;for y5 = xk5 = k5 + 1;if abs( err1(k5) - err1avg ) < 0.6745 * s1sqrtpnum = pnum + 1;%ppp = abs( err1(k5) - err1avg )elseendendpval = pnum ./ sizex;pval%p檢驗值%arr1 = x41fcast(1:6)%預測結果為區(qū)間范圍預測步長和數(shù)據(jù)長度可調(diào)整程序參數(shù)進行改進運行結果10191088132414081601ans = x1(t+1)=8908.4929exp(0.11871t)+(-7889.4929) ans = x1(t+1)=8945.2933exp(0.11871t)+(-7935.7685) x31fcast =Columns 1 through 310191122.893478570971264.43142178303Columns 4 through 61423.809872354881603.277582074421805.36675232556x41fcast =Columns 1 throug