数模竞赛-预测方法

灰色预测模型

$\operatorname{GM}(1,1)$只适合具有较强指数规律的序列，$\operatorname{GM}(2,1)$、DGM和Verhulst模型适用于非单调摆动发展序列或有饱和S形序列。

$\operatorname{GM}(1,1)$预测方法

$\operatorname{GM}(1,1)$预测模型定义：

$1$阶微分方程，只含$1$个变量。

已知参考数列$\overrightarrow{x}^{(0)}=\left(x^{(0)}(1),x^{(0)}(2),\cdots,x^{(0)}(n)\right)$，$1$次累加生成序列1-AGO定义为$\overrightarrow{x}^{(1)}=\left(x^{(1)}(1),x^{(1)}(2),\cdots,x^{(1)}(n)\right)$，式中$\displaystyle x^{(1)}(k)=\sum_{i=1}^kx^{(0)}(i),1\leqslant k\leqslant n$，且$\overrightarrow{x}^{(1)}$的均值生成序列为$\overrightarrow{z}^{(1)}=\left(z^{(1)}(2),z^{(1)}(3),\cdots,z^{(1)}(n)\right)$，其中$z^{(1)}(k)=0.5x^{(1)}(k)+0.5x^{(1)}(k+1),2\leqslant k\leqslant n$。

建立灰微分方程$x^{(0)}(k)+az^{(1)}(k)=b,2\leqslant k\leqslant n$，其对应的白化微分方程为$\dfrac{\mathrm dx^{(1)}(t)}{\mathrm dt}+ax^{(1)}(t)=b$，记：

$$
\displaylines{
\overrightarrow u=[a,b]^{\mathrm T},\overrightarrow{Y}=\left[x^{(0)}(2),x^{(0)}(3),\cdots,x^{(0)}(n)\right]^{\mathrm T},B=\left[
\begin{matrix}
-z^{(1)}(2)&1\\
-z^{(1)}(3)&1\\
\vdots&\vdots\\
-z^{(1)}(n)&1
\end{matrix}
\right]
}
$$

由最小二乘法求使得$J\left(\overrightarrow{u}\right)=\left(\overrightarrow{Y}-B\overrightarrow{u}\right)^{\mathrm T}\left(\overrightarrow{Y}-B\overrightarrow{u}\right)$达到最小估计值的$\overrightarrow u$为$\widehat u=\left[\widehat a,\widehat b\right]^{\mathrm T}=\left(B^{\mathrm T}B^{-1}\right)^{-1}B^{\mathrm T}\overrightarrow Y$，求解得：

$$
{\widehat x}^{(1)}(k+1)=\left(x^{(0)}(1)-\dfrac{\widehat b}{\widehat a}\right)e^{-\widehat ak}+\dfrac{\widehat b}{\widehat a},k\geqslant 0
$$

$\operatorname{GM}(1,1)$模型预测步骤：

当计算序列级比$\forall2\leqslant k\leqslant n,\lambda(k)=\dfrac{x^{(0)}(k-1)}{x^{(0)}(k)}\in\Theta,\Theta=\left(e^{-\frac{2}{n+1}},e^{\frac{2}{n+1}}\right)$,如果不落在范围$\Theta$中需要进行平移变换，选取充分大的$c>0$进行变换$y^{(0)}(k)=x^{(0)}(k)+c,1\leqslant k\leqslant n$，再计算级比。

检验预测值：

计算相对误差$\delta(k)=\dfrac{\left|x^{(0)}(k)-\widehat{x}^{(0)}(k)\right|}{x^{(0)}(k)},1\leqslant k\leqslant n$，当$\delta(k)<0.2$达到一般要求，$\delta(k)<0.1$达到较高要求。

计算级比偏差值$\rho(k)=\left|1-\left(\dfrac{1-0.5\widehat a}{1+0.5\widehat a}\right)\lambda(k)\right|$，当$\rho(k)<0.2$达到一般要求，$\rho(k)<0.1$达到较高要求。

例题1：

给定$L_{\mathrm{eq}}$有：$71.1,72.4,72.4,72.1,71.4,72.0,71.6$，预测下一个。

代码：

import numpy,sympy
x0=numpy.array([71.1,72.4,72.4,72.1,71.4,72.0,71.6])
n=len(x0)
lamda=x0[:-1]/x0[1:]
b1=[min(lamda),max(lamda)]
b2=[numpy.exp(-2/(n+1)),numpy.exp(2/(n+1))] #求Theta

