灰色预测模型matlab程序精确版(6页).doc
-灰色预测模型matlab程序精确版-第 6 页%x=1019,1088,1324,1408,1601;gm1(x); 测试数据 %二次拟合预测GM(1,1)模型function gmcal=gm1(x)if nargin=0 x=1019,1088,1324,1408,1601endformat long gsizex=length(x);%求数组长度k=0;for y1=x k=k+1; if k>1 x1(k)=x1(k-1)+x(k); %累加生成 z1(k-1)=-0.5*(x1(k)+x1(k-1); %z1维数减1，用于计算B yn1(k-1)=x(k); else x1(k)=x(k); endend%x1,z1,k,yn1sizez1=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 = 't'leftbra = '('rightbra = ')'%constant1,afor1,x1t1,estr,tstr,leftbra,rightbrastrcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,rightbra,'+',leftbra,num2str(ua),rightbra)%输出时间响应方程%二次拟合k2 = 0;for y2 = x1 k2 = k2 + 1; if k2 > k else ze1(k2) = exp(-(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,'+',leftbra,num2str(Bval),rightbra)%输出时间响应方程nfinal = sizex-1 + 1;（其中+1可改为+5等其他数字，即可预测更多的数字）%决定预测的步骤数5 这个步骤可以通过函数传入%nfinal = sizexd2 - 1 + 1;%预测的步骤数 1for k3=1:nfinal x3fcast(k3) = constant1*exp(afor1*k3)+ua;end%x3fcast%一次拟合累加值for k31=nfinal:-1:0 if k31>1 x31fcast(k31+1) = x3fcast(k31)-x3fcast(k31-1); else if k31>0 x31fcast(k31+1) = x3fcast(k31)-x(1); else x31fcast(k31+1) = x(1); end endendx31fcast%一次拟合预测值for k4=1:nfinal x4fcast(k4) = Aval*exp(afor1*k4)+Bval;end%x4fcastfor k41=nfinal:-1:0 if k41>1 x41fcast(k41+1) = x4fcast(k41)-x4fcast(k41-1); else if k41>0 x41fcast(k41+1) = x4fcast(k41)-x(1); else x41fcast(k41+1) = x(1); end endendx41fcast,x%二次拟合预测值%*精度检验p C*/k5 = 0;for y5 = x k5 = k5 + 1; if k5 > sizex else err1(k5) = x(k5) - x41fcast(k5); endend%err1%绝对误差xavg = mean(x);%xavg%x平均值err1avg = mean(err1);%err1avg%err1平均值k5 = 0;s1total = 0 ;for y5 = x k5 = k5 + 1; if k5 > sizex else s1total = s1total + (x(k5) - xavg)2; endends1suqare = s1total ./ sizex;s1sqrt = sqrt(s1suqare);%s1suqare,s1sqrt%s1suqare 残差数列x的方差 s1sqrt 为x方差的平方根S1k5 = 0;s2total = 0 ;for y5 = x k5 = k5 + 1; if k5 > sizex else s2total = s2total + (err1(k5) - err1avg)2; endends2suqare = s2total ./ sizex;%s2suqare 残差数列err1的方差S2Cval = sqrt(s2suqare ./ s1suqare);Cval%nnn = 0.6745 * s1sqrt%Cval C检验值k5 = 0;pnum = 0 ;for y5 = x k5 = k5 + 1; if abs( err1(k5) - err1avg ) < 0.6745 * s1sqrt pnum = pnum + 1; %ppp = abs( err1(k5) - err1avg ) else endendpval = pnum ./ sizex;pval%p检验值%arr1 = x41fcast(1:6)%预测结果为区间范围 预测步长和数据长度可调整程序参数进行改进运行结果x = 1019 1088 1324 1408 1601ans =x1(t+1)=8908.4929exp(0.11871t)+(-7889.4929)ans =x1(t+1)=8945.2933exp(0.11871t)+(-7935.7685)x31fcast = Columns 1 through 3 1019 1122.89347857097 1264.43142178303 Columns 4 through 6 1423.80987235488 1603.27758207442 1805.36675232556x41fcast = Columns 1 through 3 1019 1118.05685435129 1269.65470492098 Columns 4 through 6 1429.69153740195 1609.90061644041 1812.82460377782x = 1019 1088 1324 1408 1601Cval = 0.139501578334155pval = 1
