知然算法【2】灰色模型GM(1,1)

AnFany

已于 2023-05-08 18:01:09 修改

阅读量433

点赞数

分类专栏：知然算法文章标签： GM；灰色模型；灰色预测 GDP预测算法

于 2021-03-25 14:30:45 首次发布

本文链接：https://blog.csdn.net/qq_32882309/article/details/115207873

版权

知然算法专栏收录该内容

2 篇文章 2 订阅

订阅专栏

这篇博客介绍了灰色模型GM(1,1)，适用于短期预测，数据样本大于4即可建模。文章详细阐述了模型背景、步骤，包括级比检验、数据处理、参数计算等，并提供了Python实现代码，通过中美GDP预测对比展示了模型效果。" 81330173,7923641,解决微信客服消息发送重复回调问题,"['微信开发', '客服消息', 'PHP接口', '微信回调机制', '微信API']

摘要由CSDN通过智能技术生成

知然算法【2】灰色模型GM(1,1)

在这里插入图片描述

0、算法适用

适用于短期预测；
样本条数大于4即可建模；
不适用于较大震荡性数据；

1、符号说明

原始数据序列： $x^{0}(k), k=1,2,……,n;n$ 为数据样本条数；
一次累加数据序列： $x^{1}(k)=\sum_{j=1}^{k}x^{0}(j), k=1,2,……,n；$
累加数据预测序列： $\hat{x}^{1}(k), k=1,2,……,n$

2、模型背景

白色系统：系统的内部信息是完全已知的。黑色系统：系统的内部信息是完全未知的，只能通过它与外界的联系来加以观测研究，也就是可以知道该系统的输入输出关系，但不知道它内部是如何实现这种关系的。灰色系统：介于前面两者之间，系统内的一部分信息是已知的，另一部分信息是未知的。

灰色系统理论是由我国著名学者邓聚龙教授于1982年创立的。该理论是将具有潜在规律的时间序列数据经累加生成后，使其变为具有指数增长规律的上升形状数列，而一阶微分方程解的形式是指数增长形式，所以可对累加生成后的数列建立微分方程模型。

3、模型步骤

要求每个数据之间的时间间隔必须是一致的，并且不存在缺失值。

STEP1: 对数据序列进行级比检验，当所有级比值 $\lambda(j)$ 在区间 $e^{-2/(n+1)}, e^{2/(n+1)}]$ 内，可建立模型；否则需要对数据进行平移转
换，也就是将该数据序列同时加上一个适当的系数 $\beta$ (需要保存该系数，预测时应再减去该系数)，使得级比值全部在区间内，其中级比值：
$\lambda(j) = \frac{x^{0}(j-1)}{x^{0}(j)}, j=2,3,……,n$
此时新构成的目标数据序列为：
$x^{0}(k)=x^{0}(k)+\beta$
STEP2: 对数据序列计算一次累加数据序列 $x^{1}(k)$ ；
STEP3: 对一次累加序列计算均值序列：
$Z^{1}(k) = \frac{x^{1}(k-1)+x^{1}(k)}{2}, k=2,3,……,n$
STEP4: 计算模型参数构建的矩阵：

$Y=\begin{bmatrix} x^{0}(2) \\ x^{0}(3) \\ … \\x^{0}(n)\end{bmatrix}$

$\begin{bmatrix} -Z^{1}(2) & 1 \\ -Z^{1}(3) & 1 \\ … & …\\-Z^{1}(n) & 1 \end{bmatrix}$

计算参数

$(X^{T} X)^{-1} X^TY,其中W=\begin{bmatrix} a\\ b\end{bmatrix}$

其中 $a$ 为发展系数， $b$ 为作用量；

STEP5: 时间响应表达式：

$\hat{x}^{1}(t) = ({x}^{0}(1) - \frac{b}{a})e^{-a(t-1)} + \frac{b}{a}, t=1,2…,n,…$

目标序列的预测模型为：

$\hat{x}^{0}(1)={x}^{0}(1)$ $\hat{x}^{0}(t) = \hat{x}^{1}(t) - \hat{x}^{1}(t-1), t=2, 3, …,n,…$

