Tsmoothie：使用多種平滑技術(shù)平滑化時(shí)序數(shù)據(jù)

在處理數(shù)據(jù)的時(shí)候，我們經(jīng)常會(huì)遇到一些非連續(xù)的散點(diǎn)時(shí)間序列數(shù)據(jù)：

有些時(shí)候，這樣的散點(diǎn)數(shù)據(jù)是不利于我們進(jìn)行數(shù)據(jù)的聚類和預(yù)測(cè)的。因此我們需要把它們平滑化，如下圖所示：

如果我們將散點(diǎn)及其范圍區(qū)間都去除，平滑后的效果如下：

這樣的時(shí)序數(shù)據(jù)是不是看起來(lái)舒服多了？此外，使用平滑后的時(shí)序數(shù)據(jù)去做聚類或預(yù)測(cè)或許有令人驚艷的效果，因?yàn)樗コ艘恍┢钪挡⒓?xì)化了數(shù)據(jù)的分布范圍。

如果我們自己開(kāi)發(fā)一個(gè)這樣的平滑工具，會(huì)耗費(fèi)不少的時(shí)間。因?yàn)槠交募夹g(shù)有很多種，你需要一個(gè)個(gè)地去研究，找到最合適的技術(shù)并編寫(xiě)代碼，這是一個(gè)非常耗時(shí)的過(guò)程。平滑技術(shù)包括但不限于：

• 指數(shù)平滑
• 具有各種窗口類型（常數(shù)、漢寧、漢明、巴特利特、布萊克曼）的卷積平滑
• 傅立葉變換的頻譜平滑
• 多項(xiàng)式平滑
• 各種樣條平滑（線性、三次、自然三次）
• 高斯平滑
• 二進(jìn)制平滑

所幸，有大佬已經(jīng)為我們實(shí)現(xiàn)好了時(shí)間序列的這些平滑技術(shù)，并在GitHub上開(kāi)源了這份模塊的代碼——它就是 Tsmoothie 模塊。

1.準(zhǔn)備

開(kāi)始之前，你要確保Python和pip已經(jīng)成功安裝在電腦上，如果沒(méi)有，可以訪問(wèn)這篇文章：超詳細(xì)Python安裝指南 進(jìn)行安裝。

**(可選1) **如果你用Python的目的是數(shù)據(jù)分析，可以直接安裝Anaconda：Python數(shù)據(jù)分析與挖掘好幫手—Anaconda，它內(nèi)置了Python和pip.

**(可選2) **此外，推薦大家用VSCode編輯器，它有許多的優(yōu)點(diǎn)：Python 編程的最好搭檔—VSCode 詳細(xì)指南。

請(qǐng)選擇以下任一種方式輸入命令安裝依賴 ：

1. Windows 環(huán)境 打開(kāi) Cmd (開(kāi)始-運(yùn)行-CMD)。
2. MacOS 環(huán)境 打開(kāi) Terminal (command+空格輸入Terminal)。
3. 如果你用的是 VSCode編輯器 或 Pycharm，可以直接使用界面下方的Terminal.

pip install tsmoothie

(PS) Tsmoothie 僅支持Python 3.6 及以上的版本。

2.Tsmoothie 基本使用

為了嘗試Tsmoothie的效果，我們需要生成隨機(jī)數(shù)據(jù)：

import numpy as np
import matplotlib.pyplot as plt
from tsmoothie.utils_func import sim_randomwalk
from tsmoothie.smoother import LowessSmoother

# 生成 3 個(gè)長(zhǎng)度為200的隨機(jī)數(shù)據(jù)組
np.random.seed(123)
data = sim_randomwalk(n_series=3, timesteps=200,
                      process_noise=10, measure_noise=30)

然后使用Tsmoothie執(zhí)行平滑化：

# 平滑
smoother = LowessSmoother(smooth_fraction=0.1, iterations=1)
smoother.smooth(data)

通過(guò) smoother.smooth_data 你就可以獲取平滑后的數(shù)據(jù)：

print(smoother.smooth_data)
# [[ 5.21462928 3.07898076 0.93933646 -1.19847767 -3.32294934
# -5.40678762 -7.42425709 -9.36150892 -11.23591897 -13.05271523
# ....... ....... ....... ....... ....... ]]

繪制效果圖：

# 生成范圍區(qū)間
low, up = smoother.get_intervals('prediction_interval')

plt.figure(figsize=(18,5))

for i in range(3):
    
    plt.subplot(1,3,i+1)
    plt.plot(smoother.smooth_data[i], linewidth=3, color='blue')
    plt.plot(smoother.data[i], '.k')
    plt.title(f"timeseries {i+1}"); plt.xlabel('time')

    plt.fill_between(range(len(smoother.data[i])), low[i], up[i], alpha=0.3)

3.基于Tsmoothie的極端異常值檢測(cè)

事實(shí)上，基于smoother生成的范圍區(qū)域，我們可以進(jìn)行異常值的檢測(cè)：

可以看到，在藍(lán)色范圍以外的點(diǎn)，都屬于異常值。我們可以輕易地將這些異常值標(biāo)紅或記錄，以便后續(xù)的處理。

