灰色预测程序

灰色预测模型灰色预测是就灰色系统所做的预测. 所谓灰色系统是介于白色系统和黑箱系统之间的过渡系统，其具体的含义是：如果某一系统的全部信息已知为白色系统，全部信息未知为黑箱系统，部分信息已知，部分信息未知，那么这一系统就是灰箱系统. 一般地说，社会系统、经济系统、生态系统都是灰色系统.灰色系统理论认为对既含有已知信息又含有未知或非确定信息的系统进行预测，就是对在一定方位内变化的、与时间有关的灰色过程的预测. 尽管过程中所显示的现象是随机的、杂乱无章的，但毕竟是有序的、有界的，因此这一数据集合具备潜在的规律，灰色预测就是利用这种规律建立灰色模型对灰色系统进行预测.灰色预测模型只需要较少的观测数据即可，这和时间序列分析，多元回归分析等需要较多数据的统计模型不一样. 因此，对于只有少量观测数据的项目来说，灰色预测是一种有用的工具.一、GM(1,1)模型灰色系统理论是邓聚龙教授在1981年提出来的，是一种对含有不确定因素系统进行预测的方法. 通过鉴别系统因素之间发展趋势的相异程度，进行关联分析，并通过对原始数据进行生成处理来寻找系统的变化规律，生成较强规律性数据序列，然后建立相应微分方程模型，从而预测事物未来的发展趋势和未来状态. 目前使用最广泛的灰色预测模型是关于数列预测的一个变量、一阶微分的GM(1,1)模型.GM(1,1)模型是基于灰色系统的理论思想，将离散变量连续化，用微分方程代替差分方程，按时间累加后所形成的新的时间序列呈现的规律可用一阶线性微分方程的解来逼近，用生成数序列代替原始时间序列，弱化原始时间序列的随机性，这样可以对变化过程作较长时间的描述，进而建立微分方程形式的模型. 其建模的实质是建立微分方程的系数，将时间序列转化为微分方程，通过灰色微分方程可以建立抽象系统的发展模型. 经证明，经一阶线性微分方程的解逼近所揭示的原始时间数列呈指数变化规律时，灰色预测GM(1,1)模型的预测将是非常成功的.1.1 GM(1,1)模型的建立灰色理论认为一切随机量都是在一定范围内、一定时间段上变化的灰色量及灰色过程. 数据处理不去寻找其统计规律和概率分布, 而是对原始数据作一定处理后, 使其成为有规律的时间序列数据, 在此基础上建立数学模型.GM(1,1)模型是指一阶，一个变量的微分方案预测模型，是一阶单序列的线性动态模型，用于时间序列预测的离散形式的微分方程模型.设时间序列()0X有n 个观察值，()()()()()()(){}00001,2,,Xx x x n =，为了使其成为有规律的时间序列数据，对其作一次累加生成运算，即令()()()()101tn xt x n ==∑从而得到新的生成数列()1X，()()()()()()(){}11111,2,,Xx x x n =，新的生成数列()1X 一般近似地服从指数规律. 则生成的离散形式的微分方程具体的形式为dxax u dt+= 即表示变量对于时间的一阶微分方程是连续的. 求解上述微分方程，解为当t =1时，()(1)x t x =，即(1)c x a=-，则可根据上述公式得到离散形式微分方程的具体形式为 ()()()11a t u u x t x e a a --⎛⎫=-+ ⎪⎝⎭其中，ax 项中的x 为dxdt的背景值，也称初始值；a ，u 是待识别的灰色参数，a 为发展系数，反映x 的发展趋势；u 为灰色作用量，反映数据间的变化关系.按白化导数定义有0()()lim t dx x t t x t dt t→+-= 显然，当时间密化值定义为1时，当1t →时，则上式可记为1lim(()())t dxx t t x t dt→=+- 这表明dxdt是一次累减生成的，因此该式可以改写为 (1)(1)(1)()dxx t x t dt=+- 当t 足够小时，变量x 从()x t 到()x t t +是不会出现突变的，所以取()x t 与()x t t +的平均值作为当t 足够小时的背景值，即(1)(1)(1)1()(1)2xx t x t ⎡⎤=++⎣⎦将其值带入式子，整理得 (0)(1)(1)1(1)()(1)2x t a x t x t u ⎡⎤+=-+++⎣⎦ 由其离散形式可得到如下矩阵：(1)(1)(0)(1)(1)(0)(0)(1)(1)1(1)(2)2(2)1(2)(3)(3)2()1(1)()2x x x x x x a u x n x n x n ⎛⎫⎡⎤-+ ⎪⎣⎦⎛⎫ ⎪ ⎪ ⎪⎡⎤-+ ⎪⎣⎦ ⎪=+ ⎪ ⎪ ⎪⎪ ⎪⎝⎭ ⎪⎡⎤--+ ⎪⎣⎦⎝⎭令 (0)(0)(0)(2),(3),,()TY x x x n ⎡⎤=⎣⎦(1)(1)(1)(1)(1)(1)11(1)(2)211(2)(3)21(1)()12x x x x B x n x n ⎛⎫⎡⎤-+ ⎪⎣⎦ ⎪⎪⎡⎤-+⎣⎦ ⎪= ⎪ ⎪ ⎪⎡⎤--+ ⎪⎣⎦⎝⎭()Ta u α=称Y 为数据向量，B 为数据矩阵，α为参数向量. 则上式可简化为线性模型：Y B α=由最小二乘估计方法得()1T T a B B B Y uα-⎛⎫== ⎪⎝⎭上式即为GM(1,1)参数,a u 的矩阵辨识算式，式中()1TT B B B Y -事实上是数据矩阵B 的广义逆矩阵.将求得的a ，u 值代入微分方程的解式，则()1(1)()((1))a t u ux t x e a a--=-+其中，上式是GM(1,1)模型的时间响应函数形式，将它离散化得(1)(0)(1)ˆ()(1)a t u u xt x e a a --⎛⎫=-+ ⎪⎝⎭ 对序列()()1ˆxt 再作累减生成可进行预测. 即()(0)(1)(1)(0)(1)ˆˆˆ()()(1)(1)1a a t xt x t x t u x e ea --=--⎛⎫=-- ⎪⎝⎭ 上式便是GM(1,1)模型的预测的具体计算式. 或对()atux t cea-=+求导还原得 (0)(0)(1)ˆ()((1))a t uxt a x e a--=-- 1.2 GM(1,1)模型的检验GM(1,1)模型的检验包括残差检验、关联度检验、后验差检验三种形式.每种检验对应不同功能：残差检验属于算术检验，对模型值和实际值的误差进行逐点检验；关联度检验属于几何检验范围，通过考察模型曲线与建模序列曲线的几何相似程度进行检验，关联度越大模型越好；后验差检验属于统计检验，对残差分布的统计特性进行检验，衡量灰色模型的精度. ➢ 残差检验残差大小检验，即对模型值和实际值的残差进行逐点检验. 设模拟值的残差序列为(0)()e t ，则(0)(0)(0)ˆ()()()e t x t xt =- 令()t ε为残差相对值，即残差百分比为(0)(0)(0)ˆ()()()%()x t xt t x t ε⎡⎤-=⎢⎥⎣⎦令∆为平均残差，11()nt t n ε=∆=∑.设残差的方差为22S ，则[]22211()n t S e t e n ==-∑. 故后验差比例C 为21/C S S =，误差频率P 为{}1()0.6745P P e t e S =-<.对于,C P 检验指标如下表：检验指标好合格勉强不合格P >0.95 >0.80 >0.70 <0.70 C <0.35 <0.50 <0.65 >0.65表 1 灰色预测精确度检验等级标准一般要求()20%t ε<，最好是()10%t ε<，符合要求.➢ 关联度检验关联度是用来定量描述各变化过程之间的差别. 关联系数越大，说明预测值和实际值越接近.设 {}(0)(0)(0)(0)ˆˆˆˆ()(1),(2),,()Xt xx x n =⋯ {}(0)(0)(0)(0)()(1),(2),,()X t x x x n =⋯序列关联系数定义为(){}{}{}(0)(0)(0)(0)(0)(0)(0)(0)ˆˆmin ()()max ()(),0ˆˆ()()max ()()1,0x t x t x t x t t t x t x t x t x t t σξσ⎧-+-⎪≠⎪=⎨-+-⎪=⎪⎩ 式中,(0)(0)ˆ()()xt x t -为第t 个点(0)x 和(0)ˆx 的绝对误差，()t ξ为第t 个数据的关联系数，ρ称为分辨率，即取定的最大差百分比，0ρ<<1,一般取0.5ρ=.(0)()x t 和(0)ˆ()xt 的关联度为()11nt r t n ξ==∑精度等级 关联度均方差比值小误差概率好（1级） 0.90≥ 0.35≤ 0.95≥ 合格（2级） 0.80≥ 0.50≤ 0.80≥ 勉强（3级） 0.70≥ 0.65≤ 0.70≥ 不合格（4级）0.70< 0.65>0.70<表 2 精度检验等级关联度大于60%便满意了，原始数据与预测数据关联度越大，模型越好.➢ 后验差检验后验差检验，即对残差分布的统计特性进行检验. 检验步骤如下：1、计算原始时间数列(){}0(0)(0)(0)(1),(2),,()Xx x x n =的均值和方差()2(0)(0)2(0)11111(),()n n t t xx t S x t x n n ====-∑∑ 2、计算残差数列{}(0)(0)(0)(0)(1),(2),,()ee e e n =的均值e 和方差22s()2(0)2(0)21111(),()n n t t e e t S e t e n n ====-∑∑其中(0)(0)(0)ˆ()()(),1,2,,e t x t xt t n =-=为残差数列.3、计算后验差比值21C S S =4、计算小误差频率{}(0)1()0.6745P P e t e S =-<令0S =0.67451S ，(0)()|()|t e t e ∆=-，即{}0()P P t S =∆<.若对给定的00C >，当0C C <时，称模型为方差比合格模型；若对给定的00P >，当0P P >时，称模型为小残差概率合格模型.>0.95 <0.35 优 >0.80 <0.5 合格 >0.70 <0.65 勉强合格 <0.70>0.65不合格表 3 后验差检验判别参照表1.3 残差GM(1,1)模型当原始数据序列(0)X建立的GM(1,1)模型检验不合格时，可以用GM(1,1)残差模型来修正. 如果原始序列建立的GM(1,1)模型不够精确，也可以用GM(1,1)残差模型来提高精度.若用原始序列(0)X建立的GM(1,1)模型(1)(0)ˆ(1)[(1)]at u uxt x e a a-+=-+ 可获得生成序列(1)X 的预测值，定义残差序列(0)(1)(1)ˆ()()()e k x k x k =-. 若取k=t , t+1, …, n ，则对应的残差序列为{}(0)(0)(0)(0)()(1),(2),,()e k e e e n =计算其生成序列(1)()e k ，并据此建立相应的GM(1,1)模型(1)(0)ˆ(1)[(1)]e a k e ee eu u et e e a a -+=-+ 得修正模型(1)(0)(0)(1)(1)()()(1)e a k ak e e e u u u x t x e k t a e e a a a δ--⎡⎤⎡⎤+=-++---⎢⎥⎢⎥⎣⎦⎣⎦其中1()0k tk t k t δ≥⎧-=⎨≤⎩为修正参数.应用此模型时要考虑：1、一般不是使用全部残差数据来建立模型，而只是利用了部分残差.2、修正模型所代表的是差分微分方程，其修正作用与()k t δ-中的t 的取值有关.1.4 GM(1,1)模型的适用范围定理：当GM(1,1)发展系数||2a ≥时，GM(1,1)模型没有意义.我们通过原始序列()0i X 与模拟序列()0ˆiX 进行误差分析，随着发展系数的增大，模拟误差迅速增加. 当发展系数0.3a -≤时，模拟精度可以达到98%以上；发展系数0.5a -≤时，模拟精度可以达到95%以上；发展系数1a ->时，模拟精度低于70%；发展系数 1.5a ->时，模拟精度低于50%. 进一步对预测误差进行考虑，当发展系数0.3a -<时，1步预测精度在98%以上，2步和5步预测精度都在90%以上，10步预测精度亦高于80%；当发展系数0.8a ->时，1步预测精度已低于70%.通过以上分析，可得下述结论：1、当0.3a -<时，GM(1,1)可用于中长期预测；2、当0.30.5a <-≤时，GM(1,1)可用于短期预测，中长期预测慎用；3、当0.50.8a <-≤时，GM(1,1)作短期预测应十分谨慎；4、当0.81a <-≤时，应采用残差修正GM(1,1)模型；5、当1a ->时，不宜采用GM(1,1)模型.1.5 GM(1,1)模型实例分析例：则该学生成绩时间序列如下：()()(0)(0)(0)(0)(0)(1),(2),(3),(4)79,74.825,74.29,76.98X x x x x ==对(0)X作一次累加后的数列为()()(1)(1)(1)(1)(1)(1),(2),(3),(4)79,153.825,228.115,305.095X x x x x ==对(1)X做紧邻均值生成. 令(1)(1)(1)()0.5()0.5(1)Z k x k x k =+-，得()()(1)(1)(1)(1)(2),(3),(4)116.4125,151.47,150.1925Z z z z ==则数据矩阵B 及数据向量Y 为(1)(1)(1)(2)1116.41251(3)1151.471(4)1150.19251z B z z ⎡⎤--⎡⎤⎢⎥⎢⎥=-=-⎢⎥⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦，(0)(0)(0)(2)74.825(3)74.29(4)76.98x Y x x ⎡⎤⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦ 对参数列ˆ[,]Taa b =进行最小二乘估计，得 176.61ˆ()[,]0.0144T T T T a B B B Y B Y a u -⎡⎤====⎢⎥-⎣⎦即 0.0144a =-，76.61u = 则GM(1,1)模型为()()110.014476.61dx x dt-= 时间响应式为(1)0.0144ˆ(1)5399.13895320.1389xk e -+=- 当1k =时，我们取(1)(0)(0)ˆˆ(1)(1)(0)79xx x === 还原求出(0)X的模拟值. 由(0)(1)(1)ˆˆˆ()()(1)Xk x k x k =--，取2,3,4k =，得 ()()(0)(0)(0)(0)(0)ˆˆˆˆˆ(1),(2),(3),(4)79,74.281,74.3584,76.4513xx x x x == 通过预测，得到实际值与预测值如下表：实际值 预测值 相对误差()k ε 第一学期79 79 0 第二学期 74.825 74.2810 0.73% 第三学期 74.29 74.3584 0.0921% 第四学期76.9876.45130.7051%表 4 四学期的实际值与预测值的误差表因为()10%k ε<，那就可得学生的预测值，与现实值进行比较得出该模型精度较高，可进行预测和预报.我们对学生未来两个学期（也就是第五、六个学期）的成绩进行预测，分别为77.5602分和78.6851分.例：某大型企业1999年至2004年的产品销售额如下表，试建立GM(1,1)预测模型，并预测2005年的产品销售额。

