用Python預測「周期性時間序列」的正確姿勢

這是當初剛進公司時，leader給的一個獨立練手小項目，關於時間序列預測，情景比較簡單，整個過程實現下來代碼也僅100多行，但完成過程中踩了很多坑，覺的有必要分(tu)享(cao)一下。完整代碼和樣例數據放到了我的github上（文章僅粘貼部分)：

https://github.com/scarlettgin/cyclicalseriespredict

1、背景

公司平台上有不同的api，供內部或外部調用，這些api承擔著不同的功能，如查詢賬號、發版、搶紅包等等。日誌會記錄下每分鐘某api被訪問了多少次，即一個api每天會有1440條記錄（1440分鐘），將每天的數據連起來觀察，有點類似於股票走勢的意思。我想通過前N天的歷史數據預測出第N+1天的流量訪問情況，預測值即作為合理參考，供新一天與真實值做實時對比。當真實流量跟預測值有較大出入，則認為有異常訪問，觸發報警。

2、數據探索

我放了一份樣例數據在data文件夾下，

看一下數據大小和結構

data = pd.read_csv(filename)print(size: ,data.shape)print(data.head())

數據大小：

共10080條記錄，即10080分鐘，七天的數據。

欄位含義：

date：時間，單位分鐘

count：該分鐘該api被訪問的次數

畫圖看一下序列的走勢:（一些畫圖等探索類的方法放在了test_stationarity.py 文件中，包含時間序列圖，移動平均圖，有興趣的可以自己嘗試下）。

def draw_ts(timeseries): timeseries.plot() plt.show()data = pd.read_csv(path)data = data.set_index(date)data.index = pd.to_datetime(data.index)ts = data[count]draw_ts(ts)

看這糟心的圖，那些驟降為0的點這就是我遇到的第一個坑，我當初一拿到這份數據就開始做了。後來折騰了好久才發現，那些驟降為0的點是由於數據缺失，ETL的同學自動補零造成的，溝通晚了(TДT)。

把坑填上，用前後值的均值把缺失值補上，再看一眼：

發現這份數據有這樣幾個特點，在模型設計和數據預處理的時候要考慮到：

1、這是一個周期性的時間序列，數值有規律的以天為周期上下波動，圖中這個api，在每天下午和晚上訪問較為活躍，在早上和凌晨較為稀少。在建模之前需要做分解。

2、我的第二個坑：數據本身並不平滑，驟突驟降較多，而這樣是不利於預測的，畢竟模型需要學習好正常的序列才能對未知數據給出客觀判斷，否則會出現頻繁的誤報，令氣氛變得十分尷尬( ′Д｀)，所以必須進行平滑處理。

3、這只是一個api的序列圖，而不同的api的形態差距是很大的，畢竟承擔的功能不同，如何使模型適應不同形態的api也是需要考慮的問題。

3、預處理

3.1 劃分訓練測試集

前六天的數據做訓練，第七天做測試集。

class ModelDecomp(object): def __init__(self, file, test_size=1440): self.ts = self.read_data(file) self.test_size = test_size self.train_size = len(self.ts) - self.test_size self.train = self.ts[:len(self.ts)-test_size] self.test = self.ts[-test_size:]

3.2 對訓練數據進行平滑處理

消除數據的毛刺，可以用移動平均法，我這裡沒有採用，因為我試過發現對於我的數據來說，移動平均處理完後並不能使數據平滑，我這裡採用的方法很簡單，但效果還不錯：把每個點與上一點的變化值作為一個新的序列，對這裡邊的異常值，也就是變化比較離譜的值剃掉，用前後數據的均值填充，注意可能會連續出現變化較大的點：

def _diff_smooth(self, ts): dif = ts.diff().dropna() # 差分序列 td = dif.describe() # 描述性統計得到：min，25%，50%，75%，max值 high = td[75%] + 1.5 * (td[75%] - td[25%]) # 定義高點閾值，1.5倍四分位距之外 low = td[25%] - 1.5 * (td[75%] - td[25%]) # 定義低點閾值，同上 # 變化幅度超過閾值的點的索引 forbid_index = dif[(dif > high) | (dif < low)].index i = 0 while i < len(forbid_index) - 1: n = 1 # 發現連續多少個點變化幅度過大，大部分只有單個點 start = forbid_index[i] # 異常點的起始索引 while forbid_index[i+n] == start + timedelta(minutes=n): n += 1 i += n - 1 end = forbid_index[i] # 異常點的結束索引 # 用前後值的中間值均勻填充 value = np.linspace(ts[start - timedelta(minutes=1)], ts[end + timedelta(minutes=1)], n) ts[start: end] = value i += 1self.train = self._diff_smooth(self.train)draw_ts(self.train)

平滑後的訓練數據：

