《Python数据分析与挖掘实战》第11章——时间序列
本文是基于《Python数据分析与挖掘实战》的实战部分的第11章的数据——《应用系统负载分析与磁盘容量预测》做的分析。
旨在补充原文中的细节代码,并给出文中涉及到的内容的完整代码;
在作者所给代码的基础上增加的内容包括:
1)数据探索时画C盘/D盘已使用空间的时序图,并根据自相关和偏相关图判定平稳性,确定了所用模型是采用ARMA或者ARIMA,而不是AR或者MA;
2)模型构建构建基于ARIMA或者ARMA的模型,采用AIC/BIC/HQ信息准则对模型进行定阶,确定p,q参数,从而选择最优模型;
1 背景与目标分析
根据历史磁盘数据,采用时间序列分析法,来预测应用系统服务器磁盘已经使用空间的大小;为管理员提供定制化的预警提示。
实质:时间序列---回归
2 数据探索
2.1 数据特征分析
inputfile1 = 'discdata.xls'
data = pd.read_excel(inputfile1)
data.head()
d = data[(data['ENTITY']=='C:\\') & (data['TARGET_ID']==184 )]
import matplotlib as mpl
import matplotlib.pyplot as plt
plt.rc('figure',figsize=(9,7))
import datetime
import matplotlib.dates as mdates
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus']= False
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.set_title(u"C盘已使用空间的时序图")
# ax.set_xlabel(u'日期')
ax.set(xlabel=u'日期',ylabel=u'磁盘使用大小')
# 图上时间间隔显示为10天
ax.xaxis.set_major_locator(mdates.DayLocator(bymonthday=range(1,32), interval=10))
ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y-%m-%d"))
plt.subplots_adjust(bottom=0.13,top=0.95)
ax.plot(d['COLLECTTIME'],d['VALUE'],'ro-',)
fig.autofmt_xdate() #自动根据标签长度进行旋转
'''for label in ax.xaxis.get_ticklabels(): #此语句完成功能同上
label.set_rotation(45)
'''
plt.savefig('c.jpg')
plt.show()
同理,绘制的D盘的时序图如下:
‘
通过图中可以发现,磁盘的使用情况都不具有周期性,表现出缓慢性增长,呈现趋势性。因此,可以初步确认数据是非平稳的’
3 数据预处理
3.1 数据清洗
一般情况下默认磁盘容量是定值,所以需要剔除磁盘容量重复的数据
data.drop_duplicates(data.columns[:-1],inplace=True)
data.to_excel('dataCleaned.xlsx')
3.2 数据变换——属性构造
思路:由于每台服务器上的这三个属性值一直不变:NAME、TARGET_ID、ENTITY,将这三个属性值合并
inputfile2 = 'dataCleaned.xlsx'
data = pd.read_excel(inputfile2)
# 原书中方法一:
df = data[data['TARGET_ID'] == 184].copy() # 只选取TARGET_ID为184的数据
df_group = df.groupby('COLLECTTIME') # 以时间分组得到一个GroupBy对象
#定义属性变换函数
def attr_trans(x):
result = Series(index = ['SYS_NAME','CWXT_DB:184:C:\\','CWXT_DB:184:D:\\','COLLECTTIME'])
result['SYS_NAME'] = x['SYS_NAME'].iloc[0]
result['COLLECTTIME'] = x['COLLECTTIME'].iloc[0]
result['CWXT_DB:184:C:\\'] = x['VALUE'].iloc[0]
result['CWXT_DB:184:D:\\'] = x['VALUE'].iloc[1]
return result
data_attr_constr = df_group.apply(attr_trans)# 逐组处理
data_attr_constr.to_excel('attrsConstruction.xlsx',index=False)
data_attr_constr
部分结果如下:
#方法二,死方法,没有方法一灵活
df_g = df.groupby('COLLECTTIME').size()
indexpre = df_g.index
data_processed = DataFrame([],index = indexpre, columns=['SYS_NAME','CWXT_DB:184:C:\\','CWXT_DB:184:D:\\'])
data_processed['SYS_NAME'] = u'财务管理系统'
data_processed['CWXT_DB:184:C:\\'] = df['VALUE'][df['ENTITY']=='C:\\'].values
data_processed['CWXT_DB:184:D:\\'] = df['VALUE'][df['ENTITY']=='D:\\'].values
data_processed.head()
3.3 确定模型—— ARMA or ARIMA?
ARMA和ARIMA的区别?
网友回答:不管是ARMA还是ARIMA模型,都是对平稳数据建模。前者是直接针对平稳数据建模,无需进行差分变换;后者则需要先对数据进行差分,差分平稳后再建模。
由于ARIMA/ARMA 模型对时间序列的要求是平稳型,因此需要进行平稳性检验。(此部分参考的是某网友的博客,但是该网友的博客地址我忘记了,尝试寻找时也未找到,在此说声抱歉)
通过自相关和偏相关图判定平稳性,并确定所用模型
# 如果自相关是拖尾,偏相关截尾,则用 AR 算法
# 如果自相关截尾,偏相关拖尾,则用 MA 算法
# 如果自相关和偏相关都是拖尾,则用 ARMA 算法, ARIMA 是 ARMA 算法的扩展版,用法类似
# 查看是否序列平稳 自相关图法/时间序列的散点图
from scipy import stats
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.graphics.api import qqplot
dta = df['CWXT_DB:184:C:\\']
# 原数据的的自相关图与偏自相关图
fig = plt.figure(figsize=(12,12))
ax1=fig.add_subplot(411)# 自相关图
fig = sm.graphics.tsa.plot_acf(dta,lags=40,ax=ax1)
ax1.set_title(u'原始数据的自相关图')
ax2 = fig.add_subplot(412)# 篇自相关图
fig = sm.graphics.tsa.plot_pacf(dta,lags=40,ax=ax2)
ax2.set_title(u'原始数据的偏相关图')
# 一阶差分后的自相关图与偏自相关图
dta= dta.diff(1).dropna() # 注意一定要将查分后的空值去掉再取自相关系数
ax3=fig.add_subplot(413)# 自相关图
fig = sm.graphics.tsa.plot_acf(dta,lags=40,ax=ax3)
ax3.set_title(u'一阶差分后的自相关图')
ax4 = fig.add_subplot(414)# 篇自相关图
fig = sm.graphics.tsa.plot_pacf(dta,lags=40,ax=ax4)
ax4.set_title(u'一阶差分后的偏相关图')
plt.savefig('acf_pacf.jpg')
plt.show()
# 通过自相关图可以发现,该图自相关和偏相关都是拖尾,因此,确定是ARMA算法、ARIMA均可
4 建立模型
具体步骤在此就不重复,见下图
4.1 平稳性检验
接下来介绍几种平稳性检验方式(以C盘为例)
inputfile = 'attrsConstruction.xlsx'
#由于ARIMA 模型对时间序列的要求是平稳型,需要进行平稳性检验
data = pd.read_excel(inputfile)
df = data.iloc[:len(data)-5]# 不使用最后5个数
diff = 0
1)方法1:时间序列图法
具体见数据探索部分。通过图形可以看出该序列值不平稳(若平稳,会围绕某值上下晃动)
2)方法2:单位根检验(ADF)
# 方法:采用单位根检验(ADF)的方法或者时序图的方法(见数据探索模块)
from statsmodels.tsa.stattools import adfuller as ADF
判断D盘数据的平稳性,以及确定几次差分后平稳
x = ADF(df['CWXT_DB:184:C:\\'])
print x
3)方法3:游程检验
from statsmodels.sandbox.stats.runs import runstest_1samp as RSS
x = RSS(df['CWXT_DB:184:C:\\'])
print x
4)方法4: 自相关系数法
from statsmodels.tsa.stattools import acf as ACF
x = ACF(df['CWXT_DB:184:C:\\'])
print x
5)方法5:通过自相关和偏相关图判定
通过自相关图法/时间序列的散点图,观察平稳性,前面内容已经介绍,此处就不在赘述
4.2 构建模型
# 第 * 1 * 步--C盘---------平稳性检测
平稳性检测 :判断是否平稳,若不平稳,对其进行差分处理直至平稳
方法:采用单位根检验(ADF)的方法或者时序图的方法(见数据探索模块)
注意:其他平稳性检验方法见steadyCheck.py文件
from statsmodels.tsa.stattools import adfuller as ADF
diff = 0
# 判断D盘数据的平稳性,以及确定几次差分后平稳
adf = ADF(df['CWXT_DB:184:C:\\'])
print adf
while adf[1] >= 0.05 : # adf[1]是p值,p值小于0.05认为是平稳的
print adf[1]
diff = diff + 1
adf = ADF(df['CWXT_DB:184:C:\\'].diff(diff).dropna())#注意,差分后使用ADF检验时,必须去掉空值
print (u'原始序列经过%s阶差分后归于平稳,p值为%s') % (diff, adf[1])
df['CWXT_DB:184:C:\\_adf'] = df['CWXT_DB:184:C:\\'].diff(1)
# 第 * 2 * 步--C盘---------白噪声检验
目的:验证序列中有用信息是否已经被提取完毕,需要进行白噪声检验。若序列是白噪声序列,说明序列中有用信息已经被提取完,只剩随机扰动
方法:采用LB统计量的方法进行白噪声检验
若没有通过白噪声检验,则需要进行模型识别,识别其模型属于AR、MA还是ARMA。
inputfile2 = 'attrsConstruction.xlsx'
data1 = pd.read_excel(inputfile2)
data1 = data1.iloc[:len(data1)-5]# 不使用最后五个数据(作为预测参考)
# 白噪声检测
from statsmodels.stats.diagnostic import acorr_ljungbox
[[lb], [p]] = acorr_ljungbox(data1['CWXT_DB:184:C:\\'], lags = 1) ## lags是残差延迟个数
if p < 0.05:
print (u'原始序列为非白噪声序列,对应的p值为:%s' % p)
else:
print (u'原始序列为白噪声序列,对应的p值为:%s' % p)
[[lb], [p]] = acorr_ljungbox(data1['CWXT_DB:184:C:\\'].diff(1).dropna(), lags = 1)
if p < 0.05:
print (u'一阶差分序列为非白噪声序列,对应的p值为:%s' % p)
else:
print (u'一阶差分序列为白噪声序列,对应的p值为:%s' % p)
# 第 * 3 * 步----------模型识别
方法:
step1:采用极大似然比方法进行模型的参数估计,估计各个参数的值。
step2:然后针对各个不同模型,采用信息准则方法(有三种:BIC/AIC/HQ)对模型进行定阶,确定p,q参数,从而选择最优模型。
目前选择模型常用如下准则!!!!!
增加*参数的数目提高了拟合的优良性,AIC/BIC/HQ鼓励数据拟合的优良性但是尽量避免出现过度拟合(Overfitting)的情况。所以优先考虑的模型应是AIC/BIC/HQ值最小的那一个
# * AIC=-2 ln(L) + 2 k 中文名字:赤池信息量 akaike information criterion (AIC)
# * BIC=-2 ln(L) + ln(n)*k 中文名字:贝叶斯信息量 bayesian information criterion (BIC)
# * HQ=-2 ln(L) + ln(ln(n))*k hannan-quinn criterion (HQ)
step3:注意,进行此步时,index需要为时间序列类型
step4:确定最佳p、d、q的值
此处进行的是AIC方式定信息准则+ARMA
inputfile3 = 'attrsConstruction.xlsx'
data2 = pd.read_excel(inputfile3,index_col='COLLECTTIME')
xtest_value=data2['CWXT_DB:184:C:\\'][-5:]
data2 = data2.iloc[:len(data2)-5]# 不使用最后五个数据(作为预测参考)
xdata2 = data2['CWXT_DB:184:C:\\']
# ARIMA(p,d,q)中,AR是自回归,p为自回归项数;MA为滑动平均,q为滑动平均项数,d为使之成为平稳序列所做的差分次数(阶数),由前一步骤知d=1
# from statsmodels.tsa.arima_model import ARIMA#建立ARIMA(p,d,q)模型
from statsmodels.tsa.arima_model import ARMA #建立ARMA(p,q)模型
# AIC方式定信息准则 + ARMA --------!!!模型检验中也要对应修改!!!------------------------------
pmax = int(len(xdata2)/10) # 一般阶数不超过length/10
qmax = int(len(xdata2)/10) # 一般阶数不超过length/10
matrix = [] # aic矩阵
for p in range(pmax+1):
tmp = []
for q in range(qmax+1):
try:#存在部分为空值,会报错
# tmp.append(ARMA(xdata2, (p,q)).fit().bic) # BIC方式
tmp.append(ARMA(xdata2, (p,q)).fit().aic) # AIC方式
# tmp.append(ARMA(xdata2, (p,q)).fit().hq) # HQ方式
except:
tmp.append(None)
matrix.append(tmp)
matrix = pd.DataFrame(matrix) # 从中可以找出最小值
print matrix
print matrix.stack()
# 第 * 4 * 步--C盘---------模型检验
确定模型后,需要检验其残差序列是否是白噪声,若不是,说明,残差中还存在有用的信息,需要修改模型或者进一步提取。
若其残差不是白噪声,重新更换p,q的值,重新确定
import pandas as pd
import numpy as np
while 1:
p, q = matrix.stack().idxmin() # 先展平该表格,然后找出最小值的索引位置
print (u'当前AIC最小的p值与q值分别为:%s、%s' % (p,q))
lagnum = 12 # 残差延迟个数
arma = ARMA(xdata2, (p,q)).fit() # 建立并训练模型
xdata_pred = arma.predict() # 预测
pred_error = (xdata_pred - xdata2).dropna() # 计算残差
# 白噪声检测
from statsmodels.stats.diagnostic import acorr_ljungbox
lbx, px = acorr_ljungbox(pred_error, lags = lagnum)
h = (px < 0.05).sum() # p值小于0.05,认为是非噪声
if h > 0:
print (u'模型ARMA(%s,%s)不符合白噪声检验' % (p,q))
print ('在AIC矩阵中去掉[%s,%s]组合,重新进行计算' % (p,q))
matrix.iloc[p,q] = np.nan
arimafail = arma
continue
else:
print (p,q)
print (u'模型ARMA(%s,%s)符合白噪声检验' % (p,q))
break
arma.summary() # 当p,q值为0,0时,summary方法报错
predictdata = pd.DataFrame(xtest_value)
predictdata.insert(1,'CWXT_DB:184:C:\\_predict',forecast_values)
predictdata.rename(columns={'CWXT_DB:184:C:\\':u'实际值','CWXT_DB:184:C:\_predict':u'预测值'},inplace=True)
result_d = predictdata.applymap(lambda x: '%.2f' % x) # 将表格中各个浮点值都格式化
result_d.to_excel('pedictdata_C_AIC_ARMA.xlsx')
result_d
# 第 * 5 * 步--D盘---------模型评价
为了评价时序预测模型效果的好坏,本章采用3个衡量模型预测精度的统计量指标:平均绝对误差、均方根误差、平均绝对百分误差
# -*- coding:utf-8 -*-
import pandas as pd
inputfile4 = 'pedictdata_C_AIC_ARMA.xlsx'
result = pd.read_excel(inputfile4,index_col='COLLECTTIME')
result = result.applymap(lambda x: x/10**6)
print result
# 计算误差
abs_ = (result[u'预测值']-result[u'实际值']).abs()
mae_ = abs_.mean() # mae平均绝对误差
rmas_ = ((abs_**2).mean())**0.5 #rmas均方根误差
mape_ = (abs_/result[u'实际值']).mean() #mape平均绝对百分误差
# print abs_
print mae_
print rmas_
print mape_
errors = 1.5
print '误差阈值为%s' % errors
if (mae_ < errors) & (rmas_ < errors) & (mape_ < errors):
print (u'平均绝对误差为:%.4f, \n均方根误差为:%.4f, \n平均绝对百分误差为:%.4f' % (mae_, rmas_, mape_))
print '误差检验通过!'
else:
print '误差检验不通过!'
其他信息准则方式+ARMA/ARIMA的搭配就不赘述了,最终结果如下:
注意:
说明:由于用HQ训练模型时,都是空值,所以,本例使用HQ不合适-----ARIMA--BIC---
实际值 预测值
COLLECTTIME
2014-11-12 35.704313 35.722538
2014-11-13 35.704981 35.757104
2014-11-14 34.570385 35.791669
2014-11-15 34.673821 35.826235
2014-11-16 34.793245 35.860800
0.70232013
0.890203752645
0.0202432790493
误差阈值为1.5
BIC模型下平均绝对误差为:0.7023,
均方根误差为:0.8902,
平均绝对百分误差为:0.0202
误差检验通过!
-----ARIMA--AIC---
实际值 预测值
COLLECTTIME
2014-11-12 35.704313 35.779972
2014-11-13 35.704981 35.836938
2014-11-14 34.570385 35.889601
2014-11-15 34.673821 35.935428
2014-11-16 34.793245 35.981256
0.795290026
0.976369605661
0.0229009946085
AIC模型下平均绝对误差为:0.7953,
均方根误差为:0.9764,
平均绝对百分误差为:0.0229
误差检验通过!
通过对比AIC与BIC的结果,可以发现BIC的几个误差均较小
-----ARMA--BIC---
实际值 预测值
COLLECTTIME
2014-11-12 35.704313 35.581706
2014-11-13 35.704981 35.488223
2014-11-14 34.570385 35.405986
2014-11-15 34.673821 35.333641
2014-11-16 34.793245 35.270000
0.462308002
0.533460826783
0.0132815193493
误差阈值为1.5
平均绝对误差为:0.4623,
均方根误差为:0.5335,
平均绝对百分误差为:0.0133
误差检验通过!
综上:ARMA+BIC更优
5 预测结果可视化
inputfile = 'attrsConstruction.xlsx'
data = pd.read_excel(inputfile)
df = data.iloc[:len(data)-5]
inputfile1 = 'pedictdata_C.xlsx'# 预测值
result = pd.read_excel(inputfile1)
inputfile2 = 'pedictdata_D.xlsx'# 预测值
result1 = pd.read_excel(inputfile2)
import matplotlib as mpl
import matplotlib.pyplot as plt
plt.rc('figure',figsize=(9,9))
import datetime
import matplotlib.dates as mdates
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus']= False
fig = plt.figure()
fig.set(alpha=0.2)#设置图标透明度
ax = fig.add_subplot(2,1,1)
ax.set_title(u"C盘空间时序预测图")
ax.set(xlabel=u'日期',ylabel=u'磁盘使用大小')
# 图上时间间隔显示为10天
ax.xaxis.set_major_locator(mdates.DayLocator(bymonthday=range(1,32), interval=7))
ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y-%m-%d"))
plt.subplots_adjust(bottom=0.13,top=0.95)
ax.plot(df['COLLECTTIME'],df['CWXT_DB:184:C:\\'],'ro--',)
ax.plot(result['COLLECTTIME'],result[u'实际值'],'g+--',)
ax.plot(result['COLLECTTIME'],result[u'预测值'],'b*-',)
ax.grid(axis='y',linestyle='--')
ax.legend()
fig.autofmt_xdate() #自动根据标签长度进行旋转
'''for label in ax.xaxis.get_ticklabels(): #此语句完成功能同上
label.set_rotation(45)
'''
ax1 = fig.add_subplot(2,1,2)
ax1.set_title(u"D盘空间时序预测图")
# ax.set_xlabel(u'日期')
ax1.set(xlabel=u'日期',ylabel=u'磁盘使用大小')
# 图上时间间隔显示为10天
ax1.xaxis.set_major_locator(mdates.DayLocator(bymonthday=range(1,32), interval=7))
ax1.xaxis.set_major_formatter(mdates.DateFormatter("%Y-%m-%d"))
plt.subplots_adjust(bottom=0.13,top=0.95)
ax1.plot(df['COLLECTTIME'],df['CWXT_DB:184:D:\\'],'co--',)
ax1.plot(result1['COLLECTTIME'],result1[u'实际值'],'m+--',)
ax1.plot(result1['COLLECTTIME'],result1[u'预测值'],'y*-',)
ax1.grid(axis='y',linestyle='--')
ax1.legend()
fig.autofmt_xdate() #自动根据标签长度进行旋转
'''for label in ax.xaxis.get_ticklabels(): #此语句完成功能同上
label.set_rotation(45)
'''
plt.savefig('data_predict_pic.jpg')
plt.show()
结束语:以上所有内容即为此章全部内容。
备注:本章节完整代码详见点击打开链接
上一篇: matplotlib.pyplot库介绍
下一篇: 百度网盘卸载不掉问题