灰色预测法-编程.doc

灰色预测法一、相关知识1、灰色预测通过原始数据的处理和灰色模型的建立，发现、掌握系统发展规律，对系统的未来状态做出科学的定量预测。2、灰数简介： （1）灰数的定义：是指未明确指定的数，即处在某一范围内的数，灰数是区间数的一种推广。灰数实际上指在某一个区间或某个一般的数集内取值的不确定数，通常用记号“”表示灰数。（2）灰数的分类：()有下界而无上界的灰数或，如大树的重量必大于零，但不可能用一般手段知道其准确的重量,所以其重量为灰数。()有上界而无下界的灰数或，如一项投资工程，要有个最高投资限额，一件电器设备要有个承受电压或通过电流的最高临界值。 ()既有下界又有上界的灰数称为区间灰数，记为。如海豹的重量在20-25公斤之间，某人的身高在1.8-1.9米之间，可分别记为，()黑数：当或，即当的上、下界皆为无穷或上、下界都是灰数时，称为黑数。()白数：当且时，称为白数。(3)本征灰数是指不能或暂时还不能找到一个白数作为其“代表”的灰数，比如一般的事前预测值、宇宙的总能量、准确到秒或微妙的“年龄”等都是本征灰数。非本征灰数是指凭先验信息或某种手段，可以找到一个白数作为其“代表”的灰数。我们称此白数为相应灰数的白化值，记为，并用表示以为白化值的灰数。如托人代买一件价格100元左右的衣服，可将100作为预购衣服价格的白化数，记为。例：（1）气温不超过36，。（2）预计某地区今年夏粮产量在100万吨以上，；（3）估计某储蓄所年底居民存款总额将达7000万到9000万，；（4）如某人希望至少获得1万元科研经费，并且越多越好，；（5）有的数，从系统的高层次，即宏观层次、整体层次或认识的概括层次上看是白的，可到低层次上，即到系统的微观层次、分部层次或认识的深化层次则可能是灰的。例如，一个人的身高，以厘米度量是白的，若精确到万分之一毫米就成灰的了。3、灰数白化与灰度（1） 如今年的科研经费在5万元左右，可表示为，或，它的白化值为50000。（2） 对于一般的区间灰数，我们将白化值取为：，一般灰色系统之行为特征预测值构成的灰数，就难以给出其白化权函数。定义 形如，的白化称为等权白化。定义 在等权白化中，取而得到的白化值称为等权均值白化。当区间灰数取值的分布信息缺乏时，常采用等权均值白化。（3）灰度即为灰数的测度。灰数的灰度在一定程度上反映了人们对灰色系统之行为特征的未知程度。如果考虑一个4000左右的灰数，给出其估计值的两个灰数和，显然比更有价值，亦即比灰度小，若再考虑一个基本值为4的灰数，给出灰数，虽然与的长度都是4，但比的灰度小是显而易见的。不确定量 量化（用确定量的方法研究）灰色系统视不确定量为灰色量。提出了灰色系统建模的具体数学方法，它能利用时间序列来确定微分方程的参数。灰色预测不是把观测到的数据序列视为一个随机过程，而是看作随时间变化的灰色量或灰色过程，通过累加生成和累减生成逐步使灰色量白化，从而建立相应于微分方程解的模型并做出预报。二、灰色预测模型1、GM（1,1）模型令为GM(1,1)建模序列，为的1-AGO（一次累加生成）序列，,令为的紧邻均值（MEAN）生成序列=0.5+0.5即定义：GM(1,1)的灰微分方程模型为 式中称为发展系数，为灰色作用量。设为待估参数向量，即，则灰微分方程的最小二乘估计参数列满足 =其中=， =定义： 为灰色微分方程的白化方程，也叫影子方程。如上所述，则有1） 白化方程的解也称时间响应函数为2） GM(1,1)灰色微分方程的时间响应序列为+，3） 取，则+，4） 还原值上式即为预测方程。5）模型检验：灰色预测模型检验有残差检验，关联度检验，后验差检验。（1）残差检验：计算原始序列和原始序列的灰色预测序列之间的：绝对误差； 相对误差：；其中。相对误差越小，模型精度越高。（2） 后验差检验：首先计算原始序列的均方差：，而，。然后计算残差序列的均方差：，而，。再计算方差比，最后计算小误差概率。 根据下面的预测精度等级划分表确定模型的精度。预测精度等级划分表小误差概率值方差比值预测精度等级好合格勉强合格不合格2、GM（1,1）模型应用实例例 某大型企业1999年至2004年的产品销售额如下表，试建立GM(1,1)预测模型，并预测2005年的产品销售额。年份199920002001200220032004销售额(亿元)2.673.133.253.363.563.72解：设=2.67，3.13，3.25，3.36，3.56，3.72第1步 构造累加生成序列=2.67,5.80,9.05,12.4，15.97,19.69第2步 构造数据矩阵和数据向量，第3步 计算=第4步 得出预测模型0.043879=2.925663=69.345766.6757（=2.67；=66.6757）第5步 残差检验(1)根据预测公式，计算，得2.67，5.78，9.03，12.43，15.97，19.68，19.69(=0,1, ,6)(2)累减生成序列，=1,2, ,62.67，3.11，3.25，3.40，3.54，3.71原始序列：2.67，3.13，3.25，3.36，3.56，3.72(3)计算绝对残差和相对残差序列绝对残差序列：0，0.02，0，0.04，0.02，0.01相对残差序列：0，0.64%，0，1.19%，0.56%，0.27%相对残差不超过1.19%，模型精确度高。预测：k=7，(8)=(8)(7)=4.233、GM（1,1）模型参数估计、检验及作图及预测的matlab程序。（1）程序解读% GM（1,1）模型计算及检验、作图。文件名：fungry1.m （注释部分不能在MATLAB窗口中运行）function GM1=fungry1(x0) % 输入原始数据x0T=input(T=); % 从键盘输入最后一个历史数据算起的第T时点% x1为累加生成序列；B为，初值为0；yn 为初值为0；% Hatx1：用灰色模型+，算出的序列；% Hatx0：灰色预测累减后的预测序列；% Epsilon：绝对误差；% omega：相对误差：；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);for i=1:length(x0) % 累加生成序列 for j=1:i x1(i)=x1(i)+x0(j); endendfor i=1:length(x0)-1 % 构造数据矩阵和数据向量 B(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=for k=1:length(x0)+T Hatx1(k)=(x0(1)-HatA(2)/HatA(1)*exp(-HatA(1)*(k-1)+HatA(2)/HatA(1);% +。endHatx0(1)=Hatx1(1); % 灰色模型算出序列和累减生成序列的第一个值相等为x0(1)for k=2:length(x0)+T Hatx0(k)=Hatx1(k)-Hatx1(k-1); %累减还原得到原始序列的模拟值；endfor i=1:length(x0) epsilon(i)=x0(i)-Hatx0(i); %计算相对误差 omega(i)=(epsilon(i)/x0(i); %计算绝对误差endx0,Hatx1,Hatx0,epsilon,omega, % 显示算出的一些序列c=std(epsilon)/std(x0); % 开始模型检验。Std：均方差。std(x0)：原始序列的均方差 p=0; %即P为绝对误差序列小的概率：，。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) 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（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);for i=1:length(x0) for j=1:i x1(i)=x1(i)+x0(j); endendfor i=1:length(x0)-1 B(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)+T Hatx1(k)=(x0(1)-HatA(2)/HatA(1)*exp(-HatA(1)*(k-1)+HatA(2)/HatA(1);endHatx0(1)=Hatx1(1);for k=2:length(x0)+T Hatx0(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.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) 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 x0=2.97 3.23 3.29 3.46 3.59 3.71; fungry1(x0);提示输入：T=1HatA = % 灰色预测模型系数 -0.0365 3.0412x0 = % 原始序列 2.9700 3.2300 3.2900 3.4600 3.5900 3.7100Hatx1 = % 灰色生成序列2.9700 6.1779 9.5051 12.9560 16.5353 20.2478 24.0984Hatx0 = % 灰色累减得到原始序列的模拟