ARIMA模型構建、預測——基於Python

《服務器系統負載分析及磁盤容量預測》，附帶代碼的學習、注釋：

 

從該問題的分析思路看（有問題找方案）：建立磁盤容量使用的預警系統（避免宕機等）——>（問題背景：總容量大小基本不變，使用量根據負載情況變化）預測出某時刻的使用量——>預測使用量占比是否達到預警系統閾值——>根據閾值輸出判斷信號

從給出的數據結構及實際情況，這是tsa問題（Time Series Analysis），statsmodels模塊，ARIMA不二之選：先考慮序列平穩性(不平穩差分，確定I)/白噪聲——>從而確定可以使用ARIMA，根據BIC、AIC准則取其值最小的p,q——>最后得到ARIMA(p,I,q),模型殘差為白噪聲序列表明模型棒棒噠——>即可開始以該模型為武器預測、對比、判斷~

幾個有意思的函數摘出來留意下：

 

##數據變結構，由寬數據變長數據，即把列索引轉行的內層索引
DataFrame.stack(level=-1, dropna=True)#level用來選列索引的內外層，默認全部搬到行索引，dropna為True則當轉后的行均NaN時去除改行
p,q=bic_matrix.stack().idxmin()   #將行拆散成列，形成雙行索引結構，外層主索引為原行索引內層為列索引轉化而來，轉化過程默認去除NaN和None，idxmin取列值最小行索引
##
ARIMAResults.predict(start=None, end=None, exog=None, typ='linear', dynamic=False)
#dynamic，邏輯值，True表樣本外預測，默認False樣本預測，typ，取值'linear', 'levels'表示根據內生變量的差分做線性預測，預測原數據的水平（源數據的模型預測值）
arima.predict('2014-11-12','2014-11-16',dynamic=True,typ='levels') #預測之后5天的值

文末也給出了具體預警系統預警的判斷流程（具體和想的一致）

預警級別和預測的使用占用率：Ⅰ：>=85%，Ⅱ：>=90，Ⅲ：>=95%。

預警系統業務運行思路：

• 此ARIMA(p,I,q)模型來構建預警系統，主要系由歷史數據（每日定時采集）運行該模型（主要是P,I,Q的確定），進而預測未來5天磁盤使用量
• 由預測的磁盤使用量占總容量比率來划定預警級別，發布預警級別信息
• 實際中考慮使用容量一般不會突變，每日變化不大，遂該模型調整（主要是P,I,Q的確定）可根據業務情況半月/一個月來修正

附代碼：

# -*- coding: utf-8 -*-
"""
Created on Thu Aug 23 20:00:31 2018

@author: Luove
"""

import pandas as pd
import numpy as np
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.tsa.stattools import adfuller as ADF
from statsmodels.stats.diagnostic import acorr_ljungbox

file1='../data/discdata.xls'
#file2='../data/discdata_processed.xls'
data=pd.read_excel(file1)

data=data[data['TARGET_ID']==184]
data_group=data.groupby('COLLECTTIME')
def attribute_turn(x):
    result=pd.Series(index=['SYS_NAME','CWXT_DB:184:c\\','CWXT_DB:184:d\\','COLLECTTIME'])
    result['SYS_NAME']=x['SYS_NAME'].iloc[0]
    result['CWXT_DB:184:c\\']=x['VALUE'].iloc[0]
    result['CWXT_DB:184:d\\']=x['VALUE'].iloc[1]
    result['COLLECTTIME']=x['COLLECTTIME'].iloc[0]
    return result
data_processed=data_group.apply(attribute_turn)

data_1=data_processed[:-5]
diff=0
adf=ADF(data_processed['CWXT_DB:184:d\\'])
while adf[1]>=0.05:   #h0:存在單位根，序列不平穩
    diff+=1
    adf=ADF(data_processed['CWXT_DB:184:d\\'].diff(diff).dropna())
    
