预备知识
(1)灰色系统
白色系统是指系统内部特征是完全已知的；黑色系统是指系统内部信息完全未知的；而灰色系统是介于白色系统和黑色系统之间的一种系统，灰色系统其内部一部分信息已知，另一部分信息未知或不确定。
(2)灰色预测
灰色预测，是指对系统行为特征值的发展变化进行的预测，对既含有已知信息又含有不确定信息的系统进行的预测，也就是对在一定范围内变化的、与时间序列有关的灰过程进行预测。尽管灰过程中所显示的现象是随机的、杂乱无章的，但毕竟是有序的、有界的，因此得到的数据集合具备潜在的规律。灰色预测是利用这种规律建立灰色模型对灰色系统进行预测。
目前使用最广泛的灰色预测模型就是关于数列预测的一个变量、一阶微分的GM(1,1)模型。它是基于随机的原始时间序列，经按时间累加后所形成的新的时间序列呈现的规律可用一阶线性微分方程的解来逼近。经证明，经一阶线性微分方程的解逼近所揭示的原始时间序列呈指数变化规律。因此，当原始时间序列隐含着指数变化规律时，灰色模型GM(1,1)的预测是非常成功的。
Python代码如下:
# -*- coding: utf-8 -*-"""Created on Mon Aug 29 19:23:24 @author: DaiPW"""import numpy as npfrom pandas import Seriesfrom pandas import DataFrameimport pandas as pdimport matplotlib.pyplot as pltdef Identification_Algorithm(x): #辨识算法B = np.array([[1]*2]*5)tmp = np.cumsum(x)for i in range(len(x)-1):B[i][0] = ( tmp[i] + tmp[i+1] ) * (-1.0) / 2Y = np.transpose(x[1:])BT = np.transpose(B)a = np.linalg.inv(np.dot(BT,B))a = np.dot(a,BT)a = np.dot(a,Y)a = np.transpose(a)return a;def GM_Model(X0,a,tmp):#GM(1,1)模型A = np.ones(len(X0))for i in range(len(A)):A[i] = a[1]/a[0] + (X0[0]-a[1]/a[0])*np.exp(a[0]*(tmp[i]-1)*(-1))print ('GM(1,1)模型为:\nX(k) = ',X0[0]-a[1]/a[0],'exp(',-a[0],'(k-1))',a[1]/a[0])XK = Series(A,index=pd.period_range('2000','',freq = 'A-DEC')) print ('GM(1,1)模型计算值为:')print (XK)return XK;def Return(XK): #预测值还原tmp = np.ones(len(XK)) for i in range(len(XK)):if i == 0:tmp[i] = XK[i]else:tmp[i] = XK[i] - XK[i-1]X_Return = Series(tmp,index=pd.period_range('2000','',freq = 'A-DEC'))print ('还原值为:\n')print (X_Return)return X_Returnif __name__ == '__main__':#初始化原始数据date = pd.period_range('2000','',freq = 'A-DEC')tmp = np.array([1,2,3,4,5,6])data = np.array([132,92,118,130,187,207])X0 = Series(data,index = date)X0_copy = Series(data,index=tmp)print ('原始数据为:\n')print(X0)#对原始数据惊醒一次累加X1 = np.cumsum(X0)print ('原始数据累加为:')print(X1)#辨识算法a = Identification_Algorithm(data)print ('a矩阵为:')print (a)#GM(1,1)模型XK = GM_Model(X0,a,tmp)#预测值还原X_Return = Return(XK)#预测值即预测值精度表X_Compare1 = np.ones(len(X0))X_Compare2 = np.ones(len(X0))for i in range(len(data)):X_Compare1[i] = data[i]-X_Return[i]X_Compare2[i] = X_Compare1[i]/data[i]*100Compare = {'GM':XK,'1—AGO':np.cumsum(data),'Returnvalue':X_Return,'Realityvalue':data,'Error':X_Compare1,'RelativeError(%)':X_Compare2}X_Compare = DataFrame(Compare,index=date)print ('预测值即预测值精度表')print (X_Compare)#模型检验error_square = np.dot(X_Compare,np.transpose(X_Compare)) #残差平方和error_avg = np.mean(error_square) #平均相对误差S = 0#X0的关联度for i in range(1,len(X0)-1,1):S += X0[i]-X0[0]+(XK[-1]-XK[0])/2S = np.abs(S)SK = 0 #XK的关联度for i in range(1,len(XK)-1,1):SK += XK[i]-XK[0]+(XK[-1]-XK[0])/2SK = np.abs(SK)S_Sub = 0 #|S-SK|bfor i in range(1,len(XK)-1,1):S_Sub += X0[i]-X0[0]-(XK[i]-XK[0])+((X0[-1]-X0[0])-(XK[i]-XK[0]))/2S_Sub = np.abs(S_Sub)T = (1+S+SK)/(1+S+SK+S_Sub)if T >= 0.9:print ('精度为一级')print ('可以用GM(1,1)模型\nX(k) = ',X0[0]-a[1]/a[0],'exp(',-a[0],'(k-1))',a[1]/a[0])elif T >= 0.8:print ('精度为二级')print ('可以用GM(1,1)模型\nX(k) = ',X0[0]-a[1]/a[0],'exp(',-a[0],'(k-1))',a[1]/a[0])elif T >= 0.7:print ('精度为三级')print ('谨慎用GM(1,1)模型\nX(k) = ',X0[0]-a[1]/a[0],'exp(',-a[0],'(k-1))',a[1]/a[0])elif T >= 0.6:print ('精度为四级')print ('尽可能不用GM(1,1)模型\nX(k) = ',X0[0]-a[1]/a[0],'exp(',-a[0],'(k-1))',a[1]/a[0])X = Series(np.array([259.4489]),index=pd.period_range('','',freq = 'A-DEC'))X_Return = X_Return.append(X)print (X_Return)B = pd.DataFrame([X0,X_Return],index=['X0','X_Return'])B = np.transpose(B)B.plot()
程序运行截图如下:
实际样本曲线与灰色系统系统预测的曲线如下：
