這是當(dāng)初剛進(jìn)公司時(shí)，leader給的一個(gè)獨(dú)立練手小項(xiàng)目，關(guān)于時(shí)間序列預(yù)測，情景比較簡單，整個(gè)過程實(shí)現(xiàn)下來代碼也僅100多行，但完成過程中踩了很多坑，覺的有必要分(tu)享(cao)一下。完整代碼和樣例數(shù)據(jù)放到了我的github上（文章僅粘貼部分)：
https://github.com/scarlettgin/cyclical_series_predict

1、背景

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

2、數(shù)據(jù)探索

我放了一份樣例數(shù)據(jù)在data文件夾下，
看一下數(shù)據(jù)大小和結(jié)構(gòu)

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

data.png

數(shù)據(jù)大小：
共10080條記錄，即10080分鐘，七天的數(shù)據(jù)。
字段含義：
date：時(shí)間，單位分鐘
count：該分鐘該api被訪問的次數(shù)

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

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的點(diǎn)這就是我遇到的第一個(gè)坑，我當(dāng)初一拿到這份數(shù)據(jù)就開始做了。后來折騰了好久才發(fā)現(xiàn)，那些驟降為0的點(diǎn)是由于數(shù)據(jù)缺失，ETL的同學(xué)自動補(bǔ)零造成的，溝通晚了(TДT)。

把坑填上，用前后值的均值把缺失值補(bǔ)上，再看一眼：

填充好缺失值的序列.png

發(fā)現(xiàn)這份數(shù)據(jù)有這樣幾個(gè)特點(diǎn)，在模型設(shè)計(jì)和數(shù)據(jù)預(yù)處理的時(shí)候要考慮到：

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

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

3、這只是一個(gè)api的序列圖，而不同的api的形態(tài)差距是很大的，畢竟承擔(dān)的功能不同，如何使模型適應(yīng)不同形態(tài)的api也是需要考慮的問題。

3、預(yù)處理

3.1 劃分訓(xùn)練測試集

前六天的數(shù)據(jù)做訓(xùn)練，第七天做測試集。

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 對訓(xùn)練數(shù)據(jù)進(jìn)行平滑處理

消除數(shù)據(jù)的毛刺，可以用移動平均法，我這里沒有采用，因?yàn)槲以囘^發(fā)現(xiàn)對于我的數(shù)據(jù)來說，移動平均處理完后并不能使數(shù)據(jù)平滑，我這里采用的方法很簡單，但效果還不錯(cuò)：把每個(gè)點(diǎn)與上一點(diǎn)的變化值作為一個(gè)新的序列，對這里邊的異常值，也就是變化比較離譜的值剃掉，用前后數(shù)據(jù)的均值填充，注意可能會連續(xù)出現(xiàn)變化較大的點(diǎn)：

def _diff_smooth(self, ts):
    dif = ts.diff().dropna() # 差分序列
    td = dif.describe() # 描述性統(tǒng)計(jì)得到：min，25%，50%，75%，max值
    high = td['75%'] + 1.5 * (td['75%'] - td['25%']) # 定義高點(diǎn)閾值，1.5倍四分位距之外
    low = td['25%'] - 1.5 * (td['75%'] - td['25%']) # 定義低點(diǎn)閾值，同上

    # 變化幅度超過閾值的點(diǎn)的索引
    forbid_index = dif[(dif > high) | (dif < low)].index 
    i = 0
    while i < len(forbid_index) - 1:
        n = 1 # 發(fā)現(xiàn)連續(xù)多少個(gè)點(diǎn)變化幅度過大，大部分只有單個(gè)點(diǎn)
        start = forbid_index[i] # 異常點(diǎn)的起始索引
        while forbid_index[i+n] == start + timedelta(minutes=n):
            n += 1
        i += n - 1

        end = forbid_index[i] # 異常點(diǎn)的結(jié)束索引
        # 用前后值的中間值均勻填充
        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)

平滑后的訓(xùn)練數(shù)據(jù)：

平滑后的訓(xùn)練序列.png

3.3 將訓(xùn)練數(shù)據(jù)進(jìn)行周期性分解

采用statsmodels工具包：

from statsmodels.tsa.seasonal import seasonal_decompose

decomposition = seasonal_decompose(self.ts, freq=freq, two_sided=False)
# self.ts:時(shí)間序列，series類型; 
# freq:周期，這里為1440分鐘，即一天; 
# two_sided:觀察下圖2、4行圖，左邊空了一段，如果設(shè)為True，則會出現(xiàn)左右兩邊都空出來的情況，F(xiàn)alse保證序列在最后的時(shí)間也有數(shù)據(jù)，方便預(yù)測。

self.trend = decomposition.trend
self.seasonal = decomposition.seasonal
self.residual = decomposition.resid
decomposition.plot()
plt.show()

分解圖.png

第一行observed：原始數(shù)據(jù)；第二行trend：分解出來的趨勢部分；第三行seasonal：周期部分；最后residual：殘差部分。
我采用的是seasonal_decompose的加法模型進(jìn)行的分解，即 observed = trend + seasonal + residual，另還有乘法模型。在建模的時(shí)候，只針對trend部分學(xué)習(xí)和預(yù)測，如何將trend的預(yù)測結(jié)果加工成合理的最終結(jié)果？當(dāng)然是再做加法，后面會詳細(xì)寫。

4、模型

4.1 訓(xùn)練

對分解出來的趨勢部分單獨(dú)用arima模型做訓(xùn)練：

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

4.2 預(yù)測

預(yù)測出趨勢數(shù)據(jù)后，加上周期數(shù)據(jù)即作為最終的預(yù)測結(jié)果，但更重要的是，我們要得到的不是具體的值，而是一個(gè)合理區(qū)間，當(dāng)真實(shí)數(shù)據(jù)超過了這個(gè)區(qū)間，則觸發(fā)報(bào)警，誤差高低區(qū)間的設(shè)定來自剛剛分解出來的殘差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)

預(yù)測并完成最后的加法處理，得到第七天的預(yù)測值即高低置信區(qū)間：

    def predict_new(self):
        '''
        預(yù)測新數(shù)據(jù)
        '''
        #續(xù)接train，生成長度為n的時(shí)間索引，賦給預(yù)測序列
        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):
        '''
        為預(yù)測出的趨勢數(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í)間點(diǎn)的周期數(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)

        # 得到預(yù)測值，誤差上界和下界
        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 評估：

對第七天作出預(yù)測，評估的指標(biāo)為均方根誤差rmse，畫圖對比和真實(shí)值的差距：

    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) # 平滑過的訓(xùn)練數(shù)據(jù)加未做處理的測試數(shù)據(jù)
    plt.title(filename.split('.')[0])

    plt.subplot(212)
    pred.plot(color='blue', label='Predict') # 預(yù)測值
    test.plot(color='red', label='Original') # 真實(shí)值
    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()

預(yù)測結(jié)果.png

可以看到，均方根誤差462.8，相對于原始數(shù)據(jù)幾千的量級，還是可以的。測試數(shù)據(jù)中的兩個(gè)突變的點(diǎn)，也超過了置信區(qū)間，能準(zhǔn)確報(bào)出來。

5、結(jié)語

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