本文代码主要是基于邓聚龙教授在20实际80年代提出的灰色系统理论。



GM0.m
%该函数为GM（1，1）模型返回还原值 
function f=GM0(x0,t)  %数据数列 
[M,N]=size(x0);        %算出数据数列的大小 
x1(1)=x0(1);           %累加生成数列 
for i=2:N;     
    x1(i)=x1(i-1)+x0(i); 
end
x2=[];              %累加生成数列均值生成数列 
for j=1:(N-1);     
    x2(j)=(x1(j)+x1(j+1))/2; 
end
x=x0;              %数据数列镜像
x(1)=[];           %删除第一个数据 
Y=x';              %数据列向量 
global a; 
global b; 
B(:,1)=-x2';
B(:,2)=1; 
A=inv(B'*B)*B'*Y;   %求参量a，b组成的参数向量 
a=A(1,1);           %求参数a 
b=A(2,1);           %求参数b  
f=(1-exp(a))*(x0(1)-b/a)*exp(-a*(t-1));
f

GM1.m:
%该函数为GM（1，1）模型中数据数列进行光滑比检验 
function f=GM1(x0)   %数据数列 
N=max(size(x0));       %算出数据数列的大小  
x1=cumsum(x0);         %累加生成数列 
global J; 
global J1;
global J2; 
x0(1)=[]; 
x1(N)=[]; 
global r;
r=x0./x1;  
for j=2:(N-1);           %判断数据数列是否满足准光滑条件1   
   if(r(j)>=0.5||r(j)<0)         
       J1=0;         
       break;     
   else
       J1=1;     
   end
end
for l=1:(N-2);           %判断数据数列是否满足准光滑条件2     
    if((r(l+1)/r(l))>=1)         
        J2=0;         
        break;     
    else
        J2=1;     
    end
end
J=J1+J2;  
if(J==2)                 %判断数据数列是否为准光滑数列     
    disp('数据为准光滑数列')
else
    disp('数据不是准光滑数列') 
end

GM2.m
%该函数为GM（1，1）模型还原值参数计算 
function f=GM2(x0)  %数据数列
[M,N]=size(x0);      %算出数据数列的大小  
x1(1)=x0(1);         %累加生成数列 
for i=2:N;      
    x1(i)=x1(i-1)+x0(i); 
end
x2=[];              %累加生成数列均值生成数列 
for j=1:(N-1);      
    x2(j)=(x1(j)+x1(j+1))/2;
end
  x=x0;              %数据数列镜像      
  x(1)=[];           %删除第一个数据 
  Y=x';              %数据列向量 
  global a; 
  global b; 
  B(:,1)=-x2'; 
  B(:,2)=1;  
  A=inv(B'*B)*B'*Y;   %求参量a，b组成的参数向量 
  a=A(1,1);           %求参数a 
  disp('参数a为：') 
  a 
  b=A(2,1);           %求参数b 
  disp('参数b为：') 
  b

GM3.m
%该程序实现G(1,1)模型的精度检验 
%包括平均相对误差，绝对关联度，均方差比值，小误差概率检验
function f=GM3(x0) 
N=max(size(x0)); 
x=GM0(x0,1:N);  %利用已有程序GM2得出数据列模型估计值 
x(1)=x0(1);     %更正第一个估计值 
disp('模型模拟估计值为') 
x 
A=x-x0;         %计算绝对残差序列 
disp('模型估计值绝对残差序列为：')
A 
G=abs(A); 
Amin=min(G);    %计算最小绝对值绝对残差 
Amax=max(G);    %计算最大绝对值绝对残差 
B=A./x0;        %计算相对误差序列 
disp('模型估计值相对误差序列为：') 
B  
P=sum(abs(B))/N;     %计算平均相对误差 
disp('模型估计值平均相对误差为：')
P  
for i=1:1:N       %通过循环计算关联系数序列     
    D(i)=(Amin+0.5*Amax)/(G(i)+0.5*Amax); 
end
R=sum(D)/N;  
disp('关联度为：') 
R  
x_=sum(x0)/N;    %计算数据的均值  
S1=(sum((x0-x_).^2)/(N-1))^0.5;   %计算数据序列方均差
A_=sum(A)/N;     %计算残差平均值  
S2=(sum((A-A_).^2)/(N-1))^0.5;    %计算残差序列方均差
C=S2/S1;         %计算方均差比值 
disp('均方差比值为：') 
C  
S0=0.6745*S1;
E=A-A_; 
F=find(E<S0); 
M=max(size(F)); %计算小残差个数
p=M/N;          %计算小误差概率 
disp('小误差概率为：')
p