灰色预测算法应用

【级比判断】

clear
close all
clc

%导入数据
x0=[2.67 3.13 3.25 3.36 3.56 3.72];
n=length(x0);
%求级比
lamda=x0(1:n-1)./x0(2:n);
%级比判断
range=minmax(lamda);
range1=[exp(-2/(n+1)),exp(2/n+2)];
flag=1;

if(range(1)>range1(1)&&range(2)<range1(2))
    disp('级比检验通过');
    flag=1;
else
    flag=0;
    disp('级比检验未通过');
end

【级比检验通过】可以直接建模

%GM(1,1)建模

%AGO
x1=cumsum(x0);
%加权邻值生成
for i=2:n   
    z(i)=0.5*(x1(i)+x1(i-1));
end

B=[-z(2:n)',ones(n-1,1)];
Y=x0(2:n)';
u=B\Y;
a=u(1);b=u(2);
%预测后续数据
m=input('please input the number of the years you want to forecast:');
F=[];
F(1)=x0(1);
for i=2:(n+m)
    F(i)=(x0(1)-b/a)/exp(a*(i-1))+b/a;
end
disp('预测值为：')
F=[x0(1),diff(F)];
disp(F)

%计算残差  
epsilon=x0-F(1:n);

%计算相对误差
disp('相对误差为：');
delta=abs(epsilon./x0);
disp(delta)

%计算级比偏差
disp('级比偏差为：');
rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda ;
disp(rho)

%作图
t1=1999:2004;
t2=2000:2004+m;
plot(t1,x0,'r*');
hold on;
plot(t2,F(2:n+m),'g-');


xlabel('年份');
ylabel('销售额');
legend('实际销售额','预测销售额');
title('产品销售额变化变化情况');
grid on;

【检验不通过】

%导入数据
x0=[444,895,573,443,629,509,594,842,300,567,852,542,546,459,557,354,917,421,522,856,535,518,459,532,430,490,830,442,948,344];
n=length(x0);
%求级比
lamda=x0(1:n-1)./x0(2:n);
%级比判断
range=minmax(lamda)
range1=[exp(-2/(n+1)),exp(2/n+2)]
flag=1;
if(range(1)>range1(1)&&range(2)<range1(2))
    disp('级比检验通过');
    flag=1;
else
    flag=0;
    disp('级比检验未通过');
end

%若级比检验未通过，则进行平行变换
while(flag==0)
    C=input('请选择常数C进行平移变换：');
    for i=1:n
        y0(i)=x0(i)+C;
    end
    lamda=y0(1:n-1)./y0(2:n);
    range=minmax(lamda);
    if(range(1)>range1(1)&&range(2)<range1(2))
        flag=1;
    else flag=0;
    end
end

%GM(1,1)建模

%AGO
x1=cumsum(y0);
%加权邻值生成
for i=2:n   
    z(i)=0.5*(x1(i)+x1(i-1));
end

B=[-z(2:n)',ones(n-1,1)];
Y=y0(2:n)';
u=B\Y;
a=u(1);b=u(2);
%预测后续数据
m=input('please input the number of the years you want to forecast:');
F=[];F(1)=y0(1);
for i=2:(n+m)
    F(i)=(y0(1)-b/a)/exp(a*(i-1))+b/a;
end
disp('预测值为：')
F=[x0(1),diff(F)-C]

%计算残差  
epsilon=x0-F(1:n);

%计算相对误差
disp('相对误差为：');
delta=abs(epsilon./x0)

%计算级比偏差
disp('级比偏差为：');
rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda 

%作图
t1=1967:1996;
t2=1968:1996+m;
plot(t1,x0,'r*');
hold on;
plot(t2,F(2:n+m),'g-');
hold on;
y=ones(1,length(t2))*435;
plot(t2,y,'black');
xlabel('年份');
ylabel('降水量');
legend('实际降水量','预测降水量','阈值=435mm');
title('某地年降水量变化情况');
grid on;