python实现GM(1,1)模型

一个数值实验，实现了《统计预测和决策》第十章灰色预测法第二节的的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()

数值实验结果为

另外还有一个第一节的灰色预测理论的小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()