关联度分析、灰色预测GM(1,1)、GM(1,1)残差模型——基于Python实现

关联度分析

import numpy as np
import pandas as pd
#关联度分析
#参考序列
Y_0=[170,174,197,216.4,235.8]
#被比较序列
Y_1=[195.4,189.9,187.2,205,222.7]
Y_2=[308,310,295,346,367]

#初始化序列
X_0=np.array(Y_0)/Y_0[0]
X_1=np.array(Y_1)/Y_1[0]
X_2=np.array(Y_2)/Y_2[0]

#计算绝对差序列
deta1=np.abs(X_0-X_1)
deta2=np.abs(X_0-X_2)

#计算deta1，deta2最小值
min1=np.min([np.min(deta1),np.min(deta2)])
max1=np.max([np.max(deta1),np.max(deta2)])

#计算关联系数yita
rio=0.5
yita1 = [(min1 + rio * max1) / (deta1[i] + rio * max1) for i in range(len(deta1))]#计算关联系数yita1
yita2 = [(min1 + rio * max1) / (deta2[i] + rio * max1) for i in range(len(deta2))]#计算关联系数yita2
# #计算关联度
r1=np.mean(yita1)
r2=np.mean(yita2)
# 创建DataFrame
df = pd.DataFrame({
    'yita1': yita1,
    'yita2': yita2
})
#更改索引
df.index = ['1','2','3','4','5']

print('关联度分析：')
print(df)
#输出关联度
print('关联度：',r1,r2)

关联度分析：
      yita1     yita2
1  1.000000  1.000000
2  0.705293  0.878928
3  0.381163  0.380878
4  0.355909  0.452620
5  0.333333  0.387478
关联度： 0.5551396262321364 0.6199809489771715

灰色预测GM(1,1)模型

from math import exp
#原始序列
X_0=(6,20,40,25,45,35,21,14,18,15.5,17,15)
#累加生成序列
X_1 = np.cumsum(X_0)
#print('累加生成序列:', X_1)
#构造矩阵B和数据向量Y
B = np.zeros((len(X_1)-1,2))
Y = X_0[1:]
for i in range(len(X_1)-1):
    B[i][0] = -0.5*(X_1[i]+X_1[i+1])
    B[i][1] = 1

#print('矩阵B:', B)
#print('数据向量Y:', Y)
#计算参数a,b
A = np.dot(np.dot(np.linalg.inv(np.dot(B.T,B)),B.T),Y)
a = A[0]
miu = A[1]
#计算预测模型
X_1_predict = np.zeros(len(X_0))
X_1_predict[0] = X_0[0]
for i in range(1, len(X_0)):
    X_1_predict[i] = ((X_0[0] - miu / a) * exp(-a * i) + miu / a)
#预测方程式
print('预测方程式：', 'X_1_(k+1) = ', X_0[0]- miu/ a, '*exp(-', a, '*i) +', miu/ a)

预测方程式： X_1_(k+1) =  -476.05668934240276 *exp(- 0.0746651965600655 *i) + 482.05668934240276

模型检验

#残差检验
X_1_predict = np.zeros(len(X_0))
X_1_predict[0] = X_0[0]
for i in range(1,len(X_0)):
    X_1_predict[i] = ((X_0[0]-miu/a)*exp(-a*i)+miu/a)
    
print('预测值:', X_1_predict)
#累减生成序列X_0_predict
X_0_predict = np.zeros(len(X_0))
X_0_predict[0] = X_0[0]
for i in range(1,len(X_0)):
    X_0_predict[i] = X_1_predict[i] - X_1_predict[i-1]
    
#累减生成序列
print('累减生成序列:', X_0_predict)

预测值: [  6.          40.25030314  72.03643912 101.53569452 128.91260092
 154.31985252 177.8991578  199.78202993 220.09052023 238.93789892
 256.42928694 272.66224218]
累减生成序列: [ 6.         34.25030314 31.78613597 29.4992554  27.3769064  25.4072516
 23.57930529 21.88287213 20.30849029 18.8473787  17.49138802 16.23295524]

