2022年2022年灰色模型Matlab程序 .pdf

《2022年2022年灰色模型Matlab程序 .pdf》由会员分享，可在线阅读，更多相关《2022年2022年灰色模型Matlab程序 .pdf（9页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、%by allen 红嘴海鸥%灰色模型预测是在数据不呈现一定规律下可以采取的一种建模和预测方法，其预测数据与原始数据存在一定的规律相似性%下面程序是灰色模型GM(1,1) 程序二次拟合和等维新陈代谢改进预测程序,matlab6.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; ifk1 x1(k)=x1(k-1)+x(k); % 累加生成z1(

2、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,yn1sizez1=size(z1,2); %size(yn1); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 1 页，共 9 页 - - - - - - - - - z2 = z1; z3 = ones(1,sizez1);YN = yn1; %转置%YNB=z2 z3; au0=inv(B*B)*B

3、*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,n

4、um2str(afor1),tstr,rightbra,+,leftbra,num2str(ua),rightbra) % 输出时间响应方程%* % 二次拟合名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 2 页，共 9 页 - - - - - - - - - k2 = 0; for y2 = x1 k2 = k2 +1; if k2 k else ze1(k2) =exp(-(k2-1)*afor); end end %ze1sizeze1 = size(ze1,2); z4 =

5、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,+,leftbra,num2str(Bval),rightbra) % 输出时间响应方程nfinal = sizexd2-1 + 1; % 决定预测的步骤数5 这个步骤可以通过函数传入%n

6、final = sizexd2 - 1 + 1; % 预测的步骤数 1名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 3 页，共 9 页 - - - - - - - - - for k3=1:nfinal x3fcast(k3)= constant1*exp(afor1*k3)+ua; end %x3fcast % 一次拟合累加值for k31=nfinal:-1:0 ifk311 x31fcast(k31+1) = x3fcast(k31)-x3fcast(k31-1); else

7、 if k310 x31fcast(k31+1) = x3fcast(k31)-x(1); else x31fcast(k31+1) = x(1); end end end x31fcast % 一次拟合预测值for k4=1:nfinal x4fcast(k4)= Aval*exp(afor1*k4)+Bval; end %x4fcastfor k41=nfinal:-1:0 ifk411 x41fcast(k41+1) = x4fcast(k41)-x4fcast(k41-1); else 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - -

8、 - - - 名师精心整理 - - - - - - - 第 4 页，共 9 页 - - - - - - - - - if k410 x41fcast(k41+1) = x4fcast(k41)-x(1); else x41fcast(k41+1) = x(1); end end end x41fcast,x %二次拟合预测值%* 精度检验 pC*/ k5 = 0; for y5 = x k5 = k5 +1; if k5 sizexd2 else err1(k5) = x(k5) -x41fcast(k5); end end %err1 %绝对误差xavg = mean(x); %xavg %

9、x 平均值err1avg = mean(err1); %err1avg %err1 平均值名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 5 页，共 9 页 - - - - - - - - - k5 = 0; s1total = 0 ; for y5 = x k5 = k5 +1; if k5 sizexd2 else s1total = s1total + (x(k5) -xavg)2; end end s1suqare = s1total ./ sizexd2; s1sqrt =

10、 sqrt(s1suqare); %s1suqare,s1sqrt %s1suqare 残差数列 x 的方差s1sqrt 为 x方差的平方根S1 k5 = 0; s2total = 0 ; for y5 = x k5 = k5 +1; if k5 sizexd2 else s2total = s2total + (err1(k5) -err1avg)2; end end s2suqare = s2total ./ sizexd2; %s2suqare 残差数列 err1 的方差 S2 Cval = sqrt(s2suqare ./ s1suqare); Cval %nnn = 0.6745 *

11、 s1sqrt %Cval C 检验值名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 6 页，共 9 页 - - - - - - - - - k5 = 0; pnum = 0 ; for y5 = x k5 = k5 +1; if abs(err1(k5) - err1avg ) 0.6745 * s1sqrt pnum = pnum + 1; %ppp = abs( err1(k5) - err1avg) else end end pval = pnum ./ sizexd2; p

12、val %p 检验值%arr1 = x41fcast(1:6) %预测结果为区间范围预测步长和数据长度可调整程序参数进行改进- 程序为原创，引用请注明名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 7 页，共 9 页 - - - - - - - - - gg=x1,x2,x3,x4,x5,zeros(1,m);%矩阵 gg, 其中x1,x2,x3,x4,x5是已知数据 , 前面说到灰色理论只要四个数据以上就可以建模, 所以这里用五个比较合理 ,zeros(1,m)中的 m是你所想得到

13、的预测数据的个数,由前面用于预测的已知数据是五个, 所以 m的值最好也在五个左右 . 在下面的语句中 , 实际操作的时候不要忘了把m替换成实际的数值 . result=zeros(1,m); for j=1:m g0=gg(:,j:j+4); n=5; x0=g0; x1=zeros(n-1,1); x1(1)=x0(1); for i=2:n, x1(i)=x1(i-1)+x0(i); end for k=1:n-1 z1(k)=(x1(k)+x1(k+1)/2; end z=zeros(1,n-1); for i=1:n-1 z(i)=-z1(i); end B=z;ones(1,n-1)

14、; B=B; y=zeros(n-1,1); for i=2:n y(i-1)=x0(i); end u=inv(B*B)*B*y; a=u(1); b=u(2); for k=0:5, x2(k+1)=(x0(1)-b/a)*exp(-a*k)+b/a; end 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 8 页，共 9 页 - - - - - - - - - x3(1)=x0(1); for k=1:5, x3(k+1)=x2(k+1)-x2(k); end g=x3; result(j)=g(n+1); gg(j+5)=g(n+1); end result; r=result %r即为预测得到的数据名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 9 页，共 9 页 - - - - - - - - -