这是我曾经写过的一个灰色预测的程序：第一个文件为函数，需要在调用时输入原始数据x0和预测周期T, 第二个文件用于计算灰色关联度，使用时直接修改相关参数和原始数据。

--------------------------------------------------------------------------

第一个文件(用于灰色建模)：grymdl.m

--------------------------------------------------------------------------function GM=grymdl(x0,T)

% 输入原始数据x0

% T为从最后一个历史数据算起的第T时点

x1=zeros(1,length(x0));B=zeros(length(x0)-1,2);

yn=zeros(length(x0)-1,1);Hatx0=zeros(1,length(x0)+T);

Hatx00=zeros(1,length(x0));Hatx1=zeros(1,length(x0)+T);

epsilon=zeros(length(x0),1);omega=zeros(length(x0),1);

for i=1:length(x0)

for j=1:i

x1(i)=x1(i)+x0(j);

end

end

for i=1:length(x0)-1

B(i,1)=(-1/2)*(x1(i)+x1(i+1));

B(i,2)=1;

yn(i)=x0(i+1);

end

HatA=(inv(B'*B))*B'*yn; %GM(1,1)模型参数估计

for k=1:length(x0)+T

Hatx1(k)=(x0(1)-HatA(2)/HatA(1))*exp(-HatA(1)*(k-1))+HatA(2)/HatA(1);

end

Hatx0(1)=Hatx1(1);

for k=2:length(x0)+T

Hatx0(k)=Hatx1(k)-Hatx1(k-1);%累计还原得到历史数据的模拟值

end

for i=1:length(x0) %开始模型检验

epsilon(i)=x0(i)-Hatx0(i);

omega(i)=(epsilon(i)/x0(i))*100;

end

x0;

HatA;

Hatx0;

epsilon;

omega;

c=std(epsilon)/std(x0);

p=0;

for i=1:length(x0)

if abs(epsilon(i)-mean(epsilon))<0.6745*std(x0)

p=p+1;

end

end

p=p/length(x0);

if p>=0.95

M1=1;

elseif p>=0.8

M1=2;

elseif p>=0.7

M1=3;

else

M1=4;

end

if c<=0.35

M2=1;

elseif c<=0.5

M2=2;

elseif c<=0.65

M2=3;

else

M2=4;

end

M=max(M1,M2);

if M==1

disp('The model is good,and the forecast is:'),

disp(Hatx0(length(x0)+T))

elseif M==2

disp('The model is eligibility,and the forecast is:'),

disp(Hatx0(length(x0)+T))

elseif M==3

disp('The model is not good,and the forecast is:'),

disp(Hatx0(length(x0)+T))

else

disp('The model is bad and try again')

disp(Hatx0(length(x0)+T))

end

for i=1:length(x0)

Hatx00(i)=Hatx0(i);

end

z=1:length(x0);

plot(z,x0,'-',z,Hatx00,'*:') %将原始数据和模拟值画在一个图上帮助观察

text(2,x0(2),'History data: real line')

text(length(x0)/2,Hatx00(length(x0))/2,'Simulation data:broken line')

GM=Hatx0(length(x0)+T);

--------------------------------------------------------------------------------

第二个文件(用于计算灰色关联度)：grydgr.m

--------------------------------------------------------------------------------x=[26.4,40.83,5,18.5,27,4.9;

32.56,50.36,6.5,23.8,34.9,5.97;

41.16,63.66,7.9,30.6,45,6.83;

51.82,80.15,9.6,39.43,60.76,8.21;

70.58,109.16,11.3,52.1,81.8,9.78;

111.47,136.12,15,73.92,118.12,12.3;

127.16,166.7,18,118.17,184.51,16.75;

165.99,196.91,20.66,156.65,240.83,21.49;

246.76,257.56,24.15,218.3,333.4,27.1]; %原始数据

delta=zeros(size(x,1),size(x,2)-1); %初始化绝对差

yita=zeros(size(delta,1),size(delta,2)); %初始化关联系数

for i=1:size(x,2)

x(:,i)=x(:,i)./x(1,i);  %无量纲化处理

end

for i=1:(size(x,2)-1)

delta(:,i)=abs(x(:,1)-x(:,i+1));  %求解delta

end

delta_min=min(min(delta,[],1)); %求解最小二级差

delta_max=max(max(delta,[],1)); %求解最大二级差

rou=0.5; %设定分辨系数为0.5

for i=1:size(delta,1)

for j=1:size(delta,2)

yita(i,j)=(delta_min+rou*delta_max)/(delta(i,j)+rou*delta_max); %计算关联系数

end

end

r=sum(yita)./size(yita,1); %计算灰色关联