灰色預測模型GM(1,1)
灰色預測模型\(GM(1,1)\)是在數學建模比賽中常用的預測值方法，常用於中短期符合指數規律的預測。
其數學表達與原理分析參考文章尾部網頁與文獻資料。
- 預處理
灰色模型要求數據前后級比落入范圍 \(\displaystyle \Theta\left(e^{-\frac{2}{n+1}},e^{\frac{2}{n+2}}\right)\) ，因此做線性平移預處理使得元數據滿足要求。
該預處理保證了數據一定落入級比范圍內。
線性平移：將數據平移至不小於1，檢查級比，若不滿足要求則將數據向上平移一個最小值直到滿足要求。
可以推斷出，級比的上下界在給定數據點數越多的情況下，越趨於1。
- 代碼
經過整理，以下附上Python代碼：
import numpy as np
import matplotlib.pyplot as plt
# 線性平移預處理，確保數據級比在可容覆蓋范圍
def greyModelPreprocess(dataVec):
"Set linear-bias c for dataVec"
import numpy as np
from scipy import io, integrate, linalg, signal
from scipy.sparse.linalg import eigs
from scipy.integrate import odeint
c = 0
x0 = np.array(dataVec, float)
n = x0.shape[0]
L = np.exp(-2/(n+1))
R = np.exp(2/(n+2))
xmax = x0.max()
xmin = x0.min()
if (xmin < 1):
x0 += (1-xmin)
c += (1-xmin)
xmax = x0.max()
xmin = x0.min()
lambda_ = x0[0:-1] / x0[1:] # 計算級比
lambda_max = lambda_.max()
lambda_min = lambda_.min()
while (lambda_max > R or lambda_min < L):
x0 += xmin
c += xmin
xmax = x0.max()
xmin = x0.min()
lambda_ = x0[0:-1] / x0[1:]
lambda_max = lambda_.max()
lambda_min = lambda_.min()
return c
# 灰色預測模型
def greyModel(dataVec, predictLen):
"Grey Model for exponential prediction"
# dataVec = [1, 2, 3, 4, 5, 6]
# predictLen = 5
import numpy as np
from scipy import io, integrate, linalg, signal
from scipy.sparse.linalg import eigs
from scipy.integrate import odeint
x0 = np.array(dataVec, float)
n = x0.shape[0]
x1 = np.cumsum(x0)
B = np.array([-0.5 * (x1[0:-1] + x1[1:]), np.ones(n-1)]).T
Y = x0[1:]
u = linalg.lstsq(B, Y)[0]
def diffEqu(y, t, a, b):
return np.array(-a * y + b)
t = np.arange(n + predictLen)
sol = odeint(diffEqu, x0[0], t, args=(u[0], u[1]))
sol = sol.squeeze()
res = np.hstack((x0[0], np.diff(sol)))
return res
# 輸入數據
x = np.array([-18, 0.34, 4.68, 8.49, 29.84, 50.21, 77.65, 109.36])
c = greyModelPreprocess(x)
x_hat = greyModel(x+c, 5)-c
# 畫圖
t1 = range(x.size)
t2 = range(x_hat.size)
plt.plot(t1, x, color='r', linestyle="-", marker='*', label='True')
plt.plot(t2, x_hat, color='b', linestyle="--", marker='.', label="Predict")
plt.legend(loc='upper right')
plt.xlabel('xlabel')
plt.ylabel('ylabel')
plt.title('Prediction by Grey Model (GM(1,1))')
plt.show()
灰色預測模型GM(2,1)
根據資料，\(GM(2,1)\)適用於非單調的擺動發展序列，或有飽和的S型序列。但是從圖像上觀察，數據預測由於為指數類型，變化過於誇張、預測趨勢也有偏離的狀況，可能實用性和普適性並不如\(GM(1,1)\)。
# 灰色預測模型GM(2,1)
def greyModel2(dataVec, predictLen):
"Grey Model for exponential prediction"
# dataVec = [1, 2, 3, 4, 5, 6]
# predictLen = 5
import numpy as np
import sympy as sy
from scipy import io, integrate, linalg, signal
x0 = np.array(dataVec, float)
n = x0.shape[0]
a_x0 = np.diff(x0) # 1次差分序列
x1 = np.cumsum(x0) # 1次累加序列
z = 0.5 * (x1[0:-1] + x1[1:]) # 均值生成序列
B = np.array([-x0[1:], -z, np.ones(n-1)]).T
u = linalg.lstsq(B, a_x0)[0]
def diffEqu(x, f, a1, a2, b):
return sy.diff(f(x), x, 2) + a1*sy.diff(f(x), x) + a2*f(x) - b # f''(x)+a1*f'(x)+a2*f(x)=b 二階常系數齊次微分方程
t = np.arange(n + predictLen)
x = sy.symbols('x') # 約定變量
f = sy.Function('f') # 約定函數
eq = sy.dsolve(diffEqu(x, f, u[0], u[1], u[2]), f(x), ics={f(t[0]): x1[0], f(t[n-1]): x1[-1]})
f = sy.lambdify(x, eq.args[1], 'numpy')
sol = f(t)
res = np.hstack((x0[0], np.diff(sol)))
return res
誤差分析部分：可就絕對誤差、相對誤差、級比、殘差做數據分析，以下示例為最小二乘法線性回歸分析。
def regressionAnalysis(xVec, yVec):
import numpy as np
from scipy import linalg
x = np.array([xVec, np.ones(xVec.size)]).T
u = linalg.lstsq(x, yVec)
return u
res = regressionAnalysis(x_hat[0:x.size], x)
print(res[0]) # 回歸系數a, b
print(res[1]) # residuals 殘差平方和
參考
- [1] Improving Grey Prediction Model and Its Application in Predicting the Number of Users of a Public Road Transportation System
- [2] The analysis of GM (1, 1) grey model to predict the incidence trend of typhoid and paratyphoid fevers in Wuhan City, China
- [3] Discrete grey forecasting model and its optimization
- [4] 基於python的灰色預測模型
- [5] 司守奎,孫璽菁.數學建模算法與應用[M]北京:國防工業出版社,2011:399-402.