分数阶灰色模型的python实现

一、前言

分数阶灰色模型是在传统灰色模型的基础上引入分数阶累加，从而优化传统的灰色模型建模机制，得到更好的预测结果。
本文会首先介绍分数阶灰色模型的建模思想，再从公式到模型，一步步建立完善的FGM（1,1）模型的python代码。全手打，如有指教，请下方评论或邮件。

二、分数阶模型的建模思想

分数阶灰色模型是Wu（2013）在Commun Nonlinear Sci Numer Simulat上发表的文章，Grey system model with the fractional order accumulation。有兴趣的小伙伴可以去下载，(http://dx.doi.org/10.1016/j.cnsns.2012.11.017)
       首先，大家需要了解传统灰色模型的公式，这里不多做叙述[1-2]。灰色模型适应于小样本的单变量预测，其建模思想是将原数据一阶累加，再通过公式拟合为指数函数。那么，此中就存在一个问题，如果原数据规律不符合指数规律，预测的效果就会很不理想。很多论文也指出其模型病态性在于累加过程和最小二乘法求解参数。既然如此，可不可以找到一个方法，使累加的数据符合指数规律？解决办法自然有，比如分数阶累加，新信息优先累加，卷积累加，倒叙累加等等。有兴趣的小伙伴可以去知网上多搜搜相关进展。
       Wu在论文中提出了两条灰色模型的建模原则
       1. 信息之间存在差异，即不同时间节点的数据，对未来数据的影响是不同的。
       2. 新信息应该被赋予更大的权重。即最近时刻的事件，对我们下一步的决断有更大的影响。
       因此，分数阶累加的公式如下：
（式1）
其中

       其还原公式主要有两种方法：（1）将得到的Xr序列再经过（1-r）阶累加，得到1阶累加序列，再累减。如式2和式3。（2）对Xr序列进行（-r）阶累加（如式4），这种方法编程时可以直接用累加矩阵的逆矩阵乘以Xr的竖列形式。
（式2）
（式3）
（式4）
       注：相关模型公式可以见《灰色系统理论及其应用》和《分数阶灰色模型及其在装备费用测算中的应用》[3]。其中，式1表示序列中每个元素累加的过程，这个过程可以用矩阵来表示和计算，与式4的左面部分等价。还原过程即式2和3与式4的右面部分等价。两种方式理解任何一个都行。通俗来讲，数学方面的理解关系到您编程时，是建立for循环一个元素一个元素的进行赋值，还是直接利用矩阵进行赋值计算。下文中的代码主要利用到式1,2和3。

三、传统灰色模型的代码

在这一章，我们建立了传统灰色模型的代码，并把它定义成一个模块。当我们需要使用优化模型时，可以直接调用这个类中的相关函数，只需要改动相关变量的赋值即可。文章中有很多中文注释，相信理解起来不会费劲。需要提前加载好numpy包。

#!/usr/bin/env python
# -*- coding: UTF-8 -*-##最传统的灰色模型
import numpy as np
class Grey_model(object):def __init__(self,input_value):#初始化过程，先建立好变量空间，即累加序列，背景值序列，B和Y矩阵，拟合序列，预测值序列等。self.input_value=input_valueself.accumulation_value=np.zeros(len(input_value))self.background_value=np.zeros(len(input_value)-1)self.y_matrix_value=np.mat(np.zeros((len(input_value)-1,1)))self.b_matrix_value=np.mat(np.ones((len(input_value)-1,2)))self.accumulation()#计算累加序列def accumulation(self):for i in range(len(self.input_value)):self.accumulation_value[i]=np.sum(self.input_value[0:i+1])#计算Z值def background_values(self):for i in range(len(self.accumulation_value)-1):self.background_value[i]=(self.accumulation_value[i]+self.accumulation_value[i+1])/2#计算B矩阵def b_matrix(self):self.background_values()for i in range (self.b_matrix_value.shape[0]):self.b_matrix_value[i,0]=-self.background_value[i]#计算Y矩阵def y_matrix(self):for i in range(self.y_matrix_value.shape[0]):self.y_matrix_value[i]=self.accumulation_value[i+1]-self.accumulation_value[i]#计算参数矩阵U：U=(B^T*Y)^-1*B^T*Ydef u_matrix(self):self.y_matrix()self.b_matrix()self.u_matrix_values=(self.b_matrix_value.T*self.b_matrix_value)**-1 * self.b_matrix_value.T * self.y_matrix_value#下面把矩阵格式转化为数组格式，再转化为列表格式self.u_matrix_array=np.array(self.u_matrix_values.T)[0]self.u_matrix_value=self.u_matrix_array.tolist()#计算预测值def predict(self,number_of_forecast):self.u_matrix()self.predicted_accumulation_value=[]#使用了float("%.2f"% 来调整输出值的小数的个数self.predicted_value=[float("%.2f"%(self.input_value[0]))]for i in range(len(self.input_value)+int(number_of_forecast)):self.predicted_accumulation_value.append(((self.accumulation_value[0]-self.u_matrix_value[1]/self.u_matrix_value[0])*np.exp(-self.u_matrix_value[0]*(i)))+self.u_matrix_value[1]/self.u_matrix_value[0])#上式中e^-at,由于python中从0开始算，故t减去1for i in range(len(self.predicted_accumulation_value)-1):self.predicted_value.append(float("%.2f"%(self.predicted_accumulation_value[i+1]-self.predicted_accumulation_value[i])))#计算误差def test_error(self):MAPE_list=[]for i in range(len(self.input_value)):MAPE_list.append(abs((self.predicted_value[i]-self.input_value[i])/self.input_value[i]))self.MAPE=str(100*(np.mean(MAPE_list[1:])))[:5]+ '%'#下面是个测试实验，input_value是输入序列，预测的个数为3.
if __name__ == '__main__':input_value=[29.2,33.9,39.7,46.8,56.1,69.5,80.7,87.5,107.5]number_of_forecast=3A=Grey_model(input_value)A.predict(number_of_forecast)A.test_error()print('预测值')print(A.predicted_value)print('\nMAPE值')print(A.MAPE)

相信大家可以很轻松的看懂以上代码，这也是本青铜级编程选手第一次独立编计算程序，学艺不精敬请见谅。将代码运行一遍，可以得到以下输出结果：

预测值
[29.2, 34.72, 40.78, 47.89, 56.24, 66.06, 77.58, 91.12, 107.01, 125.68, 147.61, 173.37]MAPE值
2.64%

四、分数阶灰色模型的代码

上一章的代码文件保存为GM.py。本章代码保存为FGM.py。需要提前安装scipy包。这里编程时使用了gamma函数，见式5.
（式5）
相关代码如下：

#!/usr/bin/env python
# -*- coding: UTF-8 -*-from GM import Grey_model
import numpy as np
from scipy.special import gamma#更新累加序列（式1）
def updata_fractional_accumulation(input_value,order):accumulation_value=np.zeros(len(input_value))for i in range(len(input_value)):for j in range(i+1):tmp=gamma(i-j+order-1+1)/(gamma(i-j+1)*gamma(order-1+1))accumulation_value[i]+=tmp*input_value[j]return accumulation_value#更新还原方法（参考式2和3）
def updata_predicted_results(predicted_accumulation_value,order):if order!=1:predicted_one_accumulation_value=updata_fractional_accumulation(predicted_accumulation_value,1-order)else:predicted_one_accumulation_value=predicted_accumulation_value predicted_value=[float("%.2f"%(predicted_one_accumulation_value[0]))]for i in range(len(predicted_accumulation_value)-1):predicted_value.append(float("%.2f"%(predicted_one_accumulation_value[i+1]-predicted_one_accumulation_value[i])))return predicted_value#输入变量，order表示分数阶的阶数
input_value=[29.2,33.9,39.7,46.8,56.1,69.5,80.7,87.5,107.5]
number_of_forecast=3
order=0.1
A=Grey_model(input_value)#引入GM模型的计算模块
A.accumulation_value=updata_fractional_accumulation(input_value,order)#更新累加值
A.predict(number_of_forecast)
A.predicted_value= updata_predicted_results(A.predicted_accumulation_value,order)#更新预测值的还原方式
A.test_error()#误差用MAPE值来衡量
print('预测值')
print(A.predicted_value)
print('\nMAPE值')
print(A.MAPE)

运行以上代码（注：GM.py和FGM.py两个文件要放在一个文件夹下），可以得到以下结果：

预测值
[29.2, 33.51, 39.76, 47.29, 56.08, 66.23, 77.93, 91.37, 106.81, 124.53, 144.86, 168.2]MAPE值
1.94%

至此，我们可以发现，同样的数据，分数阶累加通过调整参数，能够得到误差更小的预测结果。当然，分数阶累加最大特点，是符合新信息优先的原则。

五、结论

灰色模型自创立以来，在很多领域都得到了很好的应用。其最大的特点在于累加的过程能过把数据的规律性体现出来。相对于其他ARIMA模型，神经网络算法等，其使用更简单，对于贫信息的处理也更为有效。当然，目前为止使用最广泛的预测方法是指数平滑，M3和M4竞赛都是这个算法赢的似乎。M3时，有人用指数平滑和线性拟合做了一个结合，并通过一个数据识别分类的框架，判断数据的规律，采用不同的theta算法。M4时，获奖的算法是指数平滑和LSTM的混合算法，即将一系列的预测模型分别预测出不同的结果，用神经网络进行数据加权，得到最佳的预测结果。
       当下，机器学习算法越来越流行，不仅参数求解可以使用算法，更在于模型间的结合也可以使用算法进行训练，能够发挥出每个模型的优点，结果自然更优。毕竟一个模型，它的规律性是几行公式定死的了，一个工具只能处理特定的一类问题，如果把几个工具一起使用，处理的结果肯定会更好。
       目前，灰色模型系列算法发展的很快。相信不久将来，会有适应范围更广，预测结果更精准的灰色模型产生。笔者正在考虑将所有的灰色模型（时序的，多参数的，verhust的，季节模型等）做一个组合框架，并通过神经训练算法将预测结果进行组合加权。
       关于分数阶灰色模型的代码就这么多了，看了其他人写的，似乎很复杂，还要考虑级比检验，我觉得没必要。如果预测的结果不好，那就换个模型就好了咯。毕竟一个模型只是一个工具，这个工具不行就换一个，或者几个工具一起用。另外，分数阶的阶数最好选择0<r<1，大于1的时候可能会拟合结果好，但是预测的结果就会很离谱。

参考文献

[1]灰色模型学习记录（有相关GM模型的公式和matlab代码，需要的自取）
[2]灰色系统理论及其应用
[3]分数阶灰色模型及其在装备费用测算中的应用