x1=numpy.cumsum(x0) #前缀和
z=(x1[:-1]+x1[1:])/2
B=numpy.vstack([-z,numpy.ones(n-1)]).T
u=numpy.linalg.pinv(B)@x0[1:]
sympy.var('t')
sympy.var('x',cls=sympy.Function)
eq=x(t).diff(t)+u[0]*x(t)-u[1]
xt0=sympy.dsolve(eq,ics={x(0):x0[0]})
xt0=xt0.args[1]
xt=sympy.lambdify(t,xt0,'numpy')
t=numpy.arange(n+1)
xh=xt(t)
x0h=numpy.hstack([x0[0],numpy.diff(xh)]) #差分：还原
x1993=x0h[-1]

cha=x0-x0h[:-1] #误差分析
delta=abs(cha/x0)*100
rho=abs(1-(1-0.5*u[0])/(1+0.5*u[0])*lamda)

print(round(x1993,4))

$\operatorname{GM}(2,1)$预测方法

设原始数列$\overrightarrow{x}^{(0)}=\left(x^{(0)}(1),x^{(0)}(2),\cdots,x^{(0)}(n)\right)$，其一次累加生成序列1-AGO为$\overrightarrow{x}^{(1)}=\left(x^{(1)}(1),x^{(1)}(2),\cdots,2^{(1)}(n)\right)$，一次累减生成序列1-IAGO为$\overrightarrow{\alpha}^{(1)}\overrightarrow{x}^{(0)}=\left(\alpha^{(1)}x^{(0)}(2),\cdots,\alpha^{(1)}x^{(0)}(n)\right)$，其中$\alpha^{(1)}x^{(0)}(k)=x^{(0)}(k)-x^{(0)}(k-1),k=2,3,\cdots,n$。设$\overrightarrow{z}^{(1)}=\left(z^{(1)}(2),z^{(1)}(3),\cdots,z^{(1)}(n)\right)$，则称$\alpha^{(1)}x^{(0)}(k)+a_1x^{(0)}(k)+a_2z^{(1)}(k)=b$为$\operatorname{GM}(2,1)$模型。称$\dfrac{\mathrm d^2x^{(1)}(t)}{\mathrm dt^2}+a_1\dfrac{\mathrm dx^{(1)}(t)}{\mathrm dt}+a_2x^{(1)}(t)=b$为$\operatorname{GM}(2,1)$模型的白化方程。

设：

$$
\displaylines{
B=\left[
\begin{matrix}
-x^{(0)}(2)&-z^{(1)}(2)&1\\
-x^{(0)}(3)&-z^{(1)}(3)&1\\
\vdots&\vdots&\vdots\\
-x^{(0)}(n)&-z^{(1)}(n)&1
\end{matrix}
\right],Y=\left[
\begin{matrix}
\alpha^{(1)}x^{(0)}(2)\\
\alpha^{(1)}x^{(0)}(3)\\
\vdots\\
\alpha^{(1)}x^{(0)}(n)
\end{matrix}
\right]
}
$$

则参数序列$\overrightarrow{u}=\left[a_1,a_2,b\right]^{\mathrm T}$最小二乘估计为$\widehat u=\left(B^{\mathrm T}B\right)^{-1}B^{\mathrm T}Y$。

例题2：

已知$\overrightarrow{x}^{(0)}=(41,49,61,78,96,104)$，建立$\operatorname{GM}(2,1)$模型。

import numpy,sympy
x0=numpy.array([41,49,61,78,96,104])
n=len(x0)
x1=numpy.cumsum(x0) #1-AGO
ax0=numpy.diff(x0) #1-IAGO
z=(x1[1:]+x1[:-1])/2
B=numpy.vstack([-x0[1:],-z,numpy.ones(n-1)]).T
u=numpy.linalg.pinv(B)@ax0
sympy.var('t')
sympy.var('x',cls=sympy.Function)
eq=x(t).diff(t,2)+u[0]*x(t).diff(t)+u[1]*x(t)-u[2]
s=sympy.dsolve(eq,ics={x(0):x1[0],x(5):x1[-1]})
xt=s.args[1] #符号表达式
x=sympy.lambdify(t,xt,'numpy') #匿名函数
xh1=x(numpy.arange(n))
xh0=numpy.hstack([x0[0],numpy.diff(xh1)])
print(numpy.round(xh0,4)) #预测值
ea=x0-xh0 #残差
er=abs(ea)/x0*100 #相对误差

$\operatorname{DGM}(2,1)$模型

一些定义同上。

称$\alpha^{(1)}x^{(0)}(k)+\alpha x^{(0)}(k)=b$为$\operatorname{DGM}(2,1)$模型，$\dfrac{\mathrm d^2x^{(1)}(t)}{\mathrm dt}+a\dfrac{\mathrm dx^{(1)}(t)}{\mathrm dt}=b$为$\operatorname{DGM}(2,1)$模型的白化方程。若$\overrightarrow{u}=[a,b]^{\mathrm T}$，更改定义：

$$
\displaylines{
B=\left[
\begin{matrix}
-x^{(0)}(2)&1\\
-x^{(0)}(3)&1\\
\vdots&\vdots\\
-x^{(0)}(n)&1
\end{matrix}
\right]
}
$$