# X_0_predict
#计算绝对误差
error = np.abs(X_0_predict-X_0)
print('绝对误差:', error)

绝对误差: [ 0.         14.25030314  8.21386403  4.4992554  17.6230936   9.5927484
  2.57930529  7.88287213  2.30849029  3.3473787   0.49138802  1.23295524]

#计算相对误差
error_rate = error/X_0
#输出相对误差,保留四位小数，输出百分比
print('相对误差:', np.round(error_rate,4)*100,'%')

相对误差: [ 0.   71.25 20.53 18.   39.16 27.41 12.28 56.31 12.82 21.6   2.89  8.22] %

#输出检验表，包括原始序列，预测序列，残差序列，绝对误差，相对误差
df = pd.DataFrame({
    '原始序列': X_0,
    '预测序列': X_1_predict,
    '残差序列': X_0_predict,
    '绝对误差': error,
    '相对误差': error_rate
})
df.index = ['1','2','3','4','5','6','7','8','9','10','11','12']
print('检验表：')
print(df)

检验表：
    原始序列        预测序列       残差序列       绝对误差      相对误差
1    6.0    6.000000   6.000000   0.000000  0.000000
2   20.0   40.250303  34.250303  14.250303  0.712515
3   40.0   72.036439  31.786136   8.213864  0.205347
4   25.0  101.535695  29.499255   4.499255  0.179970
5   45.0  128.912601  27.376906  17.623094  0.391624
6   35.0  154.319853  25.407252   9.592748  0.274079
7   21.0  177.899158  23.579305   2.579305  0.122824
8   14.0  199.782030  21.882872   7.882872  0.563062
9   18.0  220.090520  20.308490   2.308490  0.128249
10  15.5  238.937899  18.847379   3.347379  0.215960
11  17.0  256.429287  17.491388   0.491388  0.028905
12  15.0  272.662242  16.232955   1.232955  0.082197

相对误差大于于0.5%，模型检验不通过

GM(1,1)残差模型

#残差模型进行修正
#对erroe序列进行GM(1,1)模型预测
#原始序列
X_0 = error
#累加生成序列
X_1 = np.cumsum(X_0)
#构造矩阵B和数据向量Y
B = np.zeros((len(X_1)-1,2))
Y = X_0[1:]
for i in range(len(X_1)-1):
    B[i][0] = -0.5*(X_1[i]+X_1[i+1])
    B[i][1] = 1
#计算参数a,b
A = np.dot(np.dot(np.linalg.inv(np.dot(B.T,B)),B.T),Y)
a = A[0]
miu = A[1]
#计算预测模型
X_1_predict = np.zeros(len(X_0))
X_1_predict[0] = X_0[0]
for i in range(1, len(X_0)):
    X_1_predict[i] = ((X_0[0] - miu / a) * exp(-a * i) + miu / a)
# #预测方程式
print('预测方程式：', 'X_1_(k+1) = ', X_0[0]- miu/ a, '*exp(-', a, '*k) +', miu/ a)

预测方程式： X_1_(k+1) =  -47.4894117559714 *exp(- 0.21716818506240565 *k) + 47.4894117559714

#修正后的模型为原始模型加上求导后的方程式
# Differentiate the prediction equation with respect to k
X_1_predict_derivative = np.zeros(len(X_0))

for i in range(1, len(X_0)):
    deta3= 1 if (i + 1) >= 2  else 0
    X_1_predict_derivative[i] = -a * (X_0[0] - miu / a) * exp(-a * (i-1)) * deta3

# Corrected model by adding the differentiated equation to the original model
X_1_corrected = X_1_predict + X_1_predict_derivative

print('修正后的模型:', X_1_corrected)
#输出修正后的模型方程式
print('修正后的模型方程式：', 'X_1_(k+1) = ', X_0[0]- miu/ a, '*exp(-', a, '*k) +', miu/ a, ' - ', -a*(X_0[0] - miu/a), 'deta(k-1)*exp(-',a, '*k) ')
#输出修正后的模型方程式，分开为k>=2,k<2
print('修正后的模型方程式：', 'X_1_(k+1) = ', X_0[0]- miu/ a, '*exp(-', a, '*k) +', miu/ a, ' - ', -a*(X_0[0] - miu/a), 'deta(k-1)*exp(-',a, '*(k-1)) ', 'k>=2')
print('修正后的模型方程式：', 'X_1_(k+1) = ', X_0[0]- miu/ a, '*exp(-', a, '*k) +', miu/ a, 'k<2')

