一些監控指標，比如流量、連接數、請求數等，基本按天呈現周期規律，周期規律為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步，相關性依然很高，說明請求數更穩定，更容易預測。
總結
一些指標可能受很多因素影響，預測難度很大，導致很難落地。那么可以選擇一些相關性的指標進行預測，原始指標做閾值報警，新指標進行模型預測，由此間接實現對業務的相關畫像。