灰色预测代码.docxVIP

例一 function predate=greypred(ddata,lenn,stepnum) predata=[]; llen=length(ddata); sumdata(1:llen)=l; for k=1:llen sumdata(k)=sum(ddata(1:k)); end yn=ddata(2:llen); B(1:llen-l,1:2)=1; for k=l:lien-l B(k,1)=-0.5*(sumdata(k)+sumdata(k+1)); end coeff=inv(B*B)*B*yn ％用最小二乘拟合系数 for k=l:lenn+llen anser=(ddata(1)-coeff(2)/coeff(1))*exp(coeff(1)*(k-1))+coeff(2)/coeff(1); predata=[predata anser]; end predata(2:lenn+lien)=predata(2:lenn+lien)-predata(1:lenn+lien-1); step=1; X=2002:2009; Y=[2724.8 3126.1 3664.8 4193.4 4792.1]; len=8,num=5; ％已知 5 个原始数据 reg= Y(1:num); predlen=len-num; ％需要预测的数据个数 predy= greypred(reg,predlen,step); figure(1) hold on plot(x(1:5),y,-b); plot(x,predy,-ro); legend(真实值,预测值’),title(预测效果); xlabel(年),ylabel(收入增加值); hold off ep=Y-predy(1:num); ％残差 eeq=ep./y; ％相对残差 figure(2) plot(x(1:5),Y,-b.,x(1:5),predy(1:5),-ro); legend(真实值,预测值) title( 真实值与预测值的接近程度 ); xlabel(年),ylabel(收人增加值) figure(3) plot(x(1:5),ep,-b,x(1:5),eeq,-ro); legend(残差,相对误差); axis square 例2 y=[48.7 57.17 68.76 92.15] n=length(y); yy=ones(n,1); yy(1)=y(1); for i=2:n yy(i)=yy(i-1)+y(i); end B=ones(n-1,2); for i=1:(n-1) B(i,1)=-(yy(i)+yy(i+1))/2; B(i,2)=1; end BT=B; for j=1:n-1 YN(j)=y(j+1); end YN=YN; A=inv(BT*B)*BT*YN; a=A(1); u=A(2); t=u/a; t_test=10; i=1:t_test+n; yys(i+1)=(y(1)-t).*exp(-a.*i)+t; yys(1)=y(1); for j=n+t_test:-1:2 ys(j)=yys(j)-yys(j-1); end x=1:n; xs=2:n+t_test; yn=ys(2:n+t_test); plot(x,y,Ar,xs,y n, *-b); det=0; for i=2:n det=det+abs(yn(i)-y(i)); end det=det/(n-1); disp([百分绝对误差为：’,n um2str(det),%]); disp([ 预测值为： ,num2str(ys(n+1:n+t_test))]);