灰色预测模型步骤灰色预测模型是一种基于灰色理论的预测方法，其核心是建立一个数学模型来预测未来的发展趋势。

在实践中，灰色预测模型通常应用于经济、社会和环境等各个领域，以帮助决策者制定合理的规划和决策。

灰色预测模型的步骤主要包括以下5个方面：1、建立模型的数据预处理数据预处理是计算机处理向灰色预测模型输入数据的第一步。

在预处理过程中，需要对原始数据进行标准化处理，将非数值型数据转换为数值型数据，同时还需要对数据的质量进行评估，识别和剔除异常值。

2、建立灰色驱动模型在数据预处理后，需要建立一个灰色驱动模型。

该模型是一种简化的数学模型，用于描述因变量和自变量之间的离散关系。

此外，该模型还需要根据实际的情况调整参数，以提高模型的准确性。

3、对模型进行验证在灰色预测模型中，模型验证是非常重要的一步。

通过对模型进行验证，可以评估模型的预测精度，并确定预测误差的可接受范围。

如果模型的预测误差过大，则需要进一步调整模型，以获得更准确的预测结果。

4、进行预测在完成模型的验证后，需要对所建立的模型进行预测。

预测的结果通常是未来的某个时间点的数值预测。

预测结果需要基于实际情况进行解读和分析，并形成有效的决策参考。

