這是當初剛進公司時，leader給的一個獨立練手小項目，關于時間序列預測，情景比較簡單，整個過程實現(xiàn)下來代碼也僅100多行，但完成過程中踩了很多坑，覺的有必要分(tu)享(cao)一下。完整代碼和樣例數(shù)據(jù)放到了我的github上（文章僅粘貼部分)：
https://github.com/scarlettgin/cyclical_series_predict1、背景
公司平臺上有不同的api，供內(nèi)部或外部調(diào)用，這些api承擔著不同的功能，如查詢賬號、發(fā)版、搶紅包等等。日志會記錄下每分鐘某api被訪問了多少次，即一個api每天會有1440條記錄（1440分鐘），將每天的數(shù)據(jù)連起來觀察，有點類似于股票走勢的意思。我想通過前N天的歷史數(shù)據(jù)預測出第N+1天的流量訪問情況，預測值即作為合理參考，供新一天與真實值做實時對比。當真實流量跟預測值有較大出入，則認為有異常訪問，觸發(fā)報警。
2、數(shù)據(jù)探索
我放了一份樣例數(shù)據(jù)在data文件夾下，
看一下數(shù)據(jù)大小和結構
data = pd.read_csv(filename)print('size: ',data.shape)print(data.head())
data.png
數(shù)據(jù)大?。?div style="height:15px;">
共10080條記錄，即10080分鐘，七天的數(shù)據(jù)。
字段含義：
date：時間，單位分鐘
count：該分鐘該api被訪問的次數(shù)
畫圖看一下序列的走勢:（一些畫圖等探索類的方法放在了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)
序列.png
看這糟心的圖，那些驟降為0的點這就是我遇到的第一個坑，我當初一拿到這份數(shù)據(jù)就開始做了。后來折騰了好久才發(fā)現(xiàn)，那些驟降為0的點是由于數(shù)據(jù)缺失，ETL的同學自動補零造成的，溝通晚了(TДT)。
把坑填上，用前后值的均值把缺失值補上，再看一眼：
填充好缺失值的序列.png
發(fā)現(xiàn)這份數(shù)據(jù)有這樣幾個特點，在模型設計和數(shù)據(jù)預處理的時候要考慮到：
1、這是一個周期性的時間序列，數(shù)值有規(guī)律的以天為周期上下波動，圖中這個api，在每天下午和晚上訪問較為活躍，在早上和凌晨較為稀少。在建模之前需要做分解。
2、我的第二個坑：數(shù)據(jù)本身并不平滑，驟突驟降較多，而這樣是不利于預測的，畢竟模型需要學習好正常的序列才能對未知數(shù)據(jù)給出客觀判斷，否則會出現(xiàn)頻繁的誤報，令氣氛變得十分尷尬( ′Д｀)，所以必須進行平滑處理。
3、這只是一個api的序列圖，而不同的api的形態(tài)差距是很大的，畢竟承擔的功能不同，如何使模型適應不同形態(tài)的api也是需要考慮的問題。
3、預處理
3.1 劃分訓練測試集
前六天的數(shù)據(jù)做訓練，第七天做測試集。
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 對訓練數(shù)據(jù)進行平滑處理
消除數(shù)據(jù)的毛刺，可以用移動平均法，我這里沒有采用，因為我試過發(fā)現(xiàn)對于我的數(shù)據(jù)來說，移動平均處理完后并不能使數(shù)據(jù)平滑，我這里采用的方法很簡單，但效果還不錯：把每個點與上一點的變化值作為一個新的序列，對這里邊的異常值，也就是變化比較離譜的值剃掉，用前后數(shù)據(jù)的均值填充，注意可能會連續(xù)出現(xiàn)變化較大的點：
def _diff_smooth(self, ts): dif = ts.diff().dropna() # 差分序列 td = dif.describe() # 描述性統(tǒng)計得到：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 # 發(fā)現(xiàn)連續(xù)多少個點變化幅度過大，大部分只有單個點 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)
平滑后的訓練數(shù)據(jù)：
平滑后的訓練序列.png
3.3 將訓練數(shù)據(jù)進行周期性分解
采用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，則會出現(xiàn)左右兩邊都空出來的情況，F(xiàn)alse保證序列在最后的時間也有數(shù)據(jù)，方便預測。self.trend = decomposition.trendself.seasonal = decomposition.seasonalself.residual = decomposition.residdecomposition.plot()plt.show()
分解圖.png
第一行observed：原始數(shù)據(jù)；第二行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的訓練參數(shù)order =（p,d,q），具體意義查看官方文檔，調(diào)參過程略。 self.trend_model = ARIMA(train, order).fit(disp=-1, method='css')4.2 預測
預測出趨勢數(shù)據(jù)后，加上周期數(shù)據(jù)即作為最終的預測結果，但更重要的是，我們要得到的不是具體的值，而是一個合理區(qū)間，當真實數(shù)據(jù)超過了這個區(qū)間，則觸發(fā)報警，誤差高低區(qū)間的設定來自剛剛分解出來的殘差residual數(shù)據(jù)：
d = self.residual.describe()delta = d['75%'] - d['25%']self.low_error, self.high_error = (d['25%'] - 1 * delta, d['75%'] + 1 * delta)
預測并完成最后的加法處理，得到第七天的預測值即高低置信區(qū)間：
def predict_new(self): ''' 預測新數(shù)據(jù) ''' #續(xù)接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): ''' 為預測出的趨勢數(shù)據(jù)添加周期數(shù)據(jù)和殘差數(shù)據(jù) ''' 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] # 相同時間點的周期數(shù)據(jù)均值 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模型的參數(shù)order md.predict_new() pred = md.final_pred test = md.test plt.subplot(211) plt.plot(md.ts) # 平滑過的訓練數(shù)據(jù)加未做處理的測試數(shù)據(jù) 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') # 低置信區(qū)間 md.high_conf.plot(color='grey', label='high') # 高置信區(qū)間 plt.legend(loc='best') plt.title('RMSE: %.4f' % np.sqrt(sum((pred.values - test.values) ** 2) / test.size)) plt.tight_layout() plt.show()
預測結果.png
可以看到，均方根誤差462.8，相對于原始數(shù)據(jù)幾千的量級，還是可以的。測試數(shù)據(jù)中的兩個突變的點，也超過了置信區(qū)間，能準確報出來。
5、結語
前文提到不同的api形態(tài)差異巨大，本文只展示了一個，我在該項目中還接觸了其他形態(tài)的序列，有的有明顯的上升或下降趨勢；有的開始比較平緩，后面開始增長... ... ,但是都屬于典型的周期性時間序列，它的核心思想很簡單：做好分解，做好預測結果的還原，和置信區(qū)間的設置，具體操作可根據(jù)具體業(yè)務邏輯做調(diào)整，祝大家建模愉快:-D。