4、模型效果检验

4.1 残差

$\eta_{j} = {x}^{0}(j) - \hat{x}^{0}(j), j=1,2,3,…,n$

4.2 相对误差

$\bigtriangleup_{j} = \frac{\eta_{j}} {{x}^{0}(j)}, j=1,2,3,…,n$

4.3 平均相对误差

$\epsilon = \frac{\sum_{j=1}^{n}|\bigtriangleup_{j}|}{n}$

3.4 级比偏差

$\rho_{j} = 1- \frac{1-0.5a}{1+0.5a}\lambda(j), j=2,3,…n$

4.5 后验差比

$c=\frac{残差序列方差}{原始数据序列方差}=\frac{S_{\eta}}{S_{x^{0}(k)}}=\frac{\frac{1}{n}\sum_{j=1}^{n}(\eta_{j} -\bar{\eta})^{2}}{\frac{1}{n}\sum_{j=1}^{n}(x^{0}(j) -\bar{x})^{2}}$

其中 $\bar{\eta}=\frac{1}{n}\sum_{j=1}^{n}\eta_{j},\bar{x}=\frac{1}{n}\sum_{j=1}^{n}x^{0}(j)$

4.6 小概率误差

$p=P\{|\eta_{j} - \bar{\eta}|<0.6745 \sqrt{S_{\eta}}\}$

模型效果评判

模型效果评判	$p$	$c$
好	$\leq p$	$\leq 0.35$
合格	$\leq p < 0.95$	$\leq 0.5$
勉强	$\leq p < 0.8$	$\leq 0.65$
不合格	$p < 0.7$	$c > 0.65$

根据p和c值综合判断，选取两者的模型评判结果较差的作为最终的评判结果。

5、Python3实现

import pandas as pd
import numpy as np
from math import e
from IPython.display import Latex
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['FangSong']  # 指定默认字体
plt.rcParams['axes.unicode_minus'] = False

