基于MATLAB的灰色预测模型

现有给定数据（见表1），每一组编程的同学先熟悉灰色预测)1,1(GM的算法，根据算法过程（分步完成）编写相应的程序。要求：1.程序能够由编程者的需求不同预测出给定数据后任意n年的值；2.程序能够求出给定数据与预测数据的残差；3.程序能够画出初始数据与预测数据的对比图，注意图的格式，黑白；4.程序能够求出数据矩阵B及数据向量Y；5程序能够得出求预测值的方程即：cmekxnk)1(^其中，cnm、、均为需用程序求解出来的系数；6.程序中重要的程序语句用‘%’注明语义；表1给定数据年份2000200120022003200420052006数据71.172.472.472.171.472.071.6答：1:程序clc%灰色预测模型程序clearsymsab;c=[ab]';A=[71.172.472.472.171.472.071.6];%原始序列B=cumsum(A);%累加n=length(A);fori=1:(n-1)C(i)=(B(i)+B(i+1))/2;end%计算待定参数D=A;D(1)=[];D=D';E=[-C;ones(1,n-1)];c=inv(E*E')*E*D;c=c';a=c(1);b=c(2);%预测往后预测5个数据。。。。。。。此处可把5可更改为对应想得到的数据后n年的值，除了前7位数，后面的为数据后n年的值F=[];F(1)=A(1);fori=2:(n+5)F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a;endG=[];G(1)=A(1);fori=2:(n+5)G(i)=F(i)-F(i-1);endG2.3.4.5.6:程序:clcclearX0=[71.172.472.472.171.472.071.6];[m,n]=size(X0);X1=cumsum(X0);%累加X2=[];fori=1:n-1X2(i,:)=X1(i)+X1(i+1);endB=-0.5.*X2;t=ones(n-1,1);B=[B,t]%4.1.求B矩阵YN=X0(2:end)%4.2.求Y向量P_t=YN./X1(1:(length(X0)-1))%对原始数据序列X0进行准光滑性检验，%序列x0的光滑比P(t)=X0(t)/X1(t-1);A=inv(B.'*B)*B.'*YN.';a=A(1)u=A(2)c=u/a;b=X0(1)-c;X=[num2str(b),'exp','(',num2str(-a),'k',')','+',num2str(c)];strcat('X(k+1)=',X)%5.得出预测值的方程%symsk;fort=1:length(X0)k(1,t)=t-1;endkY_k_1=b*exp(-a*k)+c;forj=1:length(k)-1Y(1,j)=Y_k_1(j+1)-Y_k_1(j);endXY=[Y_k_1(1),Y]%预测值subplot(1,2,1),plot(X0,'-*k');title('原始数据');%3.绘制图像subplot(1,2,2),plot(XY,'--pk');title('预测数据');%绘制图像CA=abs(XY-X0);%残差数列Theta=CA%残差检验绝对误差序列XD_Theta=CA./X0%残差检验相对误差序列AV=mean(CA)%2.残差数列平均值R_k=(min(Theta)+0.5*max(Theta))./(Theta+0.5*max(Theta));%P=0.5R=sum(R_k)/length(R_k)%关联度Temp0=(CA-AV).^2;Temp1=sum(Temp0)/length(CA);S2=sqrt(Temp1);%绝对误差序列的标准差%----------AV_0=mean(X0)%原始序列平均值Temp_0=(X0-AV_0).^2;Temp_1=sum(Temp_0)/length(CA);S1=sqrt(Temp_1);%原始序列的标准差TempC=S2/S1*100;%方差比C=strcat(num2str(TempC),'%')%后验差检验%方差比%----------SS=0.675*S1Delta=abs(CA-AV)TempN=find(Delta=SS);N1=length(TempN);N2=length(CA);TempP=N1/N2*100;P=strcat(num2str(TempP),'%')%后验差检验%计算小误差概率