Python灰色关联度

设系统行为序列
$$\begin{gathered} X_0=(x_0(1),x_0(2),\cdots ,x_0(n)) \\ {X}_1=(x_1(1),x_1(2),\cdots ,x_1(n)) \\ \cdots \cdots \\ {X}_i=(x_i(1),x_i(2),\cdots ,x_i(n)) \\ \cdots \cdots \\ {X}_m=(x_m(1),x_m(2),\cdots ,x_m(n)) \end{gathered}$$
对于$\xi \in (0,1)$, 令
$$\begin{gathered} \gamma (x_0(k),x_i(k))=\frac{{\min }_i{\max }_k|x_0(k)-x_i(k)|+\xi {\min }_i{\max }_k|x_0(k)-x_i(k)|}{|x_0(k)-x_i(k)|+\xi {\min }_i{\max }_k|x_0(k)-x_i(k)|} \\ \gamma (X_0,X_i)=\frac{1}{n}\sum _{k=1}^n\gamma (x_0(k),x_i(k)) \end{gathered}$$
则$\gamma (X_0,X_i)$满足灰色关联公理, 其中$\xi$称为分辨系数. $\gamma (X_0,X_i)$称为$X_0$与$X_i$的灰色关联度(邓聚龙, 1985)
灰色关联度 $\gamma \left(X_0,X_i\right)$ 常简记为 ${\gamma }_{0i},k$ 点关联系数 $\gamma \left(x_0(k),x_i(k)\right)$ 简记为 ${\gamma }_{0i}(k).$
灰色关联度的计算步骤如下.
第一步: 求各序列的初值象 (或均值象). 令
$$X_i^{\prime }=\frac{X_i}{x_i(1)}=\left(x_i^{\prime }(1),x_i^{\prime }(2),\cdots ,x_i^{\prime }(n)\right),i=0,1,2,\cdots ,m$$
第二步: 求 $X_0$ 与 $X_i$ 的初值象 $($ 或均值象) 对应分量之差的绝对值序列. 记
$${\Delta }_i(k)=\left|x_0^{\prime }(k)-x_i^{\prime }(k)\right|,{\Delta }_i=\left({\Delta }_i(1),{\Delta }_i(2),\cdots ,{\Delta }_i(n)\right),i=1,2,\cdots ,m$$
第三步: 求 ${\Delta }_i(k)=\left|x_0^{\prime }(k)-x_i^{\prime }(k)\right|,k=1,2,\cdots ,n;i=1,2,\cdots ,m$ 的最大值与最小值.
记为
$$M=\underset{i}{\max }\underset{k}{\max }{\Delta }_i(k),m=\underset{i}{\min }\underset{k}{\min }{\Delta }_i(k)$$
第四步: 计算关联系数
$${\gamma }_{0i}(k)=\frac{m+\xi M}{{\Delta }_i(k)+\xi M},\xi \in (0,1),k=1,2,\cdots ,n;i=1,2,\cdots ,m$$
第五步: 求出关联系数的平均值即所要的关联度
$${\gamma }_{0i}=\frac{1}{n}\sum _{k=1}^n{\gamma }_{0i}(k),i=1,2,\cdots ,m$$