class GM:
    
    def __init__(self, data, M=10, A=None):
        
        self.data = data  # df格式，第一列时间序列，第二类是数据系列
        self.M = M        # 预测未来的期数，
        self.A = A        # 数据系列达到A值需要的期数
        
        # 打印结果标识
        self.title_sign = 1
        
        # 模型精度评判标准字典
        self.signdict ={0:'好',1:'合格',2:'勉强',3:'不合格'}
        
        self.p_rule = [[0.95,0.8,0.7],'<']
        self.c_rule = [[0.35,0.5,0.65],'>']
    
    # 根据数值和规则评判，获得最终结果
    def getmodelresult(self,n,exrule):
        nulist, sign = exrule
        if sign == '<':
            for kin, kva in enumerate(nulist):
                if n > kva:
                    return kin
            return  3
         else:
             for kin, kva in enumerate(nulist):
                if n < kva:
                    return kin
             return  3


    # 级比检验，不满足在计算平移参数
    def jibicheck(self):
        # 目标数据序列
        targetdata = self.data.values[:, 1].T
        # 计算级比值区间
        n = len(targetdata)
        min_tap = round(np.power(e,(-2/(n+1))), 6)  # 最小值
        max_tap = round(np.power(e, (2/(n+1))), 6)  # 最大值
        
        checksign = 1
        multi = 0
        start_number = np.mean(targetdata) # 平移系数的起始值
        while checksign:
            # 计算平移参数
            new_c = round(multi * start_number, 1)
            new_last_data = targetdata + new_c
            checksign = 0
            # 级比数据
            jibidata = ['%.4f' % round(f, 4) for f in new_last_data[:-1] / new_last_data[1:]]
            if not multi:
                self.jibi = jibidata
            for kk in jibidata:
                if float(kk) < min_tap or float(kk) > max_tap:
                    checksign = 1
            if not checksign:
                self.transc = new_c
                break
            multi += 0.001
        # 平移参数
        self.transpara = new_c
        # 打印原始数据
        print('(%s) 级比检验\n' % self.title_sign)
        self.title_sign += 1
        jibidatadf = pd.DataFrame()
        jibidatadf['原始数据'] = [str(f) for f in targetdata]
        if multi:
            jibidatadf['平移数据(参数%s)'% new_c] = [str(k) for k in new_last_data]
  
        jibidatadf['级比值区间[%s, %s]'%(min_tap, max_tap)] = ['_'] + [str(j) for j in jibidata]
        print(jibidatadf)
        if multi:
            return jibidatadf, '平移数据(参数%s)'% new_c
        else:
            return jibidatadf, '原始数据'
    

    # 根据参数输出预测方程
    def outputstr(self):
        base_str = r'$\hat{x}^{1}(t)=(x^{0}(t)-\frac{%s}{%s})e^{%s(t-1)} + \frac{%s}{%s},t=1,2,…k,…$' % (self.B, self.a, -self.a, self.B, self.a)
        return base_str
        
    # 建立灰色预测模型GM(1,1)
    def build_gm(self):
        gmdata, name=  self.jibicheck()
        # 累加数据
        adddata = pd.DataFrame()
        adddata['累加数据'] = np.cumsum(gmdata[name].values.astype(float))
        
        print('(%s) 累加生成序列\n' % self.title_sign)  
        self.title_sign += 1
        print(adddata)
        
        # 对目标特征计算均值序列
        print('(%s) 均值序列\n' % self.title_sign)
        self.title_sign += 1
        avgdata = pd.DataFrame()
        adddata['累加数据'] = adddata['累加数据'].astype(float)
        avgdata['均值数据'] = [str(k) for k in (adddata['累加数据'].values[:-1] + adddata['累加数据'].values[1:])/2]
        print(avgdata)
        # 构建矩阵
        print('(%s) 构建矩阵\n' % self.title_sign)

        print('(%s.1) 矩阵X\n' % self.title_sign)

        X = np.ones((len(avgdata), 1))
        cx = np.array([0-avgdata.values[:,0].astype(float)]).T.astype(str)
        X = np.hstack((cx,X))
        print(X)
        print('(%s.1) 矩阵Y\n' % self.title_sign)
     
        Y = np.array([gmdata[name].values[1:]]).T
        print(Y)
        self.title_sign +=1
        
        # 计算参数
        print('(%s)计算参数a和B\n' % self.title_sign)
        
        X= X.astype(float)
        Y = Y.astype(float)
        
        M = np.array(np.matmul(np.matmul(np.matrix(np.matmul(X.T, X)).I, X.T),Y))
        self.a = M[0][0]
        self.B= M[1][0]
        
        print('a为',self.a)
        print('B为',self.B)
        
        self.title_sign += 1
        
        return self.outputstr()
        
    # 对模型的预测效果进行分析
    def predictmodel(self):
        print('(%s) 预测效果分析\n' % self.title_sign)

        predictdf = pd.DataFrame()
        # 原始数据
        predictdf['原始数据'] = self.data.values[:,1]
        predata = [self.data.values[:,1][0]]
        # 开始预测
        fsum = predata[0]
        # 判断期数
        if self.M:
            pass
        else:
            self.M = 0
        # 判断达到的值
        if self.A:
            sign = 1
        else:
            sign = 0
            
        pbb = 0
        while pbb < len(self.data)-1+self.M or sign:

            xp = (predata[0]-self.B/self.a) * np.power(e, -self.a * (pbb+1)) + self.B/self.a
            predictnumber = xp-fsum-self.transpara
            if sign:
                if predictnumber >=self.A:
                    sign = 0
                    self.daA = pbb+1 
                    self.daAnum = predictnumber
            if pbb > 99:
                self.daA = '100'
                break
            predata.append(round(predictnumber, 4))
            fsum = xp
            pbb += 1
            
        predictdf['预测数据'] = np.array(predata[:len(self.data)])
        
        predictdf['残差'] = predictdf['原始数据'] -  predictdf['预测数据']
        predictdf['相对误差(%)'] = [round(d,2) for d in (predictdf['残差'] /  predictdf['原始数据'])*100]
        predictdf['级比偏差'] = ['-'] + [1-(1-self.a)/(1+self.a) * float(h) for h in self.jibi]
        print(predictdf)
        
        print('(%s.1) 平均相对误差(%%):%s' % (self.title_sign, np.mean(np.abs(predictdf['相对误差(%)']))))
        
        # 计算方差比
        varee = np.var(predictdf['残差'].values)
        datavar = np.var(predictdf['原始数据'].values)
        c = varee / datavar
        print('(%s.2) 后验差比:%s' % (self.title_sign,c))
        avge = predictdf['残差'].values - np.mean(predictdf['残差'].values)
        p = len([h for h in avge if abs(h)<0.6745*(datavar ** 0.5)]) / len(avge)
        print('(%s.3) 小概率误差:%s' % (self.title_sign, p))
        
        modelclass = self.signdict[max(self.getmodelresult(p, self.p_rule), self.getmodelresult(c, self.c_rule))]
        print('(%s.4) 模型结果评判【好，合格，勉强，不合格】:%s' % (self.title_sign, modelclass))
        
        
        self.title_sign += 1
        
        print('(%s) 模型预测以及图示\n' % self.title_sign)
        
        if self.M:
            print('预测期数%s\n' % self.M)
            all_data = pd.DataFrame()

            all_data[self.data.keys()[0]] = list(self.data.values[:,0]) + ['后(%s)' % (h+1) for h in range(self.M)]
            all_data[self.data.keys()[1]] = list(self.data.values[:,1]) + ['-']*self.M
            all_data['预测数据'] = np.array(predata[:(len(self.data)-1+self.M+1)]).astype(str)

            print(all_data)

            plt.figure(figsize=(10, 10))
            plt.plot(predictdf['原始数据'].astype(float),label='真实数据',marker='s')
            plt.plot(all_data['预测数据'].astype(float),label='预测数据',marker='8', **dict(linestyle=':', color='tab:orange', markersize=5))
            plt.ylabel(self.data.keys()[1])
            plt.xlabel(self.data.keys()[0])
            plt.xticks(list(range(len(all_data)))[::1], all_data[self.data.keys()[0]][::1], rotation=50)
            plt.legend()
            plt.show()
        
        if self.A:
            print('达到预测值%s\n' % self.A)
            if self.daA <= len(self.data):
                print('已经在期数%s达到，值为%s\n' % (self.daA+1, self.daAnum))
            else:
                print('根据预测所需期数为%s,值为%s\n' % (self.daA-len(self.data)+1, self.daAnum))
                
        return predata

