【python】灰色预测 GM(1,1) 模型

文章目录

前言
python代码

前言

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

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

python_9">python代码

python"># %%
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
from matplotlib.pylab import mplmpl.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 = -ku = bx0_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 = 3else:test_num = 2train_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('排污总量')