5、模型的评价和修正最后，需要对模型进行评价和修正。

因为灰色预测模型是一种逐步调整的预测方法，因此需要在应用过程中进行持续的评价和修正。

通过评价，可以确定模型的适用性和准确性，从而更好地应对未来的预测任务。

总之，灰色预测模型是一种有效的预测方法，可以在很大程度上提高预测精度和决策效率。

通过逐步的建模和修正，该方法可以为各个领域提供更好的预测和决策参考。

GM(1,1)灰色理论认为系统的行为现象尽管是朦胧的，数据是复杂的，但它毕竟是有序的，是有整体功能的。

灰数的生成，就是从杂乱中寻找出规律。

同时，灰色理论建立的是生成数据模型，不是原始数据模型，因此，灰色预测是一种对含有不确定因素的系统进行预测的方法。

灰色预测是一种对含有不确定因素的系统进行预测的方法。

灰色预测通过鉴别系统因素之间发展趋势的相异程度，即进行关联分析，并对原始数据进行生成处理来寻找系统变动的规律，生成有较强规律性的数据序列，然后建立相应的微分方程模型，从而预测事物未来发展趋势的状况。

其用等时距观测到的反应预测对象特征的一系列数量值构造灰色预测模型，预测未来某一时刻的特征量，或达到某一特征量的时间。

