function hsyc(x0)
%灰色模型GM(1,1)matlab代码在建立灰色模型前应对原始数据进行求%%%级比%%%%重要
%输入原始时序列x0，输出各种所需参数%%%%实际使用该程序时应改变图像坐标轴%%%%
%x0代表原始数据；x1代表x0经累加后；B代表所够造的矩阵B；Y代表向量Yn；au代表系数a,u构成的矩阵
n=size(x0,2); %计算x0的列数
x1=cumsum(x0,2); %将x0进行累加得到矩阵x1
for k=1:n-1;
b(k)=-0.5*(x1(k+1)+x1(k));
y(k)=x0(k+1); %计算得到向量Y的转置矩阵y
end
a=ones(n-1,1);
B=[b',a]; %计算得到矩阵B
Y=y'; %计算得到向量Y
au=inv(B'*B)*B'*Y; %计算得到系数a，u构成的矩阵au
A=au';
k=1:n+10;
X1(k)=(x0(1)-A(2)/A(1))*exp(-A(1)*(k-1))+A(2)/A(1); %计算预测累加数列的值
X0(1)=X1(1);
k=1:n+9;
X0(k+1)=X1(k+1)-X1(k); %计算预测累加数列的还原值，即预测值
for k=2:n;
e(k)=X0(k)-x0(k); %计算残差
E(k)=(X0(k)-x0(k))/x0(k); %计算预测值与实测值的差值跟实测值的比值，即残差与实测值的比值
j(k)=x0(k-1)/x0(k); %计算实测数据的级比
end
k=1:n;
p(k)=1-(1-0.5*A(1))*j(k)/(1+0.5*A(1)); %计算级比偏差e(k)、p(k)<0.1则说明模型达到较高要求，0.2>e(k)、p(k)>0.1则说明模型达到一般要求
m=min(abs(e)); %计算残差中的最小值
M=max(abs(e)); %计算残差中的最大值
k=1:n;
g(k)=(m+0.5*M)./(abs(e(k))+0.5*M); %计算关联系数
R=sum(g')/(n-1); %计算关联度
v=[1,5,120,185];
axis(v);
grid on
plot([1996:2005],x0,'o-',[1996:2015],X0,'*:'); %画出原始数据跟时间序列的图像及预测数据跟时间序列的图像
legend('原始数据','预测数据',4)
au %输出参数a,u的值
X0 %输出预测数列
e %输出残差
E %输出残差与原始数据的比值，定义的残差j %输出原始数据的级比
p %输出偏比差
R %输出关联度程度
X1
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 各年度最好成绩