灰色预测

#累加生成新序列
x0=[455,425,489,500,549,663,777]#只需将此处更换为自己的待预测序列即可
x1=[0 for i in range(len(x0))]#用于存储累加后的序列,初始化为全0

#定义累加函数
def cumsum(x0):
    n=len(x0)
    x1[0]=x0[0]
    for i in range(1,n):
        x1[i]=x0[i]+x1[i-1]
    return x1

x1=cumsum(x0)#执行累加函数,此时x1中存储了累加后的有规律序列，用于接下来建模

alpha=0.5#这里的调参很重要！！
z=[0 for i in range(len(x0)-1)]

#定义"计算x1的加权均值"的函数
def weigh_mean(x1):
    for j in range(1,len(x1)):
        #z[j-1]=alpha*x1[j]+(1-alpha)*x1[j-1]
        z[j-1]=-alpha*(x1[j]+x1[j-1])
    return z

z=weigh_mean(x1)


import numpy as np
import math

#求解灰参数a与b
def cal_params(z,x0):
    B=(np.array([[z[i] for i in range(len(z))],[1 for i in range(len(x0)-1)]])).T
    Yn=np.array([x0[j] for j in range(1,len(x0))])
    grey_params=np.array([0,0])
    grey_params=((np.linalg.inv((B.T).dot(B))).dot(B.T)).dot(Yn)
    return grey_params

#执行cal_params（），得到灰参数a与b
grey_params=cal_params(z,x0)
a=grey_params[0]
b=grey_params[1]



#建立灰色预测模型
n=len(x0)
X0 = np.zeros(n)
X0[0] = x0[0]
for i in range(1,n):
    X0[i] = (x0[0] - b/a)*(1-math.exp(a))*math.exp(-a*(i));




#模型精度的后验差检验
e = 0      #求残差平均值
for i in range(0,n):
    e += (x0[i] - X0[i])
e /= n

#求历史数据平均值
aver = 0;     
for i in range(0,n):
    aver += x0[i]
aver /= n

#求历史数据方差
s12 = 0;     
for i in range(0,n):
    s12 += (x0[i]-aver)**2;
s12 /= n

#求残差方差
s22 = 0;       
for i in range(0,n):
    s22 += ((x0[i] - X0[i]) - e)**2;
s22 /= n

#求后验差比值
C = s22 / s12   

#求小误差概率
cout = 0
for i in range(0,n):
    if abs((x0[i] - X0[i]) - e) < 0.6754*math.sqrt(s12):
        cout = cout+1
    else:
        cout = cout
P = cout / n

if (C < 0.35 and P > 0.95):
    #预测精度为一级
    m = int(input("要预测多少期的呀？"))   #请输入需要预测的年数
    #print('往后m各年负荷为：')
    f = np.zeros(m)
    for i in range(0,m):
        f[i] = (x0[0] - b/a)*(1-math.exp(a))*math.exp(-a*(i+n))#加n表示预测之后的年份
    print("各期预测值为：\n",f)
else:
    print('灰色预测法不适用，请使用其他模型')


import matplotlib.pyplot as plt
from matplotlib import font_manager
plt.rcParams['font.sans-serif']=['SimHei']
#my_font=my_font = font_manager.FontProperties(fname="C:\Windows\Fonts\STZHONGS.ttf")
#检验过去
plt.figure(figsize=(10,6))
plt.title('预测录取人数与实际录取人数比较')
plt.xlabel("年份")
plt.ylabel("录取人数（单位：人）")
plt.scatter([i for i in range(2013,2013+len(x0))],x0,color='red',label="真实值")
plt.scatter([i for i in range(2013,2013+len(X0))],X0,label="预测值")
plt.legend(loc="upper left")
plt.show()
#展望未来
plt.figure(figsize=(10,6))
plt.scatter([int(i) for i in range(2020,2020+m)],f,color='g')
plt.xlabel("年份")
plt.ylabel("录取人数（单位：人）")
plt.title("未来7年录取人数预测")
plt.show()