灰色预测--matlab&python实现
function SGrey X0 = input('请输入原始负荷数据：'); %输入原始数据 n = length(X0); %原始n年数据 %累加生成 X1 = zeros(1,n); for i = 1:n if i == 1 X1(1,i) = X0(1,i); else X1(1,i) = X0(1,i) + X1(1,i-1); end end X1 %计算数据矩阵B和数据向量Y B = zeros(n-1,2); Y = zeros(n-1,1); for i = 1:n-1 B(i,1) = -0.5*(X1(1,i) + X1(1,i+1)); B(i,2) = 1; Y(i,1) = X0(1,i+1); end B,Y %计算GM(1,1)微分方程的参数a和u A = zeros(2,1); A = inv(B'*B)*B'*Y; a = A(1,1); u = A(2,1); a,u %建立灰色预测模型 XX0(1,1) = X0(1,1); for i = 2:n XX0(1,i) = (X0(1,1) - u/a)*(1-exp(a))*exp(-a*(i-1)); end XX0 %模型精度的后验差检验 e = 0; %求残差平均值 for i =1:n e = e + (X0(1,i) - XX0(1,i)); end e = e/n; e aver = 0; %求历史数据平均值 for i = 1:n aver = aver + X0(1,i); end aver = aver / n; aver s12 = 0; %求历史数据方差 for i = 1:n s12 = s12 + (X0(1,i)-aver)^2; end s12 = s12 / n; s12 s22 = 0; %求残差方差 for i = 1:n s22 = s22 + ((X0(1,i) - XX0(1,i)) - e)^2; end s22 = s22 / n; s22 C = s22 / s12; %求后验差比值 C cout = 0; %求小误差概率 for i = 1:n if abs((X0(1,i) - XX0(1,i)) - e) < 0.6754*sqrt(s12) cout = cout+1; else cout = cout; end end P = cout / n; P if (C < 0.35 & P > 0.95) disp('预测精度为一级'); m = input('请输入需要预测的年数: m = '); %预测往后各年的负荷 disp('往后m各年负荷为：'); f = zeros(1,m); for i = 1:m f(1,i) = (X0(1,1) - u/a)*(1-exp(a))*exp(-a*(i+n-1)); end f else disp('灰色预测法不适用'); end
matlab输出
输入:[724.57 746.62 778.27 800.8 827.75 871.1 912.37 954.28 995.01 1037.2] 输出: >> SGrey 请输入原始负荷数据：[724.57 746.62 778.27 800.8 827.75 871.1 912.37 954.28 995.01 1037.2 ] X1 = 1.0e+003 * Columns 1 through 8 0.7246 1.4712 2.2495 3.0503 3.8780 4.7491 5.6615 6.6158 Columns 9 through 10 7.6108 8.6480 B = 1.0e+003 * -1.0979 0.0010 -1.8603 0.0010 -2.6499 0.0010 -3.4641 0.0010 -4.3136 0.0010 -5.2053 0.0010 -6.1386 0.0010 -7.1133 0.0010 -8.1294 0.0010 Y = 1.0e+003 * 0.7466 0.7783 0.8008 0.8277 0.8711 0.9124 0.9543 0.9950 1.0372 a = -0.0420 u = 693.9403 XX0 = 1.0e+003 * Columns 1 through 8 0.7246 0.7398 0.7715 0.8046 0.8391 0.8750 0.9125 0.9517 Columns 9 through 10 0.9925 1.0350 e = 0.1818 aver = 864.7970 s12 = 1.0357e+004 s22 = 26.8113 C = 0.0026 P = 1 预测精度为一级 请输入需要预测的年数: m = 10 往后m各年负荷为： f = 1.0e+003 * Columns 1 through 8 1.0794 1.1257 1.1739 1.2242 1.2767 1.3315 1.3885 1.4481 Columns 9 through 10 1.5101 1.5749
Python实现
# -*- 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)) else: print('灰色预测法不适用')
