备战数学建模五之灰色预测

 

 

# -*- coding: utf-8 -*-
"""
Spyder Editor
This is a temporary script file.
"""
import numpy as np
import math
from matplotlib.pyplot import plot,show,rc,legend,subplot
from scipy.optimize import curve_fit



history_data = [724.57, 746.62, 778.27, 800.8, 827.75, 871.1, 912.37, 954.28, 995.01, 1037.2]
n = len(history_data)
X0 = np.array(history_data)

n=len(X0); jibi=X0[:-1]/X0[1:]  #求级比
bd1=[jibi.min(),jibi.max()]    #求级比范围
bd2=[np.exp(-2/(n+1)),np.exp(2/(n+1))]   #q求级比的容许范围
print("级比",bd1)
print("级比容许范围",bd2)
if bd1[0]<=bd2[0] :
     print("级比下限不合适，请修改级比")
if bd1[1]>=bd2[1]:
     print("级比上限不合适，请修改级比")

# 累加生成
history_data_agg = [sum(history_data[0:i + 1]) for i in range(n)]
X1 = np.array(history_data_agg)

# 计算数据矩阵B和数据向量Y
B = np.zeros([n - 1, 2])
Y = np.zeros([n - 1, 1])
for i in range(0, n - 1):
    B[i][0] = -0.5 * (X1[i] + X1[i + 1])
    B[i][1] = 1
    Y[i][0] = X0[i + 1]

# 计算GM(1,1)微分方程的参数发展系数a和灰作用量b
# A = np.zeros([2,1])
A = np.linalg.inv(B.T.dot(B)).dot(B.T).dot(Y)#inv是求逆
a = A[0][0]
b = A[1][0]
print("最小二乘法求出的a和b：",A)
# 建立灰色预测模型
XX0 = np.zeros(n)
XX0[0] = X0[0]
for i in range(1, n):
    XX0[i] = (X0[0] - b / a) * (1 - math.exp(a)) * math.exp(-a * (i));
print("灰色预测求出的历史值",XX0)
# 模型精度的后验差检验
e = 0  # 求残差平均值
for i in range(0, n):
    e += (X0[i] - XX0[i])
e /= n
print("残差平均值：",e)
# 求历史数据平均值
aver = 0;
for i in range(0, n):
    aver += X0[i]
aver /= n
print("历史数据平均值",aver)
# 求历史数据方差
s12 = 0;
for i in range(0, n):
    s12 += (X0[i] - aver) ** 2;
s12 /= n
print("历史数据方差",s12)
# 求残差方差
s22 = 0;
for i in range(0, n):
    s22 += ((X0[i] - XX0[i]) - e) ** 2;
s22 /= n
print("残差方差",s22)
# 求后验差比值
C = s22 / s12
print("后验差比值：",C)
# 求小误差概率
cout = 0
for i in range(0, n):
    if abs((X0[i] - XX0[i]) - e) < 0.6754 * math.sqrt(s12):
        cout = cout + 1
    else:
        cout = cout
P = cout / n
print("小误差概率：",P)
if (C < 0.35 and P > 0.95):
    # 预测精度为一级
    m = 10  # 请输入需要预测的年数
    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))
    print(f)
else:
    print('灰色预测法不适用')

 一定要注意级比检验和残差检验

运行结果为

级比 [0.9502353346343703, 0.9718656343656344]
级比容许范围 [0.8337529180751806, 1.1993961020353858]
最小二乘法求出的a和b： [[-4.19774955e-02]
 [ 6.93940282e+02]]
灰色预测求出的历史值 [ 724.57        739.77422313  771.48909067  804.5636066   839.0560604
  875.02724051  912.54054142  951.66207539  992.46078901 1035.00858465]
残差平均值： 0.1817788205778811
历史数据平均值 864.7970000000001
历史数据方差 10357.277801
残差方差 26.811330299279952
后验差比值： 0.002588646439191899
小误差概率： 1.0
往后m各年负荷为：
[1079.38044724 1125.65457636 1173.91252409 1224.23933874 1276.72371473
 1331.45814889 1388.53910348 1448.0671762  1510.14727748 1574.88881535]

进程已结束,退出代码0