Matlab灰色预测模型GM(1-1)代码

functionc7fun73X0=[2.8743.2783.3073.393.679];AU=c7fun73(X0);a=AU(1);u=AU(2);m2=length(X0);fork=1:1:m2-1xx1(k+1)=(X0(1)-u/a)*exp(-a*k)+u/a;ends=0;xx0(1)=X0(1);forjj=2:1:m2;xx0(jj)=xx1(jj)-xx1(jj-1);enddisp('GM(1,1)对数列进行预测结果');xx0disp('数列1原始观测数据');X0disp('a');AU(1)disp('u');AU(2)functionau=c7fun73(X0)m=length(X0);s1=0;forjj=1:1:m;X1(jj)=s1+X0(jj);s1=X1(jj);endforii=1:1:m-1;B(ii)=-(X1(ii)+X1(ii+1))/2;endB=[B(:),ones(m-1,1)];y=X0([2:m])';au=inv((B'*B))*B'*y;%用MATLAB的灰色预测GM（1，1）模型%%程序中的变量定义；alpha是包含值的矩阵；ago是预测后累加值矩阵；var是预测值矩阵；error是残差矩阵；c是后验差比值functionmy_gm_test()clcdata=exprnd(5,1,10)%原始10个数据[agoalphavarerrorc]=gm1(data)plot(data)holdonplot(var','-r*')function[agoalphavarerrorc]=gm1(x);%定义函数gm1(x)%清屏，以使结果独立显示formatlong;%设置计算精度iflength(x(:,1))==1%对输入矩阵进行判断，如不是一维列矩阵，进行转置变换x=x';endn=length(x);%取输入数据的样本量z=0;fori=1:n%计算累加值，并将值赋予矩阵bez=z+x(i,:);be(i,:)=z;endfori=2:n%对原始数列平行移位y(i-1,:)=x(i,:);endfori=1:n-1%计算数据矩阵B的第一列数据c(i,:)=-0.5*(be(i,:)+be(i+1,:));endforj=1:n-1%计算数据矩阵B的第二列数据e(j,:)=1;endfori=1:n-1%构造数据矩阵BB(i,1)=c(i,:);B(i,2)=e(i,:);endalpha=inv(B'*B)*B'*y;%计算参数矩阵fori=1:n+1%计算数据估计值的累加数列，如改为n+1为n+m可预测后m-1个值ago(i,:)=(x(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1,:)*(i-1))+alpha(2,:)/alpha(1,:);endvar(1,:)=ago(1,:)fori=1:n%如改n为n+m-1，可预测后m-1个值var(i+1,:)=ago(i+1,:)-ago(i,:);%估计值的累加数列的还原，并计算出下一预测值endfori=1:nerror(i,:)=var(i,:)-x(i,:);%计算残差endc=std(error)/std(x);%调用统计工具箱的标准差函数计算后验差的比值cago%显示输出预测值的累加数列alpha%显示输出参数数列var%显示输出预测值error%显示输出误差c%显示后验差的比值cfunction[]=greymodel(y)%本程序主要用来计算根据灰色理论建立的模型的预测值。%应用的数学模型是GM(1,1)。%原始数据的处理方法是一次累加法。y=input('请输入数据');n=length(y);yy=ones(n,1);yy(1)=y(1);fori=2:nyy(i)=yy(i-1)+y(i);endB=ones(n-1,2);fori=1:(n-1)B(i,1)=-(yy(i)+yy(i+1))/2;B(i,2)=1;endBT=B';forj=1:n-1YN(j)=y(j+1);endYN=YN';A=inv(BT*B)*BT*YN;a=A(1);u=A(2);t=u/a;i=1:n+2;yys(i+1)=(y(1)-t).*exp(-a.*i)+t;yys(1)=y(1);forj=n+2:-1:2ys(j)=yys(j)-yys(j-1);endx=1:n;xs=2:n+2;yn=ys(2:n+2);plot(x,y,'^r',xs,yn,'*-b');det=0;fori=2:ndet=det+abs(yn(i)-y(i));enddet=det/(n-1);disp(['百分绝对误差为：',num2str(det),'%']);disp(['下个拟合值为',num2str(ys(n+1))]);disp(['再下个拟合值为',num2str(ys(n+2))]);x(0)={x(0)(1),x(0)(2),……,x(0)(N)}={724.57,746.62,778.27,800.8,827.75,871.1,912.37,954.28,995.01,1037.2}