5、数据实例(中美GDP预测比对)

# 2001-2019中美GDP数据
GDP = pd.read_excel('GDP(01-19).xlsx')
GDP

	年	中国GDP(美元)	美国GDP(美元)
0	2001	1339395718865	10581821399000
1	2002	1470550015081	10936419054000
2	2003	1660287965662	11458243878000
3	2004	1955347004963	12213729147000
4	2005	2285965892360	13036640230000
5	2006	2752131773355	13814611414000
6	2007	3550342737010	14451858656000
7	2008	4594307032660	14712844084000
8	2009	5101703073086	14448933025000
9	2010	6087163874512	14992052727000
10	2011	7551500124203	15542581104000
11	2012	8532229986993	16197007349000
12	2013	9570406235659	16784849196000
13	2014	10475682920594	17527163695000
14	2015	11061553079876	18224704440000
15	2016	11233276536737	18714960538000
16	2017	12310409370892	19519353692000
17	2018	13894817549374	20580159776000
18	2019	14342903006431	21433226000000

ChinaGDP = GDP[['年','中国GDP(美元)']]
USAGDP = GDP[['年','美国GDP(美元)']]

# 建立GM(1,1)模型
GMANFANY_China = GM(ChinaGDP,M=11,A=30000000000000)  # M为预测的期数。A为预测达到的值
# 预测表达式
prestrChina = GMANFANY_China.build_gm()
ChinaPredict = GMANFANY_China.predictmodel()
print('预测表达式为：')
Latex(r'%s' % prestrChina)