1．GM(1,1)模型预测方法已知参考数据列()(0)(0)(0)(0)(1),(2),,()x x x x n =⋅⋅⋅,1次累加生成序列(1AGO)- ()()(1)(1)(1)(1)(0)(0)(0)(0)(0)(1),(2),,()(1),(1)(2),,(1)()x x x x n x x x x x n =⋅⋅⋅=+⋅⋅⋅+⋅⋅⋅+其中：(1)(0)1()(),1,2,,ki x k x i k n ===⋅⋅⋅∑。

(1)x 的均值生成序列 ()(1)(1)(1)(1)(2),(3),,()z z z z n =⋅⋅⋅其中：(1)(1)(1)()0.5()0.5(1),2,3,,z k x k x k k n =+-=⋅⋅⋅。

建立灰微分方程(0)(1)()(),2,3,,,x k az k b k n +==⋅⋅⋅相应的白化微分方程为(1)(1)()dx ax t b dt+= 记T [,]u a b =，T (0)(0)(0)(2),(3),,()Y x x x n ⎡⎤=⋅⋅⋅⎣⎦，(1)(1)(1)(2)1(3)1()1z z B z n ⎡⎤-⎢⎥-⎢⎥=⎢⎥⎢⎥-⎣⎦,则由最小二乘法，求得使T ()()()J u Y Bu Y Bu =--达到最小值的u 的估计值为()T1T T ˆˆˆ,u a b B B B Y -⎡⎤==⎣⎦于是求解其白化微分方程得ˆ(1)(0)ˆˆ(1)(1),0,1,,1,ˆˆak b b x k x e k n a a -⎛⎫+=-+=⋅⋅⋅-⋅⋅⋅ ⎪ ⎪⎝⎭2. GM(1,1)模型预测步骤（1）数据的检验与处理首先，为了保证建模方法的可行性，需要对已知数据列作必要的检验处理。

