gm（1，1）

% x0=[3 5 1 2 7];
% n=length(x0);
% x1=cumsum(x0);
% t1=x1;
% for k=2:n
%     t1(k)=(x1(k)+x1(k-1))/2;%z1是x1的上午紧邻均值生成数列
% end
% for k=2:n
%     y(k-1)=x0(k);
% end
% for k=2:n
%     z1(k-1)=t1(k);
% end
% z1=z1';
%  fai=[-z1 ones(n-1,1)];
% theta=inv(fai'*fai)*fai*y';
%  a=theta(1);
%  b=theta(2);
% for k=0:n-1
%      x1(k+1)=(x1(1)-b/a)*exp(-a*k-1)+b/a;
%  end
%  for k=1:n-1
%      x0(k)=x1(k+1)-x1(k);
%  end
clc
clear
x0=[2 3 3 4]';
n=length(x0);
x1=cumsum(x0);
z1=zeros(n-1,1);                                %给z1预先分配空间，然后计算均值生成列
for i=1:n-1
    z1(i)=1/2*(x1(i+1)+x1(i));
end
Y=x0(2:n);
phi=[-z1(1:3) ones(3,1)];                       %计算Φ
theta=(phi'*phi)^(-1)*phi'*Y;                   %计算最小二乘参数
c=theta(2)/theta(1);           
for k=0:5                                       %计算差分后的时间响应序列
    xx1(k+1)=[x1(1)-c]*exp(-theta(1)*k)+c;
end
xx0=x0(1);                                      %还原x0的值
for k=1:5 
xx0(k+1)= xx1(k+1)-xx1(k);
end
x=2006:2011;                                    %画图对比，2009年以后的是预测值
plot(x(1:4),x0,'*',x,xx0,'ro')
disp('均方误差为：')
MSE=1/n*sum((x0'-xx0(1:4)).^2)                  %方差误差
% RMSE=sqrt(MSE);                               %均方根误差
s1=var(x0);                                     %x0的方差=  5
s2=var(x0'-xx0(1:n));                           %残差的方差
s0=sum(x0(2:n-1)-x0(1))+1/2*(x0(n)-x0(1));
ss0=sum(xx0(2:n-1)-x0(1))+1/2*(xx0(n)-x0(1));
disp('绝对关联度为：')
epsilon=(1+abs(s0)+abs(ss0))/(1+abs(s0)+abs(ss0)-abs(s0-ss0))  %绝对关联度
if epsilon>0.9
    disp('预测精度等级：好')
   elseif  epsilon>0.8
       disp('绝对关联度为：合格')
         elseif epsilon>0.7
             disp('绝对关联度为：勉强')
             else
                 disp('绝对关联度为：不合格')
end