# 建立GM(1,1)模型
GMANFANY_USA = GM(USAGDP,M=11,A=30000000000000)  # M为预测的期数。A为预测达到的值
# 预测表达式
prestrusa = GMANFANY_USA.build_gm()
USAPredict = GMANFANY_USA.predictmodel()
print('预测表达式为：')
Latex(r'%s' % prestrusa)

(1) 级比检验

              原始数据 平移数据(参数7840838422908.5) 级比值区间[0.904837, 1.105171]
0    1339395718865         9180234141773.5                         _
1    1470550015081         9311388437989.5                    0.9859
2    1660287965662         9501126388570.5                    0.9800
3    1955347004963         9796185427871.5                    0.9699
4    2285965892360        10126804315268.5                    0.9674
5    2752131773355        10592970196263.5                    0.9560
6    3550342737010        11391181159918.5                    0.9299
7    4594307032660        12435145455568.5                    0.9160
8    5101703073086        12942541495994.5                    0.9608
9    6087163874512        13928002297420.5                    0.9292
10   7551500124203        15392338547111.5                    0.9049
11   8532229986993        16373068409901.5                    0.9401
12   9570406235659        17411244658567.5                    0.9404
13  10475682920594        18316521343502.5                    0.9506
14  11061553079876        18902391502784.5                    0.9690
15  11233276536737        19074114959645.5                    0.9910
16  12310409370892        20151247793800.5                    0.9465
17  13894817549374        21735655972282.5                    0.9271
18  14342903006431        22183741429339.5                    0.9798
(2) 累加生成序列

            累加数据
0   9.180234e+12
1   1.849162e+13
2   2.799275e+13
3   3.778893e+13
4   4.791574e+13
5   5.850871e+13
6   6.989989e+13
7   8.233504e+13
8   9.527758e+13
9   1.092056e+14
10  1.245979e+14
11  1.409710e+14
12  1.583822e+14
13  1.766988e+14
14  1.956011e+14
15  2.146753e+14
16  2.348265e+14
17  2.565622e+14
18  2.787459e+14
(3) 均值序列

                  均值数据
0    13835928360768.25
1    23242185774048.25
2    32890841682269.25
3    42852336553839.25
4    53212223809605.25
5    64204299487696.25
6    76117462795439.75
7    88806306271221.25
8   102241578167928.75
9   116901748590194.75
10  132784452068701.25
11  149676608602935.75
12  167540491603970.75
13  186149948027114.25
14  205138201258329.25
15  224750882635052.25
16  245694334518093.75
17  267654033218904.75
(4) 构建矩阵

(4.1) 矩阵X

[['-13835928360768.25' '1.0']
 ['-23242185774048.25' '1.0']
 ['-32890841682269.25' '1.0']
 ['-42852336553839.25' '1.0']
 ['-53212223809605.25' '1.0']
 ['-64204299487696.25' '1.0']
 ['-76117462795439.75' '1.0']
 ['-88806306271221.25' '1.0']
 ['-102241578167928.75' '1.0']
 ['-116901748590194.75' '1.0']
 ['-132784452068701.25' '1.0']
 ['-149676608602935.75' '1.0']
 ['-167540491603970.75' '1.0']
 ['-186149948027114.25' '1.0']
 ['-205138201258329.25' '1.0']
 ['-224750882635052.25' '1.0']
 ['-245694334518093.75' '1.0']
 ['-267654033218904.75' '1.0']]
(4.1) 矩阵Y

[['9311388437989.5']
 ['9501126388570.5']
 ['9796185427871.5']
 ['10126804315268.5']
 ['10592970196263.5']
 ['11391181159918.5']
 ['12435145455568.5']
 ['12942541495994.5']
 ['13928002297420.5']
 ['15392338547111.5']
 ['16373068409901.5']
 ['17411244658567.5']
 ['18316521343502.5']
 ['18902391502784.5']
 ['19074114959645.5']
 ['20151247793800.5']
 ['21735655972282.5']
 ['22183741429339.5']]
(5)计算参数a和B