GM(1,1)模型预测方法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\Yx=dsolve('Dx+a*x=b','x(0)=x0');x=subs(x,{'a','b','x0'},{u(1),u(2),x1(1)});yuce1=subs(x,'t',[0:n-1]); %为提高预测精度，先计算预测值，再显示微分方程的解y=vpa(x,6) %其中的6 表示显示6 位数字yuce=[x0(1),diff(yuce1)] %差分运算，还原数据（预测值）epsilon=x0'-yuce %计算残差delta=abs(epsilon./x0') %计算相对误差rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda' %计算级比偏差值GM(2,1)模型预测方法clc,clearx0=[41,49,61,78,96,104];n=length(x0);x1=cumsum(x0)a_x0=diff(x0)'z=0.5*(x1(2:end)+x1(1:end-1))';B=[-x0(2:end)',-z,ones(n-1,1)];Y=a_x0;u=B\Y %最小二乘法拟合参数%下面求微分方程边值问题的解x=dsolve('D2x+a1*Dx+a2*x=b','x(0)=c1,x(5)=c2');x=subs(x,{'a1','a2','b','c1','c2'},{u(1),u(2),u(3),x1(1),x1(6)});yuce=subs(x,'t',0:n-1);digits(6),x=vpa(x)x0_hat=[yuce(1),diff(yuce)];x0_hat=round(x0_hat) %四舍五入取整数epsilon=x0-x0_hat %求残差delta=abs(epsilon./x0) %求相对误差灰色Verhulst预测模型nian=1990:2003;%原始数据x1=[1 2 3 4 5 6 7 8 9 10 11 12 13 14];%原始数据n=length(x1);plot(nian,x1,'o-');x0=diff(x1); %作累减生成x0=[x1(1),x0]z1=0.5*(x1(2:n)+x1(1:n-1)) %求紧邻均值生成序列B=[-z1',z1'.^2]Y=x0(2:end)'ab_hat=B\Y %估计参数a,b 的值x=dsolve('Dx+a*x=b*x^2','x(0)=x0'); %求解常微分方程x=subs(x,{'a','b','x0'},{ab_hat(1),ab_hat(2),x1(1)}); %代入参数值yuce=subs(x,'t',0:14) %计算预测值%下面显示微分方程的解，为了提高计算精度，把该语句放在计算预测值之后x=vpa(x,6)x1_all=[x1,9.92,10.71]; %加上2004 年的两个观测值yuce(16)=yuce(15); %2004 年有两个观测值，要对应两个相同的预测值epsilon=x1_all-yuce %计算残差delta=abs(epsilon./x1_all) %计算相对误差delta_mean=mean(delta) %计算平均相对误差x1_all_0=x1_all-x1_all(1); %观测值数据列的始点零化像yuce_0=yuce-yuce(1); %预测值数据列的始点零化像s0=abs(sum(x1_all_0(2:end-1))+0.5*x1_all_0(end));s1=abs(sum(yuce_0(2:end-1))+0.5*yuce_0(end));tt=yuce_0-x1_all_0;s1_s0=abs(sum(tt(2:end-1))+0.5*tt(end));absdegree=(1+s0+s1)/(1+s0+s1+s1_s0) %计算灰色绝对关联度c=std(epsilon,1)/std(x1_all,1) %计算标准差比值。