本文已参与「新人创作礼」活动，一起开启掘金创作之路。
上文讲到想运行程序生成多变量灰色预测结果，遂将MATLAB转成python版的。
那么优化后的代码参考如下：
#!/usr/bin/env python
# coding=utf-8
import re
import numpy as np
def summed_array(arr):
a = np.zeros(arr.shape)
for i in range(arr.shape[0]):
s = 0
for j in range(arr.shape[1]):
s += arr[i, j]
a[i, j] = s
return a
def gen_list(first, start, end, arr, index_delta, P, X0, *thetas):
a = [first]
t1, t2, t3, t4 = thetas
for k in range(start, end):
s1 = s2 = 0
for u in range(k - 1):
for i in range(1, X0.shape[0]):
s1 += t1 * t2 ** u * P[i - 1] * arr[i - index_delta, k - u]
s2 = sum(t2 ** v * ((k - v + 1) * t3 + t4) for v in range(k - 1))
s3 = t2 ** k * X0[0, 0]
a.append(s1 + s2 + s3)
return a
def matlab_string_to_array(s):
ss = s.strip().strip('[]').split(';')
a = [[float(i) for i in re.findall(r'[-0-9.]+', j)] for j in ss]
return np.array(a)
def main():
print("OBGM(1, N)模型:")
tip0 = "请输入因变量序列(第一个分号前面部分)与自变量序列(多个自变量序列之间用分号隔开)，如：[1 2 3;4 5 6]\n:"
tip00 = "请输入预测时的自变量序列(多个自变量序列之间用分号隔开)，如：[911 12 13]\n:"
try:
X0, X00 = [matlab_string_to_array(raw_input(t)) for t in (tip0, tip00)]
except NameError:
X0, X00 = [matlab_string_to_array(input(t)) for t in (tip0, tip00)]
"""
# 书上的数据
X0 = [
[287.6, 312.8, 350.4, 401.9, 480.9, 498.3, 520, 543.7],
[19978, 21989, 24725, 26738, 29073, 32903, 36469, 40321],
[3371.5, 3966.6, 3748.5, 4858.4, 5493.5, 5910.6, 6462.8, 7032.2],
[1848.5, 2053.3, 2131.7, 2303.1, 2764, 3048.8, 3294.3, 3566.4],
]
X00 = [[43910], [7562.3], [3746.8]]
X0 = np.array(X0)
X00 = np.array(X00)
"""
epsilon = 0.5
n, m = X0.shape
F = summed_array(X0)
# Calculate mode parameters
Y = X0[0][1:].T
B = np.delete(F.T[1:], 0, axis=1)
Bm = [
(epsilon - 1) * F[0, k] - epsilon * F[0, k + 1] for k in range(m - 1)
]
Li = [k + 2 for k in range(m - 1)]
Co = [1] * (m - 1)
for L in (Bm, Li, Co):
B = np.c_[B, np.array(L).T]
# Calculate b2,b3,...
P = np.linalg.inv(B.T.dot(B)).dot(B.T).dot(Y)
in_var_number = n - 1
a, c, d = P[in_var_number: in_var_number + 3]
t = 1 + epsilon * a
thetas = (1 / t, 1 - a / t, c / t, d / t)
# Calculate Simu
Simu1 = gen_list(X0[0, 0], 1, m, F, 0, P, X0, *thetas)
# Calcute g
Xf = X0[1:]
X0_F = np.c_[Xf, X00]
X0_F_1 = summed_array(X0_F)
Fore1 = gen_list(F[0, m - 1], m, X0_F.shape[1], X0_F_1, 1, P, X0, *thetas)
Fore0 = [round(a - b, 3) for a, b in zip(Fore1[1:], Fore1[:-1])]
# A
A = np.zeros((n, 5))
A[0, 0] = 1
A[0, 1] = A[0, 2] = X0[0, 0]
# mp
mp = 0
for k in range(1, len(X0)):
A[k, 0] = k + 1
A[k, 1] = X0[0, k]
A[k, 2] = round(Simu1[k] - Simu1[k - 1], 3)
A[k, 3] = round(A[k, 2] - A[k, 1], 3)
A[k, 4] = round(100 * abs(A[k, 3]) / A[k, 1], 3)
mp += 100 * abs(A[k, 3]) / A[k, 1]
mp = round(mp / (len(X0) - 1), 4)
print("-" * 20)
print("(1)模型参数b2,b3,..,bn及a,c,d分别为:")
print(np.around(P, 6))
print("\n(2)误差检验表为:")
print(" 序号 实际数据 模拟数据 残差 相对误差(%)")
print(A)
print("\n(3)平均相对模拟百分误差为(%)")
print(mp)
print("\n(4)未来[{}]步的预测值分别为:".format(len(Fore0)))
print(Fore0)
if __name__ == "__main__":
main()
