灰色预测matlab代码怎么写

这是我曾经写过的一个灰色预测的程序：第一个文件为函数，需要在调用时输入原始数据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) %计算灰色关联

自己编的，希望对你有帮助

clear

x0=[89677,99215,109655,120333,135823,159878,182321,209407,246619,300670]

pre_num=10

n=length(x0)

disp('级比检验')

lambda=x0(1:end-1)./x0(2:end)

range=minmax(lambda)

x1=cumsum(x0)

z=0.5*(x1(2:end)+x1(1:end-1))

Y=x0(2:end)'

B=[-z(1:end)' ones(n-1,1)]

u=B\Y%u=inv(B'*B)*B'*Y

a=u(1)

b=u(2)

x0_pre=[x0(1) ones(1,n+pre_num-1)]

for k=1:n-1+pre_num

x0_pre(k+1)=(x0(1)-b/a)*(exp(-a*k)-exp(-a*(k-1)))

end

err=[x0 - x0_pre(1:n)]

epsilon=abs(err)./x0(1:n).*100

disp('预测值')

disp(x0_pre)

disp('相对误差')

disp(epsilon)

t1=1999:2008

t2=1999:2018

plot(t1,x0,'d',t2,x0_pre,'LineWidth',2) %原始数据与预测数据的比较

xlabel('年份')

ylabel('利润')

---