GM(1,1)灰色预测模型
Introduction
Initial
给定原始序列：
x(0) =(x(0)(1), x(0)(2), x(0)(3)…, x(0)(n))
Step 1
一次AGO(1-AGO)生成序列，以弱化原始序列的随机性和波动性：x(1) =(x(1)(1), x(1)(2), x(1)(3)…, x(1)(n)) Matlab Program
clear
syms a b;
c=[a b]';
fid=fopen('.\Grey Model\test.txt');
x0=fscanf(fid,'%f');x0=x0';
fclose(fid);
x1=cumsum(x0); %原始数据累加
n=length(x0);
for i=1:(n-1)
z(i)=(x1(i)+x1(i+1))/2; %生成累加矩阵end
%计算待定参数的值
Y=x0;Y(1)=[];
Y=Y';
B=[-z;ones(1,n-1)];B=B';
c=inv(B'*B)*B'*Y;
c=c';
a=c(1);b=c(2);
%预测后续数据
%预测之后10个时间单位的数据
xx1=[];xx1(1)=x0(1);
for i=2:(n+10)
xx1(i)=(x0(1)-b/a)/exp(a*(i-1))+b/a; end
xx0=[];xx0(1)=x0(1);
Step 2
(1) dx (1)
dt
+ax (1)(t )=u ，式中a, u 为待定系数。

灰微分方程模型为：
x (0)(k )+az (1)(k )=u ，z 为背景值
z (1)(k )=1/2(x (1)(k )+x (1)(k −1))
(2) 构造矩阵B 和数据向量Y n
Y n =Ba ̂
Y n =[ x (0)(2)x (0)(3)⋮x (0)(n )] , B =[ −1/2(x (1)(1)+x (1)(2)),−1/2(x (1)(2)+x (1)(3)),⋮−1/2(x (1)(n −1)+x (1)(n )), 1 1 ⋮ 1]
a ̂=(a
u
)=(B T B)−1B T Y n
Step 3
模型响应函数
x ̂(1)
(k +1)=(x (0)
(1)−u a )e −ak +u a
x ̂(0)(k +1)=x ̂(1)(k +1)−x ̂(1)(k )
Step 4
检验和判断GM(1,1)模型的精度 (1) 残差检验
for i=2:(n+10)
xx0(i)=xx1(i)-xx1(i-1); end
%关联度检验 for i=1:n
e(i)=abs(x0(i)-xx0(i)); end
mmax=max(e); for i=1:n
ee(i)=0.5*mmax/(e(i)+0.5*mmax); end
r=sum(ee)/n; %后验差检验
x0bar=sum(x0)/n; s1=0; for i=1:n
s1=s1+(x0(i)-x0bar)^2; end
s1=sqrt(s1/n); s2=0;
ebar=sum(e)/n; for i=1:n
s2=s2+(e(i)-ebar)^2; end
s2=sqrt(s2/n); C=s2/s1; p=0;
for i=1:n
if abs(e(i)-ebar)<0.6745*s1
绝对误差：ε(k)=|x(0)(k)−x̂(0)(k)|
相对误差：Φ(k)=ε(k)
x(0)(k)
(2) 关联度检验
分辨率β一般取0.5，此时若关联度大于0.6则认为模型可接受(3) 后验差检验和小误差概率
原始序列标准差：S1=√∑[x(0)(i)−x̅(0)]2
n
绝对误差序列标准差：S2=√∑[ε(i)−ε̅]2
n
计算方差比：C=S2
S1
小误差概率：P=P{|ε(i)−ε̅|<0.6745S1}
p=p+1;
end
end
p=p/n;
C
p
if p>0.95&C<0.35
disp('预测精度好');
else if p>0.8&C<0.5
disp('预测合格');
else if p>0.7&C<0.65
disp('预测勉强合格'); else
disp('预测不合格'); end
end
end
%原始数据与预测数据进行比较
t1=1:n;
t2=1:(n+10);
xx0
plot(t1,x0,'o',t2,xx0)。