完整程式碼和分析

實驗介紹

實驗背景

大型企業都會有自己的辦公自動化系統，應用系統在日常執行時，會對底層軟、硬體造成負荷，顯著影響應用系統性能。底層任何一種資源負載過大，都可能會引起應用系統性能下降甚至癱瘓。因此需要關注伺服器、資料庫、中介軟體和儲存裝置的執行狀態，及時瞭解當前應用系統的負載情況，以便提前預防，確保系統安全穩定執行。

實驗目的

• 針對歷史磁碟資料，採用時間序列分析方法，預測應用系統伺服器磁碟已使用空間大小。

• 根據使用者需求設定不同的預警等級，將預測值與容量值進行比較，對其結果進行預警判斷，為系統管理員提供定製化的預警提示。

分析方法與過程

由於儲存空間隨時間變化存在很強的關聯性，且歷史資料對未來的發展存在一定的影響。故本實驗採用時間序列分析法對磁碟已使用空間進行預測分析。

模型介紹

    本次實驗我們將使用時間序列分析法來資料進行構造。

    首先來認識一下ARIMA模型。ARIMA，差分自迴歸滑動平均模型，又稱求合自迴歸滑動平均模型。是時間序列預測分析方法之一。

    ARIMA(p,d,q)中，AR是“自迴歸”，p為自迴歸項數；MA是“滑動平均”，q為滑動平均項數，d是使之成為平穩序列所做的差分次數（階數）。



ARIMA模型運用的流程：

根據時間序列的散點圖、自相關函式和偏自相關函式圖識別其平穩性。

對非平穩的時間序列資料進行平穩化處理。直到處理後的自相關函式和偏自相關函式的數值非顯著非零。

根據所識別出來的特徵建立相應的時間序列模型。平穩化處理後，若偏自相關函式是截尾的，而自相關函式是拖尾的，則建立AR模型；若偏自相關函式是拖尾的，而自相關函式是截尾的，則建立MA模型；若偏自相關函式和自相關函式均是拖尾的，則序列適合ARMA模型。

利用已通過檢驗的模型進行預測。

那麼，如何識別時間序列的平穩性？

平穩性檢驗

識別序列是否平穩，一般有三種方法：看圖法、自相關係數和偏相關係數、單位根檢驗（ADF）。

我們先來簡單介紹一下看圖法：

看圖法

這裡的“圖”指的是時序圖，即某指隨時間變化的時序圖。平穩序列的圖a是圍繞一個常數上下波動。而不平穩的圖b，則是用明顯的增長或減少的趨勢。

自相關係數和偏相關係數

這裡會涉及到兩個定義——截尾和拖尾。

平穩序列的自相關圖和偏相關圖要麼都是拖尾要麼都是截尾。截尾就是在某階之後，係數都為 0 ，怎麼理解呢，看上面偏相關的圖，當階數為 1 的時候，係數值還是很大， 0.914. 二階長的時候突然就變成了 0.050. 後面的值都很小，認為是趨於 0 ，這種狀況就是截尾。再就是拖尾，拖尾就是有一個衰減的趨勢，但是不都為 0 。自相關圖既不是拖尾也不是截尾。以上的圖的自相關是一個三角對稱的形式，這種趨勢是單調趨勢的典型圖形。

c. 單位根檢驗（ADF） 若單位根檢驗p值小於0.05則認為是平穩的。

對於不平穩的序列又該如何處理呢？ 2. 差分處理 差分即取相鄰項值差替代當前值，以此來消除一些波動，使資料趨於平穩。

資料抽取

為了抽取出磁碟資料，以屬性的標識號（TARGET_ID）與採集指標的時間（COLLECTTIME）為條件，對效能資料進行抽取。本實驗抽取2014-10-1至2014-11-16財務管理系統中某一臺資料庫伺服器的磁碟的相關資料。

資料探索分析

本實驗採用時序分析法來進行建模。為了建模的需要，需要探索資料的平穩性。

通過時序圖可以初步發現數據的平穩性。針對伺服器磁碟的使用大小，以天為單位，進行週期性分析。

新建檔案，新增如下程式碼：

import pandas as pd

dataPath = './data/discdata.xls'
data = pd.read_excel(dataPath,encoding='utf-8')

#繪製C、D盤的使用情況時序圖
import matplotlib.pyplot as plt
import matplotlib

#配置matplotlib引數
#座標軸字型
matplotlib.rc('font', **{'family': 'serif', 'serif': ['SimHei']})
#plt.rcParams['front.sans-serif'] = ['SimHei']
#座標軸負號
#plt.rcParams['axes.unicode-minus'] = False

pd.to_datetime(data['COLLECTTIME'])
data1 = data[(data['ENTITY'] == 'C:\\') & (data['TARGET_ID'] == 184)]
#設定dataframe索引，修改dataframe，不建立新物件
data1.set_index('COLLECTTIME',inplace = True)

data2 = data[(data['ENTITY'] == 'D:\\') & (data['TARGET_ID'] == 184)]
data2.set_index('COLLECTTIME',inplace = True)

