2022年2022年灰色预测法-编程 .pdf

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

1、灰色预测法一、相关知识1、灰色预测通过原始数据的处理和灰色模型的建立，发现、掌握系统发展规律，对系统的未来状态做出科学的定量预测。2、灰数简介：（1）灰数的定义：是指未明确指定的数，即处在某一范围内的数，灰数是区间数的一种推广。灰数实际上指在某一个区间或某个一般的数集内取值的不确定数，通常用记号“”表示灰数。（2）灰数的分类：( ) 有下界而无上界的灰数,a或a ，如大树的重量必大于零，但不可能用一般手段知道其准确的重量, 所以其重量为灰数,0。( ) 有上界而无下界的灰数(,a或()a，如一项投资工程，要有个最高投资限额，一件电器设备要有个承受电压或通过电流的最高临界值。( ) 既有下界 a

2、 又有上界 a 的灰数称为区间灰数，记为aa,。如海豹的重量在20-25公斤之间，某人的身高在1.8-1.9米之间，可分别记为25,201，9.1,8.12( ) 黑数：当,或21,，即当的上、下界皆为无穷或上、下界都是灰数时，称为黑数。( ) 白数：当,a a且aa时，称为白数。(3) 本征灰数是指不能或暂时还不能找到一个白数作为其“代表”的灰数，比如一般的事前预测值、宇宙的总能量、准确到秒或微妙的“年龄”等都是本征灰数。非本征灰数是指凭先验信息或某种手段，可以找到一个白数作为其“代表” 的灰数。 我们称此白数为相应灰数的白化值，记为，并用a 表示以 a 为白化值的灰数。如托人代买一件价格1

3、00 元左右的衣服，可将100 作为预购衣服价格100的白化数，记为100100。例： （1）气温不超过36，36,0。（2）预计某地区今年夏粮产量在100 万吨以上，,100；（3）估计某储蓄所年底居民存款总额将达7000 万到 9000 万，9000,7000；（4）如某人希望至少获得1 万元科研经费，并且越多越好，,10000；名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 1 页，共 11 页 - - - - - - - - - （5）有的数，从系统的高层次，即宏观层次、整体

4、层次或认识的概括层次上看是白的，可到低层次上， 即到系统的微观层次、分部层次或认识的深化层次则可能是灰的。例如，一个人的身高，以厘米度量是白的，若精确到万分之一毫米就成灰的了。3、灰数白化与灰度（1）如 今 年 的 科 研 经 费 在5 万 元 左 右 ， 可 表 示 为5000050000， 或,50000,50000，它的白化值为50000。（2）对于一般的区间灰数ba,，我们将白化值取为：ba)1(，1 ,0一般灰色系统之行为特征预测值构成的灰数，就难以给出其白化权函数。定义形如ba)1(，1 ,0的白化称为等权白化。定义在等权白化中，取21而得到的白化值称为等权均值白化。当区间灰数取值

5、的分布信息缺乏时，常采用等权均值白化。（3）灰度即为灰数的测度。灰数的灰度在一定程度上反映了人们对灰色系统之行为特征的未知程度。 如果考虑一个4000 左右的灰数， 给出其估计值的两个灰数4002,39981和4100,39002，显然1比2更有价值，亦即1比2灰度小，若再考虑一个基本值为 4 的灰数， 给出灰数6,23，虽然1与3的长度都是4，但1比3的灰度小是显而易见的。不确定量量化（用确定量的方法研究）灰色系统视不确定量为灰色量。提出了灰色系统建模的具体数学方法，它能利用时间序列来确定微分方程的参数。灰色预测不是把观测到的数据序列视为一个随机过程，而是看作随时间变化的灰色量或灰色过程，通

6、过累加生成和累减生成逐步使灰色量白化，从而建立相应于微分方程解的模型并做出预报。二、灰色预测模型1、GM（1,1）模型令( 0)X为 GM(1,1) 建模序列，( 0)(0 )( 0)(0 )(1),(2),.,()Xxxxn，(1)X为( 0)X的 1-AGO （一次累加生成）序列，(1)(1)(1)(1)(1),(2),.,()Xxxxn，名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 2 页，共 11 页 - - - - - - - - - (1)(0)1()( )kixkxi

7、,1,2,.,kn令(1)Z为(1)X的紧邻均值（ MEAN ）生成序列(1)(1)(1)(1)(2),(3),.,()Zzzzn)()1(kz=0.5)()1(kx+0.5)1()1(kx即定义： GM(1,1) 的灰微分方程模型为bkazkx)()()1()0(式中 a 称为发展系数，b为灰色作用量。设?为待估参数向量，即?(,)Ta b，则灰微分方程的最小二乘估计参数列满足=nTTYBBB1)(其中B=(1)(1)(1)(2)1(3)1.( )1zzzn，nY=( 0)( 0)( 0)(2 )(3).()xxxn定义：(1)(1)dxaxbdt为灰色微分方程bkazkx)()()1()0

