灰色预测法-编程

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

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

3、，记为 ，并用 a 表示以 a 为白化值的灰数。如托人代 买一件价格 100 元左右的衣服，可将 100 作为预购衣服价格 100 的白化数，记为100 100 。例：( 1)气温不超过 36，0,36 。2) 预计某地区今年夏粮产量在 100 万吨以上， 100, ；3) 估计某储蓄所年底居民存款总额将达7000万到 9000 万， 7000,9000 ；4) 如某人希望至少获得 1 万元科研经费，并且越多越好，10000, ；(5)有的数，从系统的高层次，即宏观层次、整体层次或认识的概括层次上看是白的，可 到低层次上，即到系统的微观层次、 分部层次或认识的深化层次则可能是灰的。例如， 一个

4、 人的身高，以厘米度量是白的，若精确到万分之一毫米就成灰的了。3、灰数白化与灰度(1) 如今年的 科研经 费在 5 万 元左右 ，可 表示为 50000 50000 ， 或50000 ,50000, ，它的白化值为 50000。(2) 对于一般的区间灰数 a,b ，我们将白化值 取为： a (1 ) b ，0,1一般灰色系统之行为特征预测值构成的灰数，就难以给出其白化权函数。定义 形如 a (1 )b ， 0,1 的白化称为等权白化。1定义 在等权白化中，取 而得到的白化值称为等权均值白化。2当区间灰数取值的分布信息缺乏时，常采用等权均值白化。( 3)灰度即为灰数的测度。灰数的灰度在一定程度上

5、反映了人们对灰色系统之行为特 征的未知程度。 如果考虑一个 4000 左右的灰数，给出其估计值的两个灰数1 3998,4002和 2 3900,4100 ，显然 1比 2更有价值，亦即 1比 2 灰度小，若再考虑一个基本 值为 4 的灰数， 给出灰数 3 2,6 ，虽然 1与 3的长度都是 4，但 1比 3 的灰度小是 显而易见的。不确定量 量化(用确定量的方法研究)灰色系统视不确定量为灰色量。提出了灰色系统建模的具体数学方法，它能利用时间序列 来确定微分方程的参数。灰色预测不是把观测到的数据序列视为一个随机过程，而是看作 随时间变化的灰色量或灰色过程，通过累加生成和累减生成逐步使灰色量白化，

6、从而建立 相应于微分方程解的模型并做出预报。二、灰色预测模型1、GM ( 1,1)模型令 X(0) 为 GM(1,1) 建模序列，X(0) (x(0) (1), x(0) (2),., x(0)(n)，X (1) 为 X (0) 的 1-AGO(一次累加生成)序列，(1) (1) (1) (1)X (x (1),x (2),., x (n) ，kx(1) ( k )x(0) (i), k 1,2,.,ni1令 Z (1) 为 X (1) 的紧邻均值( MEAN)生成序列Z(1) ( z(1) (2), z(1)(3),.,z(1)(n)z (1) ( k) =0.5 x (1) ( k ) +

7、0.5 x(1)(k 1)即定义：GM(1,1) 的灰微分方程模型为x(0) (k) az(1)(k) b式中 a称为发展系数， b为灰色作用量。设 ?为待估参数向量，即 ? (a,b)T ，则灰微分方程的最小二乘估计参数列满足=(BT B) 1BTY其中z(1)(2)x(0) (2)B=z(1)(3)(1)z (n)x(0) (3)x(0) (n)dx(1)定义： dxax(1) bdt为灰色微分方程 x(0)(k) az(1) (k) b 的白化方程，也叫影子方程。如上所述，则有dx(1)1) 白化方程 dx ax(1) b 的解也称时间响应函数为dtx?(1)(t) ( x(1) (0)

8、 b)e at b aa2) GM(1,1)灰色微分方程 x(0)(k) az (1) (k) b 的时间响应序列为x?(1) (k 1) x(1) (0) b e ak+b， k 1,2,.,n aa3) 取 x(1)(0) x(0) (1) ，则x?(1) ( k 1) x(0) (1) b e ak+b，k 1,2,., n aa4) 还原值x?(0) (k 1) x?(1) (k 1) x?(1)(k)上式即为预测方程。5) 模型检验：灰色预测模型检验有残差检验，关联度检验，后验差检验。(1)残差检验：计算原始序列和原始序列的灰色预测序列之间的：绝对误差(0) (i) x(0)(i)

9、x?(0) (i) ； i 1,2,., n ；相对误差： (i)x(0) (i)； k 1,2,.,n ；其中 x?(0) (i) x?(1) (i ) x?(1) (i 1) i 1,2,., n 。 相对误差越小，模型精度越高。2) 后验差检验：首先计算原始序列 x(0) ( i ) 的均方差：S0S02n1S0 ，而 S02x(0) (i) x(0)i1x(0)1x(0) (i) 。ni1然后计算 残差 序列 (0) ( i) 的均方差：S1S12 ，而 n1S12(0) (i) (0)i1(0)1n1 (0) (i) 。 ni1再计算方差比S1 c S0 ，最后计算小误差概率 p(0

10、) (0)0.6745 S0 。根据下面的预测精度等级划分表确定模型的精度。预测精度等级划分表小误差概率 p 值方差比 c 值预测精度等级0.950.35好0.800.5合格0.700.65勉强合格0.700.65不合格2、 GM( 1,1 )模型应用实例例 某大型企业 1999 年至 2004 年的产品销售额如下表，试建立 GM(1,1) 预测模型，并预测 2005 年的产品销售额。年份199920002001200220032004销 售额（亿 元）2.673.133.253.363.563.72解：设X (0) (k ) =2.67 ，3.13 ，3.25 ， 3.36 ，3.56 ，

