本文共 4407 字,大约阅读时间需要 14 分钟。
传统的统计学的方法:从最初的随机游走模型(RW)、历史均值(HA)、马尔科夫模型、时间序列模型和卡尔曼滤波模型。RW和HA依赖与理论假设,并未考虑交通流的波动性,以致预测结果与现实存在很大差异;而马尔科夫模型、时间序列模型和卡尔曼滤波模型则根据现有道路的历史交通流数据假定交通流符合某种概率分布,从而进行训练,估计出模型参数。今天我们介绍最经典的统计学算法——自回归滑动平均模型(ARMA)。
大家都知道,统计学处理数据,对数据的要求极为严格,需要在做分析之前,对数据进行假设检验,参数估计等等,ARMA模型也不例外。需要对时间序列的随机性和平稳性进行检验,根据检验的结果,可将序列分为不同的类型:
这里我们主要记住时间序列的任何一个时刻的序列值都是一个随机变量,然后求他们的均值和方差。
关于平稳性检验,画图是较为简单的方法,一般情况下,先画个图再再觉得要不要进行下一步的检验。
ARMA模型可以分为两个部分,AR模型和MA模型部分。
AR模型:在t时刻的随机变量 X ( t ) X_{(t)} X(t)的取值 x ( t ) x_{(t)} x(t)是关于过去p期 x ( t − 1 ) x_{(t-1)} x(t−1), x ( t − 2 ) x_{(t-2)} x(t−2), x ( t − 3 ) x_{(t-3)} x(t−3),… 多元线性回归。认为 x ( t ) x_{(t)} x(t)主要受到过去p时间的影响。误差项 ε ( t ) \varepsilon_{(t)} ε(t)是当期的随机扰动。
MA模型:在t时刻的随机变量 X ( t ) X_{(t)} X(t)的取值 x ( t ) x_{(t)} x(t)是关于过去q期随机扰动项 ε ( t − 1 ) \varepsilon_{(t-1)} ε(t−1), ε ( t − 2 ) \varepsilon_{(t-2)} ε(t−2), ε ( t − 3 ) \varepsilon_{(t-3)} ε(t−3),… 多元线性回归。认为 x ( t ) x_{(t)} x(t)主要受到过去q时间随机扰动的影响。 μ \mu μ是 X t X_{t} Xt的均值。(有限阶的MA模型一定是平稳的)
在这里解释一下随机扰动:是指 X ( t ) X_{(t)} X(t)的波动不确定的部分,我们就用无穷阶的AR模型来表示。ARMA模型:在t时刻的随机变量 X ( t ) X_{(t)} X(t)的取值 x ( t ) x_{(t)} x(t)是关于过去p期 x ( t − 1 ) x_{(t-1)} x(t−1), x ( t − 2 ) x_{(t-2)} x(t−2), x ( t − 3 ) x_{(t-3)} x(t−3),…和过去q期随机扰动项 ε ( t − 1 ) \varepsilon_{(t-1)} ε(t−1), ε ( t − 2 ) \varepsilon_{(t-2)} ε(t−2), ε ( t − 3 ) \varepsilon_{(t-3)} ε(t−3),… 的多元线性回归。
其实我们获取的数据多情况下是非平稳的数据,因而学会对非平稳数据的分析才更加重要。我们把非平稳的时间序列的分析主要分为两类:
先介绍一下python中数理统计的库==StatModels ==这里面主要包含了统计学的一些计算方法。
import pandas as pdimport matplotlib.pyplot as pltfrom statsmodels.graphics.tsaplots import plot_acfdiscfile = 'C:\\Users\\Administrator\\Desktop\\python-code\\《Python数据分析与挖掘实战(第2版)》源数据和代码\\《Python数据分析与挖掘实战(第2版)》源数据和代码-各章节\\chapter5\\demo\\data\\arima_data.xls'#forecastnum = 5#读取文件data = pd.read_excel(discfile, index_col = u'日期')#补全中文plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号data.plot()plt.show()
# 平稳性检测# 自相关图from statsmodels.graphics.tsaplots import plot_acfplot_acf(data).show()# 平稳性检测from statsmodels.tsa.stattools import adfuller as ADFprint(u'原始序列的ADF检验结果为:', ADF(data[u'销量']))# 返回值依次为adf、pvalue、usedlag、nobs、critical values、icbest、regresults、resstore
3.一阶差分、平稳性和白噪声检验
# 差分后的结果D_data = data.diff().dropna()D_data.columns = [u'销量差分']D_data.plot() # 时序图plt.show()plot_acf(D_data).show() # 自相关图from statsmodels.graphics.tsaplots import plot_pacfplot_pacf(D_data).show() # 偏自相关图print(u'差分序列的ADF检验结果为:', ADF(D_data[u'销量差分'])) # 平稳性检测# 白噪声检验from statsmodels.stats.diagnostic import acorr_ljungboxprint(u'差分序列的白噪声检验结果为:', acorr_ljungbox(D_data, lags=1)) # 返回统计量和p值
一阶差分:感觉有点像平稳的序列
自相关图:
4.ARMA模型定阶
定阶就是选择p和q的范围,这里我们选用BIC检验,组合各种p和q,获取最小BIC值的p和q。# 定阶data[u'销量'] = data[u'销量'].astype(float) pmax = int(len(D_data)/10) # 一般阶数不超过length/10qmax = int(len(D_data)/10) # 一般阶数不超过length/10bic_matrix = [] # BIC矩阵for p in range(pmax+1): tmp = [] for q in range(qmax+1): try: # 存在部分报错,所以用try来跳过报错。 tmp.append(ARIMA(data, (p,1,q)).fit().bic) except: tmp.append(None) bic_matrix.append(tmp)bic_matrix = pd.DataFrame(bic_matrix) # 从中可以找出最小值p,q = bic_matrix.stack().idxmin() # 先用stack展平,然后用idxmin找出最小值位置。print(u'BIC最小的p值和q值为:%s、%s' %(p,q))
结果:
model = ARIMA(data, (p,1,q)).fit() # 建立ARIMA(0, 1, 1)模型print('模型报告为:\n', model.summary2())print('预测未来5天,其预测结果、标准误差、置信区间如下:\n', model.forecast(5))
粉红框框的是参数检验的值,然后利用模型对2015年1月1日到2015年2月6日的销售数据进行预测。
参考文献:
书籍:《python数据分析与挖掘实战》转载地址:http://gnpv.baihongyu.com/