• 举个例子
  $2001-2005$ 年, 我国 $\mathrm{G}\mathrm{D}\mathrm{P}$ 以及第一产业、第二产业、第三产业增加值数(单
  位: 千亿元) 分别如下:
  GDP 为
  $$X_1=\left(x_1(1),x_1(2),x_1(3),x_1(4),x_1(5)\right)=(109.7,120.3,135.8,159.9,183.1)$$
  第一产业增加值为
  $$X_2=\left(x_2(1),x_2(2),x_2(3),x_2(4),x_2(5)\right)=(15.5,16.2,17.1,21.0,23.1)$$
  第二产业增加值为
  $$X_3=\left(x_3(1),x_3(2),x_3(3),x_3(4),x_3(5)\right)=(49.5,53.9,62.4,73.9,87.0)$$
  第三产业增加值为
  $$X_4=\left(x_4(1),x_4(2),x_4(3),x_4(4),x_4(5)\right)=(44.6,50.2,56.3,65.0,73.0)$$
  资料来源: 2006 中国统计年鉴.
  以 $\mathrm{G}\mathrm{D}\mathrm{P}$ 为系统行为特征序列, 计算 $X_2,X_3,X_4$ 与 $X_1$ 的灰色关联度. 解 以 $X_1$ 为系统行为特征序列计算关联度.
  第一步: 求初值象. 由 $X_i^{\prime }=\frac{X_i}{x_i(1)}=\left(x_i^{\prime }(1),x_i^{\prime }(2),x_i^{\prime }(3),x_i^{\prime }(4),x_i^{\prime }(5)\right),i=1,2,3,4,$ 得
  $$\begin{array}{rl}{X}_1^{\prime }& =(1,1.0966,1.2379,1.4576,1.6691)\\ X_2^{\prime }& =(1,1.0452,1.1032,1.3548,1.4903)\\ X_3^{\prime }& =(1,1.0889,1.2606,1.4929,1.7576)\\ X_4^{\prime }& =(1,1.1256,1.2623,1.4574,1.6368)\end{array}$$
  第二步: $x_2,X_3,X_4$ 与 $X_1$ 初值象对应分量之差的绝对值序列. 由 ${\Delta }_i(k)=$ $\left|x_1^{\prime }(k)-x_i^{\prime }(k)\right|,i=2,3,4,$ 得
  $$\begin{array}{l}{\Delta }_2=(0,0.0515,0.1347,0.1028,0.1788)\\ {\Delta }_3=(0,0.0077,0.0227,0.0353,0.0885)\\ {\Delta }_4=(0,0.0289,0.0244,0.0002,0.0323)\end{array}$$
  第三步: 求 ${\Delta }_i(k)=\left|x_1^{\prime }(k)-x_i^{\prime }(k)\right|,i=2,3,4;k=1,2,3,4,5$ 的最大值与最小值.
  $$\begin{array}{c}M={\max }_i{\max }_k{\Delta }_i(k)=0.1788\\ m={\min }_i{\min }_k{\Delta }_i(k)=0\end{array}$$
  第四步: 求关联系数. 取 $\xi =0.5,$ 有
  $${\gamma }_{1i}(k)=\frac{m+\xi M}{{\Delta }_i(k)+\xi M}=\frac{0.0894}{{\Delta }_i(k)+0.0894},i=2,3,4;k=1,2,3,4,5$$
  从而
  $$\begin{array}{llllc}{r}_{12}(1)=1,& r_{12}(2)=0.6346,& r_{12}(3)=0.3989,& r_{12}(4)=0.4652,& r_{12}(5)=0.3333\\ r_{13}(1)=1,& r_{13}(2)=0.9203,& r_{13}(3)=0.7976,& r_{13}(4)=0.7168,& r_{13}(5)=0.5026\\ r_{14}(1)=1,& r_{14}(2)=0.7555,& r_{14}(3)=0.7855,& r_{14}(4)=0.9976,& r_{14}(5)=0.7344\end{array}$$
  第五步: 计算灰色关联度
  $$\begin{array}{l}{\gamma }_{12}=\frac{1}{5}{\sum }_{k=1}^5{\gamma }_{12}(k)=0.5664\\ {\gamma }_{13}=\frac{1}{5}{\sum }_{k=1}^5{\gamma }_{13}(k)=0.7875\\ {\gamma }_{14}=\frac{1}{5}{\sum }_{k=1}^5{\gamma }_{14}(k)=0.8546\end{array}$$
  因此, 第三产业与GDP关联度最大, 第二产业关联度次之, 第一产业关联度最小.
• 代码实现

import numpy as np

X1 = [109.7, 120.3, 135.8, 159.9, 183.1]
X2 = [15.5, 16.2, 17.1, 21.0, 23.1]
X3 = [49.5, 53.9, 62.4, 73.9, 87.0]
X4 = [44.6, 50.2, 56.3 ,65.0, 73.0]

def initialize(x):
    d = []
    for i in range(len(x)):
        d.append(round(x[i]/x[0],4))
    return d

def minus(X,x):
    d = []
    for i in range(len(X)):
        d.append(round(abs(X[i]-x[i]),4))
    return d
    

X11 =  initialize(X1)
X22 =  initialize(X2)
delta2 = minus(X11,X22)
X33 =  initialize(X3)
delta3 = minus(X11,X33)
X44 =  initialize(X4)
delta4 = minus(X11,X44)

M = max(delta2 + delta3 + delta4)
m = min(delta2 + delta3 + delta4)

ksi = 0.5
gamma2 = (m + ksi*M)/(np.array(delta2) + ksi*M)
gamma3 = (m + ksi*M)/(np.array(delta3) + ksi*M)
gamma4 = (m + ksi*M)/(np.array(delta4) + ksi*M)

np.mean(gamma2)
>>> 0.5664692086702057
np.mean(gamma3)
>>> 0.787530531907082
np.mean(gamma4)
>>> 0.8546034878054456

代码结果与手动计算结果一致