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

Posted 2020-11-22

 来源公式推导连接

　　https://blog.csdn.net/qq_36387683/article/details/88554434

 

关键词：灰色预测 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和numba加速等

# coding:utf-8
import numpy as np
import numba as nb

class GM():
    def __init__(self):
        pass
    def fit(self,df):
        ‘‘‘
        :param df: 预测的数据;可以传 list 或者是 ndarray(numpy) 一维序列，等于一个向量
        :param n: 预测的个数
        :return: 预测的序列后面是预测值
        ‘‘‘
        self.x,self.consult = self.sigmod(np.array(df))

        z_1 = self.next_to_mean(np.cumsum(self.x))

        self.coefficient = self.coefficient_a_b(self.x,z_1).reshape(-1)

        del z_1

    # 归一化
    nb.jit()
    def sigmod(self,_dt):
        return _dt/(np.max(_dt) - np.min(_dt)),np.max(_dt) - np.min(_dt)

    # 计算紧邻均值数列
    nb.jit()
    def next_to_mean(self,x_1):
        n = len(x_1)
        z_1 = np.empty(n-1)
        for i in range(1,n):  # 下标从0开始，取不到最大值
            z_1[i - 1] = 0.5 * x_1[i] + 0.5 * x_1[i - 1]
        return z_1

    # 发展系数和作用量的计算 (a,b).T
    nb.jit()
    def coefficient_a_b(self,x, z_1):
        n = len(z_1)
        B = np.hstack((np.array(z_1 * -1).reshape((-1, 1)), np.ones((n, 1))))
        Y = np.array(x[1:]).reshape((-1, 1))
        # 返回的是a和b的向量转置，第一个是a 第二个是b；
        return np.dot(np.dot(np.linalg.inv(np.dot(B.T, B)), B.T), Y)

    # 预测第一个结果
    nb.jit()
    def predict(self,m):
        ‘‘‘
        :param m: 预测个数：在原数据的基础上
        :return: 返回一个一维的 ndarray 原数据 后面追加预测 个数的数据 
        ‘‘‘
        n = len(self.x)
        resut = np.empty(n + m)
        resut[:n] = self.x
        for i in range(m):
            # 伪公式： x(k+1) = f(x)*exp(-ak)
            resut[n+i] = (self.x[0] - (self.coefficient[1] / self.coefficient[0])) *                            (1 - np.exp(self.coefficient[0])) * np.exp(-1 * self.coefficient[0] * (n + i))

        # resut[-1] = (x[0]-coefficient[1]/coefficient[0])*np.exp(-1*coefficient[0]*(n))+coefficient[1]/coefficient[0]-x[-1]

        return resut*self.consult