11、3.72第 1 步 构造累加生成序列X (1) (k) =2.67,5.80,9.05,12.4， 15.97,19.69第 2 步 构造数据矩阵 B 和数据向量 Yn21 x(1) (1) x(1)(2)21 x(1) (2) x(1)(3)21 x(1) (3) x(1)(4)21 x(1) (4) x(1)(5)21 x(1) (5) x(1)(6)4.2357.42510.7314.1917.83x(0) (2)3.13x(0)(3)3.25x(0) (4)3.36x(0) (5)3.56x(0) (6)3.72Yn第 3 步 计算 ?= a =(BTB) 1BTY n bnBT B=

12、 707.4637554.4154.4150.0943191.226382(BTB)1= 0.0086670.094319?=(BTB) 1BTYn= 2.925663第 4 步 得出预测模型dx(1)dt0.043879 x(1) =2.925663(1) 0.043879kx? (k 1) =69.3457 e 66.6757( x(0) (1)=2.67 ； b =66.6757 )a第 5 步 残差检验(1) 根据预测公式，计算 X? (1) ( k) ，得 ,6)X?(1) ( k ) 2.67 ，5.78，9.03，12.43，15.97，19.68，19.69( k =0,1,(

13、2) 累减生成 X?(0)(k)序列， k =1,2, ,63.54 ，3.71 X?(0) (k) 2.67 ，3.11 ，3.25 ， 3.40 ，原始序列： X (0) (k) 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， x(0) (8)= x(1) (8) x(1) (7)=4.233、 GM( 1,1 )模型

14、参数估计、检验及作图及预测的matlab 程序。( 1)程序解读% GM(1,1 )模型计算及检验、作图。文件名：fungry1.m(注释部分不能在 MATLAB窗口中运行)function GM1=fungry1(x0)% 输入原始数据x0T=input(T=);% 从键盘输入最后一个历史数据算起的第T 时点z(1)(2)1x(0) (2)% x1 为累加生成序列； B为z(1)(3)1，初值为 0； yn 为x(0)(3)z(1)(n)1x(0)(n)初值为 0；% Hatx1 ：% Hatx0 ：用灰色模型 x?(1) (k 1) x(0) (1) b e ak + b ，aa灰色预测累

15、减后的预测序列； x?(0) (i) x?(1) (i ) x?(1)(i 1) i 1,2,., nk 1,2,., n 算出的序列；% Epsilon：绝对误差 (0)(i) x(0) (i) x?(0) (i ) ； i 1,2,., n% omega ：相对误差：(0) (i)x(0) (i) x?(0) (i )x(0)(i)k 1,2,., nx1=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,lengt

16、h(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:ix1(i)=x1(i)+x0(j);endendfor i=1:length(x0)-1 % 构造数据矩阵 B 和数据向量 YnB(i,1)=(-1/2)*(x1(i)+x1(i+1);B(i,2)=1;yn(i)=x0(i+1);endHatA=(inv(B*B)*B*yn% 参数估计， 求 a,b,HatA= ? (a,b)T(BT B) 1B

17、TYfor k=1:length(x0)+THatx1(k)=(x0(1)-HatA(2)/HatA(1)*exp(-HatA(1)*(k-1)+HatA(2)/HatA(1);% x?(1)(k) x(0) (1) b e a(k 1) + b 。aaendHatx0(1)=Hatx1(1); % 灰色模型算出序列和累减生成序列的第一个值相等为 x0(1) for k=2:length(x0)+THatx0(k)=Hatx1(k)-Hatx1(k-1);%累减还原得到原始序列的模拟值；end% 计算相对误差%计算绝对误差for i=1:length(x0)epsilon(i)=x0(i)-H

18、atx0(i);omega(i)=(epsilon(i)/x0(i);end x0,Hatx1,Hatx0,epsilon,omega, % 显示算出的一些序列c=std(epsilon)/std(x0); % 开始模型检验。 Std ：均方差。 std(x0) ：原始序列的均方差p=0;%即 P 为绝对误差序列小的概率： p (0) (0) 0.6745 S0S02 (0) (0) 2 (0) 1 (0)S00 ， S0x (i) x ， x x (i)。n 1 i 1 n i 1 for i=1:length(x0)if abs(epsilon(i)-mean(epsilon)0.95 &

19、 c0.85 & c0.70 & c0.65disp(The model is not good and the forecast is :),disp (Hatx0(length(x0)+T)else p0.65disp(The model is bad and try again )endfor i=1:length(x0)Hatx00(i)= Hatx0(i);% 得到模拟值点endz=1:length(x0);plot(z,x0,-,z,Hatx00,:) % 将原始序列和模拟值画在一个图上比较 text(2,x0(2),History data :real line)%在( 2,x0

20、(2) )位置标注text(length(x0)/2,Hatx00(length(x0)/2,Simulation data :broken line)end2) GM( 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=ze

21、ros(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);endendfor i=1:length(x0)-1B(i,1)=(-1/2)*(x1(i)+x1(i+1);B(i,2)=1;yn(i)=x0(i+1);endHatA=(inv(B*B)*B*ynfor k=1:length(x0)+THatx1(k)=(x0(1)-HatA(2)/HatA(1)*exp(-HatA(1)*(k-1)+HatA(2)/

22、HatA(1); endHatx0(1)=Hatx1(1);for k=2:length(x0)+THatx0(k)=Hatx1(k)-Hatx1(k-1);endfor 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(epsilon(i)-mean(epsilon)0.95 & c0.85 & c0.70 & c0.65disp(The model is not good and the forecast is :),disp (Hatx0(length(x0)+T)else p0.65disp(The model is bad and try again )endfor i=1:length(x0)Hatx00(i)= Hatx0(i);endz=1:le