1.生成数
在研究社会系统、经济系统等抽象系统时，往往要遇到随机干扰（即所谓“噪声”）。人们对“噪声”污染系统的研究大多基于概率统计方法。但概率统计方法有很多不足之处：要求大量数据、要求有典型的统计规律、计算工作量等。而且在某些问题中，其概率意义下的结论并不直观或信息量少。
通过对数列中的数据进行处理，产生新的数列，以此来挖掘和寻找数的规律性的方法，叫做数的生成。通过生成数，可以将原本数据间规律不明显的数列处理地更有规律性。常见的方法有累加生成、累减生成和加权累加生成等。
累加生成
把数列x 各时刻数据依次累加的过程叫做累加过程，记作AGO，累加所
得的新数列，叫做累加生成数列。
设原始数列为：
x 0 = ( x ( 0 ) ( 1 ) , x ( 0 ) ( 2 ) , ⋯ , x ( 0 ) ( n ) ) x^{0}=\left(x^{(0)}(1), x^{(0)}(2), \cdots,x^{(0)}(n)\right) x0=(x(0)(1),x(0)(2),⋯,x(0)(n))，
则累加数列记为：
x ( 1 ) = ( x ( 1 ) ( 1 ) , x ( 0 ) ( 2 ) , ⋯ , x ( 1 ) ( n ) ) x^{(1)}=\left(x^{(1)}(1), x^{(0)}(2),\cdots, x^{(1)}(n)\right) x(1)=(x(1)(1),x(0)(2),⋯,x(1)(n))
且 x ( 0 ) x^{(0)} x(0)与 x ( 1 ) x^{(1)} x(1)之间满足：
x ( 1 ) ( k ) = ∑ i = α R x ( 0 ) ( i ) , k = α , ⋯ , n x^{(1)}(k)=\sum_{i=\alpha}^{R} x^{(0)}(i), k=\alpha, \cdots, n x(1)(k)=∑i=αRx(0)(i),k=α,⋯,n
其中， α ≤ n \alpha \leq n α≤n为正整数。当 α = 1 \alpha=1 α=1时，为1次累加生成，记作 1 − A G O 1-AGO 1−AGO，当数列为 r r r次累加生成时，记作 r − A G O r-AGO r−AGO。常见的累加生成为 1 − A G O 1-AGO 1−AGO。
一般地，经济数列等实际问题的数列皆是非负数列，累加生成可使非负的摆动与非摆动的数列或任意无规律性的数列转化为非减的数列。
当然，有些实际问题的数列中有负数（例如温度等），累加时略微复杂。有时，由于出现正负抵消这种信息损失的现象，数列经过累加生成后规律性非但没得到加强，甚至可能被削弱。对于这种情形，我们可以先进行移轴，然后再做累加生成。
累减生成
累减生成一般用于对累加生成数列进行还原。
对原始数列依次做前后两数据相减的运算过程叫累减生成，记作 I A G O IAGO IAGO。若 x ( r ) x^{(r)} x(r)为 r − A G O r-AGO r−AGO数列，则称
x ( r − 1 ) ( k ) = x ( r ) ( k ) − x ( r ) ( k − 1 ) , k = 2 , 3 , ⋯ , n x^{(r-1)}(k)=x^{(r)}(k)-x^{(r)}(k-1), \quad k=2,3, \cdots, n x(r−1)(k)=x(r)(k)−x(r)(k−1),k=2,3,⋯,n
为 r r r次累减生成数列。
均值生成
设原始数列为： x 0 = ( x ( 0 ) ( 1 ) , x ( 0 ) ( 2 ) , ⋯ , x ( 0 ) ( n ) ) x^{0}=\left(x^{(0)}(1), x^{(0)}(2), \cdots,x^{(0)}(n)\right) x0=(x(0)(1),x(0)(2),⋯,x(0)(n))，则称 x ( 0 ) ( k − 1 ) x^{(0)}(k-1) x(0)(k−1)与 x ( 0 ) ( k ) x^{(0)}(k) x(0)(k)为数列 x ( 0 ) x^{(0)} x(0)的一对紧邻值，其中 x ( 0 ) ( k − 1 ) x^{(0)}(k-1) x(0)(k−1)称为前值， x ( 0 ) ( k ) x^{(0)}(k) x(0)(k)称为后值。
对于常数 α ∈ [ 0 , 1 ] \alpha \in[0,1] α∈[0,1]，则称
z ( 0 ) ( k ) = 0.5 x ( 0 ) ( k ) + 0.5 x ( 0 ) ( k − 1 ) z^{(0)}(k)=0.5 x^{(0)}(k)+0.5 x^{(0)}(k-1) z(0)(k)=0.5x(0)(k)+0.5x(0)(k−1)
为紧邻均值生成数，或等权邻值生成数。
2.灰色模型（Grey Model）
灰色系统理论是基于关联空间、光滑离散函数等概念定义灰导数与灰微分方程，进而用离散数据列建立微分方程形式的动态模型，由于这是本征灰色系统的基本模型，而且模型是近似的、非唯一的，故这种模型为灰色模型，记为GM（Grey Model），即灰色模型是利用离散随机数经过生成变为随机性被显著削弱而且较有规律的生成数，建立起的微分方程形式的模型，这样便于对其变化过程进行研究和描述。
GM（1,1）模型
设 x ( 0 ) x^{(0)} x(0)为 n n n个元素的数列： x 0 = ( x ( 0 ) ( 1 ) , x ( 0 ) ( 2 ) , ⋯ , x ( 0 ) ( n ) ) x^{0}=\left(x^{(0)}(1), x^{(0)}(2), \cdots,x^{(0)}(n)\right) x0=(x(0)(1),x(0)(2),⋯,x(0)(n))， x ( 0 ) x^{(0)} x(0)的 1 − A G O 1-AGO 1−AGO生成数列为： x ( 1 ) = ( x ( 1 ) ( 1 ) , x ( 0 ) ( 2 ) , ⋯ , x ( 1 ) ( n ) ) x^{(1)}=\left(x^{(1)}(1), x^{(0)}(2),\cdots, x^{(1)}(n)\right) x(1)=(x(1)(1),x(0)(2),⋯,x(1)(n))，则 x ( 1 ) x^{(1)} x(1)的灰导数为：
d x ( 1 ) ( k ) = x ( 0 ) ( k ) = x ( 1 ) ( k ) − x ( 1 ) ( k − 1 ) d x^{(1)}(k)=x^{(0)}(k)=x^{(1)}(k)-x^{(1)}(k-1) dx(1)(k)=x(0)(k)=x(1)(k)−x(1)(k−1)
令 z ( 0 ) z^{(0)} z(0)为 x ( 1 ) x^{(1)} x(1)的紧邻均值数列，即：
z ( 0 ) ( k ) = 0.5 x ( 0 ) ( k ) + 0.5 x ( 0 ) ( k − 1 ) z^{(0)}(k)=0.5 x^{(0)}(k)+0.5 x^{(0)}(k-1) z(0)(k)=0.5x(0)(k)+0.5x(0)(k−1)
由此定义 G M ( 1 , 1 ) GM(1,1) GM(1,1)的灰微分方程模型为：
d x ( 1 ) ( k ) + a z ( 1 ) ( k ) = b d x^{(1)}(k)+a z^{(1)}(k)=b dx(1)(k)+az(1)(k)=b
即：
x ( 0 ) ( k ) + a z ( 1 ) ( k ) = b x^{(0)}(k)+a z^{(1)}(k)=b x(0)(k)+az(1)(k)=b
将时刻 k = 2 , 3 , ⋯ , n k=2,3, \cdots, n k=2,3,⋯,n代入上式，有：
{ x ( 0 ) ( 2 ) + a z ( 1 ) ( 2 ) = b x ( 0 ) ( 3 ) + a z ( 1 ) ( 3 ) = b ⋯ ⋯ ⋯ x ( 0 ) ( n ) + a z ( 1 ) ( n ) = b \left\{\begin{array}{l}x^{(0)}(2)+a z^{(1)}(2)=b \\ x^{(0)}(3)+a z^{(1)}(3)=b \\ \cdots \cdots \cdots \\ x^{(0)}(n)+a z^{(1)}(n)=b\end{array}\right. ⎩⎪⎪⎨⎪⎪⎧x(0)(2)+az(1)(2)=bx(0)(3)+az(1)(3)=b⋯⋯⋯x(0)(n)+az(1)(n)=b
令： Y = ( x ( 0 ) ( 2 ) , x ( 0 ) ( 3 ) , ⋯ , x ( 0 ) ( n ) ) T , u = ( a , b ) T , B = [ − z ( 1 ) ( 2 ) 1 − z ( 1 ) ( 3 ) 1 ⋮ ⋮ − z ( 1 ) ( n ) 1 ] Y=\left(x^{(0)}(2), x^{(0)}(3), \cdots, x^{(0)}(n)\right)^{T}, \quad u=(a, b)^{T}, \quad B=\left[\begin{array}{cc}-z^{(1)}(2) & 1 \\ -z^{(1)}(3) & 1 \\ \vdots & \vdots \\ -z^{(1)}(n) & 1\end{array}\right] Y=(x(0)(2),x(0)(3),⋯,x(0)(n))T,u=(a,b)T,B=⎣⎢⎢⎢⎡−z(1)(2)−z(1)(3)⋮−z(1)(n)11⋮1⎦⎥⎥⎥⎤，
称为 Y Y Y数据向量， B B B数据矩阵， u u u参数向量，则 G M ( 1 , 1 ) GM(1,1) GM(1,1)模型可表示为矩阵方程： Y = B u Y=Bu Y=Bu
由最小二乘法可得：
u ^ = ( a ^ , b ^ ) T = ( B T B ) − 1 B T Y \hat{u}=(\hat{a}, \hat{b})^{T}=\left(B^{T} B\right)^{-1} B^{T} Y u^=(a^,b^)T=(BTB)−1BTY
GM（1,1）的白化型
若将时刻 k = 2 , 3 , ⋯ , n k=2,3, \cdots, n k=2,3,⋯,n视为连续的变量 t t t，则数列 x ( 1 ) x^{(1)} x(1)可视为变量 t t t的函数，即： x ( 1 ) = x ( 1 ) ( t ) x^{(1)}=x^{(1)}(t) x(1)=x(1)(t)。令灰导数 x ( 0 ) ( k ) x^{(0)}(k) x(0)(k)对应于导数 d x ( 1 ) d t \frac{d x^{(1)}}{d t} dtdx(1)，背景值 z ( 1 ) ( k ) z^{(1)}(k) z(1)(k)对应于 x ( 1 ) ( t ) x^{(1)}(t) x(1)(t)，得到上述灰微分方程对应的白微分方程：
d x ( 1 ) d t + a x ( 1 ) = b \frac{d x^{(1)}}{d t}+a x^{(1)}=b dtdx(1)+ax(1)=b
称为 G M ( 1 , 1 ) GM(1,1) GM(1,1)的白化型。
GM（1,N）模型
即1阶方程，包含有N 个变量的灰色模型。
3.灰色预测实例
北方某城市1986～1992 年道路交通噪声平均声级数据：
预测未来两年的声级数据。
（1）数据的检验与处理
首先，为了保证建模方法的可行性，需要对已知数据列做必要的检验处理。设参考数据为： x 0 = ( x ( 0 ) ( 1 ) , x ( 0 ) ( 2 ) , ⋯ , x ( 0 ) ( n ) ) x^{0}=\left(x^{(0)}(1), x^{(0)}(2), \cdots,x^{(0)}(n)\right) x0=(x(0)(1),x(0)(2),⋯,x(0)(n))，计算数列级比：
λ ( k ) = x ( 0 ) ( k − 1 ) x ( 0 ) ( k ) , k = 2 , 3 , ⋯ , n \lambda(k)=\frac{x^{(0)}(k-1)}{x^{(0)}(k)}, \quad k=2,3, \cdots, n λ(k)=x(0)(k)x(0)(k−1),k=2,3,⋯,n
如果所有的级比 λ ( k ) \lambda(k) λ(k)都落在可容覆盖 Θ = ( e − 2 n + 1 , e 2 n + 2 ) \Theta=\left(e^{-\frac{2}{n+1}}, e^{\frac{2}{n+2}}\right) Θ=(e−n+12,en+22)内，则数列可作为模型 G M ( 1 , 1 ) GM(1,1) GM(1,1)的数据进行灰色预测。否则，要对数据进行平移变换。
对于上述数据，级比的范围如下图所示：
故可以进行灰色预测。
作图代码如下：
#数据检验
lamda = []
for i in range(len(data)-1):
lamda.append(data[i+1] / data[i])
x = np.arange(len(lamda))
y1 = np.array(lamda)
y2 = np.array([math.exp(-(2/(len(data)+1)))]*len(x))
y3 = np.array([math.exp(2/(len(data)+2))]*len(x))
plt.plot(x,y1,x,y2,x,y3)
（2）建立模型
建立 G M ( 1 , 1 ) GM(1,1) GM(1,1)模型，解上述微分方程，可得到微分方程的解为：
x ^ ( 1 ) ( k + 1 ) = ( x ( 0 ) ( 1 ) − b a ) e − a k + b a , k = 0 , 1 , ⋯ , n − 1 , ⋯ \hat{x}^{(1)}(k+1)=\left(x^{(0)}(1)-\frac{b}{a}\right) e^{-a k}+\frac{b}{a}, \quad k=0,1, \cdots, n-1, \cdots x^(1)(k+1)=(x(0)(1)−ab)e−ak+ab,k=0,1,⋯,n−1,⋯
且： x ^ ( 0 ) ( k + 1 ) = x ^ ( 1 ) ( k + 1 ) − x ^ ( 1 ) ( k ) , k = 1 , 2 , ⋯ , n − 1 , ⋯ \hat{x}^{(0)}(k+1)=\hat{x}^{(1)}(k+1)-\hat{x}^{(1)}(k), \quad k=1,2, \cdots, n-1, \cdots x^(0)(k+1)=x^(1)(k+1)−x^(1)(k),k=1,2,⋯,n−1,⋯
建模代码如下：
#进行1次累加（AGO）生成
x0 = np.array(data)
x1 = np.cumsum(x0)
#构建GM（1,1）模型
n = len(x1)
z1 = -0.5*(x1[:n-1] + x1[1:n])
B = np.ones((n-1,2))
B[:,0] = z1
Y = x0[1:]
u = np.linalg.inv(B.T@B) @ B.T @ Y
a,b = u[0],u[1]
#求生成数列预测值
predict_x1 = [x0[0]]
for i in range(len(x1)):
predict_i = (x0[0] - b/a) * np.exp(-a*(i+1)) + b/a
predict_x1.append(predict_i)
predict_x0 = np.diff(predict_x1)
predict_x0 = np.insert(predict_x0,0,x0[0])
（3）预测值检验
可通过计算残差与相对误差来检验。
epsilon = x0 - predict_x0[:-1]
delta = np.abs(epsilon/x0)
xl1 = np.arange(len(predict_x0))
xl2 = np.arange(len(x0))
plt.plot(xl1,predict_x0,xl2,x0)
print('the prediction is:',predict_x0)
print('the Residual is:',epsilon)
print('the Relative Error is:',delta)
预测值与真实值结果对比如下图：
残差与相对误差：
对比通过Matlab得到的结果：
全部代码如下：
import math
import numpy as np
import matplotlib.pyplot as plt
data0 = open('voice.txt').read().replace(' ',',')
data = list(eval(data0))
''' #数据检验 lamda = [] for i in range(len(data)-1): lamda.append(data[i+1] / data[i]) x = np.arange(len(lamda)) y1 = np.array(lamda) y2 = np.array([math.exp(-(2/(len(data)+1)))]*len(x)) y3 = np.array([math.exp(2/(len(data)+2))]*len(x)) plt.plot(x,y1,x,y2,x,y3)'''
#进行1次累加（AGO）生成
x0 = np.array(data)
x1 = np.cumsum(x0)
#构建GM（1,1）模型
n = len(x1)
z1 = -0.5*(x1[:n-1] + x1[1:n])
B = np.ones((n-1,2))
B[:,0] = z1
Y = x0[1:]
u = np.linalg.inv(B.T@B) @ B.T @ Y
a,b = u[0],u[1]
#求生成数列预测值
predict_x1 = [x0[0]]
for i in range(len(x1)):
predict_i = (x0[0] - b/a) * np.exp(-a*(i+1)) + b/a
predict_x1.append(predict_i)
predict_x0 = np.diff(predict_x1)
predict_x0 = np.insert(predict_x0,0,x0[0])
#模型检验
epsilon = x0 - predict_x0[:-1]
delta = np.abs(epsilon/x0)
xl1 = np.arange(len(predict_x0))
xl2 = np.arange(len(x0))
plt.plot(xl1,predict_x0,xl2,x0)
print('the prediction is:',predict_x0)
print('the Residual is:',epsilon)
print('the Relative Error is:',delta)
参考文献
司守奎, 徐珂文, 李日华. 数学建模算法与程序[J]. 海军航空工程学院, 2007, 9: 95-98.
