灰色预测程序代码

灰色模型
灰色模型（Grey Model ），记为 GM 。

灰色建模是用原始数据列作生成后建立微分方程，如果有N 个变量，则记为GM （1，N ）。

特殊地，只有一个变量的数列，则为)1,1(GM ，用)1,1(GM 为例说明建模步骤：
6.1登录原始数列)
0(x ，并用)
0(x 累加生成数列)
1(x ：∑==
k
i i x
k x 1
)
0()
1()()(
6.2用累加生成数列构成矩阵B ，用原始数列构成数据列矩阵n y 。

⎪⎪⎪⎪⎪⎭⎫
⎝⎛+--+-+-=1))()1((5.01))3()2((5.01))2()1((5.0)1()1()
1()1()1()1(n x n x x x x x B ⎪⎪⎪⎪
⎪⎭
⎫ ⎝⎛=)()3()2()0()0()0(n x x x y n
6.3求微分方程参数u a ,。

生成数列模型的解为：
a
u e
a u x k x k
a +⋅-=+-))1(()1()0()1( 按最小二乘法，参数u a ,为：
n T
T y B B B u a ⋅⋅⋅=⎪⎪⎭
⎫
⎝⎛=-1)(ˆ 6.4建模。

将u a ,之值代入下式。

u ax dt
dx =+)1()
1(， a u e a u x k x k a +⋅-=+-))1(()1(ˆ)0()1( 为生成模型计算值。

6.5数据还原。

按下述公式还原：
)1(ˆ)(ˆ)(ˆ)1()1()
0(--=k x k x
k x
， )(ˆ)0(k x 为模型计算值。

6.6检验。

将模型计算值与实际值比较，求出残差：)(ˆ)()()0()
0(k x
k x
k q -= 根据前面的数据得出后一个数据，依次递推检验，每一个检验值对模型来说都是后验值，称为后验差检验。

记k 时刻实际值)()
0(k x
与计算值（预测值）)(ˆ)0(k x
之差为)(k q ： )(ˆ)()()0()0(k x
k x k q -= )(k q 即称为k 时刻残差。

记实际数据)()
0(k x ，1=k ，2，…，n 的平均值为x ： ∑==n k k x n x 1
)
0()(1
记残差)(k q ，1=k ，2，…，n 的平均值q ： ∑'
==n k k q n q 1
)(1，n n <'
记原始数据方差为21
S ：∑=-=n k x k x n S 1
2)0(21
))((1；记残差方差为2
2S ：
∑'
=-=n k q k q n S 1
22
2
))((1；12S S C =，称为后验差比值。

指标 C 越小越好。

C 越小，表示1
S 越大而2S 越小。

1S 大表明原始数据方差大，离散程度大。

2S 小表明残差方差小，离散程度小。

C 小表明尽管原始数据很离散，但模型所得计算值与实际值之差并不太离散。

}6745.0|)({|1S q k q P P <-=称为小频率误差。

P 越大越好，P 越大表明残差与残差
平均值之差小于给定值16745.0S 的点较多（以百分比记之）。

按C 与P 两个指标。

可综合评定预测模型的精度。

具体指标如表7一6所示：
例：
han1.txt ：
9.4 11.3 16.8 19.8 20.3 18.8 20.9 24.9 24.7 24.3 19.4 18.6 9.6 11.7 15.8 19.9 19.5 17.8 17.8 23.3 21.4 24.5 20.1 15.9 10.1 12.9 17.7 21.0 21.0 20.4 21.9 25.8 29.3 29.8 23.6 16.5 11.4 26.0 19.6 25.9 27.6 24.3 23.0 27.8 27.3 28.5 32.8 18.5 11.5 26.4 20.4 26.1 28.9 28.0 25.2 30.8 28.7 28.1 22.2 20.7 13.7 29.7 23.1 28.9 29.0 27.4 26.0 32.2 31.4 32.6 29.2 22.9 15.4 17.1 23.5 11.6 1.78 2.61 8.8 16.2 20.1 24.9 26.5 21.8
load han1.txt %把原始数据保存在纯文本文件han1.txt 中 han1(end,:)=[];m=size(han1,2); x0=mean(han1,2); x1=cumsum(x0)
alpha=0.4;n=length(x0);
z1=alpha*x1(2:n)+(1-alpha)*x1(1:n-1) Y=x0(2:n);B=[-z1,ones(n-1,1)]; ab=B\Y
k=6;
x7hat=(x0(1)-ab(2)/ab(1))*(exp(-ab(1)*k)-exp(-ab(1)*(k-1)))
z=m*x7hat
u=sum(han1)/sum(sum(han1))
v=z*u
例：（此题预测结果不理想）
X=[9.4 11.3 16.8 19.8 20.3 18.8 20.9 24.9 24.7 24.3 19.4 18.6;9.6 11.7 15.8 19.9 19.5 17.8 17.8 23.3 21.4 24.5 20.1 15.9;10.1 12.9 17.7 21.0 21.0 20.4 21.9 25.8 29.3 29.8 23.6 16.5;11.4 26.0 19.6 25.9 27.6 24.3 23.0 27.8 27.3 28.5 32.8 18.5;11.5 26.4 20.4 26.1 28.9 28.0 25.2 30.8 28.7 28.1 22.2 20.7;13.7 29.7 23.1 28.9 29.0 27.4 26.0 32.2 31.4 32.6 29.2 22.9;15.4 17.1 23.5 11.6 1.78 2.61 8.8 16.2 20.1 24.9 26.5 21.8];
%预测后三个数据
X=X';N=size(X);n1=N(1);n=N(2);
for i=1:n1
I=i
X0=X(i,:)
E=triu(ones(n)); X1=X0*E;
b1=X1; b1(1)=[]; b2=X1; b2(n)=[]; b=-0.5*(b1+b2);
B=[b;ones(1,n-1)]; B=B';
y=X0; y(1)=[];y=y';
format long
A=((inv(B'*B))*B')*y
A(2)/A(1)
X0(1)-A(2)/A(1)
k=0:n+2; X2=(X0(1)-A(2)/A(1))*exp(-k*A(1))+A(2)/A(1);
Xf=X2; Xf(n+3)=[]; Xff=[0 Xf];
X3=X2-Xff
X3(n+3)=[];X3(n+2)=[];X3(n+1)=[];
r=corrcoef(X0,X3)
Q=X0-X3;
format short
R=Q./X0
R1=max(abs(R))
s=mean(abs(R))
X0mean=mean(X0);S1=std(X0);
Qmean=mean(Q)
S2=std(Q);
C=S2/S1
D=abs(Q-Qmean);p0=0.6745*S1;m=0;
for j=1:n
d=D(j);
if (d<p0)
m=m+1;
end
end
P=m/n
E3=(Q*Q')/(n*(n-1));
E4=(1.96*sqrt(E3))/X0mean;
Pe=1-E4
D1=abs(Q);
Dmin=min(D1);Dmax=max(D1);r0=0;
for j=1:n
L=(Dmin+0.5*Dmax)/(D1(j)+0.5*Dmax);
r0=r0+L;
end
r0=r0/n
end。