python实现GM(1,1)模型

以你之姓@ 2022-10-10 13:59 169阅读 0赞

一个数值实验，实现了《统计预测和决策》第十章灰色预测法第二节的的GM(1,1)模型，  
题目是给定一个时间序列，试建立GM(1,1)模型预测模型，并预测第八期的预测值。关键步骤有代码注释，下面是代码

import numpy as np
    import pandas as pd
    import math
    np.set_printoptions(suppress=True)
    x=pd.Series([26.7,31.5,32.8,34.1,35.8,37.5])
    def multisum(df):#第一步，计算累加数列
        sum=0
        df1=[]
        for i in df:
            sum =i+sum
            df1.append(sum)
        return df1
    def getB(df):#第二三步，构造b，y，计算btb，btb1，bty
        # 这里卡了好久，一直是数组越界，最后发现是自己没有思考清楚。
        b=[]
        df=pd.Series(df)
        for i in range(5):
            t=-(df[i]+df[i+1])/2
            b.append(t)
        b=pd.Series(b)
        b1=pd.Series([1,1,1,1,1])
        b=pd.concat([b,b1],axis=1)
        # print('b is\n',b)
        y=x.drop([0])
        # print('y is \n',y)
        btb=b.T.dot(b)
        # print('btb is \n',btb)
        btb1=np.linalg.inv(btb)
        # print('btb1 is \n',btb1)
        bty=b.T.dot(y.tolist())
        # print('bty is \n',bty)
        #这里不加tolist会报错矩阵未对齐，所以尝试将数列转化为list有效。
        a=btb1.dot(bty)[0]
        n=btb1.dot(bty)[1]
        # print(a,n)
        n_a=n/a
        # print(n_a)
        return df[0]-n_a,abs(a),n_a
    def getpredict(n):#得到预测模型
        x=[]
        a=n[0]
        b=n[1]
        c=n[2]
        for k in range(6):
            xk=a*(math.exp(b*k))+c
            x.append(xk)
        print( '第八次预测为',(a * (math.exp(b * 7)) + c)-(a * (math.exp(b * 6)) + c))#这里的计算系数可以手动更改，用来预测其他结果
        return x
    def multisub(df):#计算累减数列
        df.reverse()
        b=[]
        for i in range(5):
            t=df[i]-df[i+1]#类似累加数列的思想，要点是将数列进行翻转，就可以很简单的获得结果
            b.append(t)
        b.reverse()
        b.insert(0,df[-1])
        return b
    def mysub(df1,df2):#为了计算后验差检验的前置函数
        return abs(df1-df2)#获得绝对误差序列
    def findn(df):
        min=df.min()
        max=df.max()
        return (min+max*0.5)/(df+0.5*max)#计算关联系数
    def main():
        x_add=multisum(x)
        print('累加数列的结果是 \n',x_add)
        n= getB(x_add)
        x_p=getpredict(n)
        print('预测数列的结果是 \n',x_p)
        x_sub=multisub(x_p)
        print('累减数列的结果是\n',mysub(x,x_sub).tolist())
        print('关联系数为 \n',findn(mysub(x,x_sub)).tolist())
        print('关联度为 \t',sum(findn(mysub(x,x_sub)))/6)
    if __name__ == '__main__':
        main()

数值实验结果为  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80NDAxNDgzNg_size_16_color_FFFFFF_t_70]

另外还有一个第一节的灰色预测理论的小demo，实现了两数列关联系数的计算，为GM(1,1)模型提供了部分计算支持。

import numpy as np
    import pandas as pd
    x1=pd.Series([8,8.8,16,18,24,32])
    x2=pd.Series([10,12.12,19.28,20.25,23.4,30.69])
    x3=pd.Series([6,6.35,6.57,6.98,8.35,8.75])
    def nromal(df):
        return df/df.min()
    def sub(df1,df2):
        return df1-df2
    def findn(df1,df2):
        df= df1.append(df2).abs()
        min=df.min()
        max=df.max()
        return (min+max*0.5)/(df1+0.5*max),(min+max*0.5)/(df2+0.5*max)
    def main():
        x_1=nromal(x1)
        x_2=nromal(x2)
        x_3=nromal(x3)
        x12=sub(x_1,x_2).abs()
        x13=sub(x_1,x_3).abs()
        print(x12,x13)
        print(sum(findn(x12,x13)[0])/6,
              sum(findn(x12,x13)[1])/6)
    if __name__ == '__main__':
        main()

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80NDAxNDgzNg_size_16_color_FFFFFF_t_70 1]

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80NDAxNDgzNg_size_16_color_FFFFFF_t_70]: /images/20221005/02f574ab914e40799c0cfebc40ec981e.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80NDAxNDgzNg_size_16_color_FFFFFF_t_70 1]: /images/20221005/07f6a37bc5e242b79ce07863e5541e61.png