修正后的模型: [ 0.         19.58337881 25.03078704 29.41483178 32.94308733 35.78260836
 38.06783956 39.90698129 41.38711264 42.57831436 43.53698707 44.3085217 ]
修正后的模型方程式： X_1_(k+1) =  -47.4894117559714 *exp(- 0.21716818506240565 *k) + 47.4894117559714  -  10.31318936072558 deta(k-1)*exp(- 0.21716818506240565 *k) 
修正后的模型方程式： X_1_(k+1) =  -47.4894117559714 *exp(- 0.21716818506240565 *k) + 47.4894117559714  -  10.31318936072558 deta(k-1)*exp(- 0.21716818506240565 *(k-1))  k>=2
修正后的模型方程式： X_1_(k+1) =  -47.4894117559714 *exp(- 0.21716818506240565 *k) + 47.4894117559714 k<2

#残差检验
X_1_corrected
#累减生成序列
X_0_corrected = np.zeros(len(X_1_corrected))
X_0_corrected[0] = X_1_corrected[0]
for i in range(1, len(X_1_corrected)):
    X_0_corrected[i] = X_1_corrected[i] - X_1_corrected[i-1]

print('累减生成序列:', X_0_corrected)

累减生成序列: [ 0.         19.58337881  5.44740823  4.38404474  3.52825555  2.83952103
  2.2852312   1.83914174  1.48013134  1.19120172  0.95867271  0.77153463]

# X_0_predict
#计算绝对误差
error = np.abs(X_0_corrected-X_0)
print('绝对误差:', error)

绝对误差: [0.         6.11131907 4.55574461 1.15568474 7.53661955 1.18744444
 0.14100682 2.03499833 0.386519   0.73109566 1.00028958 0.07505627]

#计算相对误差
X_0=(6,20,40,25,45,35,21,14,18,15.5,17,15)
error_rate = error/X_0
#输出相对误差,保留四位小数，输出百分比
print('相对误差:', np.round(error_rate,4)*100,'%')

相对误差: [ 0.   30.56 11.39  4.62 16.75  3.39  0.67 14.54  2.15  4.72  5.88  0.5 ] %

#输出检验表，包括原始序列，预测序列，残差序列，绝对误差，相对误差
df = pd.DataFrame({
    '原始序列': X_0,
    '预测序列': X_1_corrected,
    '残差序列': X_0_corrected,
    '绝对误差': error,
    '相对误差': error_rate
})
df.index = ['1','2','3','4','5','6','7','8','9','10','11','12']
print('检验表：')
print(df)

检验表：
    原始序列       预测序列       残差序列      绝对误差      相对误差
1    6.0   0.000000   0.000000  0.000000  0.000000
2   20.0  19.583379  19.583379  6.111319  0.305566
3   40.0  25.030787   5.447408  4.555745  0.113894
4   25.0  29.414832   4.384045  1.155685  0.046227
5   45.0  32.943087   3.528256  7.536620  0.167480
6   35.0  35.782608   2.839521  1.187444  0.033927
7   21.0  38.067840   2.285231  0.141007  0.006715
8   14.0  39.906981   1.839142  2.034998  0.145357
9   18.0  41.387113   1.480131  0.386519  0.021473
10  15.5  42.578314   1.191202  0.731096  0.047167
11  17.0  43.536987   0.958673  1.000290  0.058841
12  15.0  44.308522   0.771535  0.075056  0.005004

修正后的模型检验通过，模型预测准确

预测

#预测13月份的数据
X_13_predict = ((X_0[0]-miu/a)*exp(-a*12)+miu/a)
print('13月份的预测值:', X_13_predict)

13月份的预测值: 76.46124387795606