Matlab教程课件-灰色预测模型
灰色预测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+灰色预测模型模型GM(1,1)
GM(1,1)灰色预测模型IntroductionInitial给定原始序列：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 Programclearsyms 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; endxx0=[];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 nY 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 ̂=(au)=(B T B)−1B T Y nStep 3模型响应函数x ̂(1)(k +1)=(x (0)(1)−u a )e −ak +u ax ̂(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:ne(i)=abs(x0(i)-xx0(i)); endmmax=max(e); for i=1:nee(i)=0.5*mmax/(e(i)+0.5*mmax); endr=sum(ee)/n; %后验差检验x0bar=sum(x0)/n; s1=0; for i=1:ns1=s1+(x0(i)-x0bar)^2; ends1=sqrt(s1/n); s2=0;ebar=sum(e)/n; for i=1:ns2=s2+(e(i)-ebar)^2; ends2=sqrt(s2/n); C=s2/s1; p=0;for i=1:nif 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)]2n绝对误差序列标准差：S2=√∑[ε(i)−ε̅]2n计算方差比：C=S2S1小误差概率：P=P{|ε(i)−ε̅|<0.6745S1}p=p+1;endendp=p/n;Cpif p>0.95&C<0.35disp('预测精度好');else if p>0.8&C<0.5disp('预测合格');else if p>0.7&C<0.65disp('预测勉强合格'); elsedisp('预测不合格'); endendend%原始数据与预测数据进行比较t1=1:n;t2=1:(n+10);xx0plot(t1,x0,'o',t2,xx0)。
利用matlab进行灰色预测
(k ) 1 (
1 0.5a ) ( k ) 1 0.5a
如果 (k ) 0.2 ，则可认为达到一般要求，如果 (k ) 0.1 ，则认 为达到较高要求。 （4）预测
ˆ (0) (k ) ，根据实际需要，求出问 由模型GM（1,1）得到的预测值 x
文章地址 /logs/163249158.html 2
, n) 。
,n
求出均值数列： z (1) (k ) 0.5x(1) (k ) 0.5x(1) (k 1), k 2,3, 由此得到 z (1) ( z (1) (2), z (1) (3), 建立灰微分方程：
x(0) (k ) az (1) (k ) b, k 2,3, ,n , z (1) (n))
(0)
GM(1,1)的数据进行灰色预测。否则，需要对数列 x 入可容覆盖内。即取适当的常数c ，作平移变换
做必要的变换处理，使其落
y (0) x(0) (k ) c ,
使数列 y
(0)
k 1, 2,
,n
( y (0) (1), y (0) (2),
, y (0) (n)) 的级比 y (k ) 落入可容覆盖区内。
文章地址 /logs/163249158.html 1
, n 1
利用 Matlab 进行灰色预测
2. 灰色预测的步骤
（1）数据的检验与处理 首先，为了保证建模方法的可行性，需要对已知数据列做必要的检验处理。设参考 数据为 x
（3）检验预测值 a．残差检验：残差实际值与预测值之差，令残差为 (k ) ，计算
ˆ (0) (k ) x(0) (k ) x (k ) , x(0) (k ) k 1, 2, ,n
灰色预测模型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 %GM(1,1).m%建⽴符号变量a(发展系数)和b(灰作⽤量)syms a b;c = [a b]';%原始数列 AA = [174, 179, 183, 189, 207, 234, 220.5, 256, 270, 285];%填⼊已有的数据列！n = length(A);%对原始数列 A 做累加得到数列 BB = cumsum(A);%对数列 B 做紧邻均值⽣成for i = 2:nC(i) = (B(i) + B(i - 1))/2;endC(1) = [];%构造数据矩阵B = [-C;ones(1,n-1)];Y = A; Y(1) = []; Y = Y';%使⽤最⼩⼆乘法计算参数 a(发展系数)和b(灰作⽤量)c = inv(B*B')*B*Y;c = c';a = c(1);b = c(2);%预测后续数据F = []; F(1) = A(1);for i = 2:(n+10) %这⾥10代表向后预测的数⽬，如果只预测⼀个的话为1F(i) = (A(1)-b/a)/exp(a*(i-1))+ b/a;end%对数列 F 累减还原,得到预测出的数据G = []; G(1) = A(1);for i = 2:(n+10) %10同上G(i) = F(i) - F(i-1); %得到预测出来的数据enddisp('预测数据为：');G%模型检验H = G(1:10); %这⾥的10是已有数据的个数%计算残差序列epsilon = A - H;%法⼀：相对残差Q检验%计算相对误差序列delta = abs(epsilon./A);%计算相对误差Qdisp('相对残差Q检验：')Q = mean(delta)%法⼆：⽅差⽐C检验disp('⽅差⽐C检验：')C = std(epsilon, 1)/std(A, 1)%法三：⼩误差概率P检验S1 = std(A, 1);tmp = find(abs(epsilon - mean(epsilon))< 0.6745 * S1);disp('⼩误差概率P检验：')P = length(tmp)/n%绘制曲线图t1 = 1995:2004;%⽤⾃⼰的，如1 2 3 4 5...t2 = 1995:2014;%⽤⾃⼰的，如1 2 3 4 5... plot(t1, A,'ro'); hold on;plot(t2, G, 'g-');xlabel('年份'); ylabel('污⽔量/亿吨');legend('实际污⽔排放量','预测污⽔排放量'); title('长江污⽔排放量增长曲线'); %都⽤⾃⼰的grid on;。
第三章-灰色预测模型PPT课件
由于 x (1 )
t
涉及到累加列 x (1)
的两个时刻的值，因此，x (1) (i)
取前后两个时刻的平均代替更为合理，即将 x (1 ) ( i ) 替换为
.
16
1[x(1)(i)x(1)(i1)],(i2,3,...,N ). 2
将（3.5）写为矩阵表达式
x(0)(2)
x(0)(3)
1212[[xx((11))((32))xx((11))((21))]]
表3.7 某市2001－2005年火灾数据
.
35
x (0)(N ) a x (1)(N ) u .
.
15
把ax(1) (i) 项移到右边，并写成向量的数量积形式
x
(0)(2 )[x(1) ( ), 1 ]au
x
(0)(3 )
[
x
(1) (3 ), 1 ]
a
u
x
(0) ( N
)
[
x
(1) ( N
), 1]
a
u
(3.5)
.
25
解（1）由原始数据列计算一次累加序列 x (1) ，结
果见表3.3.
表3.3 一次累加数据
x (0)
x (1)
.
26
（2）建立矩阵：B, y
B11122212[[[[xxxx((((1111))))((((5342))))xxxx((((1111))))((((2413))))]]]]
1 4.513
dx(1) ax(1) u dt
（3.1） （3.2） （3.3）
.
13
其中，a,u分别是待定的常数，a称为发展系（灰）数；u
称为灰色作用量(内生控制灰数)。此方程满足初始条件
灰色预测模型matlab程序精确版
%x=[1019,1088,1324,1408,1601];gm1(x);测试数据%二次拟合预测GM(1,1) 模型function gmcal=gm1(x)if nargin==0x=[1019,1088,1324,1408,1601]end format long g sizex=length(x); %求数组长度k=0;for y1=x k=k+1; if k>1x1(k)=x1(k-1)+x(k);% 累加生成z1(k-1)=-0.5*(x1(k)+x1(k-1));%z1 维数减1，用于计算B yn1(k-1)=x(k);elsex1(k)=x(k);end end %x1,z1,k,yn1 sizez1=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,rightb ra,'+ ',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)%输出时间响应方程for 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);endend endx31fcast %一次拟合预测值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);endendend%二次拟合预测值%***精度检验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);%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 =运行结果x =ans =x1(t+1)=8908.4929exp(0.11871t)+(-7889.4929) ans = x1(t+1)=8945.2933exp(0.11871t)+(-7935.7685) x31fcast =Columns 1 through 3Columns 4 through 61429.691537401951609.90061644041 1812.824603777821019 1088 1324 1408 16011019 1088 13241408 160110191122.89347857097 1264.43142178303 Columns 4 through 61423.80987235488 1603.27758207442 1805.36675232556x41fcast =Columns 1 through 310191118.05685435129 1269.65470492098Cval =0.139501578334155 pval =1。
灰色预测及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)进行离散，将二者作差以便还
灰色系统预测GM(1-1)模型及其Matlab实现教学教材
灰色系统预测G M(1-1)模型及其M a t l a b实现灰色系统预测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 )＝∑=ki 1X (0)(i)＝X (1)(k －1)+ X (0)(k ) (1) 对X (1)可建立下述白化形式的微分方程：dtdX )1(十)1(aX ＝u (2)即GM(1,1)模型。
上述白化微分方程的解为(离散响应)： ∧X (1)(k +1)＝(X (0)(1)－a u )ak e -＋au(3) 或∧X (1)(k )＝(X (0)(1)－a u ))1(--k a e ＋au (4) 式中：k 为时间序列，可取年、季或月。
灰色预测及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]%产生灰色预测数据。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
2
i3 =1
k = 4, x (1) (4) = ∑ x (0) (i ) = x (1) (3) + x (0) (4) = 9.9 + 3.6 = 13.5
i =1 5
k = 5, x (5) =
(1)
(0) (1) (0) = + x i x x ( ) (4) (5) = 13.5 + 3.8 = 17.3 ∑ i =1
对于原始数据列依次做前后相邻的两个数据相减的运 算过程称为累减生成过程IAGO。如果原始数据列为
x (1) = ( x (1) (1), x (1) ( 2), , x (1) ( n))
令 x ( 0 ) ( k ) = x (1) ( k ) − x (1) ( k − 1), k = 2,3, , n,
为均值生成数，也称等权邻值生成数
9
累加生成计算示例
例：x (0)=(x (0) (k) ︱k=1，2，3，4，5) =x(0)(1)，x(0)(2)，x(0)(3)，x(0)(4)，x(0)(5) =(3.2，3.3，3.4，3.6，3.8) 求 x(1)(k) 解： (1) (0) k = 1, x = (1) x= (1) 3.2
8
3. 加权邻值生成
设原始数列为x ( 0 ) = ( x ( 0 ) (1), x ( 0 ) (2), , x ( 0 ) (n))
称为数列的邻值。 x (k − 1), x (k )
(0) (0)
x
(0)
x (0) (k − 1)为后邻值，为前邻值 x (0) (k )
对于常数，令 α ∈ [0,1]
10
累加生成的特点 一般经济数列都是非负数列。累加生成能使任意非负 数列、摆动的与非摆动的，转化为非减的、递增的
原始数列作图 1—AGO作图
2
是未知的，系统内各因素间有不确定的关系。
2. 灰色预测法 • 灰色预测法是一种对含有不确定因素的系统进行预 测的方法。 • 灰色预测是对既含有已知信息又含有不确定信息 的系统进行预则，就是对在一定范围内变化的、与时 间有关的灰色过程进行预测。
• 灰色预测通过鉴别系统因素之间发展趋势的相异程度， 即进行关联分析，并对原始数据进行生成处理来寻找系 统变动的规律，生成有较强规律性的数据序列,然后建 立相应的微分方程模型，从而预测事物未来发展趋势的 状况。 • 灰色预测法用等时距观测到的反映预测对象特征 的一系列数量值构造灰色预测模型，预测未来某一 3 时刻的特征量，或达到某一特征量的时间。
4
• 拓扑预测
将原始数据做曲线，在曲线上按定值寻找该定值发 生的所有时点，并以该定值为框架构成时点数列，然 后建立模型预测该定值所发生的时点。
5
3.2 灰色生成数列 灰色系统理论认为，尽管客观表象复杂，但总是 有整体功能的，因此必然蕴含某种内在规律。 关键在于如何选择适当的方式去挖掘和利用它。 灰色系统是通过对原始数据的整理来寻求其变化 规律的，这是一种就数据寻求数据的现实规律的途 径，即为灰色序列的生成。 一切灰色序列都能通过某种生成弱化其随机性， 显现其规律性。数据生成的常用方式有累加生成、 累减生成和加权累加生成。
3. 灰色预测的四种常见类型 • 灰色时间序列预测 即用观察到的反映预测对象特征的时间序列来构造 灰色预测模型，预测未来某一时刻的特征量，或达 到某一特征量的时间。 • 畸变预测 即通过灰色模型预测异常值出现的时刻，预测异常 值什么时候出现在特定时区内。 • 系统预测 通过对系统行为特征指标建立一组相互关联的灰 色预测模型，预测系统中众多变量间的相互协调 关系的变化。
k = 2, x (1) (2) = ∑ x (0) (i ) = x (0) (1) + x (0) (2) = 3.2 + 3.3 = 6.5
k = 3, x (1) (3) = ∑ x (0) (i ) = x (1) (2) + x (0) (3) = 6.5 + 3.4 = 9.9
i =1 4
3. 灰色预测模型
3.1 灰色预测的概念 ； 3.2 灰色生成数列； 3.3 灰色模型GM；
1
3.1 灰色预测的概念 1. 灰色系统、白色系统和黑色系统 白色系统是指一个系统的内部特征是完全已知的， 即系统的信息是完全充分的。 黑色系统是指一个系统的内部信息对外界来说是一 无所知的，只能通过它与外界联系来加以 观测研究。 灰色系统内的一部分信息是已知的，另一部分信息
称所得到的新数列为数列
x 的
k
(0)
(1) (1) (1) (1) x = x x x ( ( 1 ), ( 2 ), , (n)) 1次累加生成数列
( r −1) x ( k ) = x x 的r次累加数列 ∑ (i), k = 1,2,, n, r ≥ 1
(0)
(r )
i =1
7
2. 累减生成
z (k ) = αx (k ) + (1 − α ) x (k − 1), k = 2,3,, n,
(0) (0) (0)
由此得到的数列称为数列在权下的邻值生成数， z (0) x (0) α 权称为生成系数。 α 特别地，当生成系数 α = 0.5 时，则称
z ( 0 ) (k ) = 0.5 x ( 0 ) (k ) + 0.5 x ( 0 ) (k − 1), k = 2,3, , n,
称所得到的数列 x ( 0 ) 为 x (1) 的1次累减生成数列。
注：从这里的记号也可以看到，从原始数列得 x (0) (1) 到新数列，再通过累减生成可以还原出原始数 x ˆ (1) 列。实际运用中在数列的基础上预测出 x (1) x 通 ˆ (0)。 过累减生成得到预测数列 x
6
1. 累加生成 把数列各项（时刻）数据依次累加的过程称为累加生 成过程（AGO ）。由累加生成过程所得的数列称为 累加生成数列。 设原始数列为
令 x (1) (k ) =
k i =1
x ( 0) = ( x ( 0) (1), x ( 0) (2),, x ( 0) (n))
(0) = x ∑ (i), k 1, 2, , ,