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

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 += 1

self.train = self._diff_smooth(self.train)
draw_ts(self.train)
複製程式碼

平滑後的訓練資料：

3.3 將訓練資料進行週期性分解

採用statsmodels工具包：

from statsmodels.tsa.seasonal import seasonal_decompose

decomposition = seasonal_decompose(self.ts, freq=freq, two_sided=False)
# self.ts:時間序列，series型別; 
# freq:週期，這裡為1440分鐘，即一天; 
# two_sided:觀察下圖2、4行圖，左邊空了一段，如果設為True，則會出現左右兩邊都空出來的情況，False保證序列在最後的時間也有資料，方便預測。

self.trend = decomposition.trend
self.seasonal = decomposition.seasonal
self.residual = decomposition.resid
decomposition.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。