蔬菜数量预测¶
由1995~2001年中国蔬菜产量，具体数据见下表，建立模型预测2002年中国蔬菜产量，并对预测结果做检验。
| 年代 | 1995 | 1996 | 1997 | 1998 | 1999 | 2000 | 2001 |
|---|---|---|---|---|---|---|---|
| 产量 | 25723 | 30379 | 34473 | 38485 | 40514 | 42400 | 48337 |
数据的检验¶
记 1995 年到 2001 年的蔬菜产量分别为 $x^{(0)}(1), x^{(0)}(2), \cdots, x^{(0)}(7)$, 构造 参考序列 $x^{(0)}=\left(x^{(0)}(1), x^{(0)}(2), \cdots, x^{(0)}(7)\right)$, 经检验级比符合要求, 参考序列 $x^{(0)}$ 可以用来建立 GM(1,1)模型。
模型的建立与求解¶
构造累加序列 $$ \begin{gathered} x^{(1)}=\left(x^{(1)}(1), x^{(1)}(2), \cdots, x^{(1)}(7)\right) \\ =\left(x^{(0)}(1), x^{(0)}(1)+x^{(0)}(2), \cdots, x^{(0)}(1)+\cdots+x^{(0)}(7)\right) \end{gathered} $$ 和 $x^{(1)}$ 的均值生成序列 $$ z^{(1)}=\left(z^{(1)}(2), z^{(1)}(3), \cdots, z^{(1)}(7)\right), $$ 其中 $z^{(1)}(k)=0.5 x^{(1)}(k)+0.5 x^{(1)}(k-1), k=2,3, \cdots 7$ 。
将 $$ Y=\left[\begin{array}{c} x^{(0)}(2) \\ x^{(0)}(3) \\ \vdots \\ x^{(0)}(7) \end{array}\right], \quad B=\left[\begin{array}{cc} -z^{(1)}(2) & 1 \\ -z^{(1)}(3) & 1 \\ \vdots & \vdots \\ -z^{(1)}(7) & 1 \end{array}\right] $$ 代入式中, 得参数的估计值 $$ \hat{\boldsymbol{u}}=[\hat{\boldsymbol{a}}, \hat{\boldsymbol{b}}]^{T}=\left(\boldsymbol{B}^{T} B\right)^{-1} B^{T} Y=\left[\begin{array}{c} -0.0843 \\ 27858.4508 \end{array}\right] $$
得 GM(1,1)模型为
- (1)灰微分方程: $x^{(0)}(k)-0.0843 z^{(1)}(k)=27858.4508$.
- (2)白化方程: $\frac{d x^{(1)}(t)}{d t}-0.0843 x^{(1)}(t)=27858.4508$.
- (3)白化方程的时间响应式: $$ \hat{\boldsymbol{x}}^{(1)}(k+1)=356328.9910 e^{0.0842 k}-330605.9910 . $$
模型检验与预测¶
利用GM(1,1)式求得累加生成序列的预测值 $\hat{x}^{(1)}=\left(\hat{x}^{(1)}(1), \hat{x}^{(1)}(2), \cdots, \hat{x}^{(1)}(8)\right)$, 由 $\hat{x}^{(0)}(1)=\hat{x}^{(1)}(1)$, $\hat{x}^{(0)}(k+1)=\hat{x}^{(1)}(k+1)-\hat{x}^{(1)}(k) \quad, \quad k=1,2, \cdots, 7 \quad, \quad$ 得 $\hat{x}^{(0)}=\left(\hat{x}^{(0)}(1), \hat{x}^{(0)}(2), \cdots, \hat{x}^{(0)}(8)\right)$, 即得到 2005 年到 2002 年的预测值见表 15.2。画出预测值与实际值的变化曲线如图所示。
$$
\begin{array}{ccccc}
\hline \text { 年份 } & \text { 实际数据 } x^{(0)} & \text { 预测数据 } \hat{x}^{(0)} & \text { 残差 } x^{(0)}-\hat{\boldsymbol{x}}^{(0)} & \text { 相对误差 } \boldsymbol{\delta}(\boldsymbol{k}) / \% \\
\hline 1995 & \mathbf{2 5 7 2 3} & \mathbf{2 5 7 2 3} & 0 & 0 \\
1996 & \mathbf{3 0 3 7 9} & \mathbf{3 1 3 2 7 . 3 5 6 7} & \mathbf{- 9 4 8 . 3 5 6 7} & 3.12 \\
1997 & 34473 & 34081.5622 & 391.4378 & 1.14 \\
1998 & 38485 & \mathbf{3 7 0 7 7 . 9 0 9 1} & 1407.0909 & 3.66 \\
1999 & 40514 & \mathbf{4 0 3 3 7 . 6 8 5 7} & \mathbf{1 7 6 . 3 1 4 3} & 0.44 \\
2000 & 42400 & \mathbf{4 3 8 8 4 . 0 5 1 8} & -1484.0518 & 3.50 \\
2001 & 48337 & \mathbf{4 7 7 4 2 . 2 0 3 6} & \mathbf{5 9 4 . 7 9 6 4} & 1.23 \\
\mathbf{2 0 0 2} & & \mathbf{5 1 9 3 9 . 5 5 2 4} & & \\
\hline
\end{array}
$$
所有的相对误差 $\delta(k)<0.10(k=1,2, \cdots, 7)$, 所以模型精度为优, 可以用 来预测, 2002 年的预测值为 51939.5524。
代码
import numpy as np
import sympy as sp
from matplotlib.pyplot import plot,show,rc,legend,xticks
rc('font',size=16); rc('font',family='SimHei')
x0=np.array([25723,30379,34473,38485,40514,42400,48337])
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求级比的容许范围
x1=np.cumsum(x0) #求累加序列
z=(x1[:-1]+x1[1:])/2.0
B=np.vstack([-z,np.ones(n-1)]).T
u=np.linalg.pinv(B)@x0[1:] #最小二乘法拟合参数
sp.var('t'); sp.var('x',cls=sp.Function) #定义符号变量和函数
eq=x(t).diff(t)+u[0]*x(t)-u[1] #定义符号微分方程
xt=sp.dsolve(eq,ics={x(0):x0[0]}) #求解符号微分方程
xt=xt.args[1] #提取方程中的符号解
xt=sp.lambdify(t,xt,'numpy') #转换为匿名函数
t=np.arange(n+1)
xt1=xt(t) #求模型的预测值
x0_pred=np.hstack([x0[0],np.diff(xt1)]) #还原数据
x2002=x0_pred[-1] #提取2002年的预测值
cha=x0-x0_pred[:-1]; delta=np.abs(cha/x0)*100
print('1995~2002的预测值：',x0_pred)
print('\n-------------------\n','相对误差',delta)
t0=np.arange(1995,2002); plot(t0,x0,'*--')
plot(t0,x0_pred[:-1],'s-'); legend(('实际值','预测值'));
xticks(np.arange(1995,2002)); show()
