数模竞赛-预测方法

灰色预测模型

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

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

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

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

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

建立灰微分方程$x^{(0)}(k)+a{z}^{(1)}(k)=b,2⩽k⩽n$，其对应的白化微分方程为${\displaystyle \frac{\mathrm{d}{x}^{(1)}(t)}{\mathrm{d}t}}+a{x}^{(1)}(t)=b$，记：

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

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

$${\hat{x}}^{(1)}(k+1)=(x^{(0)}(1)-{\displaystyle \frac{\hat{b}}{\hat{a}}})e^{-\hat{a}k}+{\displaystyle \frac{\hat{b}}{\hat{a}}},k⩾0$$

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

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

检验预测值：

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

计算级比偏差值$\rho (k)=|1-\left({\displaystyle \frac{1-0.5\hat{a}}{1+0.5\hat{a}}}\right)\lambda (k)|$，当$\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))

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

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

设：

$$\begin{array}{c}B=\left[\begin{array}{ccc}-x^{(0)}(2)& -z^{(1)}(2)& 1\\ -x^{(0)}(3)& -z^{(1)}(3)& 1\\ ⋮& ⋮& ⋮\\ -x^{(0)}(n)& -z^{(1)}(n)& 1\end{array}\right],Y=\left[\begin{array}{c}{\alpha }^{(1)}{x}^{(0)}(2)\\ {\alpha }^{(1)}{x}^{(0)}(3)\\ ⋮\\ {\alpha }^{(1)}{x}^{(0)}(n)\end{array}\right]\end{array}$$

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

例题2：

已知${\overrightarrow{x}}^{(0)}=(41,49,61,78,96,104)$，建立$\mathrm{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 #相对误差

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

一些定义同上。

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

$$\begin{array}{c}B=\left[\begin{array}{cc}-x^{(0)}(2)& 1\\ -x^{(0)}(3)& 1\\ ⋮& ⋮\\ -x^{(0)}(n)& 1\end{array}\right]\end{array}$$

该模型参数最小二乘估计满足$\hat{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)$，建立$\mathrm{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模型

一些定义同$\mathrm{GM}(1,1)$模型，设$a,b$为参数，称$x^{(0)}(k)+a{z}^{(1)}(k)=b{[z^{(1)}(k)]}^2$，${\displaystyle \frac{\mathrm{d}{x}^{(1)}(t)}{\mathrm{d}t}}+a{x}^{(1)}(t)=b{[x^{(1)}(t)]}^2$为白化方程，$t$为时间。

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

$$\begin{array}{c}B=\left[\begin{array}{cc}-z^{(1)}(2)& {(z^{(1)}(2))}^2\\ -z^{(1)}(3)& {(z^{(1)}(3))}^2\\ ⋮& ⋮\\ -z^{(1)}(n)& {(z^{(1)}(n))}^2\end{array}\right],Y=\left[\begin{array}{c}{x}^{(0)}(2)\\ x^{(0)}(3)\\ ⋮\\ x^{(0)}(n)\end{array}\right]\end{array}$$

最小二乘估计同上，白化方程的解为$x^{(1)}(t)={\displaystyle \frac{\hat{a}{x}^{(0)}(1)}{\hat{b}{x}^{(0)}(1)+(\hat{a}-\hat{b}{x}^{(0)}(1))e^{\hat{a}t}}}$，相应的时间相应序列为$x^{(1)}(k+1)={\displaystyle \frac{\hat{a}{x}^{(0)}(1)}{\hat{b}{x}^{(0)}(1)+(\hat{a}-\hat{b}{x}^{(0)}(1))e^{\hat{a}k}}}$，同样累减还原式为${\hat{x}}^{(0)}(k+1)={\hat{x}}^{(1)}(k+1)-{\hat{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}={[x_1,x_2,\cdots ,x_m]}^{\mathrm{T}},\overrightarrow{W}={[w_1,w_2,\cdots ,2_m]}^{\mathrm{T}}$，于是网络输入为${\displaystyle u=\sum _{i=1}^m{w}_i{x}_i={\overrightarrow{W}}^{\mathrm{T}}\overrightarrow{X}}$。

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

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

$$\begin{array}{c}v=\sum _{i=1}^m{w}_i{x}_i+b,y=\{\begin{array}{ll}1,& v⩾0,\\ 0,& v<0\end{array}\end{array}$$

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

$$\begin{array}{c}P=\left[\begin{array}{cccc}-0.5& -0.5& 0.3& 0.0\\ -0.5& 0.5& -0.5& 1.0\end{array}\right]\end{array}$$

目标分类向量$\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函数时使用$\tilde{x}={\displaystyle \frac{x-x_{min}}{x_{max}-x_{min}}}$。当激活函数为tanh函数时使用$\tilde{x}={\displaystyle \frac{2(x-x_{min})}{x_{max}-x_{min}}}-1$。也可以用一般标准化变换$\tilde{x}={\displaystyle \frac{x-\bar{x}}{\sigma }}$。

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

对效益型指标$x_1,x_3,x_4,x_5,x_6,x_7,x_8$用公式$\widetilde{x_i}={\displaystyle \frac{x_i-x_i^{min}}{x_i^{max}-x_i^{min}}}$，成本型指标$x_2$用公式$\widetilde{x_2}={\displaystyle \frac{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}={\displaystyle \frac{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()