灰色預測法 —— python

目錄

1.簡介

2.演算法詳解

2.1 生成累加資料

2.2 累加后的資料運算式

2.3 求解2.2的未知引數

3.實體分析

3.1 匯入資料

3.2 進行累加資料

3.3 求解系數

3.4 預測資料及對比

完整代碼

1.簡介

灰色系統理論認為對既含有已知資訊又含有未知或非確定資訊的系統進行預測，就是對在一定方位內變化的、與時間有關的灰色程序的預測，盡管程序中所顯示的現象是隨機的、雜亂無章的，但畢竟是有序的、有界的，因此這一資料集合具備潛在的規律，灰色預測就是利用這種規律建立灰色模型對灰色系統進行預測，

灰色預測通過鑒別系統因素之間發展趨勢的相異程度，即進行關聯分析，并對原始資料進行生成處理來尋找系統變動的規律，生成有較強規律性的資料序列，然后建立相應的微分方程模型，從而預測事物未來發展趨勢的狀況，其用等時距觀測到的反應預測物件特征的一系列數量值構造灰色預測模型，預測未來某一時刻的特征量，或達到某一特征量的時間，

2.演算法詳解

2.1 生成累加資料

2.2 累加后的資料運算式

2.3 求解2.2的未知引數

3.實體分析

現有1997—2002年各項指標相關統計資料如下表：

用灰色預測方法預測2003—2009年各項指標的資料，且已知實際的預測資料如下：將預測資料與實際資料進行比較

3.1 匯入資料

#原資料
data=np.array([[72.03,241.2,1592.74],[73.84,241.2,1855.36],[74.49,244.8,2129.60],[76.68,250.9,2486.86],[78.00,250.9,2728.94],[79.68,252.2,3038.90]])
#要預測資料的真實值
data_T=np.array([[81.21,256.5,3458.05],[82.84,259.4,3900.27],[84.5,262.4,4399.06],[86.19,265.3,4961.62],[87.92,268.3,5596.1],[89.69,271.4, 6311.79],[91.49,274.5,7118.96]])

回傳結果，請自行列印查看

3.2 進行累加資料

#累加資料
data1=np.cumsum(data.T,1) #按列相加
print(data1)

回傳：

3.3 求解系數

[m,n]=data1.shape #得到行數和列數 m=3,n=6
#對這三列分別進行預測
X=[i for i in range(1997,2003)]#已知年份資料
X=np.array(X)
X_p=[i for i in range(2003,2010)]#預測年份資料
X_p=np.array(X_p)
X_sta=X[0]-1#最開始參考資料
#求解未知數
for j in range(3):
    B=np.zeros((n-1,2))
    for i in range(n-1):
        B[i,0]=-1/2*(data1[j,i]+data1[j,i+1])
        B[i,1]=1
    Y=data.T[j,1:7]
    a_u=np.dot(np.dot(np.linalg.inv(np.dot(B.T,B)),B.T),Y.T)
    print(a_u)
    #進行資料預測
    a=a_u[0]
    u=a_u[1]

回傳：

得到3對 a和u

3.4 預測資料及對比

需在3.3的基礎上進行預測

[m,n]=data1.shape #得到行數和列數 m=3,n=6
#對這三列分別進行預測
X=[i for i in range(1997,2003)]#已知年份資料
X=np.array(X)
X_p=[i for i in range(2003,2010)]#預測年份資料
X_p=np.array(X_p)
X_sta=X[0]-1#最開始參考資料
#求解未知數
for j in range(3):
    B=np.zeros((n-1,2))
    for i in range(n-1):
        B[i,0]=-1/2*(data1[j,i]+data1[j,i+1])
        B[i,1]=1
    Y=data.T[j,1:7]
    a_u=np.dot(np.dot(np.linalg.inv(np.dot(B.T,B)),B.T),Y.T)
#     print(a_u)
    #進行資料預測
    a=a_u[0]
    u=a_u[1]
    T=[i for i in range(1997,2010)]
    T=np.array(T)
    data_p=(data1[0,j]-u/a)*np.exp(-a*(T-X_sta-1))+u/a #累加資料