print('原序列經過%d階差分平穩，此時p值為%f'%(diff,adf[1]))
[[lb],p]=acorr_ljungbox(data_1['CWXT_DB:184:d\\'],lags=1) #LB時間序列是否存在滯后相關的一種統計檢驗，選一階滯后
if p<0.05:  #h0序列是白噪聲序列
    print('原序列是非白噪聲序列，此時p值為%f'%p)
else :
    print('原序列是白噪聲序列，此時p值為%f'%p)
[[lb],p]=acorr_ljungbox(data_1['CWXT_DB:184:d\\'].diff().dropna(),lags=1) 
if p<0.05:  #h0序列是白噪聲序列
    print('一階差分序列是非白噪聲序列，此時p值為%f'%p)
else :
    print('一階差分序列是白噪聲序列，此時p值為%f'%p)

xdata=data_1['CWXT_DB:184:d\\']
pmax=int(len(xdata)/10)
qmax=int(len(xdata)/10)
bic_matrix=[]

for p in range(pmax+1):
    tmp=[]
    for q in range(qmax+1):  # 生成一個在p=p時，q=0,1,2,3,4時的ARIMA模型的BIC，q為之后的列索引
        try :
            tmp.append(ARIMA(xdata,(p,1,q)).fit().bic)
        except :
            tmp.append(None)
    bic_matrix.append(tmp)  #此時的行列索引分別對應p,q
bic_matrix=pd.DataFrame(bic_matrix)
p,q=bic_matrix.stack().idxmin()   #將行拆散成列，形成雙行索引結構，外層主索引為原行索引內層為列索引轉化而來，轉化過程默認去除NaN和None
print('BIC最小的p和q分別是%d和%d'%(p,q))  # 確定了p,i,q，模型設定完畢可以開始回歸和預測

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

[lb],p=acorr_ljungbox(pred_error,lags=1) #LB時間序列是否存在滯后相關的一種統計檢驗，選一階滯后
if p<0.05: # 即拒絕原假設，h0是白噪聲序列
    print('ARIMA(0,1,2)殘差序列不是白噪聲序列，p值為%f'%p)
else:
    print('ARIMA(0,1,2)殘差序列是白噪聲序列，p值為%f'%p)

#arima.predict('2014-11-12','2014-11-16',dynamic=True,typ='levels') #預測之后5天的值
#ARIMAResults.predict(start=None, end=None, exog=None, typ='linear', dynamic=False)#dynamic，邏輯值，True表樣本外預測，默認False樣本預測，typ，取值'linear', 'levels'表示根據內生變量的差分做線性預測，預測原數據的水平（源數據的模型預測值）
compare_data=pd.concat([data_processed[-5:]['CWXT_DB:184:d\\'],arima.predict('2014-11-12','2014-11-16',dynamic=True,typ='levels')],axis=1)#合並列
compare_data.columns=['real','pred']

abs_=(compare_data['real']-compare_data['pred']).abs()
mae=abs_/mean()#Mean Absolute Error ，平均絕對誤差
rmse=(abs_**2.mean())**0.5#Root Mean Square Error,均方根誤差
mape=(abs_/compare_data['real']).mean()# mean absolute percentage error，平均絕對百分比誤差

print('平均絕對誤差MAE=%f；\n均方根誤差RMSE=%f；\n平均絕對百分比誤差MAPE=%f。'%(mae,rmse,mape))

#至此，模型構建、計算、評價已完成；評價即誤差在誤差設定閾值之內，說明模型可以使用
#由此模型來構建預警系統，主要系由歷史數據（每日定時采集）運行該模型（主要是P,I,Q的確定），進而預測未來5天磁盤使用量
#由預測的磁盤使用量占總容量比率來划定預警級別，發布預警級別信息
#實際中考慮使用容量一般不會突變，每日變化不大，遂該模型調整（主要是P,I,Q的確定）可根據業務情況半月/一個月來修正

 

 

REF：

《數據分析與挖掘實戰》

源代碼及數據需要可自取：https://github.com/Luove/Data