数学建模1：灰色预测模型

灰色预测模型：

灰色预测模型是一种用于处理不完全信息的预测方法，适用于数据样本少、系统不确定性高的情况。

灰色预测模型的核心思想是利用已有的部分信息对未来进行预测。它的基本形式是灰色模型GM(1,1)，其中“1”表示一个方程、一个变量的模型。这个模型通过对原始数据序列进行累加生成序列，然后建立一个差分方程进行建模，最后利用这个模型来预测未来的数据。本文按照最基础的单序列，一阶模型

解题：

假设我们有一个时间序列数据（比如月销售量）如下：

月销量	30	35	40	45	50
第几月	1	2	3	4	5

求解：使用灰色模型预测第6个月的销量

1:计算累加生成序列 $x_{1}(t)$

$x_{1}(1)=x(1)=30$

$x_{1}(2)=x(1)+x(2)=30+35=65$

$x_{1}(3)=x(1)+x(2)+x(3)=30+35+40=105$

$x_{1}(4)=x(1)+x(2)+x(3)+x(4)=30+35+40+45=150$

$x_{1}(5)=x(1)+x(2)+x(3)+x(4)+x(5)=30+35+40+45+50=200$

解得累加生成序列：

$x_{1}(t)=[30,65,105,150,200]$

2.建立差分方程： $z(t)$

$z(2) = \frac{x_{1}(1) + x_{1}(2)}{2} = \frac{30 + 65}{2} = 47.5$

$z(3) = \frac{x_{1}(2) + x_{1}(3)}{2} = \frac{65 + 105}{2} = 85$

$z(4) = \frac{x_{1}(3) + x_{1}(4)}{2} = \frac{105 + 150}{2} = 127.5$

$z(5) = \frac{x_{1}(4) + x_{1}(5)}{2} = \frac{150 + 200}{2} = 175$

解均值序列为：

$z(t) = [47.5,85,127.5,175]$

3.构造矩阵：

$B = \begin{pmatrix} -47.5\\ -85\\ -127.5\\ -175\\ \end{pmatrix}$

目标矩阵：

$Y = \begin{pmatrix} 35\\ 40\\45\\50\\ \end{pmatrix}$

4.使用最小二乘法求解模型参数

求解目标方程为：

$x_{0}(t) = u - az(t)$

化为矩阵形式：

$Y = BU$ 也就是如下形式

$\begin{pmatrix} 35\\ 40\\45\\50\\ \end{pmatrix} = \begin{pmatrix} -47.5& 1\\-85 &1 \\ -127.5&1 \\ -175&1 \end{pmatrix}\begin{pmatrix} a\\ b \end{pmatrix}$

最小二乘法就是求 $(Y-BU)^{T}(Y-BU)$ 取最小值时的U

求解U的估计值为： $\hat{U} = [\hat{a},\hat{b}]^{T} = (B^{T}B)^{-1}B^{T}Y$

求出a和b的值带入到微分方程中：

$x_{1}(t + 1) = (x_{0}(1) - \frac{b}{a})e^{-at} + \frac{b}{a}$

计算得：

a = -0.1173

b = 29.7412

带入得

$x_{1}(6) = 256.2$

$x_{0}(6) = 56.2$

python代码：

import numpy as np
import matplotlib.pyplot as pltclass GM11:def __init__(self, X, T):self.X = Xself.T = Tself.X1= []self.Z = []self.U = Nonedef caculateXSum(self):#计算X1矩阵和Z矩阵count = 0X1 = []Z  = []for i in range(0, len(self.X),1):count = count + self.X[i]X1.append(count)for j in range(0, len(X1) - 1,1):Z.append((X1[j] + X1[j + 1])/2)self.X1 = X1self.Z = Zdef caculateAB(self):#计算U矩阵的值Y = (np.array(self.X[1:])).reshape(-1,1)B = np.array(self.Z)*-1ones_column = np.ones((B.shape[0], 1))B = np.hstack((B.reshape(-1,1), ones_column))#print(B)#print(Y)self.U = np.linalg.inv(B.T @ B) @ B.T @ Y#self.U = np.dot(np.linalg.inv(np.dot(B.T,B)),np.dot((B.T),Y))def caculateX(self, t):#计算最终的Xif self.U is None:self.caculateXSum()self.caculateAB()t = t - 1a = self.U[0]b = self.U[1]result = (self.X[0] - b/a)*(np.e ** (-1*a*(t))) + b/areturn  result[0]X = np.array([30, 35, 40, 45, 50])
T = np.array([1, 2, 3, 4, 5])
GM11 = GM11(X, T)
n = 6#预测第n个值
result = GM11.caculateX(6)
print(f'解得u矩阵为：{GM11.U[0]},{GM11.U[1]}')
print(f'预测得x({n})的值为：{result}')
#print(X)#绘图
X1 =  GM11.X1
for i in range(1,4,1):X1.append(GM11.caculateX(5 + i))
#print(X1)
X = []
X.append(X1[0])
for i in range(1, len(X1), 1):X.append(X1[i] - X1[i - 1])
plt.plot(range(len(X) - 3), X[0:-3], 'o-', label='Original')
plt.plot(range(len(X) - 4, len(X)), X[-4:], 'x--', label='Future Prediction')
plt.legend()
plt.show()