灰色关联分析/灰色预测

Halcom

       灰色理论认为系统的行为现象尽管是朦胧的，数据是复杂的，但它毕竟是有序的，是有整体功能的。灰数的生成，就是从杂乱中寻找出规律。同时，灰色理论建立的是生成数据模型，不是原始数据模型，因此，灰色预测的数据是通过生成数据的GM(1，1)模型所得到的预测值的逆处理结果。
       所谓灰色系统是介于白色系统和黑箱系统之间的过渡系统，其具体的含义是：如果某一系统的全部信息已知为白色系统，全部信息未知为黑箱系统，部分信息已知，部分信息未知，那么这一系统就是灰色系统。一般地说，社会系统、经济系统、生态系统都是灰色系统。例如物价系统，导致物价上涨的因素很多，但已知的却不多，因此对物价这一灰色系统的预测可以用灰色预测方法。
       灰色系统理论认为对既含有已知信息又含有未知或非确定信息的系统进行预测，就是对在一定方位内变化的、与时间有关的灰色过程的预测。尽管过程中所显示的现象是随机的、杂乱无章的，但毕竟是有序的、有界的，因此这一数据集合具备潜在的规律，灰色预测就是利用这种规律建立灰色模型对灰色系统进行预测。
      灰色预测通过鉴别系统因素之间发展趋势的相异程度，即进行关联分析，并对原始数据进行生成处理来寻找系统变动的规律，生成有较强规律性的数据序列，然后建立相应的微分方程模型，从而预测事物未来发展趋势的状况。其用等时距观测到的反应预测对象特征的一系列数量值构造灰色预测模型，预测未来某一时刻的特征量，或达到某一特征量的时间。
      灰色预测模型求解步骤，
      Step 1： 设有原始数据列： （x1,x2,x3,……xn）；
      Step 2： 通过累加得到数列：
                                                   （y1,y2,y3,y4……yn）=(x1,x1+x2, x1+x2+x3, ……， x1+x2+x3+……xn )
      Step 3： 再求均值数列：
                                                       (z1,z2,z3,……，zn)=(0.5y1+0.5y2, 0.5y2+0.5y3, ……，0.5yn-1 + 0.5yn)
设有：
                                                          x2+a z2=b
                                                          x3+a z3=b     
                                                          x4+a z4=b
                                                          x5+a z5=b
                                                         …………………
                                                          xn+a zn=b     
Step 4：用最小二乘法寻找其中常数 a，b 。于是可得如下的白化微分方程：
                                                            dy/dt + a y=b
Step 5：解此微分方程可得该题拟合的函数模型为：
                                                           y(k+1)= (x1-b/a)*exp(-a*k) + b/a
其中， xk+1 = yk+1-yk，于是 xk = (x1-b/a)*( exp(-a*(k-1)) -exp( -a*(k-2) ) )。
灰色预测模型称为GM模型，G为grey的第一个字母，M为model的第一个字母。GM（1，1）表示一阶的，一个变量的微分方程型预测模型。GM（1，1）是一阶单序列的线性动态模型，主要用于时间序列预测。

灰色预测建模
上述已经对灰色预测的原理以及理论分析做了系统的阐述，灰色预测中GM(1，1)预测模型的求解过程如下。
Step 1：设有数列x(0) 共有 n 个观察值，
                                                             x(0)(1), x(0)(1), x(0)(2), ……x(0)(n)
Step 2：对 x(0) 作累加生成，得到新的数列 x(1)，其元素
                                                             x(1)(i) = sum( x(0)(m) ),  m = 1,2,……i
                                                                   i = 1,2,3, ……,n
有：
                                                                x(1)(1) = x(0)(1)
                                                                x(1)(2) = x(0)(1)+x(0)(2) = x(1)(1)+x(0)(2)
                                                                x(1)(3) = x(0)(1)+x(0)(2)++x(0)(3) = x(1)(2) + x(0)(3)
                                                                ……………………
                                                                x(1)(n) = x(1)(n-1)+x(0)(n)
Step:3：对数列x(1) ，可建立预测模型的白化形式方程，
                                                                dx(1)/dt + ax(1) = u
式中：a, u —为待估计参数。分别称为发展灰数和内生控制灰数。设 aa为待估计参数向量。
则
                                                                 aa = [a; u]