该模型参数最小二乘估计满足$\widehat u=\left(B^{\mathrm T}B\right)^{-1}B^{\mathrm T}Y$。

例题3：

已知$\overrightarrow{x}^{(0)}=(2.874,3.278,3.39,3.679,3.77,3.8)$，建立$\operatorname{DGM}(2,1)$模型。

import numpy,sympy
x0=numpy.array([2.874,3.278,3.39,3.679,3.77,3.8])
n=len(x0)
ax0=numpy.diff(x0)
B=numpy.vstack([-x0[1:],numpy.ones(n-1)]).T
u=numpy.linalg.pinv(B)@ax0
sympy.var('t')
sympy.var('x',cls=sympy.Function)
eq=x(t).diff(t,2)+u[0]*x(t).diff(t)-u[1]
s=sympy.dsolve(eq,ics={x(0):x0[0],x(t).diff(t).subs(t,0):x0[0]})
xt=s.args[1]
x=sympy.lambdify(t,xt,'numpy')
xh1=x(numpy.arange(n))
xh0=numpy.hstack([x0[0],numpy.diff(xh1)])
print(numpy.round(xh0,4))
ea=x0-xh0
er=abs(ea)/x0*100

Verhulst模型

一些定义同$\operatorname{GM}(1,1)$模型，设$a,b$为参数，称$x^{(0)}(k)+az^{(1)}(k)=b\left[z^{(1)}(k)\right]^2$，$\dfrac{\mathrm dx^{(1)}(t)}{\mathrm dt}+ax^{(1)}(t)=b\left[x^{(1)}(t)\right]^2$为白化方程，$t$为时间。

参数序列为$\overrightarrow{u}=[a,b]^{\mathrm T}$更改定义：

$$
\displaylines{
    B=\left[
     \begin{matrix}
     -z^{(1)}(2)&\left(z^{(1)}(2)\right)^2\\
     -z^{(1)}(3)&\left(z^{(1)}(3)\right)^2\\
     \vdots&\vdots\\
     -z^{(1)}(n)&\left(z^{(1)}(n)\right)^2
     \end{matrix}
    \right],Y=\left[
     \begin{matrix}
     x^{(0)}(2)\\
     x^{(0)}(3)\\
     \vdots\\
     x^{(0)}(n)
     \end{matrix}
    \right]
}
$$

最小二乘估计同上，白化方程的解为$x^{(1)}(t)=\dfrac{\widehat ax^{(0)}(1)}{\widehat bx^{(0)}(1)+\left(\widehat a-\widehat bx^{(0)}(1)\right)e^{\widehat at}}$，相应的时间相应序列为$x^{(1)}(k+1)=\dfrac{\widehat ax^{(0)}(1)}{\widehat bx^{(0)}(1)+\left(\widehat a-\widehat bx^{(0)}(1)\right)e^{\widehat ak}}$，同样累减还原式为${\widehat x}^{(0)}(k+1)={\widehat x}^{(1)}(k+1)-{\widehat x}^{(1)}(k)$。

例题4：

已知$\overrightarrow{x}^{(0)}=(4.93,2.33,3.87,4.35,6.63,7.15,5.37,6.39,7.81,8.35)$，建立Verhulst模型。

import numpy
x0=numpy.array([4.93,2.33,3.87,4.35,6.63,7.15,5.37,6.39,7.81,8.35])
n=len(x0)
x1=numpy.cumsum(x0)
z=(x1[1:]+x1[:-1])/2
B=numpy.vstack([-z,z**2]).T
u=numpy.linalg.pinv(B)@x0[1:]
x=lambda t:u[0]*x0[0]/(u[1]*x0[0]+(u[0]-u[1]*x0[0])*numpy.exp(u[0]*t)) #别求解了直接上解吧
xh1=x(numpy.arange(n))
xh0=numpy.hstack([x0[0],numpy.diff(xh1)])
print(numpy.round(xh0,4))
ea=x0-xh0
er=abs(ea)/x0*100

神经元网络

单层感知器模型

设连接权为$\overrightarrow X=\left[x_1,x_2,\cdots,x_m\right]^{\mathrm T},\overrightarrow W=\left[w_1,w_2,\cdots,2_m\right]^{\mathrm T}$，于是网络输入为$\displaystyle u=\sum_{i=1}^mw_ix_i={\overrightarrow W}^{\mathrm T}\overrightarrow X$。

激活/激励/活化函数有四种：线性函数$f(u)=ku+c$；非线性斜面函数；阈值/阶跃函数；sigmoid函数$f(u)=\dfrac{1}{1+e^{-u}}$，将$(-\infty,+\infty)$映射到$(0,1)$；tanh函数$f(u)=\dfrac{e^u-e^{-u}}{e^u+e^{-u}}$，将$(-\infty,+\infty)$映射到$(-1,1)$。

