python版多变量灰色预测

淩亂°似流年 2022-05-08 14:14 277阅读 0赞

应同学邀请，运行程序生成多变量灰色预测结果。因为很久没用过MATLAB，对它不熟，遂花了点时间，把该程序转成python版的。

原版MATLAB程序（省略部分注释）：

clc,clear all;
    format long g;
    fprintf(2, 'OBGM(1,N)模型: \n');
    
    X0=input('请输入因变量序列(第一个分号前面部分)与自变量序列(多个自变量序列之间用分号隔开)，如：[1 2 3;4 5 6]\n:')
    X00=input('请输入预测时的自变量序列(多个自变量序列之间用分号隔开)，如：[911 12 13]\n:')
    % 书上的数据
    %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];
    
    % epsilon:
    epsilon=0.5;
    [n,m]=size(X0);
    F = zeros(n, m);
    for u=1:n
        ago=0;
        for v=1:m
            ago=ago+X0(u,v);
            F(u,v)=ago;
        end
    end
    
    % Calculate mode parameters
    Y=X0';
    Y=Y(:,1);
    Y(1,:)=[];
    
    B=F';
    B(1,:)=[];B(:,1)=[];
    Bm=[];Li=[];Co=[];
    for k=1:m-1
        Bm=[Bm,(epsilon-1)*F(1,k)-epsilon*F(1,k+1)];
        Li=[Li,k+1];
        Co=[Co,1];
    end
    B=[B,Bm'];
    B=[B,Li'];
    B=[B,Co'];
    disp('B:');
    disp(B);
    
    % Calculate b2,b3,...
    P=(B'*B)^(-1)*B'*Y;
    disp('B.T*B:');
    disp(B'*B);
    disp('B.T*B**-1:');
    disp((B'*B)^(-1));
    disp('_*B.T:');
    disp((B'*B)^(-1)*B');
    disp('P:');
    disp(P);
    
    in_var_number=n-1;
    a = P(in_var_number+1);
    c = P(in_var_number+2);
    d = P(in_var_number+3);
    
    theta_1 = 1/(1+epsilon*a);
    theta_2 = 1-a/(1+epsilon*a);
    theta_3 = c/(1+epsilon*a);
    theta_4 = d/(1+epsilon*a);
    
    % Calcute Simu1
    Simu1 = [];
    x01=X0(1,1);
    Simu1=[Simu1,x01];
    for k=2:m
        s1=0;s2=0;
        for u=1:k-1
            for i=2:n
                s1=s1+theta_1*power(theta_2,u-1)*P(i-1)*F(i,k-u+1);
            end
        end
        for v=0:k-2
            s2=s2+power(theta_2,v)*(k-v)*theta_3+theta_4;
        end
        s3=power(theta_2,k-1)*x01;
        Simu1=[Simu1,s1+s2+s3];
    end
    
    % Calculate Fore0
    Xf=X0;
    Xf(1,:)=[];
    X0_F=[Xf,X00];
    
    [g, h] =size(X0_F);
    
    for u=1:g
        ago=0;
        for v=1:h
            ago=ago+X0_F(u,v);
            X0_F_1(u,v)=ago;
        end
    end
    
    Fore1=[];
    for k=m+1:h
        s1=0;s2=0;
        for u=1:k-1
            for i=2:n
                s1=s1+theta_1*power(theta_2,u-1)*P(i-1)*X0_F_1(i-1,k-u+1);
            end
        end
        for v=0:k-2
            s2=s2+power(theta_2,v)*((k-v)*theta_3+theta_4);
        end
        s3=power(theta_2,k-1)*x01;
        Fore1=[Fore1,s1+s2+s3];
    end
    Fore1=[F(1,m),Fore1];
    Fore0=[];
    for fo=1:length(Fore1)-1
        ss=Fore1(fo+1)-Fore1(fo);
        Fore0=[Fore0,roundn(ss,-3)];
    end
    
    % A;
    A=zeros(n,5);
    A(1,1)=1;
    A(1,2)=X0(1,1);
    A(1,3)=X0(1,1);
    
    % mp:
    mp=0;
    disp('length X0');
    disp(length(X0));
    disp(' X0 size:');
    disp(size(X0));
    for k=2:length(X0)
        A(k,1)=k;
        A(k,2)=X0(1,k);
        A(k,3)=roundn(Simu1(k)-Simu1(k-1),-3);
        A(k,4)=roundn(A(k,3)-A(k,2),-3);
        A(k,5)=roundn(100*abs(A(k,4))/A(k,2),-3);
        mp=mp+100*abs(A(k,4))/A(k,2);
    end
    mp=roundn(mp/(length(X0)-1),-4);
    
    disp('---------------------------------------');
    disp('(1)模型参数b2,b3,..,bn及a,c,d分别为:');
    disp(roundn(P,-6));
    disp('(2)误差检验表为:');
    disp('   序号    实际数据    模拟数据    残差    相对误差(%)');
    disp(A);
    
    disp('(3)平均相对模拟百分误差为(%)');
    disp(mp);
    
    str1='(4)未来[';
    str2=int2str(length(Fore0));
    str3=']步的预测值分别为:';
    str=[str1, str2, str3];
    disp(str);
    disp(Fore0);

直接转化为python版（需要有numpy库，可通过pip install numpy安装）：

#!/usr/bin/env python3
    import math
    import numpy as np
    
    print("OBGM(1, N)模型:")
    X0 = input('请输入因变量序列:')
    X00 = input('请输入预测时的自变量序列:')
    # 书上的数据
    #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 = np.zeros((n, m))
    for u in range(n):
        ago = 0
        for v in range(m):
            ago += X0[u, v]
            F[u, v] = ago
    
    # Calculate mode parameters
    Y = np.transpose(X0[0][1:])
    B = np.delete(F.T[1:], 0, axis=1)
    
    Bm = []
    Li = []
    Co = []
    for k in range(m - 1):
        Bm.append((epsilon - 1) * F[0, k] - epsilon * F[0, k + 1])
        Li.append(k + 2)
        Co.append(1)
    for L in (Bm, Li, Co):
        B = np.c_[B, np.transpose(np.array(L))]
    
    # Calculate b2,b3,...
    P = np.linalg.inv(B.T.dot(B)).dot(B.T).dot(Y)
    in_var_number = n - 1
    a = P[in_var_number]
    c = P[in_var_number + 1]
    d = P[in_var_number + 2]
    
    _t = 1 + epsilon * a
    theta_1 = 1 / _t
    theta_2 = 1 - a / _t
    theta_3 = c / _t
    theta_4 = d / _t
    
    # Calculate Simu
    Simu1 = [X0[0, 0]]
    for k in range(1, m):
        s1 = s2 = s3 = 0
        for u in range(k - 1):
            for i in range(1, n):
                s1 += theta_1 * math.pow(theta_2, u) * P[i - 1] * F[i, k - u]
        for v in range(k - 1):
            s2 += math.pow(theta_2, v) * ((k - v + 1) * theta_3 + theta_4)
        s3 = math.pow(theta_2, k) * Simu1[0]
        Simu1.append(s1 + s2 + s3)
    
    # % Calcute g
    Xf = X0[1:]
    X0_F = np.c_[Xf, np.array(X00)]
    X0_F_1 = np.zeros(X0_F.shape)
    g, h = X0_F.shape
    for u in range(g):
        ago = 0
        for v in range(h):
            ago += X0_F[u, v]
            X0_F_1[u, v] = ago
    
    Fore1 = [F[0, m - 1]]
    for k in range(m, h):
        s1 = s2 = s3 = 0
        for u in range(k - 1):
            for i in range(1, n):
                s1 += (
                    theta_1
                    * math.pow(theta_2, u)
                    * P[i - 1]
                    * X0_F_1[i - 1, k - u]
                )
        for v in range(k - 1):
            s2 += math.pow(theta_2, v) * ((k - v + 1) * theta_3 + theta_4)
        s3 = math.pow(theta_2, k) * X0[0, 0]
        Fore1.append(s1 + s2 + s3)
    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("(2)误差检验表为:")
    print("   序号    实际数据    模拟数据    残差    相对误差(%)")
    print(A)
    
    print("(3)平均相对模拟百分误差为(%)")
    print(mp)
    
    print(f"(4)未来[{len(Fore0)]步的预测值分别为:")
    print(Fore0)

优化后：

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