GM(1,N)模型程序

function GM1_1(X0)

%format long

[m,n]=size(X0)

X1=cumsum(X0) %累加

X2=[]

for i=1:n-1

X2(i,:)=X1(i)+X1(i+1)

end

B=-0.5.*X2

t=ones(n-1,1)

B=[B,t] % 求B矩阵

YN=X0(2:end)

P_t=YN./X1(1:(length(X0)-1)) %对原始数据序列X0进行准光滑性检验，

%序列x0的光滑比P(t)=X0(t)/X1(t-1)

A=inv(B.'*B)*B.'*YN.'

a=A(1)

u=A(2)

c=u/a

b=X0(1)-c

X=[num2str(b),'exp','(',num2str(-a),'k',')',num2str(c)]

strcat('X(k+1)=',X)

%syms k

for t=1:length(X0)

k(1,t)=t-1

end

k

Y_k_1=b*exp(-a*k)+c

for j=1:length(k)-1

Y(1,j)=Y_k_1(j+1)-Y_k_1(j)

end

XY=[Y_k_1(1),Y]%预测值

CA=abs(XY-X0)%残差数列

Theta=CA %残差检验 绝对误差序列

XD_Theta= CA ./ X0 %残差检验 相对误差序列

AV=mean(CA) % 残差数列平均值

R_k=(min(Theta)+0.5*max(Theta))./(Theta+0.5*max(Theta)) % P=0.5

R=sum(R_k)/length(R_k) %关联度

Temp0=(CA-AV).^2

Temp1=sum(Temp0)/length(CA)

S2=sqrt(Temp1)%绝对误差序列的标准差

%----------

AV_0=mean(X0)% 原始序列平均值

Temp_0=(X0-AV_0).^2

Temp_1=sum(Temp_0)/length(CA)

S1=sqrt(Temp_1) %原始序列的标准差

TempC=S2/S1*100 %方差比

C=strcat(num2str(TempC),'%') %后验差检验 %方差比

%----------

SS=0.675*S1

Delta=abs(CA-AV)

TempN=find(Delta<=SS)

N1=length(TempN)

N2=length(CA)

TempP=N1/N2*100

P=strcat(num2str(TempP),'%') %后验差检验%计算小误差概率

% GM(1,1)模型计算及检验、作图。文件名fungry1.m

function GM1=fungry1(x0) %输入原始数据x0

T=input('T=')

x1=zeros(1,length(x0))

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

yn=zeros(length(x0)-1,1)

Hatx0=zeros(1,length(x0)-1,2)

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)模型参数估计

a=HatA(1)

b=HatA(2)

for k=1:length(x0)+T

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

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

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

p=p/length(x0)

if p>0.95 &c<0.35

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

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

elseif p>0.85 &c<0.5

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

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

elseif p>0.70 &c<0.65

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

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

else p<=0.70 &c>0.65

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

end

for i=1:length(x0)

Hatx00(i)=Hatx0(i)

end

z=1:length(x0)

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

end

---