a为 -0.05493103593203573
B为 8281321853114.156
(6) 预测效果分析

              原始数据          预测数据            残差  相对误差(%)        级比偏差
0    1339395718865  1.339396e+12  0.000000e+00     0.00           -
1    1470550015081  7.477894e+11  7.227606e+11    49.15  -0.0166784
2    1660287965662  1.232770e+12  4.275180e+11    25.75   0.0113394
3    1955347004963  1.745136e+12  2.102108e+11    10.75   0.0521941
4    2285965892360  2.286435e+12 -4.687932e+08    -0.02   0.0451617
5    2752131773355  2.858299e+12 -1.061673e+11    -3.86   0.0728447
6    3550342737010  3.462455e+12  8.788743e+10     2.48    0.134685
7    4594307032660  4.100727e+12  4.935801e+11    10.74    0.137364
8    5101703073086  4.775040e+12  3.266628e+11     6.40 -0.00518103
9    6087163874512  5.487431e+12  5.997332e+11     9.85   0.0644728
10   7551500124203  6.240048e+12  1.311452e+12    17.37    0.100193
11   8532229986993  7.035164e+12  1.497066e+12    17.55   0.0120092
12   9570406235659  7.875179e+12  1.695228e+12    17.71   0.0048652
13  10475682920594  8.762627e+12  1.713056e+12    16.35  -0.0198039
14  11061553079876  9.700187e+12  1.361366e+12    12.31  -0.0570866
15  11233276536737  1.069069e+13  5.425871e+11     4.83  -0.0991691
16  12310409370892  1.173712e+13  5.732862e+11     4.66   -0.018576
17  13894817549374  1.284265e+13  1.052171e+12     7.57   0.0110046
18  14342903006431  1.401060e+13  3.323064e+11     2.32  -0.0814207
(6.1) 平均相对误差(%):11.561578947368421
(6.2) 后验差比:0.0172173146519471
(6.3) 小概率误差:1.0
(6.4) 模型结果评判【好，合格，勉强，不合格】:好
(7) 模型预测以及图示

预测期数11

        年       中国GDP(美元)                预测数据
0    2001   1339395718865     1339395718865.0
1    2002   1470550015081   747789421673.4062
2    2003   1660287965662  1232769944703.2812
3    2004   1955347004963  1745136228529.1562
4    2005   2285965892360  2286434685550.8125
5    2006   2752131773355     2858299050602.5
6    2007   3550342737010  3462455311854.9375
7    2008   4594307032660  4100726920154.2817
8    2009   5101703073086   4775040292519.969
9    2010   6087163874512   5487430626413.531
10   2011   7551500124203     6240048042326.0
11   2012   8532229986993   7035164073222.781
12   2013   9570406235659   7875178520433.656
13   2014  10475682920594   8762626696679.501
14   2015  11061553079876   9700187078095.438
15   2016  11233276536737  10690689388347.875
16   2017  12310409370892  11737123139242.688
17   2018  13894817549374    12842646653602.5
18   2019  14342903006431  14010596597645.562
19   后(1)               -  15244498051636.938
20   后(2)               -  16548075149206.623
21   后(3)               -    17925262317445.0
22   后(4)               -  19380216151703.938
23   后(5)               -  20917327960937.562
24   后(6)               -    22541237021453.5
25   后(7)               -  24256844579072.375
26   后(8)               -   26069328641956.25
27   后(9)               -    27984159608759.0
28  后(10)               -  30007116779257.375
29  后(11)               -  32144305797307.316

在这里插入图片描述

达到预测值30000000000000

根据预测所需期数为10,值为30007116779257.375

预测表达式为：
(1) 级比检验

              原始数据 级比值区间[0.904837, 1.105171]
0   10581821399000                         _
1   10936419054000                    0.9676
2   11458243878000                    0.9545
3   12213729147000                    0.9381
4   13036640230000                    0.9369
5   13814611414000                    0.9437
6   14451858656000                    0.9559
7   14712844084000                    0.9823
8   14448933025000                    1.0183
9   14992052727000                    0.9638
10  15542581104000                    0.9646
11  16197007349000                    0.9596
12  16784849196000                    0.9650
13  17527163695000                    0.9576
14  18224704440000                    0.9617
15  18714960538000                    0.9738
16  19519353692000                    0.9588
17  20580159776000                    0.9485
18  21433226000000                    0.9602
(2) 累加生成序列

            累加数据