#     print(data_p)
    data_p1=data_p
    data_p1[1:len(data_p)]=data_p1[1:len(data_p)]-data_p1[0:len(data_p)-1]
#     print(data_p1)
    title_str=['第一產業GDP預測','居民消費價格指數預測','第三產業GDP預測']
    plt.subplot(221+j)
    data_n=data_p1
    plt.scatter(range(1997,2003),data[:,j])
    plt.plot(range(1997,2003),data_n[X-X_sta])
    plt.scatter(range(2003,2010),data_T[:,j])
    plt. plot(range(2003,2010),data_n[X_p-X_sta-1])
#     plt.title(title_str[j])
    plt.legend(['實際原資料','擬合資料','預測參考資料','預測資料'])
    y_n=data_n[X_p-X_sta-1].T
    y=data_T[:,j]
    wucha=sum(abs(y_n-y)/y)/len(y)
    titlestr1=[title_str[j],'預測相對誤差：',wucha]
    plt.title(titlestr1)
    plt.show()

回傳：

完整代碼

import numpy as np
import matplotlib.pyplot as plt
import math

# 解決圖示題中文亂碼問題
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默認字體
mpl.rcParams['axes.unicode_minus'] = False  # 解決保存影像是負號'-'顯示為方塊的問題

#原資料
data=np.array([[72.03,241.2,1592.74],[73.84,241.2,1855.36],[74.49,244.8,2129.60],[76.68,250.9,2486.86],[78.00,250.9,2728.94],[79.68,252.2,3038.90]])
#要預測資料的真實值
data_T=np.array([[81.21,256.5,3458.05],[82.84,259.4,3900.27],[84.5,262.4,4399.06],[86.19,265.3,4961.62],[87.92,268.3,5596.1],[89.69,271.4, 6311.79],[91.49,274.5,7118.96]])

#累加資料
data1=np.cumsum(data.T,1)
print(data1)

[m,n]=data1.shape #得到行數和列數 m=3,n=6
#對這三列分別進行預測
X=[i for i in range(1997,2003)]#已知年份資料
X=np.array(X)
X_p=[i for i in range(2003,2010)]#預測年份資料
X_p=np.array(X_p)
X_sta=X[0]-1#最開始參考資料
#求解未知數
for j in range(3):
    B=np.zeros((n-1,2))
    for i in range(n-1):
        B[i,0]=-1/2*(data1[j,i]+data1[j,i+1])
        B[i,1]=1
    Y=data.T[j,1:7]
    a_u=np.dot(np.dot(np.linalg.inv(np.dot(B.T,B)),B.T),Y.T)
#     print(a_u)
    #進行資料預測
    a=a_u[0]
    u=a_u[1]
    T=[i for i in range(1997,2010)]
    T=np.array(T)
    data_p=(data1[0,j]-u/a)*np.exp(-a*(T-X_sta-1))+u/a #累加資料
#     print(data_p)
    data_p1=data_p
    data_p1[1:len(data_p)]=data_p1[1:len(data_p)]-data_p1[0:len(data_p)-1]
#     print(data_p1)
    title_str=['第一產業GDP預測','居民消費價格指數預測','第三產業GDP預測']
    plt.subplot(221+j)
    data_n=data_p1
    plt.scatter(range(1997,2003),data[:,j])
    plt.plot(range(1997,2003),data_n[X-X_sta])
    plt.scatter(range(2003,2010),data_T[:,j])
    plt. plot(range(2003,2010),data_n[X_p-X_sta-1])
#     plt.title(title_str[j])
    plt.legend(['實際原資料','擬合資料','預測參考資料','預測資料'])
    y_n=data_n[X_p-X_sta-1].T
    y=data_T[:,j]
    wucha=sum(abs(y_n-y)/y)/len(y)
    titlestr1=[title_str[j],'預測相對誤差：',wucha]
    plt.title(titlestr1)
    plt.show()