GM（1，1）和GM（1，n）预测模型

%建立符号变量a(发展系数)和b(灰作用量)
syms a b;
c &＃61; [a b]&＃39;;%原始数列 A
A &＃61; [353.12,354.39,355.61,356.45,357.1,358.83,360.82,362.61,363.73,366.7,368.38,369.55,371.14,373.28,375.8,377.52,379.8,381.9,383.79 385.6,387.43,389.9,391.65,393.85,396.52,398.65,400.83,404.24,406.55,407.11;];
n &＃61; length(A);%对原始数列 A 做累加得到数列 B
B &＃61; cumsum(A);%对数列 B 做紧邻均值生成
for i &＃61; 2:nC(i) &＃61; (B(i) &＃43; B(i - 1))/2;
end
C(1) &＃61; [];%构造数据矩阵
B &＃61; [-C;ones(1,n-1)];
Y &＃61; A; Y(1) &＃61; []; Y &＃61; Y&＃39;;%使用最小二乘法计算参数 a(发展系数)和b(灰作用量)
c &＃61; inv(B*B&＃39;)*B*Y;
c &＃61; c&＃39;;
a &＃61; c(1); b &＃61; c(2);%预测后续数据
F &＃61; []; F(1) &＃61; A(1);
for i &＃61; 2:(n&＃43;25)F(i) &＃61; (A(1)-b/a)/exp(a*(i-1))&＃43; b/a;
end%对数列 F 累减还原,得到预测出的数据
G &＃61; []; G(1) &＃61; A(1);
for i &＃61; 2:(n&＃43;25)G(i) &＃61; F(i) - F(i-1); %得到预测出来的数据
enddisp(&＃39;预测数据为&＃xff1a;&＃39;);
G%模型检验H &＃61; G(1:30);
%计算残差序列
epsilon &＃61; A - H;%法一&＃xff1a;相对残差Q检验
%计算相对误差序列
delta &＃61; abs(epsilon./A);
%计算相对误差Q
disp(&＃39;相对残差Q检验&＃xff1a;&＃39;)
Q &＃61; mean(delta)%法二&＃xff1a;方差比C检验
disp(&＃39;方差比C检验&＃xff1a;&＃39;)
C &＃61; std(epsilon, 1)/std(A, 1)%法三&＃xff1a;小误差概率P检验
S1 &＃61; std(A, 1);
tmp &＃61; find(abs(epsilon - mean(epsilon))<0.6745 * S1);
disp(&＃39;小误差概率P检验&＃xff1a;&＃39;)
P &＃61; length(tmp)/n%绘制曲线图
t1 &＃61; 1989:2018;
t2 &＃61; 1989:2043;plot(t1, A,&＃39;ro&＃39;);
hold on;
plot(t2, G, &＃39;g-&＃39;);
xlabel(&＃39;年份&＃39;);
ylabel(&＃39;全球冰浓度年平均值/L&＃39;);
legend(&＃39;实际全球冰浓度年平均值&＃39;,&＃39;预测全球冰浓度年平均值&＃39;);
title(&＃39;全球冰浓度年平均值增长曲线&＃39;);
grid on;

预测数据为&＃xff1a;

G &＃61;

  列 1 至 15

  353.1200  351.7591  353.5883  355.4271  357.2755  359.1334  361.0010  362.8784  364.7654  366.6623  368.5691  370.4858  372.4125  374.3491  376.2959

  列 16 至 30

  378.2527  380.2198  382.1971  384.1846  386.1825  388.1908  390.2095  392.2387  394.2785  396.3289  398.3899  400.4617  402.5442  404.6376  406.7419

  列 31 至 45

  408.8571  410.9833  413.1205  415.2689  417.4284  419.5992  421.7812  423.9747  426.1795  428.3957  430.6236  432.8629  435.1140  437.3767  439.6512

  列 46 至 55

  441.9376  444.2358  446.5460  448.8682  451.2024  453.5488  455.9074  458.2783  460.6615  463.0571

相对残差Q检验&＃xff1a;

Q &＃61;

    0.0019

方差比C检验&＃xff1a;

C &＃61;

    0.0561

小误差概率P检验&＃xff1a;

P &＃61;

     1

#!/usr/bin/env python
# coding&＃61;utf-8
import re
import numpy as npdef summed_array(arr):a &＃61; np.zeros(arr.shape)for i in range(arr.shape[0]):s &＃61; 0for j in range(arr.shape[1]):s &＃43;&＃61; arr[i, j]a[i, j] &＃61; sreturn adef gen_list(first, start, end, arr, index_delta, P, X0, *thetas):a &＃61; [first]t1, t2, t3, t4 &＃61; thetasfor k in range(start, end):s1 &＃61; s2 &＃61; 0for u in range(k - 1):for i in range(1, X0.shape[0]):s1 &＃43;&＃61; t1 * t2 ** u * P[i - 1] * arr[i - index_delta, k - u]s2 &＃61; sum(t2 ** v * ((k - v &＃43; 1) * t3 &＃43; t4) for v in range(k - 1))s3 &＃61; t2 ** k * X0[0, 0]a.append(s1 &＃43; s2 &＃43; s3)return adef matlab_string_to_array(s):ss &＃61; s.strip().strip(&＃39;[]&＃39;).split(&＃39;;&＃39;)a &＃61; [[float(i) for i in re.findall(r&＃39;[-0-9.]&＃43;&＃39;, j)] for j in ss]return np.array(a)def main():print("OBGM(1, N)模型:")tip0 &＃61; "请输入因变量序列(第一个分号前面部分)与自变量序列(多个自变量序列之间用分号隔开)&＃xff0c;如&＃xff1a;[1 2 3;4 5 6]\n:"tip00 &＃61; "请输入预测时的自变量序列(多个自变量序列之间用分号隔开)&＃xff0c;如&＃xff1a;[911 12 13]\n:"try:X0, X00 &＃61; [matlab_string_to_array(raw_input(t)) for t in (tip0, tip00)]except NameError:X0, X00 &＃61; [matlab_string_to_array(input(t)) for t in (tip0, tip00)]"""# 书上的数据X0 &＃61; [[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 &＃61; [[43910], [7562.3], [3746.8]]X0 &＃61; np.array(X0)X00 &＃61; np.array(X00)"""epsilon &＃61; 0.5n, m &＃61; X0.shapeF &＃61; summed_array(X0)# Calculate mode parametersY &＃61; X0[0][1:].TB &＃61; np.delete(F.T[1:], 0, axis&＃61;1)Bm &＃61; [(epsilon - 1) * F[0, k] - epsilon * F[0, k &＃43; 1] for k in range(m - 1)]Li &＃61; [k &＃43; 2 for k in range(m - 1)]Co &＃61; [1] * (m - 1)for L in (Bm, Li, Co):B &＃61; np.c_[B, np.array(L).T]# Calculate b2,b3,...P &＃61; np.linalg.inv(B.T.dot(B)).dot(B.T).dot(Y)in_var_number &＃61; n - 1a, c, d &＃61; P[in_var_number: in_var_number &＃43; 3]t &＃61; 1 &＃43; epsilon * athetas &＃61; (1 / t, 1 - a / t, c / t, d / t)# Calculate SimuSimu1 &＃61; gen_list(X0[0, 0], 1, m, F, 0, P, X0, *thetas)# Calcute gXf &＃61; X0[1:]X0_F &＃61; np.c_[Xf, X00]X0_F_1 &＃61; summed_array(X0_F)Fore1 &＃61; gen_list(F[0, m - 1], m, X0_F.shape[1], X0_F_1, 1, P, X0, *thetas)Fore0 &＃61; [round(a - b, 3) for a, b in zip(Fore1[1:], Fore1[:-1])]# AA &＃61; np.zeros((n, 5))A[0, 0] &＃61; 1A[0, 1] &＃61; A[0, 2] &＃61; X0[0, 0]# mpmp &＃61; 0for k in range(1, len(X0)):A[k, 0] &＃61; k &＃43; 1A[k, 1] &＃61; X0[0, k]A[k, 2] &＃61; round(Simu1[k] - Simu1[k - 1], 3)A[k, 3] &＃61; round(A[k, 2] - A[k, 1], 3)A[k, 4] &＃61; round(100 * abs(A[k, 3]) / A[k, 1], 3)mp &＃43;&＃61; 100 * abs(A[k, 3]) / A[k, 1]mp &＃61; 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__ &＃61;&＃61; "__main__":main()

runfile(&＃39;/home/rosfun/Downloads/PyTorch-YOLOv3-master/GM.py&＃39;, wdir&＃61;&＃39;/home/rosfun/Downloads/PyTorch-YOLOv3-master&＃39;)
OBGM(1, N)模型:
请输入因变量序列(第一个分号前面部分)与自变量序列(多个自变量序列之间用分号隔开)&＃xff0c;如&＃xff1a;[1 2 3;4 5 6]
:>? [14.24,14.39,14.38,14.19,14.21,14.29,14.43,14.33,14.46,14.62,14.40,14.41,14.53,14.62,14.60,14.51,14.66,14.59,14.63,14.49,14.60,14.67,14.55,14.57,14.60,14.68,14.73,14.81,14.92,15.0;353.12    354.39    355.61    356.45    357.1    358.83    360.82    362.61    363.73    366.7    368.38    369.55    371.14    373.28    375.8    377.52    379.8    381.9    383.79    385.6    387.43    389.9    391.65    393.85    396.52    398.65    400.83    404.24    406.55    407.11;13.281    13.3149    13.2771    13.2418    13.2517    13.3213    13.2737    13.319    13.5334    13.3334    13.2779    13.3559    13.447    13.4413    13.452    13.4793    13.4416    13.49    13.3087    13.4371    13.575    13.4271    13.3958    13.508    13.4925    13.5685    13.7189    13.6053    13.5419    13.6738;221.6520151    222.1802287    221.2058581    219.3821231    219.6195571    219.95827    221.5673472    220.3310659    221.136813    221.7323415    219.6845982    219.3498937    220.8361833    220.7716359    222.3663721    222.0425624    222.9544563    222.3201405    223.3460201    222.1782235    222.8916138    221.7041123    222.9582143    222.0703228    221.2904285    221.0724525    221.0537916    221.1392098    221.2785377    222.0068288;11067.38245    11067.15991    11067.20213    11067.2973    11067.21165    11067.41693    11067.23402    11067.38291    11067.15333    11067.27023    11080.3715    11080.54704    11080.26253    11076.14917    11076.36529    10726.12941    9982.330009    9910.526507    10024.66611    10025.18875    10025.17692    10025.15234    9798.937467    9663.557062    9615.226197    9925.141851    10024.83751    10024.08732    10023.65689    10024.14915]
请输入预测时的自变量序列(多个自变量序列之间用分号隔开)&＃xff0c;如&＃xff1a;[911 12 13]
:>? [408.8570584 410.983255    413.1205086    415.2688767    417.4284169    419.5991876    421.781247    423.9746538    426.1794671    428.3957462    430.6235507    432.8629406    435.113976    437.3767176    439.6512262    441.9375631    444.2357897    446.5459678    448.8681597    451.2024278    453.5488348    455.907444    458.2783187    460.6615228    463.0571204;13.61691872    13.62941986    13.64193248    13.65445659    13.66699219    13.67953931    13.69209794    13.7046681    13.71724981    13.72984306    13.74244788    13.75506426    13.76769223    13.78033179    13.79298296    13.80564574    13.81832015    13.83100619    13.84370388    13.85641323    13.86913424    13.88186693    13.89461131    13.9073674    13.92013519;222.2003166    222.2541984    222.3080933    222.3620013    222.4159223    222.4698564    222.5238036    222.5777639    222.6317373    222.6857237    222.7397233    222.7937359    222.8477616    222.9018005    222.9558524    223.0099174    223.0639956    223.1180869    223.1721913    223.2263088    223.2804394    223.3345832    223.38874    223.4429101    223.4970932;9657.134207    9603.439593    9550.043526    9496.944346    9444.140403    9391.630055    9339.41167    9287.483624    9235.844303    9184.492102    9133.425424    9082.642682    9032.142297    8981.922699    8931.982327    8882.319628    8832.933059    8783.821084    8734.982177    8686.414818    8638.117499    8590.088717    8542.326981    8494.830804    8447.598711]
--------------------
(1)模型参数b2,b3,..,bn及a,c,d分别为:
[ 1.376000e-03  5.531890e-01  5.450000e-04  6.000000e-05  2.486720e-01
 -5.069032e&＃43;00  1.255887e&＃43;01]
(2)误差检验表为:
   序号    实际数据    模拟数据    残差    相对误差(%)
[[  1.     14.24   14.24    0.      0.   ]
 [  2.     14.39   -3.149 -17.539 121.883]
 [  3.     14.38   18.21    3.83   26.634]
 [  4.     14.19   17.323   3.133  22.079]
 [  5.     14.21   16.639   2.429  17.094]]
(3)平均相对模拟百分误差为(%)
46.9225
(4)未来[25]步的预测值分别为:
[14.908, 14.966, 14.984, 15.003, 15.025, 15.048, 15.071, 15.096, 15.121, 15.147, 15.174, 15.201, 15.228, 15.255, 15.283, 15.311, 15.339, 15.368, 15.397, 15.425, 15.454, 15.484, 15.513, 15.543, 15.572]