0   1.058182e+13
1   2.151824e+13
2   3.297648e+13
3   4.519021e+13
4   5.822685e+13
5   7.204147e+13
6   8.649332e+13
7   1.012062e+14
8   1.156551e+14
9   1.306472e+14
10  1.461897e+14
11  1.623867e+14
12  1.791716e+14
13  1.966988e+14
14  2.149235e+14
15  2.336384e+14
16  2.531578e+14
17  2.737379e+14
18  2.951712e+14
(3) 均值序列

                 均值数据
0    16050030926000.0
1    27247362392000.0
2    39083348904500.0
3    51708533593000.0
4    65134159415000.0
5    79267394450000.0
6    93849745820000.0
7   108430634374500.0
8   123151127250500.0
9   138418444166000.0
10  154288238392500.0
11  170779166665000.0
12  187935173110500.0
13  205811107178000.0
14  224280939667000.0
15  243398096782000.0
16  263447853516000.0
17  284454546404000.0
(4) 构建矩阵

(4.1) 矩阵X

[['-16050030926000.0' '1.0']
 ['-27247362392000.0' '1.0']
 ['-39083348904500.0' '1.0']
 ['-51708533593000.0' '1.0']
 ['-65134159415000.0' '1.0']
 ['-79267394450000.0' '1.0']
 ['-93849745820000.0' '1.0']
 ['-108430634374500.0' '1.0']
 ['-123151127250500.0' '1.0']
 ['-138418444166000.0' '1.0']
 ['-154288238392500.0' '1.0']
 ['-170779166665000.0' '1.0']
 ['-187935173110500.0' '1.0']
 ['-205811107178000.0' '1.0']
 ['-224280939667000.0' '1.0']
 ['-243398096782000.0' '1.0']
 ['-263447853516000.0' '1.0']
 ['-284454546404000.0' '1.0']]
(4.1) 矩阵Y

[['10936419054000']
 ['11458243878000']
 ['12213729147000']
 ['13036640230000']
 ['13814611414000']
 ['14451858656000']
 ['14712844084000']
 ['14448933025000']
 ['14992052727000']
 ['15542581104000']
 ['16197007349000']
 ['16784849196000']
 ['17527163695000']
 ['18224704440000']
 ['18714960538000']
 ['19519353692000']
 ['20580159776000']
 ['21433226000000']]
(5)计算参数a和B

a为 -0.036135494166086246
B为 10838403457276.238
(6) 预测效果分析

              原始数据          预测数据            残差  相对误差(%)        级比偏差
0   10581821399000  1.058182e+13  0.000000e+00     0.00           -
1   10936419054000  1.142598e+13 -4.895622e+11    -4.48  -0.0401511
2   11458243878000  1.184642e+13 -3.881714e+11    -3.39  -0.0260688
3   12213729147000  1.228232e+13 -6.859059e+10    -0.56 -0.00843915
4   13036640230000  1.273426e+13  3.023764e+11     2.32 -0.00714918
5   13814611414000  1.320284e+13  6.117735e+11     4.43   -0.014459
6   14451858656000  1.368865e+13  7.632050e+11     5.28  -0.0275738
7   14712844084000  1.419235e+13  5.204984e+11     3.54  -0.0559533
8   14448933025000  1.471457e+13 -2.656387e+11    -1.84  -0.0946526
9   14992052727000  1.525601e+13 -2.639611e+11    -1.76  -0.0360662
10  15542581104000  1.581738e+13 -2.747978e+11    -1.77  -0.0369261
11  16197007349000  1.639940e+13 -2.023929e+11    -1.25  -0.0315512
12  16784849196000  1.700284e+13 -2.179885e+11    -1.30  -0.0373561
13  17527163695000  1.762848e+13 -1.013158e+11    -0.58  -0.0294013
14  18224704440000  1.827714e+13 -5.243822e+10    -0.29  -0.0338087
15  18714960538000  1.894967e+13 -2.347137e+11    -1.25   -0.046816
16  19519353692000  1.964695e+13 -1.275987e+11    -0.65  -0.0306912
17  20580159776000  2.036989e+13  2.102719e+11     1.02   -0.019619
18  21433226000000  2.111942e+13  3.138013e+11     1.46  -0.0321962
(6.1) 平均相对误差(%):1.9563157894736842
(6.2) 后验差比:0.012118255490748613
(6.3) 小概率误差:1.0
(6.4) 模型结果评判【好，合格，勉强，不合格】:好
(7) 模型预测以及图示

