python 实现GM21
GM11和GM21虽然都可进行预测,但是GM11是对单调性的数据进行预测(单调递增或递减),而GM21是对周期性的数据进行预测。
GM21代码如下:
# -*- coding=utf-8 -*- import numpy as np from sympy import * from sympy.abc import x, y import re import math init_printing() # 定义符号常量x 与 f(x) g(x)。这里的f g还可以用其他字母替换,用于表示函数 f, g = symbols('f g', cls=Function) def solving_equation(x1, equation_parameter, n): # 二阶齐次微分方程的根分为两个不同的实根、两个相同的实根以及两个虚根 # 不同情况方程的形式不同,根据predict函数中求得的带参数的方程来决定使用的策略 # 下面以两个不同的实根为例 parameter = solve( [x * math.exp(equation_parameter[0] * 0) + y * math.exp(equation_parameter[1] * 0) + equation_parameter[2] - x1[0], x * math.exp(equation_parameter[0] * (len(x1)-1)) + y * math.exp(equation_parameter[1] * (len(x1)-1)) + equation_parameter[2] - x1[len(x1) - 1]]) print("parameter", parameter) # 返回X1的预测值 return [parameter[x] * math.exp(equation_parameter[0] * i) + parameter[y] * math.exp(equation_parameter[1] * i) + equation_parameter[2] for i in range(len(x1) + n)] # n代表需要预测的个数,默认为0 def predict(data,n=0): x1 = data.cumsum() a_x0 = np.ediff1d(data).T z = (x1[:len(x1) - 1] + x1[1:]) / 2.0 B = np.array([-data[1:], -z, np.ones([len(data) - 1])]).T Y = a_x0 u = np.dot(np.dot(np.linalg.inv(np.dot(B.T, B)), B.T), Y) a1, a2, b = u[0], u[1], u[2] # 用diffeq代表微分方程 diffeq = Eq(f(x).diff(x, x) + a1 * f(x).diff(x) + a2 * f(x), b) # 调用dsolve函数,返回一个Eq对象,并提取带参数方程 differential_equation = str(dsolve(diffeq, f(x)).args[1]) # 使用正则表达式提取齐次微分方程的根与非齐次微分方程的特解 equation_parameter = re.findall("-?\d+.?\d+", differential_equation.replace(' ', '')) # str转为float for i in range(len(equation_parameter)): equation_parameter[i] = float(equation_parameter[i]) # 利用边界条件,取X1中第一个数和最后一个数,构造方程组,求参数C1和C2,并返回预测值 return solving_equation(x1, equation_parameter, n) if __name__ == '__main__': data = np.array([41, 49, 61, 78, 96, 104]) predict_data = predict(data) # 预测x1 result = np.ediff1d(predict_data) # 递减,求出最终结果 print('原数据:', data[1:]) print('预测结果:', result)
结果如下: