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

最新推荐文章于 2024-08-22 19:54:12 发布

Tine Chan

最新推荐文章于 2024-08-22 19:54:12 发布

阅读量7.6k

点赞数 6

本文链接：https://blog.csdn.net/chanleoo/article/details/101170956

版权

%建立符号变量a(发展系数)和b(灰作用量)
syms a b;
c = [a b]';

%原始数列 A
A = [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 = length(A);

%对原始数列 A 做累加得到数列 B
B = cumsum(A);

%对数列 B 做紧邻均值生成
for i = 2:n
    C(i) = (B(i) + B(i - 1))/2; 
end
C(1) = [];

%构造数据矩阵 
B = [-C;ones(1,n-1)];
Y = A; Y(1) = []; Y = Y';

%使用最小二乘法计算参数 a(发展系数)和b(灰作用量)
c = inv(B*B')*B*Y;
c = c';
a = c(1); b = c(2);

%预测后续数据
F = []; F(1) = A(1);
for i = 2:(n+25)
    F(i) = (A(1)-b/a)/exp(a*(i-1))+ b/a;
end

%对数列 F 累减还原,得到预测出的数据
G = []; G(1) = A(1);
for i = 2:(n+25)
    G(i) = F(i) - F(i-1); %得到预测出来的数据
end

disp('预测数据为：');
G

%模型检验

H = G(1:30);
%计算残差序列
epsilon = A - H;

%法一：相对残差Q检验
%计算相对误差序列
delta = abs(epsilon./A);
%计算相对误差Q
disp('相对残差Q检验：')
Q = mean(delta)

%法二：方差比C检验
disp('方差比C检验：')
C = std(epsilon, 1)/std(A, 1)

%法三：小误差概率P检验
S1 = std(A, 1);
tmp = find(abs(epsilon - mean(epsilon))< 0.6745 * S1);
disp('小误差概率P检验：')
P = length(tmp)/n

%绘制曲线图
t1 = 1989:2018;
t2 = 1989:2043;

plot(t1, A,'ro');
hold on;
plot(t2, G, 'g-');
xlabel('年份'); 
ylabel('全球冰浓度年平均值/L');
legend('实际全球冰浓度年平均值','预测全球冰浓度年平均值');
title('全球冰浓度年平均值增长曲线');
grid on;

预测数据为：

G =

列 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检验：

Q =

0.0019

方差比C检验：

C =

0.0561

小误差概率P检验：

P =

#!/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()

runfile('/home/rosfun/Downloads/PyTorch-YOLOv3-master/GM.py', wdir='/home/rosfun/Downloads/PyTorch-YOLOv3-master')
OBGM(1, N)模型:
请输入因变量序列(第一个分号前面部分)与自变量序列(多个自变量序列之间用分号隔开)，如：[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]
请输入预测时的自变量序列(多个自变量序列之间用分号隔开)，如：[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+00 1.255887e+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]