预测期数11

        年       美国GDP(美元)                预测数据
0    2001  10581821399000    10581821399000.0
1    2002  10936419054000    11425981256910.0
2    2003  11458243878000  11846415283294.125
3    2004  12213729147000  12282319733317.688
4    2005  13036640230000    12734263861590.5
5    2006  13814611414000  13202837869195.125
6    2007  14451858656000  13688653674440.312
7    2008  14712844084000  14192345711974.625
8    2009  14448933025000  14714571761305.252
9    2010  14992052727000  15256013805801.936
10   2011  15542581104000  15817378923311.252
11   2012  16197007349000   16399400209539.81
12   2013  16784849196000  17002837735416.314
13   2014  17527163695000  17628479539680.125
14   2015  18224704440000  18277142657993.312
15   2016  18714960538000  18949674189921.938
16   2017  19519353692000   19646952405175.25
17   2018  20580159776000   20369887890553.25
18   2019  21433226000000    21119424739096.5
19   后(1)               -  21896541782991.625
20   后(2)               -  22702253871845.875
21   后(3)               -  23537613197992.125
22   后(4)               -    24403710670567.5
23   后(5)               -  25301677340146.375
24   后(6)               -    26232685875799.0
25   后(7)               -  27197952096497.625
26   后(8)               -    28198736558875.0
27   后(9)               -  29236346203405.375
28  后(10)               -  30312136061158.125
29  后(11)               -  31427511023355.812

在这里插入图片描述

达到预测值30000000000000

根据预测所需期数为10,值为30312136061158.125

预测表达式为：

$\hat{x}^{1}(t)=(x^{0}(t)-\frac{10838403457276.238}{-0.036135494166086246})e^{0.036135494166086246(t-1)} + \frac{10838403457276.238}{-0.036135494166086246},t=1,2,…k,…$

6、绘制中美GDP预测对比

# 01-19年的真实数据
ChinaReal = GDP[['中国GDP(美元)']].values
USAReal = GDP[['美国GDP(美元)']].values

# 计算超越时刻
cc = 0
mm= 0
for kin, kva in enumerate(ChinaPredict):
    if kva > USAPredict[kin]:
        cc = kin
        mm = kva
        break
        
plt.figure(figsize=(10, 6.5))
plt.scatter(range(len(ChinaReal)), ChinaReal,label='中国实际', s=45, color='k')
plt.plot(ChinaPredict,label='中国预测-GM(1,1)',marker='D',color='tab:red',ms=5,linestyle='--',lw=2)
plt.scatter(range(len(USAReal)), USAReal,label='美国实际', s=45, color='orange')
plt.plot(USAPredict,label='美国预测-GM(1,1)',marker='s',color='tab:green',ms=5,linestyle='-.',lw=2)
plt.legend(fontsize=15)
plt.title('中美GDP(美元)对比', fontsize=17)
plt.xlabel('年', fontsize=15)
plt.ylabel('GDP(美元)', fontsize=15)
ctick = [int(2000+h) for h in list(range(1,len(ChinaPredict)+1))]
if cc:
    plt.text(cc, mm, '赶超年份：%s' % ctick[cc],fontsize=16,color='red',horizontalalignment='center')
plt.xticks(list(range(len(ChinaPredict))), ctick, rotation=-40)
plt.show()