如何自动判断流量是否异常？| SARIMAX实战
一些监控指标，比如流量、连接数、请求数等，基本按天呈现周期规律，周期规律为1天，直观上符合季节时间序列，本次采用 SARIMAX 模型进行实战
SARIMAX 组成名词
AR(Auto Regressive Model)自回归模型
AR是线性时间序列分析模型中最简单的模型。通过自身前面部分的数据与后面部分的数据之间的相关关系（自相关）来建立回归方程，从而可以进行预测或者分析。下图中展示了一个时间如果可以表示成如下结构，那么就说明它服从p阶的自回归过程，表示为AR(p)。其中，\(u_t\)表示白噪声，是时间序列中的数值的随机波动，但是这些波动会相互抵消，最终是0。\(\phi\)表示自回归系数。
\(x_t = \phi_1x_{t-1} + \phi_2x_{t-2} + ... + \phi_px_{t-p} + u_t\)
所以当只有一个时间记录点时，称为一阶自回归过程，即AR(1)。
\(x_t = \phi_1x_{t-1} + u_t\)
MA(Moving Average Model)移动平均模型
通过将一段时间序列中白噪声序列进行加权和，可以得到移动平均方程。如下图所示为q阶移动平均过程，表示为MA(q)。\(\thetasym\)表示移动回归系数。\(u_t\)表示不同时间点的白噪声。
\(x_t = u_t + \thetasym_1u_{t-1} + \thetasym_2u_{t-2} + ... + \thetasym_qu_{t-q}\)
ARMA(Auto Regressive and Moving Average Model)自回归移动平均模型
自回归移动平均模型是与自回归和移动平均模型两部分组成。所以可以表示为ARMA(p, q)。p是自回归阶数，q是移动平均阶数。
\(x_t = u_t + \phi_1x_{t-1} + \phi_2x_{t-2} + ... + \phi_px_{t-p} + \thetasym_1u_{t-1} + \thetasym_2u_{t-2} + ... + \thetasym_qu_{t-q} \)
从式子中就可以看出，自回归模型结合了两个模型的特点，其中，AR可以解决当前数据与后期数据之间的关系，MA则可以解决随机变动也就是噪声的问题。
ARIMA(Auto Regressive Integrate Moving Average Model)差分自回归移动平均模型
同前面的三种模型，ARIMA模型也是基于平稳的时间序列的或者差分化后是稳定的，另外前面的几种模型都可以看作ARIMA的某种特殊形式。表示为ARIMA(p, d, q)。p为自回归阶数，q为移动平均阶数，d为时间成为平稳时所做的差分次数，也就是Integrate单词的在这里的意思。
SARIMA 季节性差分自回归移动平均模型
季节性自回归整合移动平均线，SARIMA 或季节性 ARIMA，是 ARIMA 的扩展，明确支持具有季节性成分的单变量时间序列数据。
它增加了三个新的超参数来指定系列季节性成分的自回归（AR），差分（I）和移动平均（MA），以及季节性周期的附加参数。
实战(集群请求数预测)
数据清洗
- 获取数据转化格式
- 先df绘图，查看异常点，然后移除
- 周期太大计算过慢，可以压缩返回
def read_json():
indexs = []
count = []
with open("./cn-hb-1-10m.json", 'r') as fd:
data = json.load(fd)
for item in data["aggregations"]["dates"]["buckets"]:
ts = int(item['key'] / 1000)
indexs.append(pd.datetime.fromtimestamp(ts))
count.append(item['doc_count'])
df = pd.DataFrame(index=indexs, data=count, columns=['count'])
df = df.fillna(df.bfill())
df2 = df[df['count'] <= 200000]
df3 = df2.resample(rule='1H').mean()
return df3
数据可视化与平稳检测
ADF 直接检测
def pre_check(df):
decompostion = seasonal_decompose(df, period=24)
decompostion.plot()
plt.show()
dftest = adfuller(df, autolag='AIC')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic', 'p-value', '#Lags Used', 'Number of Observations Used'])
for key, value in dftest[4].items():
dfoutput['Critical Value (%s)' % key] = value
# 查看 p-value ，看数据是否平稳
print(dfoutput)
参考诊断绘图如下
执行获取结果
Test Statistic -4.806309
p-value 0.000053
#Lags Used 15.000000
Number of Observations Used 177.000000
Critical Value (1%) -3.467845
Critical Value (5%) -2.878012
Critical Value (10%) -2.575551
dtype: float64
如果小于0.05，则不需要进行序列平稳处理，否则参考如下操作，进行序列化平稳操作
序列平稳化
消除时间序列的趋势效应和季节性效应最常用方法之一是差分法。对原始的序列作1阶12步差分来提取原序列的趋势效应和季节效应，差分后的序列图及检验结果如下所示。
#消除消除趋势和季节性
mte_first_difference = mte - mte.shift(1) #一阶差分
mte_seasonal_first_difference = mte_first_difference - mte_first_difference.shift(12) #12步差分
# 备注
df.diff(1) = df - df.shift(1)
最终目的就是 p-value 小于 0.05，实现时间序列平稳
建立SARIMA模型
图解法识别
绘制差分后序列的自相关和偏自相关图，通过观察自相关和偏自相关图找出模型的最优参数。通过Python绘制的差分后序列的自相关和偏自相关图如下所示。
def autocorrelation(timeseries, lags):
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(211)
sm.graphics.tsa.plot_acf(timeseries, lags=lags, ax=ax1)
ax2 = fig.add_subplot(212)
sm.graphics.tsa.plot_pacf(timeseries, lags=lags, ax=ax2)
plt.show()
看起来有点晕，只做引申，有兴趣可自行实战。
网格搜寻
用图解法求解ARIMA模型的最优参数并非易事，主观性很大，而且耗费时间。所以本文进一步考虑利用网格搜索的方法系统地选择最优的参数值。
网格搜索可以遍历探索参数的不同组合。对于每个参数组合，可以利用Python中statsmodels模块的SARIMAX()函数拟合一个新的季节性ARIMA模型，并评估其整体质量。当网格搜索遍历完整个参数环境是，我们可以依据评价时间序列模型的准则从参数集中选出最佳性能的参数。
在评估和比较具有不同参数的统计模型时，可以根据数据的拟合程度或准确预测未来数据点的能力对每个模型进行排序
AIC
AIC信息准则即Akaike information criterion，是衡量统计模型拟合优良性(Goodness of fit)的一种标准，由于它为日本统计学家赤池弘次创立和发展的，因此又称赤池信息量准则。它建立在熵的概念基础上，可以权衡所估计模型的复杂度和此模型拟合数据的优良性。
AIC = 2k - 2ln(L)
- k为参数数量
- L是似然函数
增加自由参数的数目提高了拟合的优良性，AIC鼓励数据拟合的优良性但是尽量避免出现过拟合的情况。所以优先考虑的模型应是AIC值最小的那一个，假设在n个模型中作出选择，可一次算出n个模型的AIC值，并找出最小AIC值对应的模型作为选择对象。
一般而言，当模型复杂度提高（k）增大时，似然函数L也会增大，从而使AIC变小，但是k过大时，似然函数增速减缓，导致AIC增大，模型过于复杂容易造成过拟合现象。
BIC
BIC= Bayesian Information Criterions，贝叶斯信息准则。
BIC=ln(n)k-2ln(L)
- L是似然函数
- n是样本大小
- K是参数数量
pk
BIC的惩罚项比AIC大，考虑了样本个数，样本数量多，可以防止模型精度过高造成的模型复杂度过高。
AIC和BIC前半部分是一样的，BIC考虑了样本数量，样本数量过多时，可有效防止模型精度过高造成的模型复杂度过高。
所以我们实战选择 BIC
代码实战
def p_q_choice(df):
p = d = q = range(0, 2)
pdq = list(itertools.product(p, d, q))
seasonal_pdq = [(x[0], x[1], x[2], 24) for x in list(itertools.product(p, d, q))]
warnings.filterwarnings("ignore")
pre_result = float('inf')
pre_param = None
pre_param_seasonal = None
for param in pdq:
for param_seasonal in seasonal_pdq:
try:
mod = sm.tsa.statespace.SARIMAX(df,
order=param,
seasonal_order=param_seasonal,
enforce_stationarity=False,
enforce_invertibility=False)
results = mod.fit()
if results.bic < pre_result:
pre_result = results.bic
pre_param = param
pre_param_seasonal = param_seasonal
# print('ARIMA{}x{} - BIC:{}'.format(param, param_seasonal, results.aic))
except:
continue
print('ARIMA{}x{} - BIC:{}'.format(pre_param, pre_param_seasonal, pre_result))
return pre_param, pre_param_seasonal
网络搜索不易选择过大值，满足需求即可，可以减少时间消耗，参考当前样本获取模型如下
CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
ARIMA(0, 1, 1)x(0, 1, 1, 24)12 - BIC:2905.7000757306228
模型诊断 - 可视化
def train(df, param, reasonal_param):
mod = sm.tsa.statespace.SARIMAX(df,
order=param,
seasonal_order=reasonal_param,
enforce_stationarity=False,
enforce_invertibility=False)
results = mod.fit()
# 初期使用，后期注釋掉
print(results.summary())
results.plot_diagnostics(figsize=(10,5))
plt.show()
return results
模型训练完，统计结果如下
SARIMAX Results
==========================================================================================
Dep. Variable: count No. Observations: 193
Model: SARIMAX(0, 1, 1)x(0, 1, 1, 24) Log Likelihood -1445.416
Date: Tue, 19 Oct 2021 AIC 2896.833
Time: 20:50:12 BIC 2905.700
Sample: 05-17-2021 HQIC 2900.436
- 05-25-2021
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ma.L1 -0.3530 0.054 -6.512 0.000 -0.459 -0.247
ma.S.L24 -0.4468 0.044 -10.069 0.000 -0.534 -0.360
sigma2 4.696e+07 1.94e-10 2.42e+17 0.000 4.7e+07 4.7e+07
===================================================================================
Ljung-Box (L1) (Q): 0.16 Jarque-Bera (JB): 171.85
Prob(Q): 0.69 Prob(JB): 0.00
Heteroskedasticity (H): 0.23 Skew: 0.31
Prob(H) (two-sided): 0.00 Kurtosis: 8.35
===================================================================================
由上述结果可知，coef列显示每个变量的权重(即重要性)，以及每个变量如何影响时间序列。P>|z|列是对每个变量系数的检验。每个变量的P值均小于0.01，所以在0.01的显著性水平下，拒绝加假设，模型中每个变量的系数通过显著性检验，可以认为拟合的模型中包含这些变量是合理的。
模型诊断展示如下，
残差序列延迟1-12阶时，Q统计量的P值均大于0.05，所以在0.05的显著性水平下，不拒绝原假设，即残差为白噪声序列，说明残差序列的波动已经没有任何统计规律。因此，可以认为拟合的模型已经充分提取了时间序列中的信息。
的模型诊断结果可知，残差的的时序图基本稳定，残差随着时间的波动并没有很大的波动。
由右上角的正态分布图可知，模型的残差是正态分布的。红色的KDE线紧随着N(0，1)线。其中，N(0,1)是均值为0，标准差为1的标准正态分布。这很好地说明了残差是正态分布的。而图中左下角的正太Q—Q图也说明了残差服从正态分布。（合理的是这样的，不过此次正态分布有点歪，可能不太行）
图中右下角残差的自相关图显示残差不存在自相关，说明残差序列是白噪声序列。
上述是很专业的判定，感觉图不是很顺畅，SARIMAX(0，1，1)x(0，1，1，24)模型的拟合效果可能有点问题，下面用r2来诊断下，辅助判定。
验证预测 (Validating the Forecasts)
def test(y, results, steps):
# future
# Get forecast 20% steps ahead in future
train_length = len(results.data.dates)
y = y[:train_length + steps]
pred_uc = results.get_forecast(steps=steps)
# Get confidence intervals of forecasts
pred_ci = pred_uc.conf_int(alpha=0.2)
y_forecasted = pred_uc.predicted_mean
y_truth = y.iloc[-steps:]
r2 = r2_score(y_truth, y_forecasted)
ax = y.plot(label='观测值')
y_forecasted.plot(ax=ax, label="动态预测", alpha=.7)
ax.fill_between(pred_ci.index,
pred_ci.iloc[:, 0],
pred_ci.iloc[:, 1], color='k', alpha=.2)
ax.set_xlabel('时间')
ax.set_ylabel('GET数量')
plt.title("SARIMAX 预测相关系数r2 <%0.4f>" % (r2))
plt.legend()
plt.show()
主要看相关系数，以及在观测数据是否都在可信区间内，都在，说明可以用来区间来诊断观测值是否异常，用来作为预警。
一般大于90%认为有效，说明此次样本并没有成功，搞了这么久，心理不服，硬着头皮继续往下走。
模型有效期
使用动态预测，与测试数据进行校验，可以发现时间步数的预测的相关系数并不相同。也就是说预测准确率与步数有关系，或者说模型在一定步数内有效。
def predict_step(y, results, max_steps):
r2_list = []
train_length = len(results.data.dates)
step_range = range(int(max_steps*0.1), max_steps + 1)
step_length = len(step_range)
for step in step_range:
pred_uc = results.get_forecast(steps=step)
y_truth = y.iloc[train_length:train_length+step]
y_forecasted = pred_uc.predicted_mean
r2 = r2_score(y_truth, y_forecasted)
r2_list.append(round(r2*100, 4))
df = pd.DataFrame(data=r2_list, index=step_range)
ax = df.plot()
ax.fill_between(df.index,
[90]*step_length,
[100]*step_length, color='k', alpha=.2)
plt.show()
残酷的事实胜于雄辩，说明r2还是很有参考意义的，这下死心了。
如何实现业务指标画像
预测流量翻车了，我思考下能够影响流量的因素很多，比如有时是个大文件，有时是个小文件，但是访问的次数应该是变化不大的。是否可以用其他指标代替，这次数据换为请求数，方法不变，我们直接看最后总结、验证和预测。
右上角的正态分布似乎更集中了一点，感觉比流量靠谱些，接下来继续看
这次相关系数符合，那我们再看看预测曲线
预测30步，相关性依然很高，说明请求数更稳定，更容易预测。
总结
一些指标可能受很多因素影响，预测难度很大，导致很难落地。那么可以选择一些相关性的指标进行预测，原始指标做阈值报警，新指标进行模型预测，由此间接实现对业务的相关画像。