3.3 將訓練數據進行周期性分解

採用statsmodels工具包：

from statsmodels.tsa.seasonal import seasonal_decomposedecomposition = seasonal_decompose(self.ts, freq=freq, two_sided=False)# self.ts:時間序列，series類型; # freq:周期，這裡為1440分鐘，即一天; # two_sided:觀察下圖2、4行圖，左邊空了一段，如果設為True，則會出現左右兩邊都空出來的情況，False保證序列在最後的時間也有數據，方便預測。self.trend = decomposition.trendself.seasonal = decomposition.seasonalself.residual = decomposition.residdecomposition.plot()plt.show()

第一行observed：原始數據；

第二行trend：分解出來的趨勢部分；

第三行seasonal：周期部分；最後residual：殘差部分。

我採用的是seasonal_decompose的加法模型進行的分解，即

observed = trend + seasonal + residual，另還有乘法模型。在建模的時候，只針對trend部分學習和預測，如何將trend的預測結果加工成合理的最終結果？當然是再做加法，後面會詳細寫。

4、模型

4.1 訓練

對分解出來的趨勢部分單獨用arima模型做訓練：

def trend_model(self, order): self.trend.dropna(inplace=True) train = self.trend[:len(self.trend)-self.test_size] #arima的訓練參數order =（p,d,q），具體意義查看官方文檔，調參過程略。 self.trend_model = ARIMA(train, order).fit(disp=-1, method=css)

4.2 預測

預測出趨勢數據後，加上周期數據即作為最終的預測結果，但更重要的是，我們要得到的不是具體的值，而是一個合理區間，當真實數據超過了這個區間，則觸發報警，誤差高低區間的設定來自剛剛分解出來的殘差residual數據：

d = self.residual.describe()delta = d[75%] - d[25%]self.low_error, self.high_error = (d[25%] - 1 * delta, d[75%] + 1 * delta)

預測並完成最後的加法處理，得到第七天的預測值即高低置信區間：

def predict_new(self): 預測新數據 #續接train，生成長度為n的時間索引，賦給預測序列 n = self.test_size self.pred_time_index= pd.date_range(start=self.train.index[-1], periods=n+1, freq=1min)[1:] self.trend_pred= self.trend_model.forecast(n)[0] self.add_season() def add_season(self): 為預測出的趨勢數據添加周期數據和殘差數據 self.train_season = self.seasonal[:self.train_size] values = [] low_conf_values = [] high_conf_values = [] for i, t in enumerate(self.pred_time_index): trend_part = self.trend_pred[i] # 相同時間點的周期數據均值 season_part = self.train_season[ self.train_season.index.time == t.time() ].mean() # 趨勢 + 周期 + 誤差界限 predict = trend_part + season_part low_bound = trend_part + season_part + self.low_error high_bound = trend_part + season_part + self.high_error values.append(predict) low_conf_values.append(low_bound) high_conf_values.append(high_bound) # 得到預測值，誤差上界和下界 self.final_pred = pd.Series(values, index=self.pred_time_index, name=predict) self.low_conf = pd.Series(low_conf_values, index=self.pred_time_index, name=low_conf) self.high_conf = pd.Series(high_conf_values, index=self.pred_time_index, name=high_conf)

4.3 評估：

對第七天作出預測，評估的指標為均方根誤差rmse，畫圖對比和真實值的差距：

md = ModelDecomp(file=filename, test_size=1440) md.decomp(freq=1440) md.trend_model(order=(1, 1, 3)) # arima模型的參數order md.predict_new() pred = md.final_pred test = md.test plt.subplot(211) plt.plot(md.ts) # 平滑過的訓練數據加未做處理的測試數據 plt.title(filename.split(.)[0]) plt.subplot(212) pred.plot(color=blue, label=Predict) # 預測值 test.plot(color=red, label=Original) # 真實值 md.low_conf.plot(color=grey, label=low) # 低置信區間 md.high_conf.plot(color=grey, label=high) # 高置信區間 plt.legend(loc=best) plt.title(RMSE: %.4f % np.sqrt(sum((pred.values - test.values) ** 2) / test.size)) plt.tight_layout() plt.show()

可以看到，均方根誤差462.8，相對於原始數據幾千的量級，還是可以的。測試數據中的兩個突變的點，也超過了置信區間，能準確報出來。

5、結語

前文提到不同的api形態差異巨大，本文只展示了一個，我在該項目中還接觸了其他形態的序列，有的有明顯的上升或下降趨勢；有的開始比較平緩，後面開始增長... ... ,但是都屬於典型的周期性時間序列，它的核心思想很簡單：做好分解，做好預測結果的還原，和置信區間的設置，具體操作可根據具體業務邏輯做調整，祝大家建模愉快:-D。