print(data.head())
print(data1.head())
print(data2.head())

執行後，可以看到我們將在這些資料上進行視覺化操作：

在原始檔中新增如下程式碼：

plt.plot(data1.index, data1['VALUE'], 'ro-')
plt.title(u"C盤已使用空間的時序圖")
plt.xlabel(u'日期')
plt.ylabel(u'磁碟使用大小')
#刻度位置
plt.xticks(rotation = 30)
plt.show()

執行後可以看到C盤容量的時序圖：

繼續新增如下程式碼：

plt.plot(data2.index,data2['VALUE'],'ko-')
plt.title(u'D盤已使用空間的時序圖')
plt.xlabel(u'日期')
plt.ylabel(u'磁碟使用大小')
plt.xticks(rotation = 30)
plt.show()

執行後可以看到D盤容量的時序圖：

根據上圖可以得知，磁碟的使用情況都不具備有周期性，它們表現出緩慢性增長，呈現趨勢線。

初步判斷，資料是非平穩的。

資料清洗

在實際業務中，監控系統會每天定時對磁碟的資訊進行收集。但是一般情況下，磁碟的容量屬性是一個定值，因此磁碟原始資料中會存在磁碟容量的重複資料。

在資料清洗過程中，剔除磁碟容量的重複資料，並且將所有伺服器的磁碟容量作為一個固定值，方便模型預警。

開啟檔案，新增如下程式碼：

data.drop_duplicates(data.columns[:-1],inplace=True)

將原始資料中重複資料剔除。

屬性構造

在資料儲存中，磁碟容量以KB為單位。因為每臺伺服器的磁碟資訊可以通過表中NAME、TARGET_ID、ENTITY三個屬性進行區分，且每臺伺服器的上述三個屬性是不變的，所以可以將這個三個屬性值合併，構造新的屬性。

開啟檔案，新增如下程式碼：

data3 = pd.DataFrame(index = data1.index,columns=['SYS_NAME','CWXT_DB:184:C:\\','CWXT_DB:184:D:\\','COLLECTTIME'])
data3['SYS_NAME'] = data1['SYS_NAME']
data3['CWXT_DB:184:C:\\'] = data1['VALUE']
data3['CWXT_DB:184:D:\\'] = data2['VALUE']
data3['COLLECTTIME'] = data1.index
data3.to_excel('./data/data_processed.xls')
print(data3.head())

執行後可以看到，屬性構造後的資料內容為：

容量預測模型

將處理後的資料分為兩部分，一部分是建模樣本資料，一部分是模型驗證資料。選取建模資料的最後5條資料為驗證資料，其他資料為建模樣本資料。

（1）我們首先要做的是對觀測值進行平穩性檢驗，如果不平穩需要進行差分運算直到差分後的資料平穩

（2）在資料平穩後進行白噪聲檢驗，如果沒有通過白噪聲檢驗，就進行模型識別，識別模型屬於AR,MA和ARMA中的哪一種模型，並通過貝葉斯資訊準則對模型進行定階，確定ARIMA模型的p、q引數 。

（3）在模型識別後需要進行模型檢驗，檢測模型殘差序列是否為白噪聲序列，如果沒有通過檢驗，需要對其進行重新識別，對已經通過檢驗的模型採用極大似然估計方法進行模型引數估計。

（4）應用模型進行預測，將實際值與預測值進行誤差分析。如果誤差比較小，表明模型擬合效果較好，則模型可以結束，反之需要重新估計引數。

平穩性檢驗

為了確定原始資料序列中沒有隨機趨勢或確定趨勢，需要對資料進行平穩性檢驗，否則將會產生“偽迴歸”的現象。本節採用ADF方法來進行平穩性檢驗。

需要注意的是，本次實驗我們只使用C盤容量來進行模型的構建和預測。學生可自行使用D盤容量進行模型構建和預測。

#-*- coding:utf-8 -*-
import pandas as pd
#ADF檢驗
#引數初始化
discfile = './data/data_processed.xls'

data4 = pd.read_excel(discfile)
#不使用最後5個數據
data = data4.iloc[:len(data)-5]

#平穩性測試
from statsmodels.tsa.stattools import adfuller as ADF
diff = 0 
adf = ADF(data['CWXT_DB:184:C:\\'])
#print(adf)
#adf[1]為p值，p值小於0.05認為是平穩的
while adf[1] >= 0.05:
    diff = diff + 1
    adf = ADF(data['CWXT_DB:184:C:\\'].diff(diff).dropna())
    #print(adf)
    
print(u'原始序列經過%s階差分後歸於平穩，p值為%s'%(diff,adf[1]))

執行後可以看到：

此時將d值確定為1。

白噪聲檢驗

為了驗證序列中有用的資訊是否已被提取完畢，需要對序列進行白噪聲檢驗。如果序列檢驗為白噪聲序列，就說明序列中有用的資訊已經被提取完畢了，剩下的全是隨機擾動，無法進行預測和使用。本實驗採用LB統計量的方法進行白噪聲檢驗。

