% GM(1,1)模型
x0=load('data.txt');% 載入原始數據
n=length(x0);
% 一次累加1-GAO

x1(1)=x0(1);
for ii=1:n-1
    x1(ii+1)=x1(ii)+x0(ii+1);
end

% 求Y 陣
y=ones(n-1,1);
for ii=2:n
    y(ii-1)=x0(ii);
end

% 求B陣

b=ones(n-1,2);
for ii=1:n-1
    b(ii,1)=-0.5*(x1(ii)+x1(ii+1));
end

% 求參數a,u
w=inv(b'*b)*b'*y;
a=w(1,1)
u=w(2,1)
% 計算模型x11
m=u/a;
g=x0(1)-m;
x11=ones(n,1);
x11(1)=x0(1);
for ii=1:n-1
    x11(ii+1)=exp(-a*ii)*g+m;
end
%還原計算
x00=ones(n,1);
x00(1)=x0(1);
for ii=1:n-1
    x00(ii+1)=x11(ii+1)-x11(ii);
end
%模型檢驗
e=abs(x00-x0) %絕對誤差
e1=(e./x0)*100 %相對誤差
%關聯度檢驗
mine=min(e);
maxe=max(e);
for jj=1:n
    gg(jj)=(mine+0.5*maxe)/(e(jj)+0.5*maxe);
end
r=mean(g)
%后驗差檢驗
s1=std(x0);%原始序列標準差
s2=std(e);%絕對誤差序列標準差
c=s2/s1
me=ones(n,1)*mean(e);
ex=abs(e-me);
h=0.674*s1;
k=0;
for kk=1:n
    if ex(kk)<h
        k=k+1;
    end
end
p=k/n
%模型預測
ycn=6;%預測時段長度
yx=zeros(ycn,1);
for dd=1:ycn
    yx(dd)=exp((-a*(dd+n-1)))*(1-exp(a))*g;
end