Step 4：按最小二乘法求解有：
                                                                 aa = inv(B'*B)*B'*yn
式中：
                                                                 B=[ -0.5*(x(1)(1)+x(1)(2)),  1;
                                                                        -0.5*(x(1)(2)+x(1)(3)),  1;
                                                                        -0.5*(x(1)(3)+x(1)(4)),  1;
                                                                        ……………………
                                                                        -0.5*(x(1)(n-1)+x(1)(n)),  1;]
                                                                yn = [x(0)(2),  x(0)(3),  x(0)(4),  ……,x(0)(n)]'
Step 5：转到Step 3，将 aa 代入其微分方程式，并解微分方程，有 GM(1, 1)预测模型为：
                                                             x(1)(i+1) = ( x(0)(1) - u/a )*exp( -a*i ) +u/a  

1. clc,clear,close all
   
2. y=[ 102.4        102.8        103.1        102.9        103.3        103.5        103.6        104.4        105.1        104.6        104.9        104.9        105.4
   
3. 105.2        105.9        106.1        105.7        106.8        107.5        108        110.1        117.1        109.6        110.3        111        111.7];%%
   
4. plot(1:13,y');figure(gcf);
   
5. grid on
   
6. axis([0,14,90,120])
   
7. legend('居民消费价格指数（以上年同期价格为100）','城市食品零售价格指数（上年同月=100）')
   
8. gtext('资料来源：中国国家统计局网站')
   

复制代码

Halcom

1. clc,clear,close all
   
2. y=[ 102.4   102.8   103.1   102.9   103.3   103.5   103.6   104.4   105.1   104.6   104.9   104.9   105.4
   
3. 105.2   105.9   106.1   105.7   106.8   107.5   108 110.1   117.1   109.6   110.3   111 111.7
   
4. 101.7   101.7   101.7   101.7   101.6   101.5   101.4   101.5   101.6   101.8   101.8   101.9   102.1
   
5. 98.9    98.7    98.8    99  99.2    98.8    98.5    98.7    99.3    100.1   99.8    100.4   100.8
   
6. 99.3    99.5    99.7    100 100.2   100.4   100.4   100.5   100.7   101.2   101.4   101.4   101.9
   
7. 102.5   102.8   103.2   103.2   103.3   103.3   103.4   103.7   104 104 103.2   103 103.2
   
8. 100 100 100.1   99.7    99.3    99.4    99.3    99.5    99.3    99.3    99.9    99.7    100.1
   
9. 100.3   100.4   100.6   100.9   101.1   101.2   101.2   100.9   100.6   100.7   101 100.3   100.5
   
10. 100.3   104.5   105 105 104.8   104.4   104.3   104.9   105.8   106 106.8   106.1   106.6];
    
11. y1=mean(y');
    
12. y1=y1';
    
13. for i=1:9
    
14.     for j=1:13
    
15.         y2(i,j)=y(i,j)/y1(i);
    
16.     end
    
17. end
    
18. for i=2:9
    
19.     for j=1:13
    
20.         y3(i-1,j)=abs(y2(i,j)-y2((i-1) ,j));
    
21.     end
    
22. end
    
23. a=1;b=0;
    
24. for i=1:8
    
25.        for j=1:13
    
26.             if(y3(i,j)<=a)  a=y3(i,j);
    
27.                 elseif(y3(i,j)>=b) b=y3(i,j);
    
28.             end
    
29.         end
    
30. end
    
31. for i=1:8
    
32.     for j=1:13
    
33.             y4(i,j)=(a+0.5*b)/(y3(i,j)+0.5*b);
    
34.     end
    
35. end
    
36. y5=sum(y4')/12
    
37. y5 =
    
38.   Columns 1 through 6
    
39.     0.8223    0.7225    0.9665    0.9951    0.9544    0.9469
    
40.   Columns 7 through 8
    
41.     0.9641    0.8760
    
42. >>
    

复制代码

具体参考书籍《MATLAB优化算法案例分析与应用》、《MATLAB数学建模案例应用与实战》，余胜威著

Halcom

灰色关联分析的计算步骤如下:
Step 1，建立各相关指标的原始数据矩阵
                                                            xi = ( xi(1), xi(2), xi(3)，……xi(k)，…… )
式中 xi(k) 表示 i 因素在第 k 年的原始数据。
Step 2，求初值像矩阵 xi'
                                                            xi' = (xi(1)/xi(1), xi(2)/xi(1),xi(3)/xi(1), ……， xi(k)/xi(1) ) = (xi(1)' , xi(2)' , xi(3)' ,xi(4)' ,……xi(k)',…… )
Step 3，求差序列
                                                            delta(i, k) = abs( xi(k)' - xi(k-1)' )
Step 4，计算关联系数 ksai 和灰色关联度  delta(k)
                                                           ksai = ( minmin(delta(i, k)) + fai* maxmax(delta(i, k)) )/ ( delta(i, k) +fai*maxmax(delta(i, k)) )
其中 fai 为分辨系数，其作用在于提高关联系数间的差异显著性，0<fai<1  ，一般取 fai = 0.5。
则灰色关联度gamma ,
                                                           gamma = sum( ksai  ) / (n-1)