#白噪聲檢驗
#LB統計量
from statsmodels.stats.diagnostic import acorr_ljungbox

[[lb],[p]] = acorr_ljungbox(data['CWXT_DB:184:C:\\'], lags=1)
if p < 0.05:
    print(u'原始序列為非白噪聲序列，對應的p值為：%s'%p)
else:
    print(u'原始序列為白噪聲序列，對應的p值為：%s'%p)
    
[[lb],[p]] = acorr_ljungbox(data['CWXT_DB:184:C:\\'].diff(1).dropna(),lags=1)
if p < 0.05:
    print(u'一階差分序列為非白噪聲序列，對應的p值為：%s'%p)
else:
    print(u'一階差分為白噪聲序列，對應的p值為：%s'%p)

執行後可以看到：

同樣需要一階差分將原始資料的非平穩性消除。

模型識別

本節採用極大似然法進行模型的引數估計，估計各個引數的值。針對各個不同模型，採用BIC資訊準則對模型進行定階，確定p，q引數，從而選擇最優模型。

目前我們已經確定ARIMA模型中d=1。

開啟model.py檔案，新增如下程式碼：

#模型識別
#確定最佳p、d、q值
#xdata = data['CWXT_DB:184:D:\\']
xdata = data['CWXT_DB:184:C:\\']
from statsmodels.tsa.arima_model import ARIMA

#定階
pmax = int(len(xdata)/10) #一般階數不超過length/10
qmax = int(len(xdata)/10)
bic_matrix = [] #bic矩陣
for p in range(pmax+1):
    tmp = []
    for q in range(qmax+1):
        try:
            tmp.append(ARIMA(xdata,(p,1,q)).fit().bic)
        except:
            tmp.append(None)
    bic_matrix.append(tmp)
bic_matrix = pd.DataFrame(bic_matrix) #取值區域
#stack()將資料from columns to indexs
p,q = bic_matrix.stack().astype('float64').idxmin()
print(u'BIC最小的p值和q值為:%s、%s'%(p,q))
#D:1,1
#C:0,0

執行後，可以看到針對C盤容量的最優模型為ARIMA(0,1,0)。 針對D盤容量的最優模型為ARIMA(1,1,1)。

模型檢驗

模型確定後，檢驗其殘差序列是否為白噪聲。如果不是白噪聲，說明殘差中還存在有用資訊，需要修改模型或者進一步提取。

#模型檢驗
lagnum = 12
from statsmodels.tsa.arima_model import ARIMA

arima = ARIMA(xdata,(p,1,q)).fit()
xdata_pred = arima.predict(typ = 'levels')#predict
#print(xdata_pred)
pred_error = (xdata_pred - xdata).dropna()#殘差

from statsmodels.stats.diagnostic import acorr_ljungbox

lb,p_l = acorr_ljungbox(pred_error, lags = lagnum)
h = (p_l < 0.05).sum()#p值小於0.05，認為是非白噪聲
if h > 0:
    print(u'模型ARIMA（%s,1,%s）不符合白噪聲檢驗'%(p,q))
else:
    print(u'模型ARIMA（%s,1,%s）符合白噪聲檢驗'%(p,q))

執行後，針對C盤容量的ARIMA(0,1,0)符合白噪聲檢驗。

模型預測

應用通過檢驗的模型進行預測，獲取未來5天的預測值，並且與實際值比較。

我們在建模時未使用最後5個數據，我們用這些資料來做預測驗證。並將單位換算為GB。

#模型預測
#forecast向前預測5個值
test_predict = arima.forecast(5)[0]

#預測對比
test_data = pd.DataFrame(columns = [u'實際容量',u'預測容量'])
test_data[u'實際容量'] = data4[(len(data4)-5):]['CWXT_DB:184:C:\\']
test_data[u'預測容量'] = test_predict
test_data = test_data.applymap(lambda x :'%.2f'%x)
print(test_data)

執行後，可以看到預測值：

模型評價

至此，我們的模型已經構建完成。為了評價時序預測模型效果的好壞，本實驗採用3個衡量模型預測精度的統計量指標：平均絕對誤差、均方根誤差和平均絕對百分誤差。這3個指標從不同側面反映了演算法的預測精度。

結合實際業務分析，將誤差閾值設為1.5。

#計算誤差
#列操作
test_data[u'預測容量'] = test_data[u'預測容量'].astype(float)
test_data[u'實際容量'] = test_data[u'實際容量'].astype(float)
#10**6單位換算
abs_ = (test_data[u'預測容量'] - test_data[u'實際容量']).abs()/10**6
mae_ = abs_.mean()
rmse_ = ((abs_**2).mean())**0.05
mape_ = (abs_/test_data[u'實際容量']).mean()

print(u'平均絕對誤差為：%0.4f,\n均方根誤差為：%0.4f,\n平均絕對百分誤差為：%0.6f。'%(mae_,rmse_,mape_))

執行後可以看到：

實際值與預測值之間的誤差全都小於誤差閾值。因此模型的預測效果在實際業務可接受的範圍內，可以採用此模型進行預測。

參考：