灰色预测模型的Matlab程序及检验程序(精)
灰色系统G(1,1)预测步骤【模板带代码】
=3499.075e -0.1062t
3641.075
编写程序
u=alpha(2)/alpha(1) v=X0(1)-u v=3499.075 u=—3641.075
(5)进行参差检验
1）根据预测公式，计算
v=3499.075 u=—3641.075
Xˆ
1
k
1
X
0
1
for n=0:10
X2(n+1)=v*exp(-alpha(1)*n)+u
end
X2
2.0690
2）累减生成序列
Xˆ X3 =1.0e+003 * (0) 0.1420 0.4079 0.4536
0.5044
0.7713 0.8577 0.9537 1.0605
源程序：X3(1)=X2(1)
for m=1:10
kesi =
4.4388 339.0664 176.2445 203.6132
0 0.1998 1.2682 0.2130 0.8524 0.1330
0.0089
0.3767
0.2203
0.4155
{0%，19.98%，126.82%，0.89%，37.67% ，22.03% ，41.55% ，21.30%，85.24%，13.30%}
e=
179.4592 111.5134 74.1747 175.0204 159.6072 29.2461 215.2168 33.1910
3.2147 24.1540
源程序：S0=0.6745*X0std e=abs(daita0-daita0mean) 对所有的 e 都小于 S0 ，故小参差概率 P(k S0) 1 0.95
灰色预测MATLAB程序
灰色预测作用：求累加数列、求a b的值、求预测方程、求残差clc %清屏，以使结果独立显示x=[71.1 72.4 72.4 72.1 71.4 72.0 71.6];format long; %设置计算精度if length(x(:,1))==1 %对输入矩阵进行判断，如不是一维列矩阵，进行转置变换 x=x';endn=length(x); %取输入数据的样本量z=0;for i=1:n %计算累加值，并将值赋予矩阵bez=z+x(i,:);be(i,:)=z;endfor i=2:n %对原始数列平行移位y(i-1,:)=x(i,:);endfor i=1:n-1 %计算数据矩阵B的第一列数据c(i,:)=-0.5*(be(i,:)+be(i+1,:));endfor j=1:n-1 %计算数据矩阵B的第二列数据e(j,:)=1;endfor i=1:n-1 %构造数据矩阵BB(i,1)=c(i,:);B(i,2)=e(i,:);endalpha=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,: );%显示输出预测值的累加数列endvar(1,:)=ago(1,:) %显示输出预测值for i=1:n %如改n为n+m-1，可预测后m-1个值var(i+1,:)=ago(i+1,:)-ago(i,:); %估计值的累加数列的还原，并计算出下一预测值endfor i=1:nerror(i,:)=x(i,:)-var(i,:); %计算残差endc=std(error)/std(x); %调用统计工具箱的标准差函数计算后验差的比值c ago %显示输出预测值的累加数列alpha %显示输出参数数列var %显示输出预测值error %显示输出误差c %显示后验差的比值作用：数据处理判断是否可以用灰色预测、求级比、求累加数列、求a b的值、求预测方程clc,clearx0=[71.1 72.4 72.4 72.1 71.4 72.0 71.6]'; %注意这里为列向量n=length(x0);lamda=x0(1:n-1)./x0(2:n) %计算级比range=minmax(lamda') %计算级比的范围x1=cumsum(x0) %累加运算B=[-0.5*(x1(1:n-1)+x1(2:n)),ones(n-1,1)];Y=x0(2:n);u=B\Y %拟合参数u(1)=a,u(2)=bx=dsolve('Dx+a*x=b','x(0)=x0'); %求微分方程的符号解x=subs(x,{'a','b','x0'},{u(1),u(2),x0(1)}) %代入估计参数值和初始值yuce1=subs(x,'t',[0:n-1]); %求已知数据的预测值y=vpa(x,6) %其中的6表示显示6位数字yuce=[x0(1),diff(yuce1)] %差分运算，还原数据。
灰色预测MATLAB程序
作用：求累加数列、求a b的值、求预测方程、求残差clc %清屏，以使结果独立显示x=[ ];format long; %设置计算精度if length(x(:,1))==1 %对输入矩阵进行判断，如不是一维列矩阵，进行转置变换 x=x';endn=length(x); %取输入数据的样本量z=0;for i=1:n %计算累加值，并将值赋予矩阵bez=z+x(i,:);be(i,:)=z;endfor i=2:n %对原始数列平行移位y(i-1,:)=x(i,:);endfor i=1:n-1 %计算数据矩阵B的第一列数据c(i,:)=*(be(i,:)+be(i+1,:));endfor j=1:n-1 %计算数据矩阵B的第二列数据e(j,:)=1;endfor i=1:n-1 %构造数据矩阵BB(i,1)=c(i,:);B(i,2)=e(i,:);endalpha=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,: );%显示输出预测值的累加数列endvar(1,:)=ago(1,:) %显示输出预测值for i=1:n %如改n为n+m-1，可预测后m-1个值var(i+1,:)=ago(i+1,:)-ago(i,:); %估计值的累加数列的还原，并计算出下一预测值endfor i=1:nerror(i,:)=x(i,:)-var(i,:); %计算残差endc=std(error)/std(x); %调用统计工具箱的标准差函数计算后验差的比值c ago %显示输出预测值的累加数列alpha %显示输出参数数列var %显示输出预测值error %显示输出误差c %显示后验差的比值作用：数据处理判断是否可以用灰色预测、求级比、求累加数列、求a b的值、求预测方程clc,clearx0=[ ]'; %注意这里为列向量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)=bx=dsolve('Dx+a*x=b','x(0)=x0'); %求微分方程的符号解x=subs(x,{'a','b','x0'},{u(1),u(2),x0(1)}) %代入估计参数值和初始值yuce1=subs(x,'t',[0:n-1]); %求已知数据的预测值y=vpa(x,6) %其中的6表示显示6位数字yuce=[x0(1),diff(yuce1)] %差分运算，还原数据。
MATLAB_+_灰色预测程序,数学建模
MATLAB实现灰色预测程序灰色预测很好的东西呐，······~~··`~··~~~~~~~~~~~~~~~~~~~~~````````````` fon [feval,au,ec,C,P]=GM1_1(x, r)if nrgin<2myar=0;end[mx,nx]=size(x);if mx==1x=x';endn=length(x);for i=2:nz(i-1)=0.5*x1(i)+0.5*x1(i-1);endY=x(2:end);B(:,1)=-z;2)/au(1));yc(1)=x(1);for k=1:n+myear-1y1(k+1)=pm*exp(-au*k)+a(2)/au(1);yc(k+1)=y1(k+1)-y1(k);endfeval=yc';ex=ec./x;r=0;rou=0.5;for k=1:nr=r+rou* s(ec(k))+rou*max(a (ec))); endr=r/n;%%====%原始序列的标准差s1=std(x);%计算残差的标准差s2=std(ec);%计算CC=s2/s1;%计算后验概率deta=ec-mean(ec);index=fineta)<0.6745*s1);P=length(index)/n;%%if C<0.35&P>0.95disp('预测精度为一级')elsP>0.8disp('预测精度为二级')elseif >0.7disp('预测精度为三级')elsedisp('预测精度过低，需要对模型进行修正') endif r>0.6disp('关联度符合检验要求')end%%%%=========t1=1:length(x);t2=1:lengt);plot(t1,x,'b--+',t2,feval,'r-o')legend('原始数据','预测数据')另一个程序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.3437.0],5)if nargout>3;r('Too maoutput argument.');enif nargin==1,k=1;x_orig=X;elseif ==0|nargin>2errr('Wrong nu arguments.');endx_rig=X;predict=k; %AGO 处理，即是对初始数列进行一阶累加x=cumsum(x_orig); %计算系数(a 和u)------------------------n=leh(x_orig); %生成矩阵Bfor i=1:(n-1);B(i)=-(x(i)+x(i+1))/2;enB=[B' ones(n-1,1)]; %生成矩阵Yfor i=1:(n-1);y(i)=x_ori(i+1);edY=y'; %计算系数a=au(1) u=au(2) au=(inv(B'*B))*(B'*Y); %--------------------------------------------------------%把huise模型公式转换成符号coef1=au(2)/au(1);coef2=x_or (1)-coef1;co3=0-au(1);costr1=nm2str(coef1);costr2=numstr(abs(coef2));costr3=ntr(coef3);eq=strcat(ctr1,'+',costr2,'e^',costr3,'*(t-1))'); %计算每一个值for t=1:(n+predict)mcv(t)=co1+coef2*exp(coef3*(t-1));endx_mcv0=diff(mcv);x_mcve=[x_orig(1) x_mcv0] %输出图形中的各点x_c_error=x_orig_n-x_mcv;x_errr=mn(abs(x_c_error./x_orig_n));if x_error>0.2 %相对误差的均值disp('del disqualification!');elseif x_error>0.1dip('model check out');disp('model is perfect!');endplot(1:n,x_orig,'o',1:n+predict,x_mcve);p=x_mcve(end-predict+1:end); %画出预测模型和初始数列的点xlabel('年份(从第一个数据年份起)');ylabel('产水量（万吨）');tie('灰度模型GM(1,1)');grid ony=eq;e=x_error;p=x_mcve(end-predict+1:end);。
灰色预测模型matlab程序精确版
%x=[1019,1088,1324,1408,1601];gm1(x); 测试数据%二次拟合预测GM(1,1)模型function gmcal=gm1(x)if nargin==0x=[1019,1088,1324,1408,1601]endformat long gsizex=length(x);%求数组长度k=0;for y1=xk=k+1;if k>1x1(k)=x1(k-1)+x(k);%累加生成z1(k-1)=-0.5*(x1(k)+x1(k-1));%z1维数减1，用于计算Byn1(k-1)=x(k);elsex1(k)=x(k);endend%x1,z1,k,yn1sizez1=length(z1);%size(yn1);z2 = z1';z3 = ones(1,sizez1)';YN = yn1'; %转置%YNB=[z2 z3];au0=inv(B'*B)*B'*YN;au = au0';%B,au0,auafor = 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,rightbrastrcat(x1t1,'=',num2str(constant1),estr,leftbra,num2str(afor1),tstr,rightbra,'+ ',leftbra,num2str(ua),rightbra)%输出时间响应方程%******************************************************%二次拟合k2 = 0;for y2 = x1k2 = k2 + 1;if k2 > kelseze1(k2) = exp(-(k2-1)*afor);endend%ze1sizeze1=length(ze1);z4 = ones(1,sizeze1)';G=[ze1' z4];X1 = x1';au20=inv(G'*G)*G'*X1;au2 = au20';%z4,X1,G,au20Aval = 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;%预测的步骤数 1for k3=1:nfinalx3fcast(k3) = constant1*exp(afor1*k3)+ua;end%x3fcast%一次拟合累加值for k31=nfinal:-1:0if k31>1x31fcast(k31+1) = x3fcast(k31)-x3fcast(k31-1);elseif k31>0x31fcast(k31+1) = x3fcast(k31)-x(1);elsex31fcast(k31+1) = x(1);endendendx31fcast%一次拟合预测值for k4=1:nfinalx4fcast(k4) = Aval*exp(afor1*k4)+Bval;end%x4fcastfor k41=nfinal:-1:0if k41>1x41fcast(k41+1) = x4fcast(k41)-x4fcast(k41-1);elseif k41>0x41fcast(k41+1) = x4fcast(k41)-x(1);elsex41fcast(k41+1) = x(1);endendendx41fcast,x%二次拟合预测值%***精度检验p C************////////////////////////////////// k5 = 0;for y5 = xk5 = k5 + 1;if k5 > sizexelseerr1(k5) = x(k5) - x41fcast(k5);endend%err1%绝对误差xavg = mean(x);%xavg%x平均值err1avg = mean(err1);%err1avg%err1平均值k5 = 0;s1total = 0 ;for y5 = xk5 = k5 + 1;if k5 > sizexelses1total = s1total + (x(k5) - xavg)^2;endends1suqare = s1total ./ sizex;s1sqrt = sqrt(s1suqare);%s1suqare,s1sqrt%s1suqare 残差数列x的方差 s1sqrt 为x方差的平方根S1 k5 = 0;s2total = 0 ;for y5 = xk5 = k5 + 1;if k5 > sizexelses2total = s2total + (err1(k5) - err1avg)^2; endends2suqare = s2total ./ sizex;%s2suqare 残差数列err1的方差S2Cval = sqrt(s2suqare ./ s1suqare);Cval%nnn = 0.6745 * s1sqrt%Cval C检验值k5 = 0;pnum = 0 ;for y5 = xk5 = k5 + 1;if abs( err1(k5) - err1avg ) < 0.6745 * s1sqrtpnum = pnum + 1;%ppp = abs( err1(k5) - err1avg )elseendendpval = pnum ./ sizex;pval%p检验值%arr1 = x41fcast(1:6)%预测结果为区间范围预测步长和数据长度可调整程序参数进行改进运行结果x =1019 1088 1324 1408 1601ans =x1(t+1)=8908.4929exp(0.11871t)+(-7889.4929)ans =x1(t+1)=8945.2933exp(0.11871t)+(-7935.7685)x31fcast =Columns 1 through 31019 1122.89347857097 1264.43142178303 Columns 4 through 61423.80987235488 1603.27758207442 1805.36675232556 x41fcast =Columns 1 through 31019 1118.05685435129 1269.65470492098 Columns 4 through 61429.69153740195 1609.90061644041 1812.82460377782 x =1019 1088 1324 1408 1601Cval =0.139501578334155 pval =1。
灰色预测MATLAB程序
灰色预测专设工⑼他QA—叫吋)为原始数列.其1次累❖加生成数列为恥=妙①曲⑵，…卅何),其中X° 仇)二工* ° (0.址=1=2= -:n5-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「—ay=代⑶ B =Ib*- :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"・如果対所有的|p 伙)|<0・1 -则认为达到较高的要求,否则 若旳所有的|。
灰色预测MATLAB程序
灰色预测心设尹曲⑴#为原始数列，其1次累＜加生成数列为炉=(孝①宀2\S,其中©=2^°:⑺卫=12…止i-1尋定文沙的灰导数为d(Jt)=玄㈣(Jt)=尤⑴的-工⑴(*-1).令尹为数列壬⑴的邻值生成数列，即尹)(町=加小(町十(1—a)x山(k-1).于是定文GM(1T1)的灰微分方程模型为d(k)+az①(上)=&_即或.严⑹+盘⑴懐)=乩⑴在式(1)中口①的称为灰导数’熬称为发展系数'弧称为白化背景值，b称为灰作用量。
将时刻表庄=23…用代入(O式有j<0)(2)-az⑴(2)=工®⑶—俺叫巧=»于是GMIL)樫型可表示为r=现在问题归结为求巧h在值。
用一元绒性回归，即最小二垂進求它们的估计值住=[]卜护跖护F奕厢上回归分析中求诂计值是用软件计算的，有标淮程博求解，如山訥甜等。
GM(1.1)的白化型对于的(1-1)的获微分方程⑴，如果将解导教矿悶的时報=%…屮观対连续叢里"则工⑴衩为时间i函敕卅®,于是-<'W耐应于导敕重级必％),白化背杲值刃(時对应于导數申⑴。
于是GM(1，1)的换微分方嗨对应于的白微分方程为写®4曲％「)=也⑵GAI(1>1)换色预刪的步叢1-數堀的椅噓弓处理为了保证©M(B1)屋複方达的可行性・需要対已却皴堀锁必要的检峻处Ho 设療皓数攥列为了-计算埶列的级比如果所有的级比都落在可容覆盖区间盂-內・则數摒列X糾可咲建立G*ICL-1)複型且可以避行页色预测。
否则，丙軌据懺适当的叢换处理，如平移銮换：取C使得敕培列严⑹二工蚀盘)+匚用二12…”的级比都落在可啓禎盖内。
(1)残差檢验：计算相对薙差Z 建立GM (L T 1)複型不妬设少弋以m 叫唠霸足上面的要求，以它芮議堀列建立GM(1>1)型蛊(仍(i)+血C1\A)=b ・用回归分祈求得目上的估计值"于是相应的白化模型为 气^十小卄工解为工叱)=0)①—勺中1-色-⑶ 应Q于是停到预测值壬⑴(上+1)=0叫1)一勺>加+仝血二12…卫一1=aa伙而相应地得到预«=x co \t +1)=x 0)(t+l)-x a)(i)3i =1,2,-?n-l ?如果对所有的^<0.1・则认为达到鞭嵩的要求：否则，若耐所有的|^)1<0^，则认対达到一般要求©(2)级比偏差値桧验：计算能)=1-呂学©如果对所有的|，则认为达列较高的要求孑吾则若对斫有的，则认为达到一般要求O灰色预测计算实例^…;=:=-■■■■昏例北方某城市1986—1992年道路交通噪声平均声级数据见表6序号年吶寺表拆市近年来交通噪声数据[眶(应)]二諾;二319S872.4第—爭:级比检验建立丢通噪屛均声级数锯时间序列如下：4198972.1j 1990?1.4 619?17201199771.6艸=(•严①卫购(2)厂卅⑺) =(711,72.4.71.4,72.1.71.4,7UQ.71.6)些(1)求级比k(k)忠防护住T)2=(几⑵山⑶.…也⑺)g=(0.982JJ.0042J.0098-0.9917J.0056)(2)级比判断由于所有的X.(10e[0.982J.009S],k=2,3.6故可以用双0)作满意的GM(1,1)建模’第二步：GM(1,1)建模(1)对原始数据X®作一次累加，即卞⑴=(71.L143.5215.9.288359.4.431.4,503)(2)构造数据矩阵B及数据向量Y-2)—H 弋3/>1⑶讦算1T心求解得F'⑴=(工倒〔1〉_-)e 弋Q f+-1*^+1)=0＜l,U)--)£-t +-=-3092^--^+31000⑶求生咸数列值歸型齊看:n令“is 那血由上面的碉醯数可甲得，其中取菱由龙⑴(i}=恥壮曲5加得丁I —"炉閃=进悶-进德-尊(71儿72.4.72.2:72.1:71.9:71.7,71.6)^}=(s"a ＞亍⑴⑵，…，网⑺A＜第三步;模型检验•＞模型的各种检验指标值的计算结果见表工 •t*表7GM(1检验表＜序号年俯原始值模型值残差相对误差级比偏差•>1 19S6 71.1 71.1<219S7 72.4 72.4 -0.0057 0.01%0.0023 <3 19S S 72.4 72.2 0.163S 0.23%0.0203 •>4 19S9 72.1 72.1 0.0329 0.05%-O.(K H8 •>5199071.4 71.9 -0-49S4 0.7%-0.0074 <61991 72.0 71.7 0.21599 037%0.0107<71992 71.6 71.6 0.037S0.05%-0.0032于是得到目=山的餡,立欖型7-B)'1B TY=(dt0.0023 72.6573dt+0.002ix (1>=72.657^心经验证・该模型的精度较高.可进行预测和预报计算的Matlab 程序如下:仃坝测和预报n=length(x); z=0;%取输入数据的样本量for i=1:nz=z+x(i,:)be(i,:)=z; %计算累加值，并将值赋予矩阵beend for i=2:n %对y(i-1,:)=x(i,:)%对原始数列平行移位 endfor i=1：n-1%计算数据矩阵B 的第一列数据c(i,:)=-0.5*(be(i,:)+be(i+1,:)); clCjdearxO=[71H 72.472A 72J71477m c n.lengthtxO)；*'b%注意这里为列帖lamda =xD(l ：n-1),A0(2：n)%计算级比range =minmaxflamda f )%计算级比的范阖 X1=cumsum(xO)；%累加运算B=['0,5*(xl(l ；n ^l)+xl(2:n))t ones(n -1,1)]TY 二甸(2：町；口=B\Y%拟合参数u(l>=a .u(2)=bx=dsolve (+a 'x =b\f x(0)-xO^J ；%求徴分方程的特号解x =subs(xJ*a\,b r /xO ，Mu(l)P u(2)t xO(l)|)i%代入荷计痹擞值和初蜡值yucel =subs %求巳知数擁的扳测位y-vpa(x,6)奄其中的石表示显不白位数字yuce=[x0(l)T diff(yucel)]%羔分运算，还原数据 epsiIon=-yuce%计算战羞作用：求累加数列、求ab 的值、求预测方程、求残差clc %清屏，以使结果独立显示x=[71.172.472.472.171.472.071.6]; format long ;%设置计算精度if length(x(:,1))==1%对输入矩阵进行判断，如不是一维列矩阵，进行转置变换x=x endM.I-JTVorhlllst 模型endfor j=1:n-1%计算数据矩阵B的第二列数据e(j,:)=1;endfor i=1:n-1%构造数据矩阵BB(i,1)=c(i,:);B(i,2)=e(i,:);endalpha=inv(B'*B)*B'*y;%计算参数矩阵即ab的值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,:);%显示输出预测值的累加数列endvar(1,:)=ago(1,: )for i=1:n%显示输出预测值%如改n为n+m-1，可预测后m-1个值var(i+1,:)=ago(i+1,:)-ago(i,:);%估计值的累加数列的还原，并计算出下一预测值endfor i=1:nerror(i,:)=x(i,:)-var(i,:);%计算残差endc=std(error)/std(x);%调用统计工具箱的标准差函数计算后验差的比值cago alpha var%显示输出预测值的累加数列%显示输出参数数列%显示输出预测值error %显示输出误差c %显示后验差的比值作用：数据处理判断是否可以用灰色预测、求级比、求累加数列、求ab的值、求预测方程clc,clearx0=[71.172.472.472.171.472.071.6]';%注意这里为列向量n=length(x0);lamda=x0(1:n-1)./x0(2:n)%计算级比range=minmax(lamda')%计算级比的范围x1=cumsum(x0)%累加运算B=[-0.5*(x1(1:n-1)+x1(2:n)),ones(n-1,1)];Y=x0(2:n);u=B\Y%拟合参数u(1)=a,u(2)=bx=dsolve('Dx+a*x=b','x(0)=x0');%求微分方程的符号解x=subs(x,{'a','b','x0'},{u(1),u(2),x0(1)})%代入估计参数值和初始值yuce1=subs(x,'t',[0:n-1]);%求已知数据的预测值y=vpa(x,6)%其中的6表示显示6位数字yuce=[x0(1),diff(yuce1)]%差分运算，还原数据。
灰色系统预测GM(1-1)模型及其Matlab实现
灰色系统预测GM(1,1)模型及其Matlab 实现预备知识（1）灰色系统白色系统是指系统内部特征是完全已知的;黑色系统是指系统内部信息完全未知的;而灰色系统是介于白色系统和黑色系统之间的一种系统，灰色系统其内部一部分信息已知,另一部分信息未知或不确定。
(2)灰色预测 灰色预测，是指对系统行为特征值的发展变化进行的预测，对既含有已知信息又含有不确定信息的系统进行的预测,也就是对在一定范围内变化的、与时间序列有关的灰过程进行 预测。
尽管灰过程中所显示的现象是随机的、杂乱无章的，但毕竟是有序的、有界的，因此得到的数据集合具备潜在的规律。
灰色预测是利用这种规律建立灰色模型对灰色系统进行预测。
目前使用最广泛的灰色预测模型就是关于数列预测的一个变量、一阶微分的ＧＭ（1，1)模型。
它是基于随机的原始时间序列，经按时间累加后所形成的新的时间序列呈现的规律可用一阶线性微分方程的解来逼近。
经证明，经一阶线性微分方程的解逼近所揭示的原始时间序列呈指数变化规律。
因此，当原始时间序列隐含着指数变化规律时，灰色模型GM(1,1)的预测是非常成功的。
1 灰色系统的模型GM(1,１）１．1 GM(1，1)的一般形式设有变量X (0)={X (0)(i ）,i ＝1,2,．..，n}为某一预测对象的非负单调原始数据列，为建立灰色预测模型:首先对X（0)进行一次累加（１—ＡGO ， Acum ｕl ａt ｅｄ Ge ｎeｒa ｔｉng Opera ｔo ｒ)生成一次累加序列: X (1)={X（1）(k ),k =1，2,…,n}其中X (1）(k )＝∑=ki 1X (0）(i)＝X (1)(ｋ-１）+ X （0）(ｋ) (1)对X（1)可建立下述白化形式的微分方程：dtdX )1(十)1(aX =u (2)即G Ｍ（１，1)模型。
上述白化微分方程的解为(离散响应): ∧X（1)(k +１)＝(X (０)(1)－a u ）ake -+au (3）或∧X （1）（k )=(X (0）(1)-a u ))1(--k a e +au （４)式中：k 为时间序列，可取年、季或月。
灰色理论预测模型及GM(1,1)matlab程序
灰色理论预测模型及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:nyy(i)=yy(i-1)+y(i);endB=ones(n-1,2);for i=1:(n-1)B(i,1)=-(yy(i)+yy(i+1))/2;B(i,2)=1;endBT=B';for j=1:n-1YN(j)=y(j+1);endYN=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:2ys(j)=yys(j)-yys(j-1);endx=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:ndet=det+abs(yn(i)-y(i));enddet=det/(n-1);disp(['百分绝对误差为：',num2str(det),'%']); disp(['预测值为：',num2str(ys(n+1:n+t_test))]);。
灰色系统预测GM1,1模型及其Matlab实现
灰色系统预测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)灰色预测灰色预测，是指对系统行为特征值的发展变化进行的预测，对既含有已知信息又含有不确定信息的系统进行的预测，也就是对在一定范围内变化的、与时间序列有关的灰过程进行预测。
灰色预测的matlab实现,代码大全
disp('确定的模型为:x1(k+1)=(x(1)-u_a)*exp(-k*a)+u_a') 改进灰色模型代码： for i=1:inf if i==1
e(i)=input('输入数据:'); else if e(i-1)==0 break; else end
e(i)=input('输入数据:');
6
1
x2(k+1)=(x0(1)-u_a)*exp(-k*a)+u_a; end %求出新数列的模拟值 x2 x3=x2; x3(m2+k0)=[]; x4=[0 x3]; x5=x2-x4;%利用累减生成法求出原始数据的模拟值 x5 disp('3.一次累加（1-AGO）生成的数据的模拟值:') Y=''; for z=1:m2+k0 Y=strcat(Y,'(',num2str(x2(z)),')'); end disp(Y); disp('4.原始数据的模拟值:') Y=''; for z=1:m2+k0 Y=strcat(Y,'(',num2str(x5(z)),')'); end disp(Y); % 4. 模型检验（算出的值到等级参照表中检查其精度等级） %计算后验差比 C 和残差序列 Q x6=x5(1:m2); Q=x0-x6;%Q 为残差序列 s1=std(Q);%s1 为残差序列 Q 的标准差 s2=std(x0);%s2 为初始序列 x0 的标准差 C=s1/s2;%后验差比 C(越小越好) w1=1:m2; w1=[ones(m2,1) w1']; w2=Q'; [bb,bint,r1,rint,stats]=regress(w2,w1); rcoplot(r1,rint) C1=strcat('5.后验差比(均方差比值): C=',num2str(C));disp(C1); if C<=0.35 disp(' 由于 C<=0.35，则此模型精度等级为 1 级（好） 。'); else if C<=0.5 disp(' 由于 0.35<C<=0.5，则此模型精度等级为 2 级（合格） 。'); else if C<=0.65 disp(' 由于 0.5<C<=0.65，则此模型精度等级为 3 级（勉强） 。'); else disp(' 由于 C>0.65，则此模型精度等级为 4 级（不合格） 。'); end end end
灰色预测模型实验以及例题分析实验报告
1.0691 1.2371
预测数据与实际数据的比较图如下：
2.P3-2（P34）
结果如下：
a=-0.0624
b =156.6162
预测数据与实际数据的比较图如下：
3.预测与会代表人数”（P35）
结果如下：
G=
Columns 1 through 9
0.1810 0.1980 0.1735 0.1520 0.1331 0.1166 0.1021 0.0895 0.0784
实验名称
灰色预测模型实验
实验目的
用MATLAB编程验证程序P3-1（P32）及P3-2（P34），并用其解决“预测与会代表人数”（P35）问题。
学会使用matlab编程做灰色预测实验。
实验内容（算法、程序、步骤和方法）
实验内容（算法、程序、步骤和方法）
实验内容（算பைடு நூலகம்、程序、步骤和方法）
1.P3-1（P32）
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); %预测后续数据
F=[];F(1)=A(1);
for i=2:(n+10)
F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a;
end
G=[];G(1)=A(1);
for i=2:(n+10)
G(i)=F(i)-F(i-1); %得到预测出来的数据
end
t1=1:4;
t2=1:14;
G
plot(t1,A,'o',t2,G) %原始数据与预测数据的比较
灰色预测及MATLAB实现
（3）对累加生成数据做均值生成 B 矩阵与常数项向量Yn ，即
0.5(x(1) (1) x(1) (2))
B
0.5(
x(1)
(2)
x(1)
(3))
，Yn
(x(0) (2), x(0) (3),
, x(0) (n))T
0.5(x(1) (n 1) x(1) (n))
3.2 灰色预测的MATLAB程序
3.2.1 典型程序结构
（1）对原始数据进行累加。
矩阵处理， MATLAB的长
处
（2）构造累加矩阵B与常数向量。
（3）求解灰参数。
（4）将参数代入预测模型进行数据预测。
【例】某公司1999-2008年利润为（元/年）：[89677 99215 109655 120333 135823 159878 182321 209407 246619 300670]， 预测该公司未来几年的利润情况。
已知本届会议的回执情况（表1），往几届会议代表回执和 与会情况（表2），根据这些数据预测本届与会代表。
表1 回执中对住房的要求
要求 男 女
合住1 154 78
合住2 104 48
合住3 32 17
独住1 107 59
独住2 68 28
独住3 41 19
表2 以往几届代表的回执参会情况表
届次
第一届
第二届
dx
由于aˆ 是通过最小二乘法求出的近似值，因此 xˆ(1) (t 1)事近似表达
式，与原序列区分，多了一个“帽子”。
（6）对函数表达式 xˆ(1) (t 1)及 xˆ(1) (t)进行离散，将二者作差以便还
灰色预测及MATLAB实现
3.1灰色预测基础知识
什么是灰色预测？
灰色预测是就灰色系统所做的预测。所谓灰色系统是介于白 色系统和黑箱系统之间的过渡系统，其具体的含义是：如果某一 系统的全部信息已知为白色系统，全部信息未知为黑箱系统，部 分信息已知，部分信息未知，那么这一系统就是灰色系统。一般 地说，社会系统、经济系统、生态系统都是灰色系统。例如物价 系统，导致物价上涨的因素很多，但已知的却不多，因此对物价 这一灰色系统的预测可以用灰色预测方法。 灰色系统理论认为对既含有已知信息又含有未知或非确定信 息的系统进行预测，就是对在一定方位内变化的、与时间有关的 灰色过程的预测。尽管过程中所显示的现象是随机的、杂乱无章 的，但毕竟是有序的、有界的，因此这一数据集合具备潜在的规 律，灰色预测就是利用这种规律建立灰色模型对灰色系统进行预 测。
ˆ （4）用最小二乘法求解灰参数 a ( a, ) ( B B )
T T
1
B Yn 。
T
ˆ （5）将灰参数 a 代入
ˆ x dy
(1)
dx
(1)
ax
(1)
，求解得
(1)
dt
(t 1) ( x
(0)
a
)e
at
a
dx ˆ ˆ (1) 由于 a 是通过最小二乘法求出的近似值，因此 x (t 1) 事近似表达
(1)
ˆ 序列，得到近似数据序列 x ˆ x
(0)
(0)
ˆ (t 1) x
(1)
ˆ (t 1) x
(t )
（7）建立灰色预测模型进行检验，步骤如下：
① 计算 x
(0)
与x
(0)
(t ) 之间的残差和相对误差 (t ) x
灰色理论灰色预测模型和灰色关联度分析matlab通用代码(精)
%该程序用于灰色关联分析，其中原始数据的第一行为参考序列，1至15行为正相关序列，16至17为负相关序列clc,clearload x.txt %把原始数据存放在纯文本文件x.txt 中%如果全为正相关序列，则将两个循环替换为下列代码%for i=1:size(x,1%x(i,=x(i,/x(i,1;%endfor i=1:15x(i,=x(i,:/x(i,1; %标准化数据endfor i=16:17x(i,:=x(i,1./x(i,:; %标准化数据enddata=x;n=size(data,1;ck=data(1,:;%分离参考序列bj=data(2:n,:;m1=size(bj,1;for j=1:m1t(j,:=bj(j,:-ck;endjc1=min(min(abs(t';jc2=max(max(abs(t';rho=0.5;%灰色关联度为0.5ksi=(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,kn=length(x0;x1=cumsum(x0;for i=2:nz1(i=alpha*x1(i+(1-alpha*x1(i-1;endz1=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]%产生灰色预测数据。
灰色预测例题matlab求解(精)
建立GM（1，1）模型对产品销售额预测祁诗阳冯晓凯申静某大型企业1999年至2004年的产品销售额如下表，试建立GM(1,1预测模型，并预测2005年的产品销售额。
年份1999 2000 2001 2002 2003 2004销售额2.673.13 3.25 3.36 3.56 3.72(亿元有题目知构造累加生成序列对作紧邻均值生成编程如下：x=[2.67 5.8 9.05 12.41 15.97 19.69];z(1=x(1;for i=2:6z(i=0.5*(x(i+x(i-1;endformat long gz结果如下：z =Columns 1 through 42.67 4.235 7.425 10.73Columns 5 through 614.19 17.83因此于是构造B矩阵和Y矩阵如下：对参数进行最小二乘估计，采用matlab编程完成解答如下：B=[[-4.235 -7.425 -10.73 -14.19 -17.83]',ones(5,1];Y=[3.13 3.25 3.36 3.56 3.72]';format long ga=inv(B'*B*B'*Y结果如下：a =-0.04396098154759662.92561659879905即=-0.044，u=2.96 =-66.55则GM(1,1白化方程为预测模型为：1、关联度检验法：采用matlab编程得到模拟序列for i=1:6X(i=69.22*exp(0.044*(i-1-66.55;endformat long gx(1=X(1;for i=2:6x(i=X(i-X(i-1;endX结果如下：x =Columns 1 through 42.673.11367860537808 3.25373920141375 3.40010005288617 Columns 5 through 63.55304456012134 3.71286887145915因此模拟序列为求模拟序列和原始序列的相关度初始化，即将该序列所有数据分别除以第一个数据。