文章目录
- 前言
- python代码
前言
用 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('排污总量')
运行结果：
