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