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

 来源公式推导连接

　　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)模型适合具有较强的指数规律的数列，只能描述单调的变化过程。已知元素序列数据：x^{(0)}=(X^{(0)}(1),X^{(0)}(2),...,X^{(0)}(n))

做一次累加生成（1-AGO）序列：

x^{(1)}=(X^{(1)}(1),X^{(1)}(2),...,X^{(1)}(n))
其中

，x^{(1)}(k)=\sum_{i=1}^kx^{(0)}(i), \quad k=1,...,n

令Z^{(1)}为X^{(1)}的紧邻均值生成序列：

Z^{(1)}=(z^{(1)}(2), z^{(1)}(3), ...,z^{(1)}(n))

其中，

z^{(1)}(k)=0.5x^{(1)}(k)+0.5x^{(1)}(k-1)

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

x^{(0)}(k)+az^{(1)}(k)=b

其中，a为发展系数，b为灰色作用量。设\hat{a}为待估参数向量，即\hat{a}=(a,b)^T，则灰微分方程的最小二乘估计参数列满足\hat{a}=(B^TB)^{-1}B^TY_n

其中

B=\left[ \begin{matrix} -z^{(1)}(2) , 1 \\ -z^{(1)}(3) , 1 \\ ... , ... \\ -z^{(1)}(n) , 1 \end{matrix}\right]， Y_n=\left[ \begin{matrix} x^{(0)}(2)\\ x^{(0)}(3) \\ ... , \\ x^{(0)}(n) \end{matrix}\right]
再建立灰色微分方程的白化方程（也叫影子方程）：

\frac{dx^{(1)}}{dt}+ax^{(1)}=b

白化方程的解（也叫时间响应函数）为

\hat{x}^{(1)}(t)=(x^{(1)}(0)-\frac{b}{a})e^{-at}+\frac{b}{a}

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

\hat{x}^{(1)}(k+1) = [x^{(1)}(0)-\frac{b}{a}]e^{-at}+\frac{b}{a}取x^{(1)}(0)=x^{(0)}(1)，

则\hat{x}^{(1)}(k+1)=[x^{(0)}(1)-\frac{b}{a}]e^{-ak}+\frac{b}{a},\quad k=1,....n-1

再做累减还原可得

\hat{x}^{(0)}(k+1)=\hat{x}^{(1)}(k+1)-\hat{x}^{(1)}(k)=[x^{(0)}(1)-\frac{b}{a}](1-e^a)e^{-ak}, \quad k=1,...,n-1即为预测方程。

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

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

3.算法步骤

(1) 数据的级比检验
  为了保证灰色预测的可行性，需要对原始序列数据进行级比检验。
  对原始数据列X^{(0)}=(x^{(0)}(1),x^{(0)}(2),...,x^{(0)}(3))，

计算序列的级比：\lambda(k)=\frac{x^{0}(k-1)}{x^{(0)}(k)}, \quad k=2,...,n  

若所有的级比\lambda(k)都落在可容覆盖\Theta=(e^{-2/(n+1)},e^{2/(n+2)})内，则可进行灰色预测；否则需要对X^{(0)}做平移变换，Y^{(0)}=X^{(0)}+c，使得Y^{(0)}满足级比要求。

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

(3) 检验预测值：

① 相对残差检验，计算

\varepsilon(k)=\frac{x^{(0)}(k-1)}{x^{(0)}(k)}, \quad k=2,...,n  

若\varepsilon(k)<0.2 ，则认为达到一般要求，若\varepsilon(k)<0.1 ，则认为达到较高要求；

② 级比偏差值检验

  根据前面计算出来的级比\lambda(k), 和发展系数a, 计算相应的级比偏差：

\rho(k)=1-(\frac{1-0.5a}{1+0.5a})]\lambda(k)  

若\rho(k)<0.2, 则认为达到一般要求，若\rho(k)<0.1, 则认为达到较高要求。

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

三、程序实现 python实现

　　需要安装numpy和Torch加速等

 

?
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
# 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))

　　

posted @   算法浪客  阅读(13209)  评论(0)  编辑  收藏  举报
编辑推荐：
· .NET Core 中如何实现缓存的预热？
· 从 HTTP 原因短语缺失研究 HTTP/2 和 HTTP/3 的设计差异
· AI与.NET技术实操系列：向量存储与相似性搜索在 .NET 中的实现
· 基于Microsoft.Extensions.AI核心库实现RAG应用
· Linux系列：如何用heaptrack跟踪.NET程序的非托管内存泄露
阅读排行：
· TypeScript + Deepseek 打造卜卦网站：技术与玄学的结合
· Manus的开源复刻OpenManus初探
· 三行代码完成国际化适配，妙~啊~
· .NET Core 中如何实现缓存的预热？
· 如何调用 DeepSeek 的自然语言处理 API 接口并集成到在线客服系统
点击右上角即可分享
微信分享提示
SQL AI 助手