来源公式推导连接

关键词：灰色预测 python 实现 灰色预测 GM(1,1)模型 灰色系统 预测 灰色预测公式推导

一、前言

本文的目的是用Python和类对灰色预测进行封装

二、原理简述

1.灰色预测概述

灰色预测是用灰色模型GM(1,1)来进行定量分析的，通常分为以下几类：

(1) 灰色时间序列预测。用等时距观测到的反映预测对象特征的一系列数量(如产量、销量、人口数量、存款数量、利率等)构造灰色预测模型，预测未来某一时刻的特征量，或者达到某特征量的时间。

(2) 畸变预测(灾变预测)。通过模型预测异常值出现的时刻，预测异常值什么时候出现在特定时区内。

(3) 波形预测，或称为拓扑预测，它是通过灰色模型预测事物未来变动的轨迹。

(4) 系统预测，对系统行为特征指标建立一族相互关联的灰色预测理论模型，在预测系统整体变化的同时，预测系统各个环节的变化。

上述灰色预测方法的共同特点是：

(1)允许少数据预测；

(2)允许对灰因果律事件进行预测，例如：

灰因白果律事件：在粮食生产预测中，影响粮食生产的因子很多，多到无法枚举，故为灰因，然而粮食产量却是具体的，故为白果。粮食预测即为灰因白果律事件预测。

白因灰果律事件：在开发项目前景预测时，开发项目的投入是具体的，为白因，而项目的效益暂时不很清楚，为灰果。项目前景预测即为灰因白果律事件预测。

(3)具有可检验性，包括：建模可行性的级比检验(事前检验)，建模精度检验(模型检验)，预测的滚动检验(预测检验)。

2.GM(1,1)模型理论

GM(1,1)模型适合具有较强的指数规律的数列，只能描述单调的变化过程。已知元素序列数据：

做一次累加生成(1-AGO)序列：

其中

，

令

为

的紧邻均值生成序列：

其中，

建立GM(1,1)的灰微分方程模型为：

其中，

为发展系数，

为灰色作用量。设

为待估参数向量，即

，则灰微分方程的最小二乘估计参数列满足

其中

再建立灰色微分方程的白化方程(也叫影子方程)：

白化方程的解(也叫时间响应函数)为

那么相应的GM(1,1)灰色微分方程的时间响应序列为：

取

，

则

再做累减还原可得

即为预测方程。

注1：原始序列数据不一定要全部使用，相应建立的模型也会不同，即

和

不同；

注2：原始序列数据必须要等时间间隔、不间断。

3.算法步骤

(1) 数据的级比检验

为了保证灰色预测的可行性，需要对原始序列数据进行级比检验。

对原始数据列

，

计算序列的级比：

若所有的级比

都落在可容覆盖

内，则可进行灰色预测；否则需要对

做平移变换，

，使得

满足级比要求。

(2) 建立GM(1,1)模型，计算出预测值列。

(3) 检验预测值：

① 相对残差检验，计算

若

，则认为达到一般要求，若

，则认为达到较高要求；

② 级比偏差值检验

根据前面计算出来的级比

, 和发展系数

, 计算相应的级比偏差：

若

, 则认为达到一般要求，若

, 则认为达到较高要求。

(4) 利用模型进行预测。

三、程序实现 python实现

需要安装numpy和Torch加速等

# condig:utf-8

import torch as th

import numpy as np

class GM():

def __init__(self):

# 判断是否可用 gpu 编程 , 大量级计算使用GPU

self._is_gpu = False # th.cuda.is_available()

def fit(self,dt:list or np.ndarray):

self._df :th.Tensor = th.from_numpy(np.array(dt,dtype=np.float32))

if self._is_gpu:

self._df.cuda()

self._n:int = len(self._df)

self._x,self._max_value = self._sigmod(self._df)

z:th.Tensor = self._next_to_mean(th.cumsum(self._x,dim=0))

self.coef:th.Tensor = self._coefficient(self._x, z)

del z

self._x0:th.Tensor = self._x[0]

self._pre:th.Tensor = self._pred()

# 归一化

def _sigmod(self,x:th.Tensor):

_maxv:th.Tensor = th.max(x)

return th.div(x,_maxv),_maxv

# 计算紧邻均值数列

def _next_to_mean(self, x_1:th.Tensor):

z:th.Tensor = th.zeros(self._n-1)

if self._is_gpu:

z.cuda()

for i in range(1,self._n): # 下标从0开始，取不到最大值

z[i - 1] = 0.5 * x_1[i] + 0.5 * x_1[i - 1]

return z

# 计算系数 a,b

def _coefficient(self,x:th.Tensor,z:th.Tensor):

B:th.Tensor = th.stack((-1*z, th.ones(self._n-1)),dim=1)

Y:th.Tensor = th.tensor(x[1:],dtype=th.float32).reshape((-1,1))

if self._is_gpu:

B.cuda()

Y.cuda()

# 返回的是a和b的向量转置，第一个是a 第二个是b；

return th.matmul(th.matmul(th.inverse(th.matmul(B.t(), B)), B.t()),Y)

def _pred(self,start:int=1,end:int=0):

les:int = self._n+end

resut:th.Tensor = th.zeros(les)

if self._is_gpu:

resut.cuda()

resut[0] = self._x0

for i in range(start,les):

resut[i] = (self._x0 - (self.coef[1] / self.coef[0])) * \

(1 - th.exp(self.coef[0])) * th.exp(-1 * self.coef[0] * (i))

del les

return resut

# 计算绝对误差

def confidence(self):

return round((th.sum(th.abs(th.div((self._x-self._pre),self._x)))/self._n).item(),4)

# 预测个数，默认个数大于等于0，

def predict(self,m:int=1,decimals:int=4):

y_pred:th.Tensor = th.mul(self._pre,self._max_value)

y_pred_ = th.zeros(1)

if m<0:

return "预测个数需大于等于0"

elif m>0:

y_pred_:th.Tensor = self._pred(self._n,m)[-m:].mul(self._max_value)

else:

if self._is_gpu:

return list(map(lambda _: round(_, decimals), y_pred.cpu().numpy().tolist()))

else:

return list(map(lambda _:round(_,decimals),y_pred.numpy().tolist()))

# cat 拼接 0 x水平拼接，1y垂直拼接

result:th.Tensor = th.cat((y_pred,y_pred_),dim=0)

del y_pred,y_pred_

if self._is_gpu:

return list(map(lambda _: round(_, decimals), result.cpu().numpy().tolist()))

return list(map(lambda _:round(_,decimals),result.numpy().tolist()))

if __name__=="__main__":

ls = np.arange(91,100,2)

print(type(ls))

# ls = list(range(91, 100, 2))

gm = GM()

gm.fit(ls)

print(gm.confidence())

print(ls)

print(gm.predict(m=2))