原文地址:灰色预测模型matlab程序作者:红嘴海鸥%by allen @ 红嘴海鸥 %灰色模型预测是在数据不呈现一定规律下可以采取的一种建模和预测方法,其预测数据与原始数据存在一定的规律相似性 %下面程序是灰色模型gm(1,1)程序二次拟合和等维新陈代谢改进预测程序,matlab6.5 ,使用本程序请注明,程序存储为gm1.m %x = [5999,5903,5848,5700,7884];gm1(x); 测试数据 %二次拟合预测gm(1,1)模型 f nction gmcal=gm1(x) sizexd2 = size(x,2); %求数组长度 k=0; for y1=x k=k+1; if k1 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); end end %x1,z1,k,yn1 sizez1=size(z1,2); %size(yn1); z2 = z1 ; z3 = ones(1,sizez1) ; yn = yn1 ; %转置 %yn b=[z2 z3]; a 0=inv(b *b)*b *yn; a = a 0 ; %b,a 0,a afor = a (1); for = a (2); a = a (2).a (1); %afor, for, a %输出预测的 a 和 a的值 constant1 = x(1)- a; afor1 = -afor; x1t1 = x1(t+1) ; estr = exp ; tstr = t ; leftbra = ( ; rightbra = ) ; %constant1,afor1,x1t1,estr,tstr,leftbra,rightbra strcat(x1t1, = ,n m2str(constant1),estr,leftbra,n m2str(afor1),tstr,rightbra, + ,leftbra,n m2str( a),rightbra) %输出时间响应方程 % %二次拟合 k2 = 0; for y2 = x1 k2 = k2 + 1; if k2  k else ze1(k2) = exp(-(k2-1)*afor); end end %ze1 sizeze1 = size(ze1,2); z4 = ones(1,sizeze1) ; g=[ze1 z4]; x1 = x1 ; a 20=inv(g *g)*g *x1; a 2 = a 20 ; %z4,x1,g,a 20 aval = a 2(1); bval = a 2(2); %aval,bval %输出预测的 a,b的值 strcat(x1t1, = ,n m2str(aval),estr,leftbra,n m2str(afor1),tstr,rightbra, + ,leftbra,n m2str(bval),rightbra) %输出时间响应方程 nfinal = sizexd2-1 + 1; %决定预测的步骤数5 这个步骤可以通过函数传入 %nfinal = sizexd2 - 1 + 1; %预测的步骤数 1 for k3=1:nfinal x3fcast(k3) = constant1*exp(afor1*k3)+ a; end %x3fcast %一次拟合累加值 for k31=nfinal:-1:0 if k311 x31fcast(k31+1) = x3fcast(k31)-x3fcast(k31-1); else if k310 x31fcast(k31+1) = x3fcast(k31)-x(1); else x31fcast(k31+1) = x(1); end end end x31fcast %一次拟合预测值 for k4=1:nfinal x4fcast(k4) = aval*exp(afor1*k4)+bval; end %x4fcast for k41=nfinal:-1:0 if k411 x41fcast(k41+1) = x4fcast(k41)-x4fcast(k41-1); else if k410 x41fcast(k41+1) = x4fcast(k41)-x(1); else x41fcast(k41+1) = x(1); end end end x41fcast,x %二次拟合预测值 %精度检验p c k5 = 0; for y5 = x k5 = k5 + 1; if k5  sizexd2 else err1(k5) = x(k5) - x41fcast(k5); end end %err1 %绝对误差 xavg = mean(x); %xavg %x平均值 err1avg = mean(err1); %err1avg %err1平均值 k5 = 0; s1total = 0 ; for y5 = x k5 = k5 + 1; if k5  sizexd2 else s1total = s1total + (x(k5) - xavg)^2; end end s1s qare = s1total . sizexd2; s1sqrt = sqrt(s1s qare); %s1s qare,s1sqrt %s1s qare 残差数列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 s2s qare = s2total . sizexd2; %s2s qare 残差数列err1的方差s2 cval = sqrt(s2s qare . s1s qare); cval %nnn = 0.6745 * s1sqrt %cval c检验值 k5 = 0; pn m = 0 ; for y5 = x k5 = k5 + 1; if abs( err1(k5) - err1avg )  0.6745 * s1sqrt pn m = pn m + 1; %ppp = abs( err1(k5) - err1avg ) else end end pval = pn m . sizexd2;