感知器模型，其中$w_i$为系数，$b$为常数项，$y$为分类预测值：

$$
\displaylines{
v=\sum_{i=1}^mw_ix_i+b,y=\begin{cases}
1,&v\geqslant0,\\
0,&v<0
\end{cases}
}
$$

例题：解决简单分类问题，将四个输入列向量分成两类，输入向量为：

$$
\displaylines{
P=\left[
\begin{matrix}
-0.5&-0.5&0.3&0.0\\
-0.5&0.5&-0.5&1.0
\end{matrix}
\right]
}
$$

目标分类向量$\overrightarrow T=[1,1,0,0]$，试预测新输入向量$\overrightarrow p=[-0.5,0.2]^{\mathrm T}$目标值。

代码：

from sklearn.linear_model import Perceptron
import numpy
x0=numpy.array([[-0.5,-0.5,0.3,0.0],[-0.5,0.5,-0.5,1.0]]).T
y0=numpy.array([1,1,0,0])
md=Perceptron().fit(x0,y0)
print(md.coef_,md.intercept_,md.score(x0,y0),md.predict([[-0.5,0.2]])) #系数 常数项 精度 预测值

BP神经网络

数据预处理时，将效益型数据映射到相应需要的$[0,1]$或$[1,1]$区间上，有三种方法。当激活函数为sigmoid函数时使用$\widetilde x=\dfrac{x-x_{\min}}{x_{\max}-x_{\min}}$。当激活函数为tanh函数时使用$\widetilde x=\dfrac{2(x-x_{\min})}{x_{\max}-x_{\min}}-1$。也可以用一般标准化变换$\widetilde x=\dfrac{x-\overline x}{\sigma}$。

例题：根据前十所公司信用，判断后两所公司是否为ST公司。

对效益型指标$x_1,x_3,x_4,x_5,x_6,x_7,x_8$用公式$\widetilde{x_i}=\dfrac{x_i-x_i^{\min}}{x_i^{\max}-x_i^{\min}}$，成本型指标$x_2$用公式$\widetilde{x_2}=\dfrac{x_2^{\max}-x_2}{x_2^{\max}-x_2^{\min}}$。有一个隐层，隐层神经元个数用$30$，激活函数用sigmoid函数。

代码：

from sklearn.neural_network import MLPClassifier
import numpy
a=numpy.loadtxt("xxx.txt")
x0=a[:10,:]
x=a[10:,:]
m1=x0.max(axis=0)
m2=x0.min(axis=0)
bx0=(x0-m2)/(m1-m2)
bx0[:,1]=(m1[1]-x0[:,1])/(m1[1]-m2[1])
y0=numpy.hstack([numpy.zeros(5),numpy.ones(5)])
md=MLPClassifier(solver="lbfgs",activation="logistic",hidden_layer_sizes=30).fit(bx0,y0)
bx=(x-m2)/(m1-m2)
bx[:,1]=(m1[1]-x[:,1])/(m1[1]-m2[1])
yh=md.predict(bx)
print(yh,md.predict_proba(bx),md.score(bx0,y0)) #待测样本类别 属于各类别概率 回代准确率

例题：给出1990-2009减公路客运人口数量、机动车数量、公路面积、客运量，已知2010和2011年人口数量、机动车数量、公路面积，预测客运量。

处理人口数量、机动车数量、公路面积分别为$x_1,x_2,x_3$，选择公式$\widetilde{x_i}=\dfrac{2(x_i-x_i^{\min})}{x_i^{\max}-x_i^{\min}}-1$。设置隐层神经元为$10$个，激活函数取线性函数。

代码：

from sklearn.neural_network import MLPRegressor
import numpy,pylab
a=numpy.loadtxt('xxx.txt')
x0=a[:,:3]
y0=a[:,3]
m1=x0.max(axis=0)
m2=x0.min(axis=0)
bx0=2*(x0-m2)/(m1-m2)-1
md=MLPRegressor(solver='lbfgs',activation='identity',hidden_layer_sizes=10).fit(bx0,y0)
x=numpy.array([[73.39,75.55],[3.9635,4.0975],[0.9880,1.0268]]).T
bx=2*(x-m2)/(m1-m2)-1
yh=md.predict(bx)
print(numpy.round(yh,4)) #预测值
yh0=md.predict(bx0)
delta=abs(yh0-y0)/y0*100
print(numpy.round(delta,4)) #相对误差
t=numpy.arange(1990,2010)
pylab.rc('font',sizeof=15)
pylab.rc('font',family='SimHei')
pylab.plot(t,y0,'--o',label='原始数据')
pylab.plot(t,yh0,'-*',label='预测数据')
pylab.xticks(t,rotation=55)
pylab.legend()
pylab.show()