灰色系统预测GM(1,1)模型及其Matlab 实现
预备知识
(1)灰色系统
白色系统是指系统内部特征是完全已知的;黑色系统是指系统内部信息完全未知的;而灰色系统是介于白色系统和黑色系统之间的一种系统,灰色系统其内部一部分信息已知,另一部分信息未知或不确定。
(2)灰色预测 灰色预测,是指对系统行为特征值的发展变化进行的预测,对既含有已知信息又含有不确定信息的系统进行的预测,也就是对在一定范围内变化的、与时间序列有关的灰过程进行 预测。尽管灰过程中所显示的现象是随机的、杂乱无章的,但毕竟是有序的、有界的,因此得到的数据集合具备潜在的规律。灰色预测是利用这种规律建立灰色模型对灰色系统进行预测。
目前使用最广泛的灰色预测模型就是关于数列预测的一个变量、一阶微分的GM(1,1)模型。它是基于随机的原始时间序列,经按时间累加后所形成的新的时间序列呈现的规律可用一阶线性微分方程的解来逼近。经证明,经一阶线性微分方程的解逼近所揭示的原始时间序列呈指数变化规律。因此,当原始时间序列隐含着指数变化规律时,灰色模型GM(1,1)的预测是非常成功的。
1 灰色系统的模型GM(1,1)
1.1 GM(1,1)的一般形式
设有变量X (0)={X (0)(i),i=1,2,...,n}为某一预测对象的非负单调原始数据列,为建立灰色预测模型:首先对X (0)进行一次累加(1—AGO, Acumulated Generating Operator)生成一次累加序列:
X (1)={X (1)(k ),k =1,2,…,n}
其中
X (1)(k )=
∑
=k
i 1
X (0)(i)
=X (1)(k -1)+ X (0)(k ) (1) 对X (1)可建立下述白化形式的微分方程:
dt
dX )1(十)
1(aX =u (2)
即GM(1,1)模型。
上述白化微分方程的解为(离散响应): ∧
X (1)(k +1)=(X (0)(1)-
a u )ak e -+a
u
(3)
或
∧
X (1)(k )=(X (0)(1)-
a u ))1(--k a e +a
u (4)
式中:k 为时间序列,可取年、季或月。 1.2 辩识算法
记参数序列为∧a , ∧a =[a,u]T
,∧
a 可用下式求解:
∧
a =(B T B)-1B T Y n (5)
式中:B —数据阵;Y n —数据列
B =?????????
???????????+++- 1 (n))X 1)-(n (X 21 ... 1 (3))X (2)X (211 (2))X (1)X (21(1)1(1)
(1)(1)(1))(-- (6) Y n =(X (0)(2), X (0)(3),…, X (0)(n))T (7)
1.3 预测值的还原
由于GM 模型得到的是一次累加量,k ∈{n+1,n+2,…}时刻的预测值,必须将GM 模型所得数据∧
X
(1)(k +1)(或
∧
X
(1)(k ))经过逆生成,即累减生成(I —AGO)还原为∧
X (0)(k +1)(或
∧
X (0)(k ))才能用。
∧
X
(1)(k )=
∑
=k
i 1
∧
X (0)(i)
=
∑
-=1
1
k i ∧
X (0)(i)+∧
X (0)(k)
∧
X
(0)(k)=∧
X
(1)(k )-
∑
-=1
1
k i ∧
X (0)(i)
因为∧
X
(1)(k -1)=
∑
-=1
1
k i ∧
X
(0)(i),所以∧
X (0)(k)=∧
X (1)(k )-∧
X (1)(k -1)。
2 应用举例
取某高校1998年~2003年的某专业招生数据建模,见表1。
以表1中的数据构造原始数据列X (0),即
X (0)={X (0)(1),X (0)(2),X (0)(3),X (0)(4),X (0)(5),X (0)(6)} ={132,92,118,130,187,207}
对X (0)进行一次累加(1—AGO),生成数列:
X (1)
(k)=
∑
=k
i 1
X (0)(i)即
X (1)={X (1)(1),X (1)(2),X (1)(3),X (1)(4),X (1)(5),X (1)(6)}
={132,224,342,472,659,866} 和数据阵B 、数据列Y n
B =???
????
??????
???----- 1 5.762 1
5.5651 4071 2831
178,Y n =(92,118,130,187,207)T 由式(5)得
∧
a =[a,u]T
=??
????-7878184.56205
.0 由式(4)得灰色预测模型GM(1,1)为 ∧
X (1)(k )=(X (0)(1)-
a u ))1(--k a e +a
u
=(132+277.0137483))
1(205.0-k e -277.0137483
=409.0137483)
1(205.0-k e -277.0137483
预测值及预测精度见表2。
表2 某高校专业招生预测值及预测精度表
由表2知预测精度较高。2006年某专业招生人数预测值为259人。由于人数为整数,所以结果取整数部分。
GM(1,1)是一种长期预测模型,在没有大的市场波动及政策性变化的前提下,该预
测值应是可信的。如前所述,影响招生人数的因素很多且难以预测。因此,在采用灰色系统理论进行定量预测时,如果存在对预测对象影响较大的因素,就要在定性分析的基础上,寻找原始数据信息的突变点的量化值,然后再对预测值进行必要的修正,使预测值更接近实际情况,提高预测值的可信度,为科学决策提供可靠的数据。另外,若作长期预测,要考虑对上限值的约束条件。
应用灰色预测模型GM(1,1),对某专业招生人数进行了预测,具有较高的预测精度。应用灰色模型进行预测较之其它常规的时间序列预测法有以下显著的特点。
(1)灰色模型是一种长期预测模型,将预测系统中的随机元素作为灰色数据进行处理,而找出数据的内在规律。进行预测所需原始数据量小,预测精度较高,无须像其它预测法要么需要数据量大且规律性强,要么需要凭经验给出系数。
(2)理论性强,计算方便,籍助计算机及其程序设计语言,使得数据处理简便、快速、准确性好。
(3)用有限的表征系统行为特征的外部元素,分析系统的内在规律。灰色系统理论采用对系统的行为特征数据进行生成的方法,对杂乱无章的系统的行为特征数据进行处理,从杂乱无章的现像中发现系统的内在规律,这是该方法的独特之处。
(4) 适用性强。用灰色模型既可对周期性变化的系统行为进行预测,亦可对非周期性变化的系统行为进行预测;既可进行宏观长期的预测,亦可用于微观短期的预测。 2灰色系统模型的检验 定义1.
设原始序列
{})(,),2(),1()0()0()0()0(n x x x X Λ=
相应的模型模拟序列为
{
}
)(?,),2(?),1(??)0()0()0()0(n x x x X
Λ= 残差序列
{})(),2(),1()0(n εεεεΛ=
{
}
)(?)(,),2(?)2(),1(?)1()0()0()0()0()
0()0(n x n x x
x x x ---=Λ 相对误差序列
???
???=?)()(,,)
2()2(,)1()1()0()0()0(n x n x x εεεΛ
{}n
k 1?=
1.对于k <n,称)
()
()0(k x k k ε=
?为k 点模拟相对误差,称)
()
()0(n x n n ε=
?为滤波相对误差,
称∑=?=?n
k k n 1
1为平均模拟相对误差;
2.称?-1为平均相对精度,n ?-1为滤波精度;
3.给定α,当α,且α
设)0(X 为原始序列,)0(?X
为相应的模拟误差序列,ε为)0(X 与)0(?X 的绝对关联度,若对于给定的00,0εεε>>,则称模型为关联合格模型。
定义3
设)0(X 为原始序列,)
0(?X
为相应的模拟误差序列,)
0(ε
为残差序列。
∑==n k k x n x 1
)
0()(1为)0(X 的均值,
21)0(2
1))((1x k x n s n k -=∑=为)0(x 的方差,
∑==n
k k n 1)(1εε为残差均值,
∑=-=n k k n s 122
2))((1εε为残差方差,
1. 称1
2s s
c =为均方差比值;对于给定的00>c ,当0c c <时,称模型为均方差比合格模
型。 2. 称()16745.0)(s k P
p <-=εε为小误差概率,
对于给定的00>p ,当0p p >时,称模型为小误差概率合格模型。
应用举例2、设原始序列
{}
)5(),4(),3(),2(),1()0()0()0()0()0()0(x x x x x X =
()679.3,390.3,337.3,278.3,874.2=
建立GM(1,1)模型,并进行检验。 解:1)对)0(X
作1-AGO ,得
[D 为)
0(X 的一次累加生成算子,记为1-AGO]
{
}
)5(),4(),3(),2(),1()1()1()1()1()1()
1(x x x x x X
= ()558.16,579.12,489.9,152.6,874.2=
2)对)
1(X
作紧邻均值生成,令
)1(5.0)(5.0)()1()1()
1(-+=k x k x k Z
{}
)5(),4(),3(),2(),1()1()1()1()1()1()1(z z z z z Z =
()718.14,84.11,820.7,513.4,874.2=
于是,
??????????----=
??????
??????----=1718.14184
.111820.71513.41)5(1)
4(1)3(1)2()1()
1()1()1(z z z z B ,??????????=?????
???????=679.3390.3337.3278.3)5()4()3()2()0()0()0()
0(x x x x Y
????
?
?????----???????----=T 1718
.14184
.111820.71513
.41111718.14184.11820.7513.4B B ??
????--=4235.38235.38221.423
??
????==??????--=--T
832371.11665542.0165542.0017318.04235.38235.38221.423)(1
1
B B
??
???
???????????
????=??????-?=221.423235.38235.384969.2301221.423235.38235.384235.384221.42312
??
??
?????????????----???????==T -T 679.3390.3337.3278.31111718.14184.11820.7513.4832371.11665542.0165542.0017318.0)(?1Y B B B a
??
??
?
????????????---=679.3390.3337.3278.3604076.10019051.0537833.0085280.1089344.0028143.0030115.0087386.0 ??
????-=065318.3037156.0 3)确定模型
065318.3037156.0)1()
1(=-x dt
dx 及时间响应式
a b
e a b x k x
ak +-=+-))1(()1(?)0()1( 4986.823728.85037156.0-=k
e
4)求)
1(X 的模拟值
{})5(?),4(?),3(?),2(?),1(??)1()1()1()1()1()1(x x x x x X
= =(2.8740,6.1058,9.4599,12.9410,16.5538)
5)还原出)
0(X 的模拟值,由
)(?)1(?)1(?)1()1()0(k x k x k x
-+=+ 得 {})5(?),4(?),3(?),2(?),1(??)0()0()0()0()0()0(x x x x x X
=
=(2.8740,3.2318,3.3541,3.4811,3.6128)
残差平方和
[]???
?
???
????==T )5()4()3()2()5()4()3()2(εεεεεεεεεεs []????
?
?????--
?--=0662.00911.00171.00462.00662.00911.00171.00462.0
=0.0151085
平均相对误差
%)80.1%69.2%51.0%41.1(4
1
4151+++=?=?∑=k k
=1.0625%
计算X 与X
?的灰色关联度 ))1()5((21)1()((4
2
x x x k x S k -+-=
∑= =)874.2679.3(2
1
)874.2390.3()874.2337.3()874.2278.3(-+
-+-+- 0.40250.5160.4630.404+++=
=1.7855
)1(?)5(?(21
)1(?)(?(?4
2
x x x k x S
k -+-=∑= )874.26128.3(2
1
)874.24811.3()874.23541.3()874.22318.3(-+-+-+-=
3694.06071.04801.03578.0+++=
=1.8144
[][]∑=---+---=-4
2
))1(?)5(?())1()5((2
1
))1(?)(?())1()((?k x x x x x k x
x k x S S
)4025.03694.0(2
1
)516.06071.0()463.04801.0()404.03578.0(-+
-+-+-=01655.0091.00171.00462.0-++-=
=0.04535
64525
.45999
.404535.08144.17855.118144.17855.11??1?1=+++++=
-+++++=
S S S
S S S ε
=0.9902>0.90
精度为一级,可以用
4986.823728.85)1(?037156.0)1(-=+k e k x
)(?)1(?)1(?)1()1()0(k x k x k x
-+=+预测。
function [y,p,e]=gm_1_1(X,k)
%gray model: GM(1,1)
%Example [y,p]=gm_1_1([200 250 300 350],2)
if nargout>3,error('Too many output argument.');end
if nargin==1,k=1;x_orig=X;
elseif nargin==0|nargin>2
error('Wrong number of input arguments.');
end
x_orig=X;
predict=k;
%AGO process
x=cumsum(x_orig);
%compute the coefficient(a and u)------------------------ n=length(x_orig);
%first generate the matrix B
for i=1:(n-1);
B(i)=-(x(i)+x(i+1))/2;
end
B=[B' ones(n-1,1)];
%then generate the matrix Y
for i=1:(n-1);
y(i)=x_orig(i+1);
end
Y=y';
%get the coefficient. a=au(1) u=au(2)
au=(inv(B'*B))*(B'*Y);
%-------------------------------------------------------- %change the grey model to symbolic expression
coef1=au(2)/au(1);
coef2=x_orig(1)-coef1;
coef3=0-au(1);
costr1=num2str(coef1);
costr2=num2str(abs(coef2));
costr3=num2str(coef3);
eq=strcat(costr1,'+',costr2,'e^(',costr3,'*(t-1))');
%comparison of calculated and observed value
for t=1:n+predict
mcv(t)=coef1+coef2*exp(coef3*(t-1));
end
x_mcv0=diff(mcv);
x_mcve=[x_orig(1) x_mcv0];
x_mcv=diff(mcv(1:end-predict));
x_orig_n=x_orig(2:end);
x_c_error=x_orig_n-x_mcv;
x_error=mean(abs(x_c_error./x_orig_n));
if x_error>0.2
disp('model disqualification!');
elseif x_error>0.1
disp('model check out');
else
disp('model is perfect!');
end
%predicting model and plot gragh
plot(1:n,x_orig,'diamond',1:n+predict,x_mcve);
p=x_mcve(end-predict+1:end);
xlabel('CURVE OF GREY MODEL ANALYSIS');
title('GM(1,1)');
grid on
y=eq;
e=x_error;
p=x_mcve(end-predict+1:end);
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 +ax1t=u,式中a, u为待定系数。 灰微分方程模型为: x0k+az1k=u,z为背景值 z1k=1/2(x1k+x1k?1) (2) 构造矩阵B和数据向量Y n Y n=Ba Y n=x02 x03 ? x0n , B= ?1/2(x11+x12), ?1/2(x12+x13), ? ?1/2(x1n?1+x1n), 1 1 ? 1 a=a u=(B T B)?1B T Y n Step 3 模型响应函数 x1k+1= x01?u e?ak+ u x0k+1=x1k+1?x1k 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
灰色预测模型的Matlab 程序及检验程序 %灰色预测模型程序 clear syms a b; c=[a b]'; A=[46.2 32.6 26.7 23.0 20.0 18.9 17.5 16.3];% 原始序列 B=cumsum(A);%累加n=length(A); for i=1:(n-1) C(i)=(B(i)+B(i+1))/2; end %计算待定参数 D=A; D(1)=[]; D=D'; E=[-C; ones(1,n-1)]; c=inv(E*E')*E*D; c=c'; a=c(1); b=c(2); %预测往后预测5个数据 F=[];F(1)=A(1); for i=2:(n+5) F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a; end G=[];G(1)=A(1); for i=2:(n+5) G(i)=F(i)-F(i-1); end t1=2002:2009; t2=2002:2014; G plot(t1,A,'o',t2,G) %灰色预测模型检验程序 function [ q,c,p ] = checkgm( x0,x1 ) %GM 检验函数 %x0 原始序列
%x1 预测序列 %·返回值 % q –- 相对误差 % c -- ·方差比 % p -- 小误差概率 e0=x0-x1; q=e0/x0; s1=var(x0); %qpa=mean(e0); s2=var(e0); c=s2/s1; len=length(e0); p=0; for i=1:len if(abs(e0(i)) < 0.6745*s1) p=p+1; end end p=p/len; end
作用:求累加数列、求a b的值、求预测方程、求残差 clc %清屏,以使结果独立显示 x=[ ]; format long; %设置计算精度 if length(x(:,1))==1 %对输入矩阵进行判断,如不是一维列矩阵,进行转置变换 x=x'; end n=length(x); %取输入数据的样本量 z=0; for i=1:n %计算累加值,并将值赋予矩阵be z=z+x(i,:); be(i,:)=z; end for i=2:n %对原始数列平行移位 y(i-1,:)=x(i,:); end
for i=1:n-1 %计算数据矩阵B的第一列数据 c(i,:)=*(be(i,:)+be(i+1,:)); end for j=1:n-1 %计算数据矩阵B的第二列数据 e(j,:)=1; end for i=1:n-1 %构造数据矩阵B B(i,1)=c(i,:); B(i,2)=e(i,:); end alpha=inv(B'*B)*B'*y; %计算参数矩阵即a b的值 for i=1:n+1 %计算数据估计值的累加数列,如改为n+1为n+m可预测后m-1个值 ago(i,:)=(x(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1,:)*(i-1))+alpha(2,:)/alpha(1,: );%显示输出预测值的累加数列 end var(1,:)=ago(1,:) %显示输出预测值 for i=1:n %如改n为n+m-1,可预测后m-1个值 var(i+1,:)=ago(i+1,:)-ago(i,:); %估计值的累加数列的还原,并计算出下一预测值end for i=1:n error(i,:)=x(i,:)-var(i,:); %计算残差 end c=std(error)/std(x); %调用统计工具箱的标准差函数计算后验差的比值c ago %显示输出预测值的累加数列 alpha %显示输出参数数列 var %显示输出预测值 error %显示输出误差 c %显示后验差的比值 作用:数据处理判断是否可以用灰色预测、求级比、求累加数列、求a b的值、求预测方程 clc,clear x0=[ ]'; %注意这里为列向量 n=length(x0); lamda=x0(1:n-1)./x0(2:n) %计算级比 range=minmax(lamda') %计算级比的范围 x1=cumsum(x0) %累加运算 B=[*(x1(1:n-1)+x1(2:n)),ones(n-1,1)]; Y=x0(2:n); u=B\Y %拟合参数u(1)=a,u(2)=b x=dsolve('Dx+a*x=b','x(0)=x0'); %求微分方程的符号解
function g(x); %定义函数gm(x) for j=1:10 clc; %清屏,以使计算结果独立显示 format long; %设置计算精度 if length(x(:,1))==1 %对输入矩阵进行判断,如不是一维列矩阵,进行转置变换x=x'; end; n=length(x); %取输入数据的样本量 z=0; for i=1:n %计算累加值,并将值赋与矩阵be z=z+x(i,:); be(i,:)=z; end for i=2:n %对原始数列平行移位 y(i-1,:)=x(i,:); end for i=1:n-1 %计算数据矩阵B的第一列数据 c(i,:)=-0.5*(be(i,:)+be(i+1,:)); end for j=1:n-1 %计算数据矩阵B的第二列数据 e(j,:)=1; end for i=1:n-1 %构造数据矩阵B B(i,1)=c(i,:); B(i,2)=e(i,:); end alpha=inv(B.'*B)*B.'*y; %计算参数α、μ矩阵 for i=1:n+1 %计算数据估计值的累加数列,如改n+1为n+m可预测后m-1个值ago(i,:)=(x(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1,:)*(i-1))+alpha(2,:)/alpha(1,:); end var(1,:)=ago(1,:) for i=1:n %如改n为n+m-1,可预测后m-1个值 var(i+1,:)=ago(i+1,:)-ago(i,:); %估计值的累加数列的还原,并计算出下一预测值end for i=1:n error(i,:)=var(i,:)-x(i,:); %计算残差 end c=std(error)/std(x); %调用统计工具箱的标准差函数计算后验差的比值c ago %显示输出预测值的累加数列 alpha %显示输出参数α、μ数列 var %显示输出预测值 error %显示输出误差 c %显示后验差的比值c var=fix(var);
Matlab的灰色预测程序: y=input('请输入数据'); n=length(y); yy=ones(n,1); yy(1)=y(1); for i=2:n yy(i)=yy(i-1)+y(i) end B=ones(n-1,2); for i=1:(n-1) B(i,1)=-(yy(i)+yy(i+1))/2; B(i,2)=1; end BT=B'; for j=1:(n-1) YN(j)=y(j+1); end YN=YN'; A=inv(BT*B)*BT*YN; a=A(1); u=A(2); t=u/a; t_test=input('输入需要预测的个数'); i=1:t_test+n; yys(i+1)=(y(1)-t).*exp(-a.*i)+t; yys(1)=y(1); for j=n+t_test:-1:2 ys(j)=yys(j)-yys(j-1); end x=1:n; xs=2:n+t_test; yn=ys(2:n+t_test); plot(x,y,'^r',xs,yn,'*-b'); det=0; for i=2:n det=det+abs(yn(i)-y(i)); end det=det/(n-1); disp(['百分绝对误差为:',num2str(det),'%']); disp(['预测值为:',num2str(ys(n+1:n+t_test))]); 请输入数据[29.8 30.11 41.05 70.12 77.79 77.79 104.82 65.22 82.7 100.79] 输入需要预测的个数4
百分绝对误差为:14.5128% 预测值为:110.5718 120.8171 132.0116 144.2434
灰色预测模型matlab程序 %下面程序是灰色模型GM(1,1)程序二次拟合和等维新陈代谢改进预测程序,mat lab6.5 ,使用本程序请注明,程序存储为gm1.m %x = [5999,5903,5848,5700,7884];gm1(x); 测试数据 %二次拟合预测GM(1,1)模型 function gmcal=gm1(x) sizexd2 = size(x,2); %求数组长度 k=0; for y1=x k=k+1; if k>1 x1(k)=x1(k-1)+x(k); %累加生成 z1(k-1)=-0.5*(x1(k)+x1(k-1)); %z1维数减1,用于计算B yn1(k-1)=x(k); else x1(k)=x(k); end end %x1,z1,k,yn1 sizez1=size(z1,2); %size(yn1); z2 = z1'; z3 = ones(1,sizez1)'; YN = yn1'; %转置 %YN B=[z2 z3]; au0=inv(B'*B)*B'*YN; au = au0'; %B,au0,au
ufor = au(2); ua = au(2)./au(1); %afor,ufor,ua %输出预测的 a u 和 u/a的值 constant1 = x(1)-ua; afor1 = -afor; x1t1 = 'x1(t+1)'; estr = 'exp'; tstr = 't'; leftbra = '('; rightbra = ')'; %constant1,afor1,x1t1,estr,tstr,leftbra,rightbra strcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,r ightbra,'+',leftbra,num2str(ua),rightbra) %输出时间响应方程 %****************************************************** %二次拟合 k2 = 0; for y2 = x1 k2 = k2 + 1; if k2 > k else ze1(k2) = exp(-(k2-1)*afor); end end %ze1 sizeze1 = size(ze1,2); z4 = ones(1,sizeze1)'; G=[ze1' z4]; X1 = x1'; au20=inv(G'*G)*G'*X1; au2 = au20'; %z4,X1,G,au20
%x=[1019,1088,1324,1408,1601];gm1(x); 测试数据%二次拟合预测GM(1,1)模型 function gmcal=gm1(x) if nargin==0 x=[1019,1088,1324,1408,1601] end format long g sizex=length(x); %求数组长度 k=0; for y1=x k=k+1; if k>1 x1(k)=x1(k-1)+x(k); %累加生成 z1(k-1)=-0.5*(x1(k)+x1(k-1)); %z1维数减1,用于计算B yn1(k-1)=x(k); else x1(k)=x(k); end end %x1,z1,k,yn1 sizez1=length(z1); %size(yn1); z2 = z1'; z3 = ones(1,sizez1)'; YN = yn1'; %转置 %YN B=[z2 z3]; au0=inv(B'*B)*B'*YN; au = au0'; %B,au0,au afor = au(1); ufor = au(2); ua = au(2)./au(1); %afor,ufor,ua %输出预测的 a u 和 u/a的值 constant1 = x(1)-ua; afor1 = -afor; x1t1 = 'x1(t+1)'; estr = 'exp'; tstr = 't'; leftbra = '(';
rightbra = ')'; %constant1,afor1,x1t1,estr,tstr,leftbra,rightbra strcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,rightbra,'+ ',leftbra,num2str(ua),rightbra) %输出时间响应方程 %****************************************************** %二次拟合 k2 = 0; for y2 = x1 k2 = k2 + 1; if k2 > k else ze1(k2) = exp(-(k2-1)*afor); end end %ze1 sizeze1=length(ze1); z4 = ones(1,sizeze1)'; G=[ze1' z4]; X1 = x1'; au20=inv(G'*G)*G'*X1; au2 = au20'; %z4,X1,G,au20 Aval = au2(1); Bval = au2(2); %Aval,Bval %输出预测的 A,B的值 strcat(x1t1,'=',num2str(Aval),estr,leftbra,num2str(afor1),tstr,rightbra,'+',lef tbra,num2str(Bval),rightbra) %输出时间响应方程 nfinal = sizex-1 + 1;(其中+1可改为+5等其他数字,即可预测更多的数字) %决定预测的步骤数5 这个步骤可以通过函数传入 %nfinal = sizexd2 - 1 + 1; %预测的步骤数 1 for k3=1:nfinal x3fcast(k3) = constant1*exp(afor1*k3)+ua; end %x3fcast %一次拟合累加值 for k31=nfinal:-1:0 if k31>1 x31fcast(k31+1) = x3fcast(k31)-x3fcast(k31-1); else if k31>0
灰色预测 专设工⑼他QA—叫吋)为原始数列.其1次累 ?加生成数列为恥=妙①曲⑵,…卅何),其中 X° 仇)二工* ° (0.址=1=2= -:n 5-1 卷定义卫的决导数为 d(k) = *町(上)=x 叫咼-x cl)(Jt-l). 令为数列工①的邻值生成数列.即 却(去)=^(*) + (1- a)x0)(t-lX 于是定义GM (L 1)的灰微分方程模型为 d(k)-血⑴住)=K 即或严>(£) + “尹⑻=人⑴ 在式(1)中』。>(灼称为灰导数,我称为发展系数, 弧称为白化背景值,b称为灰作用量乜 将时刻表殳二2「3「/代入(1)式有 V!1「— a y=代⑶ B =I b*- : X闵0)-Z,:](K)1
于是G\I <1?1)複至可表示为Y = Bu. 現在问题归结为求sb 在值。用一元线性回归?即最小二秦法求它们的活计值 为 注二实陌上回归分析中求估计值是用软件计尊的?有标准程序求解,iOmaClab 等。 GM <1? 1>的白化晏 対于G\I <1> 1)的灰微分方程(1) >如果将灰导数打(Q 的时刻 视为连绫变里"则x°)视为时问(函数卅⑺,于是*〉(Q 対血于导数里级 心2 >白化背臬值申的对应于导数卅⑴。于是G\I (1,1)的坝徽 分方樂対应于的白微分 方程为 内?则数堀列X?可以塗互G\I <19 1) 且可以进行页 色预测。否朋,対数摄做适当的克换处理■如平移叢换: 取C 使得鞍据列 严伙)=工⑴伙)+ G 上=1,2,…, 的级比都華住可吝禎盖内。 心⑴⑴ + o?i> (r)二dr <2) GM mi )质色预测的步骤 1 ?教摇的枪绘与处連 为了ftilGAl (1,1)建複方法的可行性,亲要为已知期S 做必要的检蛉 处理。 设原始教据列为了 逛=(乂°(1)*6(2)严炉00; >计算数列的级比 如果所有的级比都落在可容覆盖区间 ? fc = A- 2,3"?
灰色预测模型的Matlab程序及检验程序%灰色预测模型程序 clear syms a b; c=[a b]'; A=[46.232.626.723.020.018.917.516.3];%原始序列B=cumsum(A);%累加 n=length(A); for i=1:(n-1) C(i)=(B(i)+B(i+1))/2; end %计算待定参数 D=A; D(1)=[]; D=D'; E=[-C;ones(1,n-1)]; c=inv(E*E')*E*D; c=c'; a=c(1); b=c(2); %预测往后预测5个数据 F=[];F(1)=A(1); for i=2:(n+5) F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a; end G=[];G(1)=A(1); for i=2:(n+5) G(i)=F(i)-F(i-1); end t1=2002:2009; t2=2002:2014; G plot(t1,A,'o',t2,G) %灰色预测模型检验程序 function[q,c,p]=checkgm(x0,x1) %GM检验函数 %x0原始序列 %x1预测序列 %·返回值
%q–-相对误差 %c--·方差比 %p--小误差概率 e0=x0-x1; q=e0/x0; s1=var(x0); %qpa=mean(e0); s2=var(e0); c=s2/s1; len=length(e0); p=0; for i=1:len if(abs(e0(i))<0.6745*s1) p=p+1; end end p=p/len; end 等级相对误差q方差比C小误差概论P I级<0.01<0.35>0.95 II级<0.05<0.50<0.80 III级<0.10<0.65<0.70 IV级>0.20>0.80<0.60
灰色预测M A T L A B程 序 文档编制序号:[KKIDT-LLE0828-LLETD298-POI08]
灰色预测 作用:求累加数列、求a b的值、求预测方程、求残差 clc %清屏,以使结果独立显示 x=[ ]; format long; %设置计算精度 if length(x(:,1))==1 %对输入矩阵进行判断,如不是一维列矩阵,进行转置变换 x=x'; end n=length(x); %取输入数据的样本量 z=0; for i=1:n %计算累加值,并将值赋予矩阵be z=z+x(i,:); be(i,:)=z; end for i=2:n %对原始数列平行移位 y(i-1,:)=x(i,:); end for i=1:n-1 %计算数据矩阵B的第一列数据 c(i,:)=*(be(i,:)+be(i+1,:)); end for j=1:n-1 %计算数据矩阵B的第二列数据 e(j,:)=1; end for i=1:n-1 %构造数据矩阵B B(i,1)=c(i,:); B(i,2)=e(i,:); end alpha=inv(B'*B)*B'*y; %计算参数矩阵即a b的值 for i=1:n+1 %计算数据估计值的累加数列,如改为n+1为n+m可预测后m-1个值 ago(i,:)=(x(1,:)-alpha(2,:)/alpha(1,:))*exp(-alpha(1,:)*(i- 1))+alpha(2,:)/alpha(1,:);%显示输出预测值的累加数列 end var(1,:)=ago(1,:) %显示输出预测值 for i=1:n %如改n为n+m-1,可预测后m-1个值 var(i+1,:)=ago(i+1,:)-ago(i,:); %估计值的累加数列的还原,并计算出下一预测值end for i=1:n error(i,:)=x(i,:)-var(i,:); %计算残差 end c=std(error)/std(x); %调用统计工具箱的标准差函数计算后验差的比值c ago %显示输出预测值的累加数列
灰色预测模型matlab程序 灰色模型预测是在数据不呈现一定规律下可以采取的一种建模和预测方法,其预测数据与原始数据存在一定的规律相似性 %下面程序是灰色模型GM(1,1)程序二次拟合和等维新陈代谢改进预测程序,mat lab6.5 ,使用本程序请注明,程序存储为gm1.m %x = [5999,5903,5848,5700,7884];gm1(x); 测试数据 %二次拟合预测GM(1,1)模型 function gmcal=gm1(x) sizexd2 = size(x,2); %求数组长度 k=0; for y1=x k=k+1; if k>1 x1(k)=x1(k-1)+x(k); %累加生成 z1(k-1)=-0.5*(x1(k)+x1(k-1)); %z1维数减1,用于计算B yn1(k-1)=x(k); else x1(k)=x(k); end end %x1,z1,k,yn1 sizez1=size(z1,2); %size(yn1); z2 = z1'; z3 = ones(1,sizez1)'; YN = yn1'; %转置 %YN B=[z2 z3]; au0=inv(B'*B)*B'*YN;
%B,au0,au afor = au(1); ufor = au(2); ua = au(2)./au(1); %afor,ufor,ua %输出预测的 a u 和 u/a的值 constant1 = x(1)-ua; afor1 = -afor; x1t1 = 'x1(t+1)'; estr = 'exp'; tstr = 't'; leftbra = '('; rightbra = ')'; %constant1,afor1,x1t1,estr,tstr,leftbra,rightbra strcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,r ightbra,'+',leftbra,num2str(ua),rightbra) %输出时间响应方程 %****************************************************** %二次拟合 k2 = 0; for y2 = x1 k2 = k2 + 1; if k2 > k else ze1(k2) = exp(-(k2-1)*afor); end end %ze1 sizeze1 = size(ze1,2); z4 = ones(1,sizeze1)'; G=[ze1' z4]; X1 = x1'; au20=inv(G'*G)*G'*X1;
灰色理论预测模型及GM(1,1)matlab程序灰色预测方法简介 灰色预测是一种对含有不确定因素的系统进行预测的方法。灰色预测通过鉴别系统因素之间发展趋势的相异程度,即进行关联分析,并对原始数据进行生成处理来寻找系统变动的规律,生成有较强规律性的数据序列,然后建立相应的微分方程模型,从而预测事物未来发展趋势的状况。其用等时距观测到的反应预测对象特征的一系列数量值构造灰色预测模型,预测未来某一时刻的特征量,或达到某一特征量的时间。通过对原始数据的整理寻找数的规律,分为三类: a、累加生成:通过数列间各时刻数据的依个累加得到新的数据与数列。累加前数列为原始数列,累加后为生成数列。 b、累减生成:前后两个数据之差,累加生成的逆运算。累减生成可将累加生成还原成非生成数列。 c、映射生成:累加、累减以外的生成方式。 建模步骤 a、建模机理 b、把原始数据加工成生成数; c、对残差(模型计算值与实际值之差)修订后,建立差分微分方程模型; d、基于关联度收敛的分析; e、gm模型所得数据须经过逆生成还原后才能用。 f、采用“五步建模(系统定性分析、因素分析、初步量化、动态量化、优化)”法,建立一种差分微分方程模型gm(1,1)预测模型。 GM(1,1)程序: % 本程序主要用来计算根据灰色理论建立的模型的预测值。 % 应用的数学模型是GM(1,1)。 % 原始数据的处理方法是一次累加法。 clear;clc; % load ('data.txt');
% y=data'; y=[3 4 5 4 7 7]; n=length(y); yy=ones(n,1); yy(1)=y(1); for i=2:n yy(i)=yy(i-1)+y(i); end B=ones(n-1,2); for i=1:(n-1) B(i,1)=-(yy(i)+yy(i+1))/2; B(i,2)=1; end BT=B'; for j=1:n-1 YN(j)=y(j+1); end YN=YN'; A=inv(BT*B)*BT*YN; a=A(1); u=A(2); t=u/a; t_test=input('请输入需要预测个数:'); i=1:t_test+n; yys(i+1)=(y(1)-t).*exp(-a.*i)+t; yys(1)=y(1); for j=n+t_test:-1:2 ys(j)=yys(j)-yys(j-1); end x=1:n; xs=2:n+t_test; yn=ys(2:n+t_test); plot(x,y,'^r',xs,yn,'*-b'); det=0; for i=2:n det=det+abs(yn(i)-y(i)); end det=det/(n-1); disp(['百分绝对误差为:',num2str(det),'%']); disp(['预测值为:',num2str(ys(n+1:n+t_test))]);
灰色系统预测GM(1,1)模型及其Matlab 实现 预备知识 (1)灰色系统 白色系统是指系统内部特征是完全已知的;黑色系统是指系统内部信息完全未知的;而灰色系统是介于白色系统和黑色系统之间的一种系统,灰色系统其内部一部分信息已知,另一部分信息未知或不确定。 (2)灰色预测 灰色预测,是指对系统行为特征值的发展变化进行的预测,对既含有已知信息又含有不确定信息的系统进行的预测,也就是对在一定范围内变化的、与时间序列有关的灰过程进行 预测。尽管灰过程中所显示的现象是随机的、杂乱无章的,但毕竟是有序的、有界的,因此得到的数据集合具备潜在的规律。灰色预测是利用这种规律建立灰色模型对灰色系统进行预测。 目前使用最广泛的灰色预测模型就是关于数列预测的一个变量、一阶微分的GM(1,1)模型。它是基于随机的原始时间序列,经按时间累加后所形成的新的时间序列呈现的规律可用一阶线性微分方程的解来逼近。经证明,经一阶线性微分方程的解逼近所揭示的原始时间序列呈指数变化规律。因此,当原始时间序列隐含着指数变化规律时,灰色模型GM(1,1)的预测是非常成功的。 1 灰色系统的模型GM(1,1) 1.1 GM(1,1)的一般形式 设有变量X (0)={X (0)(i),i=1,2,...,n}为某一预测对象的非负单调原始数据列,为建立灰色预测模型:首先对X (0)进行一次累加(1—AGO, Acumulated Generating Operator)生成一次累加序列: X (1)={X (1)(k ),k =1,2,…,n} 其中 X (1)(k )= ∑ =k i 1 X (0)(i) =X (1)(k -1)+ X (0)(k ) (1) 对X (1)可建立下述白化形式的微分方程: dt dX )1(十) 1(aX =u (2) 即GM(1,1)模型。 上述白化微分方程的解为(离散响应): ∧ X (1)(k +1)=(X (0)(1)- a u )ak e -+a u (3) 或 ∧ X (1)(k )=(X (0)(1)- a u ))1(--k a e +a u (4)
MATLAB实现灰色预测程序 function [y,p,e]=huise_1_1(X,k) %灰色模型的malab程序 %Example [y,p]=gm_1_1([200 250 300 350],2) %接口描述:X的预测的初始数列,|X|>4,K是指向后进行预测的个数 %命令格式:程序保存的文件名,eg:huise.m 则命令是: huise([579.8 547.5 527.0 492.3 437.0],5) if nargout>3; error('Too many output argument.'); end if nargin==1,k=1;x_orig=X; elseif nargin==0|nargin>2 error('Wrong number of input arguments.'); end x_orig=X; predict=k; %AGO 处理,即是对初始数列进行一阶累加x=cumsum(x_orig); %计算系数(a 和u)------------------------ n=length(x_orig); %生成矩阵B for i=1:(n-1); B(i)=-(x(i)+x(i+1))/2; end B=[B' ones(n-1,1)]; %生成矩阵Y for i=1:(n-1); y(i)=x_orig(i+1); end Y=y'; %计算系数a=au(1) u=au(2) au=(inv(B'*B))*(B'*Y); %-------------------------------------------------------- %把huise模型公式转换成符号 coef1=au(2)/au(1); coef2=x_orig(1)-coef1; coef3=0-au(1); costr1=num2str(coef1); costr2=num2str(abs(coef2)); costr3=num2str(coef3); eq=strcat(costr1,'+',costr2,'e^',costr3,'*(t-1))'); %计算每一个值 for t=1:(n+predict) mcv(t)=coef1+coef2*exp(coef3*(t-1)); end x_mcv0=diff(mcv); x_mcve=[x_orig(1) x_mcv0] %输出图形中的各点
%该程序用于灰色关联分析,其中原始数据的第一行为参考序列,1至15行为正相关序列,16至17为负相关序列 clc,clear load x.txt %把原始数据存放在纯文本文件x.txt 中 %如果全为正相关序列,则将两个循环替换为下列代码 %for i=1:size(x,1) %x(i,=x(i,/x(i,1); %end for i=1:15 x(i,=x(i,:)/x(i,1); %标准化数据 end for i=16:17 x(i,:)=x(i,1)./x(i,:); %标准化数据 end data=x; n=size(data,1); ck=data(1,:);%分离参考序列 bj=data(2:n,:);m1=size(bj,1); for j=1:m1 t(j,:)=bj(j,:)-ck; end jc1=min(min(abs(t')));jc2=max(max(abs(t'))); rho=0.5;%灰色关联度为0.5 ksi=(jc1+rho*jc2)./(abs(t)+rho*jc2); r=sum(ksi')/size(ksi,2); r %灰色关联度向量 [rs,rind]=sort(r,'descend') %对关联度进行降序排序 %该函数用于灰色预测模型,其中x0为列向量,alpha一般取0.5,将第一个数据视为序号为0,k从0开始的序号矩阵 function y=huiseyuce(x0,alpha,k) n=length(x0); x1=cumsum(x0); for i=2:n z1(i)=alpha*x1(i)+(1-alpha)*x1(i-1); end z1=z1'; B=[-z1(2:n),ones(n-1,1)]; Y=x0(2:n); ab=B\Y; y1=(x0(1)-ab(2)/ab(1))*exp(-ab(1)*k)+ab(2)/ab(1);%产生预测累加生成序列 y=[x0(1) diff(y1)]%产生灰色预测数据 1 / 1
灰色系统预测GM(1,1)模型及其Matlab实现 三天三夜72小时: 读懂题目-》查找文献资料-》选择题目-》重查找文献资料-》精读其中几篇-》查找资料的资料。。。。 在数学建模中常常会遇到数据的预测问题,有些赛题中,预测占主导地位,例如: 2003年A题 SARS的传播问题; 2005年A题长江水质的评价和预测问题; 2006年B题艾滋病疗法的评价及疗效的预测问题; 2007年A题中国人口增长预测问题。 有些问题则是需要在求解的过程中进行预测,如2009年D题“会议筹备”对与会人数的确定等。 参考资料: 《灰色系统理论及其应用第五版》作者:刘思峰,党耀国等著出版时间:2010.05 校超星数字图书馆可阅读。 灰色模型(Gray Model)有严格的理论基础,最大优点是实用。用灰色模型预测的结果比较稳定,不仅适用于大数据量的预测,在数据量较少时(>3)预测结果依然较准确。 预备知识 (1)灰色系统 白色系统是指系统内部特征是完全已知的,即人们不仅知道该系统的输入——输出关系,而且知道实现输入——输出关系的结构与过
程;黑色系统是指系统内部信息完全未知的,即人们只知道该系统输入——输出关系,但不知道实现输入——输出关系的结构与过程;而灰色系统是介于白色系统和黑色系统之间的一种系统,灰色系统其内部一部分信息已知,另一部分信息未知或不确定。 例如,一个加有电压的电阻,也是一个系统,根据欧姆定律,I=U/R,当电阻的大小知道后,便可由多大电压算出能得到多大电流。电压与电流之间有明确的关系或函数,这便是白色系统。因此,这样的系统要求有明确的作用原理,一个有明确作用原理的系统必定是具有确定结构的,必定是有物理原型的。然而许多社会经济系统都没有物理原型,虽然知道影响系统的某些因素,但很难明确全部因素,更不可能确定因素之间的映射关系。这种没有确定的映射关系(函数关系)的系统是灰色系统。 (2)灰色预测 灰色预测,是指对系统行为特征值的发展变化进行的预测,对既含有已知信息又含有不确定信息的系统进行的预测,也就是对在一定范围内变化的、与时间序列有关的灰过程进行预测。尽管灰过程中所显示的现象是随机的、杂乱无章的,但毕竟是有序的、有界的,因此得到的数据集合具备潜在的规律。灰色预测是利用这种规律建立灰色模型对灰色系统进行预测。 目前使用最广泛的灰色预测模型就是关于数列预测的一个变量、一阶微分的GM(1,1)模型。它是基于随机的原始时间序列,经按时间累加后所形成的新的时间序列呈现的规律可用一阶线性微分方程的解来逼近。经证明,经一阶线性微分方程的解逼近所揭示的原始
灰色预测模型Matlab代码 sr(1) = 11.985; sr(2) = 12.1121; sr(3) = 12.2389; sr(4) = 12.3626; sr(5) = 12.481; x1=[0 0 0 0 0]; for j = 1 : 5 for i = 1 : j x1(j) = x1(j)+sr(i); end end for k = 2 : 5 z1(k) = 0.5 * (x1(k) + x1(k - 1)); end B=[-z1(2) 1;-z1(3) 1;-z1(4) 1;-z1(5) 1]; Yn=[sr(2);sr(3);sr(4);sr(5)]; ajg=inv(B'*B)*B'*Yn; a=ajg(1,1); b=ajg(2,1); for i = 0 : 100 sc1(i+1) = (sr(1) - b / a) * exp(-(a*i)) + b / a; end for j = 1 : 100 sc0(j + 1) = sc1(j + 1) - sc1(j); end sc0 sc1; 等维递补灰色预测模型Matlab代码 format compact sr(1) = 11.985; sr(2) = 12.1121; sr(3) = 12.2389; sr(4) = 12.3626; sr(5) = 12.481; p=0; while p<=100 x1=[0 0 0 0 0]; for j = 1 : 5 for i = 1 : j x1(j) = x1(j)+sr(i); end end for k = 2 : 5
§12.5 灰色预测 我们通常所说的系统是指:由客观世界中相同或相似的事物和因素按一定的秩序相互关联、相互制约而构成的一个整体.例如:工程技术系统、社会系统、经济系统等.如果一个系统中具有充足的信息量,其发展变化的规律明显、定量描述方便、结构与参数具体,则这种系统通常称为白色系统.如果一个系统的内部特征全部是未知的,则称此系统为黑色系统.如果系统内部信息和特征是部分已知的,另一部分是未知的,这种系统称为灰色系统.例如:社会系统、农业系统、经济系统、气象系统、生物系统等.对于这类系统,内部因素难以辨识,相互之间的关系较为隐蔽,人们难以准确了解这类系统的行为特征.因此,对于这类问题进行定量描述,即建立模型难度较大.区别白色系统与灰色系统的重要标志是系统内各因素之间是否具有确定的关系. 灰色系统分析方法主要是根据具体灰色系统的行为特征数据,充分利用数量不多的数据和信息寻求相关因素自身与各因素之间的数学关系,建立相应的数学模型.目前,灰色系统理论在实际中已得到了广泛的应用,例如:在工程技术、经济管理、气象预报以及政治、社会、工业、农业等领域都取得了一定的应用成果. 我们往往要对农业问题、商业问题等做未来的预测工作,另外,进行军事战争以及治理生态环境也需对未来的发展情形做一可靠的分析,这就产生了灰色预测.灰色预测是对灰色系统问题进行未来的预测,实际问题中,应用最多的灰色预测模型是以GM(1,1)(即GM(1,N)当N=1时的特例)模型为基础的. 12.5.1 GM(1,1)模型的建立 设X(0)=(X(0)(1),X(0)(2),…,X(0)(n)),做1-AGO,得
(1)(1)(1)(1)((1),(2),,()) X X X X n =L (1)(1)(0)(1)(0)((1),(1)(2),,(1)())X X X X n X n =+-+L 则GM(1,1)模型相应的微分方程为: (1) (1)dX aX u dt += (1) 式中:a 称为发展灰数;μ称为内生控制灰数. 设?α =(a ,μ)T ,按最小二乘法得到 11?()T T B B B Y α -= (2) 其中 (1)(1) (1)(1)(1) (1) 1((1)(2))121((2)(3))12 1((1)())12X X X X B X n X n ??-+ ? ? ?-+ ?= ? ? ?--+ ??? M M (0)(0)1(0)(2)(3)()X X Y X n ?? ? ? = ? ? ??? M 易求得,方程(1)的解为 (1)(0)?(1)((1))ak u u X k X e a a -+=-+ (3) 例4 100m 成绩预测 1983~1990年世界男子和中国女子100m 最好成绩如表6. 表6 各年度最好成绩