数学建模之灰色预测模型代码（matlab版）

一、灰色理论中GM（1，1）模型

function []=greymodel(y)
% 本程序主要用来计算根据灰色理论建立的模型的预测值。
% 应用的数学模型是 GM(1,1)。
% 原始数据的处理方法是一次累加法。
y=input('请输入数据 ');
n=length(y);
yy=ones(n,1);
yy(1)=y(1);
for i=2:n
    yy(i)=yy(i-1)+y(i);
end
B=ones(n-1,2);
for i=1:(n-1)
    B(i,1)=-(yy(i)+yy(i+1))/2;
    B(i,2)=1;
end
BT=B';
for j=1:n-1
    YN(j)=y(j+1);
end
YN=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);
for j=n+2:-1:2
    ys(j)=yys(j)-yys(j-1);
end
x=1:n;
xs=2:n+2;
yn=ys(2:n+2);
plot(x,y,'^r',xs,yn,'*-b');
det=0;

sum1=0;
sumpe=0;
for i=1:n
    sumpe=sumpe+y(i);
end
pe=sumpe/n;
for i=1:n;
    sum1=sum1+(y(i)-pe).^2;
end
s1=sqrt(sum1/n);
sumce=0;
for i=2:n
    sumce=sumce+(y(i)-yn(i));
end
ce=sumce/(n-1);
sum2=0;
for i=2:n;
    sum2=sum2+(y(i)-yn(i)-ce).^2;
end
s2=sqrt(sum2/(n-1));
c=(s2)/(s1);
disp(['后验差比值为：',num2str(c)]);
if c<0.35
    disp('系统预测精度好')
else if c<0.5
        disp('系统预测精度合格')
    else if c<0.65
            disp('系统预测精度勉强')
        else
            disp('系统预测精度不合格')
        end
    end
end
            
disp(['下个拟合值为 ',num2str(ys(n+1))]);
disp(['再下个拟合值为',num2str(ys(n+2))]);

二、灰色关联分析

灰色关联分析步骤
【1】确定比较对象（评价对象）(就是数据，并且需要进行规范化处理，就是标准化处理，见下面例题的表格数据)和参考数列（评价标准，一般该列数列都是1，就是最优的的情况）
【2】确定各个指标权重，可用层次分析确定
【3】计算灰色关联系数
【4】计算灰色加权关联度
【5】评价分析

x1=[1.14 1.49 1.69 2.12 2.43 4.32 5.92 6.07 7.85;3.30 3.47 3.61 3.80 4.00 4.19 4.42 4.61 4.80;6.00 6.00 6.00 7.50 7.50 7.50 9.00 9.00 9.00;1.20 1.20 1.80 1.80 1.80 2.40 2.70 3.60 4.00;4.87 5.89 6.76 7.97 8.84 10.05 11.31 12.25 11.64]%原始数据5行9列
x=x1;
for i=1:5
    for j=1:9
  x(i,j)=x(i,j)/x1(1,j)
 end
end
x1=x
for i=1:5
    for j=1:9
  x(i,j)=abs(x(i,j)-x1(i,1))
 end
end
max=x(1,1)
min=x(1,1)
for i=1:5
    for j=1:9
 if x(i,j)>=max
     max=x(i,j)
 end
  end
end
for i=1:5
    for j=1:9
 if x(i,j)<=min
     min=x(i,j)
  end
  end
end
k=0.5 %分辨系数取值
l=(min+k*max)./(x+k*max)%求关联系数矩阵
guanliandu=sum(l')/n
[rs,rind]=sort(guanliandu,'descend') %对关联度进行排序

 三、灰色预测模型

灰色预测步骤
（1）输入前期的小样本数据
（2）输入预测个数
（3）运行

y=input('请输入数据');
n=length(y);
yy=ones(n,1);
yy(1)=y(1);
for i=2:n
    yy(i)=yy(i-1)+y(i)
end
B=ones(n-1,2);
for i=1:(n-1)
    B(i,1)=-(yy(i)+yy(i+1))/2;
    B(i,2)=1;
end
BT=B';
for j=1:(n-1)
    YN(j)=y(j+1);
end
YN=YN';
A=inv(BT*B)*BT*YN;
a=A(1);
u=A(2);
t=u/a;
t_test=input('输入需要预测的个数');
i=1:t_test+n;
yys(i+1)=(y(1)-t).*exp(-a.*i)+t;
yys(1)=y(1);
for j=n+t_test:-1:2
    ys(j)=yys(j)-yys(j-1);
end
x=1:n;
xs=2:n+t_test;
yn=ys(2:n+t_test);
plot(x,y,'^r',xs,yn,'*-b');
det=0;
for i=2:n
    det=det+abs(yn(i)-y(i));
end
det=det/(n-1);
disp(['百分绝对误差为：',num2str(det),'%']);
    disp(['预测值为：',num2str(ys(n+1:n+t_test))]);

 第一个代码可以检验是否预测良好，最后一个可以输入预测多少个