灰色预测模型：通过少量的、不完全的信息来预测。（如：极端气候预测）
针对：时间序列短、统计数据少、信息不完全
常用的灰色预测有四种：
（1）数列预测：如1 3 5 7 9 ？
（2）灾变与异常值预测：通过灰色模型预测异常值出现的时刻，预测异常值什么时候出现在特定时区内
（3）拓扑预测：将原始数据作曲线，在曲线上按定值寻找该定值发生的所有时点，并以该定值为框架构成时点数列，然后建立模型预测该定值所发生的时点
（4）系统预测：通过对系统行为特征指标建立一组相互关联的灰色预测模型，预测系统中众多变量间的相互协调关系的变化
建立一个基于模型的灰色预测
1.数据预处理
举个例子：
原始数据序列如下
对数据累加
得一个新数据序列
归纳上面的式子写为
此式简称为一次累加生成
现将两个序列分别作图
可见原始数据的起伏已显著弱化
现考虑后减运算
归纳式子写为
2.建模原理
设满足一阶常微分方程
a为常数，称为发展灰数
u为内生控制灰数，是对系统的常定输入
此方程满足初始条件如下
解得
解中有a和u两个未知数，需要两组数据即可求解，但初始数据多于两组
现使用最小二乘法找最合适的解
可推导得
现想要求a和u，故令矩阵形式为
该方程组的最小二乘估计为
时间响应方程如下
用后减运算还原即可
3.精度检验
由于模型是基于一阶常微分方程建立的，故称为一阶一元灰色模型，记为GM(1,1)
须指出的是， 建模时先要作一次累加，因此要求原始数据均为非负数
4.GM(1,1)的建模步骤
# -*- coding: utf-8 -*-
"""
Spyder Editor
This is a temporary script file.
"""
import numpy as np
import math
history_data = [724.57,746.62,778.27,800.8,827.75,871.1,912.37,954.28,995.01,1037.2]
n = len(history_data)
X0 = np.array(history_data)
#累加生成
history_data_agg = [sum(history_data[0:i+1]) for i in range(n)]
X1 = np.array(history_data_agg)
#计算数据矩阵B和数据向量Y
B = np.zeros([n-1,2])
Y = np.zeros([n-1,1])
for i in range(0,n-1):
B[i][0] = -0.5*(X1[i] + X1[i+1])
B[i][1] = 1
Y[i][0] = X0[i+1]
#计算GM(1,1)微分方程的参数a和u
#A = np.zeros([2,1])
A = np.linalg.inv(B.T.dot(B)).dot(B.T).dot(Y)
a = A[0][0]
u = A[1][0]
#建立灰色预测模型
XX0 = np.zeros(n)
XX0[0] = X0[0]
for i in range(1,n):
XX0[i] = (X0[0] - u/a)*(1-math.exp(a))*math.exp(-a*(i));
#模型精度的后验差检验
e = 0 #求残差平均值
for i in range(0,n):
e += (X0[i] - XX0[i])
e /= n
#求历史数据平均值
aver = 0;
for i in range(0,n):
aver += X0[i]
aver /= n
#求历史数据方差
s12 = 0;
for i in range(0,n):
s12 += (X0[i]-aver)**2;
s12 /= n
#求残差方差
s22 = 0;
for i in range(0,n):
s22 += ((X0[i] - XX0[i]) - e)**2;
s22 /= n
#求后验差比值
C = s22 / s12
#求小误差概率
cout = 0
for i in range(0,n):
if abs((X0[i] - XX0[i]) - e) < 0.6754*math.sqrt(s12):
cout = cout+1
else:
cout = cout
P = cout / n
if (C < 0.35 and P > 0.95):
#预测精度为一级
m = 10 #请输入需要预测的年数
print('往后m各年负荷为：')
f = np.zeros(m)
for i in range(0,m):
f[i] = (X0[0] - u/a)*(1-math.exp(a))*math.exp(-a*(i+n))
print(f)
else:
print('灰色预测法不适用')