8、(的白化方程，也叫影子方程。如上所述，则有1） 白化方程(1)(1)dxaxbdt的解也称时间响应函数为(1)(1)?( )(0)atbbxtxeaa2） GM(1,1) 灰色微分方程bkazkx)()()1()0(的时间响应序列为?(1)xk(1)(0)bxaake+ba，1,2,.,kn3） 取(1)( 0)(0)(1)xx，则(1)?(1)xk( 0)(1)bxaake+ba，1,2,.,kn4） 还原值名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 3 页，共 11 页 -

9、- - - - - - - - ( 0)?(1)xk(1)?(1)xk(1)?()xk上式即为预测方程。5）模型检验：灰色预测模型检验有残差检验，关联度检验，后验差检验。（1）残差检验：计算原始序列和原始序列的灰色预测序列之间的：绝对误差( 0)( 0)(0 )?( )( )( )ixixi；1,2,.,in；相对误差：( 0)(0 )(0)(0)?( )( )( )( )xixiixi；1, 2,.,kn；其中(0 )?( )xi(1)?( )xi(1)?(1)xi1,2,.,in。相对误差越小，模型精度越高。（2）后验差检验：首先计算原始序列( 0)( )xi的均方差：2001SSn，而2

10、2( 0)(0 )01( )niSxix，(0)( 0)11( )nixxin。然后计算 残差 序列(0 )( ) i的均方差：2111SSn，而22( 0)(0)11( )niSi，( 0)(0 )11( )niin。再计算方差比10ScS，最后计算小误差概率(0 )(0 )00.6745pS。根据下面的预测精度等级划分表确定模型的精度。预测精度等级划分表小误差概率p值方差比 c 值预测精度等级0.950.35好0.800.5合格0.700.65勉强合格0.700.65不合格2、GM （1,1 ）模型应用实例例 某大型企业1999 年至 2004 年的产品销售额如下表，试建立GM(1,1)

11、预测模型，并预测2005 年的产品销售额。名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 4 页，共 11 页 - - - - - - - - - 年份1999 2000 2001 2002 2003 2004 销售 额 ( 亿元) 2.67 3.13 3.25 3.36 3.56 3.72 解：设( 0)()Xk=2.67 ，3.13 ， 3.25 ，3.36 ，3.56 ，3.72 第 1 步构造累加生成序列(1)()Xk=2.67,5.80,9.05,12.4，15.97,19

12、.69 第 2 步构造数据矩阵B和数据向量nY(1)(1)(1)(1)(1)(1)(1)(1)(1)(1)1(1)(2)1214.2351(2)(3)127.4251110.731(3)(4)1214.1911(4)(5)117.83121(5)(6)12xxxxBxxxxxx，(0)( 0)(0)( 0)(0)3.13(2)3.25(3)3.36(4)3.56(5)3.72(6)nxxYxxx第 3 步计算 ? =ab=nTTYBBB1)(BBT=541.5441.5446375.7071)(BBT=226382.1094319.0094319.0008667.0名师资料总结 - - -精品

13、资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 5 页，共 11 页 - - - - - - - - - ? =nTTYBBB1)(=0.0438792.925663第 4 步得出预测模型(1)dxdt0.043879(1)x=2.925663 (1)?(1)xk=69.34570.043879ke66.6757 （(0 )(1)x=2.67 ；ba=66.6757 ）第 5 步残差检验(1) 根据预测公式，计算(1)?()Xk，得(1)?()Xk 2.67 ，5.78 ，9.03 ，12.43 ，15.97

14、，19.68 ，19.69 (k=0,1, , ,6) (2) 累减生成( 0)?()Xk序列，k=1,2, , ,6 (0 )?()Xk 2.67 ， 3.11 ，3.25 ，3.40 ，3.54 ，3.71 原始序列：( 0)()Xk 2.67 ，3.13 ，3.25 ，3.36 ，3.56 ，3.72 (3) 计算绝对残差和相对残差序列绝对残差序列：( 0) 0，0.02 ，0， 0.04 ，0.02 ，0.01 相对残差序列： 0，0.64%，0，1.19%，0.56%，0.27%相对残差不超过1.19%，模型精确度高。预测：k=7，(0)x(8)=(1)x(8)(1)x(7)=4.2

15、3 3、GM （1,1 ）模型参数估计、检验及作图及预测的matlab 程序。（1）程序解读% GM （1,1 ）模型计算及检验、作图。文件名：fungry1.m （注释部分不能在MATLAB 窗口中运行）function GM1=fungry1(x0) % 输入原始数据x0 T=input(T=); % 从键盘输入最后一个历史数据算起的第T 时点名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 6 页，共 11 页 - - - - - - - - - % x1 为累加生成序列；B为(

16、1)(1)(1)(2)1(3)1.()1zzzn，初值为0；yn 为(0 )(0 )(0 )(2)(3).()xxxn初值为 0；% Hatx1：用灰色模型(1)?(1)xk(0 )(1)bxaake+ba，1,2,.,kn算出的序列；% Hatx0：灰色预测累减后的预测序列；( 0)?( )xi(1)?( )xi(1)?(1)xi1,2,.,in% Epsilon ：绝对误差(0 )(0)( 0)?( )( )( )ixixi；1,2,.,in% omega：相对误差：(0)( 0 )(0 )(0)?( )( )( )( )xixiixi；1,2,.,knx1=zeros(1,length(

17、x0);B=zeros(length(x0)-1,2); yn=zeros(length(x0)-1,1);Hatx0=zeros(1,length(x0)+T); Hatx00=zeros(1,length(x0);Hatx1=zeros(1,length(x0)+T); epsilon=zeros(length(x0),1);omega=zeros(length(x0),1); for i=1:length(x0) % 累加生成序列 for j=1:i x1(i)=x1(i)+x0(j); end end for i=1:length(x0)-1 % 构造数据矩阵B和数据向量nY B(i,

18、1)=(-1/2)*(x1(i)+x1(i+1); B(i,2)=1; yn(i)=x0(i+1); end HatA=(inv(B*B)*B*yn % 参数估计，求 a,b,HatA=1?(,)()TTTna bBBB Yfor k=1:length(x0)+T Hatx1(k)=(x0(1)-HatA(2)/HatA(1)*exp(-HatA(1)*(k-1)+HatA(2)/HatA(1); % (1)?()xk( 0)(1)bxa(1)ake+ba。end Hatx0(1)=Hatx1(1); % 灰色模型算出序列和累减生成序列的第一个值相等为x0(1) for k=2:length(

19、x0)+T Hatx0(k)=Hatx1(k)-Hatx1(k-1); % 累减还原得到原始序列的模拟值；end for i=1:length(x0) epsilon(i)=x0(i)-Hatx0(i); %计算相对误差 omega(i)=(epsilon(i)/x0(i); %计算绝对误差end 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 7 页，共 11 页 - - - - - - - - - x0,Hatx1,Hatx0,epsilon,omega, %显示算出的一些序列c

20、=std(epsilon)/std(x0); %开始模型检验。Std：均方差。 std(x0)：原始序列的均方差 p=0; % 即 P 为绝对误差序列小的概率：( 0)( 0)00.6745pS2001SSn，22(0)(0 )01( )niSxix，(0)( 0)11( )nixxin。for i=1:length(x0) if abs(epsilon(i)-mean(epsilon)0.95 & c0.85 & c0.70 & c0.65 disp(The model is not good and the forecast is :), disp (Hatx0(length(x0)+T)

21、 else p0.65 disp(The model is bad and try again ) end for i=1:length(x0) Hatx00(i)= Hatx0(i); % 得到模拟值点 end z=1:length(x0); plot(z,x0,-,z,Hatx00,:) % 将原始序列和模拟值画在一个图上比较 text(2,x0(2),History data :real line) %在（ 2,x0(2)）位置标注 text(length(x0)/2,Hatx00(length(x0)/2,Simulation data :broken line) end （2）GM

22、（1,1 ）模型参数估计、检验及作图及预测的matlab 程序。文件名：fungry1.mfunction GM1=fungry1(x0) T=input(T=); x1=zeros(1,length(x0);B=zeros(length(x0)-1,2); yn=zeros(length(x0)-1,1);Hatx0=zeros(1,length(x0)+T); Hatx00=zeros(1,length(x0);Hatx1=zeros(1,length(x0)+T); epsilon=zeros(length(x0),1);omega=zeros(length(x0),1); 名师资料总结

23、 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 8 页，共 11 页 - - - - - - - - - for i=1:length(x0) for j=1:i x1(i)=x1(i)+x0(j); end end for i=1:length(x0)-1 B(i,1)=(-1/2)*(x1(i)+x1(i+1); B(i,2)=1; yn(i)=x0(i+1); end HatA=(inv(B*B)*B*yn for k=1:length(x0)+T Hatx1(k)=(x0(1)-HatA

24、(2)/HatA(1)*exp(-HatA(1)*(k-1)+HatA(2)/HatA(1); end Hatx0(1)=Hatx1(1); for k=2:length(x0)+T Hatx0(k)=Hatx1(k)-Hatx1(k-1); end for i=1:length(x0) epsilon(i)=x0(i)-Hatx0(i); omega(i)=(epsilon(i)/x0(i)*100; end %x0,Hatx1,Hatx0,epsilon,omega, c=std(epsilon)/std(x0); p=0; for i=1:length(x0) if abs(epsilo

25、n(i)-mean(epsilon)0.95 & c0.85 & c0.70 & c0.65 disp(The model is not good and the forecast is :), disp (Hatx0(length(x0)+T) else p0.65 disp(The model is bad and try again ) end for i=1:length(x0) 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 9 页，共 11 页 - - - - - - -

26、 - - 11.522.533.544.555.5633.13.23.33.43.53.63.73.8History data :real lineSimulation data :broken line Hatx00(i)= Hatx0(i); end z=1:length(x0); plot(z,x0,-,z,Hatx00,:) text(2,x0(2),History data :real line) text(length(x0)/2,Hatx00(length(x0)/2,Simulation data :broken line) end 在 MATLAB 窗口中输入： clear

27、x0=2.97 3.23 3.29 3.46 3.59 3.71; fungry1(x0); 提示输入： T=1 HatA = % 灰色预测模型系数 -0.0365 3.0412 x0 = % 原始序列 2.9700 3.2300 3.2900 3.4600 3.5900 3.7100 Hatx1 = % 灰色生成序列2.9700 6.1779 9.5051 12.9560 16.5353 20.2478 24.0984 Hatx0 = % 灰色累减得到原始序列的模拟值序列 2.9700 3.2079 3.3272 3.4510 3.5793 3.7125 3.8506 epsilon = %

28、 绝对误差序列 0.0000 0.0221 -0.0372 0.0090 0.0107 -0.0025 omega = % 相对误差序列 0.0000 0.0069 -0.0113 0.0026 0.0030 -0.0007 p = 1 % 绝对误差小误差概率为1The model is good and the forecast is : 3.8506 % 这个模型很好，第7 年人口预测值为3.8506 。名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 10 页，共 11 页 - - - - - - - - - 3.8506 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 11 页，共 11 页 - - - - - - - - -