python庫Tsmoothie模塊數(shù)據(jù)平滑化異常點抓取

更新時間：2022年06月10日 15:20:31 作者：Python編程學(xué)習(xí)圈

這篇文章主要為大家介紹了python庫Tsmoothie模塊數(shù)據(jù)平滑化技術(shù)實現(xiàn)異常點抓取，有需要的朋友可以借鑒參考下，希望能夠有所幫助，祝大家多多進步，早日升職加薪

前言

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

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

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

這樣的時序數(shù)據(jù)是不是看起來舒服多了？此外，使用平滑后的時序數(shù)據(jù)去做聚類或預(yù)測或許有令人驚艷的效果，因為它去除了一些偏差值并細化了數(shù)據(jù)的分布范圍。

如果我們自己開發(fā)一個這樣的平滑工具，會耗費不少的時間。因為平滑的技術(shù)有很多種，你需要一個個地去研究，找到最合適的技術(shù)并編寫代碼，這是一個非常耗時的過程。平滑技術(shù)包括但不限于：

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

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

1.準備

開始之前，你要確保Python和pip已經(jīng)成功安裝在電腦上。

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

(可選2) 此外，推薦大家用VSCode編輯器，它有許多的優(yōu)點。

請選擇以下任一種方式輸入命令安裝依賴：

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

pip install tsmoothie

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

2.Tsmoothie 基本使用

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

import numpy as np
import matplotlib.pyplot as plt
from tsmoothie.utils_func import sim_randomwalk
from tsmoothie.smoother import LowessSmoother
# 生成 3 個長度為200的隨機數(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)

通過 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
# ....... ....... ....... ....... ....... ]]

繪制效果圖：

3.基于Tsmoothie的極端異常值檢測

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

可以看到，在藍色范圍以外的點，都屬于異常值。我們可以輕易地將這些異常值標紅或記錄，以便后續(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è)藍色范圍interval的最大值為up、最小值為low，如果存在 data > up 或 data < low 則表明此數(shù)據(jù)是異常點。

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

上滑查看更多代碼

# Origin: https://github.com/cerlymarco/MEDIUM_NoteBook/blob/master/Anomaly_Detection_RealTime/Anomaly_Detection_RealTime.ipynb
import numpy as np
import matplotlib.pyplot as plt
from celluloid import Camera
from collections import defaultdict
from functools import partial
from tqdm import tqdm
from tsmoothie.utils_func import sim_randomwalk, sim_seasonal_data
from tsmoothie.smoother import *
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 ax in range(n_series)]
series = defaultdict(partial(np.ndarray, shape=(n_series,1), dtype='float32'))
for i in 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 s in range(n_series):
            pltargs = {k:v[s,:] for k,v in 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')

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