前言

用 python 复刻上一篇博客的 Matlab 代码。

【学习笔记】灰色预测 GM(1,1) 模型 —— Matlab

python代码

# %%
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl

mpl.rcParams['font.sans-serif'] = ['SimHei']   #设置字体
mpl.rcParams['axes.unicode_minus'] = False     # - 号设置

year =np.array([1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004]).T  # 横坐标表示年份
x0 = np.array([174,179,183,189,207,234,220.5,256,270,285]).T # 原始数据序列

# 创建第一个图形
plt.figure(1)
n = x0.shape[0]
x1 = np.cumsum(x0)
plt.plot(year,x0,'o-')
plt.plot(year,x1,'r-')
plt.legend('x(0)','x(1)')
plt.grid(True)
plt.xlabel('年份')  
plt.ylabel('排污总量')

# %%
# 级比检验
rho = np.zeros((n,))
# 计算 rho
for i in range(1, n):
    rho[i] = x0[i] / x1[i-1]

# 创建图表
plt.figure(2)
plt.plot(year[1:], rho[1:], 'o-', label='rho')
plt.plot([year[1], year[-1]], [0.5, 0.5], '-', label='临界线')
plt.grid(True)

# 在指定坐标添加文本
plt.text(year[-2] + 0.2, 0.55, '临界线')

# 设置x轴刻度
plt.xticks(year[1:])

# 添加标签
plt.xlabel('年份')
plt.ylabel('原始数据的光滑度')

# 显示图表
plt.legend()
plt.show()

# %%
# 指标1:光滑比小于0.5的数据占比
num1 = np.sum(rho<0.5)/(n-1)
# 指标2:除去前两个时期外，光滑比小于0.5的数据占比
num2 = np.sum(rho[2:]<0.5)/(n-3)

print("指标一:",num1*100,"%")
print("指标二:",num2*100,"%")

# %%
def gm11(x0, predict_num):
    n = x0.shape[0]
    x1 = np.cumsum(x0)
    z1 = 0.5 * x1[1:n] + 0.5 * x1[0:n-1]
    y = x0[1:]
    x = z1
    # 最小二乘法求解
    k = ((n-1)*np.sum(x*y)-np.sum(x)*np.sum(y))/((n-1)*np.sum(x*x)-np.sum(x)*np.sum(x))
    b = (np.sum(x*x)*np.sum(y)-np.sum(x)*np.sum(x*y))/((n-1)*np.sum(x*x)-np.sum(x)*np.sum(x))
    a = -k
    u = b
    x0_hat = np.zeros((n,))
    x0_hat[0] = x0[0]
    for m in range(n-1):
        x0_hat[m+1] = (1-np.exp(a))*(x0[0]-u/a)*np.exp(-a*(m+1))
    result = np.zeros((predict_num,))
    for i in range(predict_num):
        result[i] = (1-np.exp(a))*(x0[0]-u/a)*np.exp(-a*(n+i))

    # 计算绝对残差和相对残差
    absolute_residuals = x0[1:] - x0_hat[1:]   # 从第二项开始计算绝对残差，因为第一项是相同的
    relative_residuals = np.abs(absolute_residuals) / x0[1:]  # 计算相对残差
    # 计算级比和级比偏差
    class_ratio = x0[1:] / x0[0:n-1]
    eta = np.abs(1-(1-0.5*a)/(1+0.5*a)*(1./class_ratio)) # 计算级比偏差
    return result, x0_hat, relative_residuals, eta


# %%
if num1 > 0.6 and num2 > 0.9:
    if n > 7:    # 将数据分为训练组和试验组(根据原数据量大小n来取，n小于7则取最后两年为试验组，n大于7则取最后三年为试验组)
        test_num = 3
    else:
        test_num = 2
    train_x0 = x0[0:n-test_num]   # 训练数据
    print('训练数据是: ',train_x0)

    test_x0 =  x0[n-test_num:]  # 试验数据
    print('试验数据是: ',test_x0)
    
    # 使用GM(1,1)模型对训练数据进行训练，返回的result就是往后预测test_num期的数据
    print('GM(1,1)模型预测')
    result1,_,_,_ = gm11(train_x0, test_num) # 使用传统的GM(1,1)模型对训练数据，并预测后test_num期的结果
    # 绘制对试验数据进行预测的图形
    test_year = year[n-test_num:]  # 试验组对应的年份
    plt.figure(3)
    plt.plot(test_year,test_x0,'o-',label='试验组的真实数据')
    plt.plot(test_year,result1,'*-',label='预测值')
    plt.grid(True) 
    # 设置x轴刻度
    plt.xticks(year[n-test_num:])
    plt.legend() 
    plt.xlabel('年份')
    plt.ylabel('排污总量')
 
    predict_num = int(input('请输入你要往后面预测的期数：'))
    # 计算使用传统GM模型的结果
    result, x0_hat, relative_residuals, eta = gm11(x0, predict_num)
    
    ## 绘制相对残差和级比偏差的图形
    plt.figure(4)
    # 创建一个2行1列的子图布局
    plt.subplot(2, 1, 1)  # 第1个子图
    plt.plot(year[1:], relative_residuals,'*-',label='相对残差')
    plt.xticks(year[1:])
    plt.legend()
    plt.grid(True) 

    plt.subplot(2, 1, 2)  # 第2个子图
    plt.plot(year[1:], eta,'*-',label='级比偏差')
    plt.xticks(year[1:])
    plt.legend()
    plt.grid(True) 
    plt.xlabel('年份')
    
    ## 残差检验
    average_relative_residuals = np.mean(relative_residuals)  # 计算平均相对残差 mean函数用来均值
    print('平均相对残差为',average_relative_residuals)
    
    ## 级比偏差检验
    average_eta = np.mean(eta)  # 计算平均级比偏差
    print('平均级比偏差为',average_eta)
    
    ## 绘制最终的预测效果图
    plt.figure(5)
    plt.plot(year, x0, '-o', label='原始数据')
    plt.plot(year, x0_hat, '-*m', label='拟合数据')

    year_predict = np.arange(year[n-1], year[n-1] + predict_num + 1)
    res = np.append(x0[n-1],result)
    plt.plot(year_predict, res, '-*k', label='预测数据' )
    plt.grid(True) 

    plt.legend()  
    plt.xticks(year[1:] + predict_num)
    plt.xlabel('年份')
    plt.ylabel('排污总量')

运行结果：