_low, _up = smoother.get_intervals('sigma_interval', n_sigma=2)
series['low'] = np.hstack([series['low'], _low[:,[-1]]])
series['up'] = np.hstack([series['up'], _up[:,[-1]]])
is_anomaly = np.logical_or(
    series['original'][:,-1] > series['up'][:,-1],
    series['original'][:,-1] < series['low'][:,-1]
).reshape(-1,1)

假設(shè)藍(lán)色范圍interval的最大值為up、最小值為low，如果存在 data > up 或 data < low 則表明此數(shù)據(jù)是異常點(diǎn)。

使用以下代碼通過(guò)滾動(dòng)數(shù)據(jù)點(diǎn)進(jìn)行平滑化和異常檢測(cè)，就能保存得到上方的GIF動(dòng)圖。

上滑查看更多代碼

# Origin: https://github.com/cerlymarco/MEDIUM_NoteBook/blob/master/Anomaly_Detection_RealTime/Anomaly_Detection_RealTime.ipynb
import numpyas np
import matplotlib.pyplotas plt
from celluloidimport Camera
from collectionsimport defaultdict
from functoolsimport partial
from tqdmimport tqdm

from tsmoothie.utils_funcimport sim_randomwalk, sim_seasonal_data
from tsmoothie.smootherimport *


def plot_history(ax, i, is_anomaly, window_len, color='blue', **pltargs):

    posrange = np.arange(0,i)

    ax.fill_between(posrange[window_len:],
                    pltargs['low'][1:], pltargs['up'][1:],
                    color=color, alpha=0.2)
if is_anomaly:
        ax.scatter(i-1, pltargs['original'][-1], c='red')
else:
        ax.scatter(i-1, pltargs['original'][-1], c='black')
    ax.scatter(i-1, pltargs['smooth'][-1], c=color)

    ax.plot(posrange, pltargs['original'][1:],'.k')
    ax.plot(posrange[window_len:],
            pltargs['smooth'][1:], color=color, linewidth=3)

if 'ano_id' in pltargs.keys():
if pltargs['ano_id'].sum() >0:
            not_zeros = pltargs['ano_id'][pltargs['ano_id']!=0]-1
            ax.scatter(not_zeros, pltargs['original'][1:][not_zeros],
                       c='red', alpha=1.)

np.random.seed(42)

n_series, timesteps =3,200

data = sim_randomwalk(n_series=n_series, timesteps=timesteps,
                      process_noise=10, measure_noise=30)

window_len =20

fig = plt.figure(figsize=(18,10))
camera = Camera(fig)

axes = [plt.subplot(n_series,1,ax+1)for axin range(n_series)]
series = defaultdict(partial(np.ndarray, shape=(n_series,1), dtype='float32'))

for iin tqdm(range(timesteps+1), total=(timesteps+1)):

if i >window_len:

        smoother = ConvolutionSmoother(window_len=window_len, window_type='ones')
        smoother.smooth(series['original'][:,-window_len:])

        series['smooth'] = np.hstack([series['smooth'], smoother.smooth_data[:,[-1]]])

        _low, _up = smoother.get_intervals('sigma_interval', n_sigma=2)
        series['low'] = np.hstack([series['low'], _low[:,[-1]]])
        series['up'] = np.hstack([series['up'], _up[:,[-1]]])

        is_anomaly = np.logical_or(
            series['original'][:,-1] > series['up'][:,-1],
            series['original'][:,-1] < series['low'][:,-1]
        ).reshape(-1,1)

if is_anomaly.any():
            series['ano_id'] = np.hstack([series['ano_id'], is_anomaly*i]).astype(int)

for sin range(n_series):
            pltargs = {k:v[s,:]for k,vin series.items()}
            plot_history(axes[s], i, is_anomaly[s], window_len,
                         **pltargs)

        camera.snap()

if i >=timesteps:
continue

    series['original'] = np.hstack([series['original'], data[:,[i]]])


print('CREATING GIF...')# it may take a few seconds
camera._photos = [camera._photos[-1]] + camera._photos
animation = camera.animate()
animation.save('animation1.gif', codec="gif", writer='imagemagick')
plt.close(fig)
print('DONE')

注意，異常點(diǎn)并非都是負(fù)面作用，在不同的應(yīng)用場(chǎng)景下，它們可能代表了不同的意義。

比如在股票中，它或許可以代表著震蕩行情中某種趨勢(shì)反轉(zhuǎn)的信號(hào)。

或者在家庭用電量分析中，它可能代表著某個(gè)時(shí)刻的用電峰值，根據(jù)這個(gè)峰值我們可以此時(shí)此刻開(kāi)啟了什么樣的電器。

所以異常點(diǎn)的作用需要根據(jù)不同應(yīng)用場(chǎng)景進(jìn)行不同的分析，才能找到它真正的價(jià)值。

總而言之，Tsmoothie 不僅可以使用多種平滑技術(shù)平滑化我們的時(shí)序數(shù)據(jù)，讓我們的模型訓(xùn)練更加有效，還可以根據(jù)平滑結(jié)果找出數(shù)據(jù)中的離群點(diǎn)，是我們做數(shù)據(jù)分析和研究的一個(gè)好幫手，非常有價(jià)值。