MSTICPy timeseries_anomalies_stlを用いた異常検出 (Splunkのanomaly系コマンドとの比較付き)

# 伝えたいこと

MSTICPyの`timeseries_anomalies_stl`関数による異常検出および`display_timeseries_anomolies`による可視化が周期・季節変動を考慮した異常検出としてSplunkよりもサクッと利用しやすかったので紹介です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/646242/ee852928-0619-6029-3c64-c61f3dcfa028.png)

**異常検出** と一概に言っても **何に基づいた異常か** 次第で、たとえ同じデータであっても結果が変わります。
いろんなメソッドを適材適所で選択したり、組み合わせたりしましょう。

※ 「MSTICPyとは？」については以前投稿しましたこちらをご参考ください。
https://qiita.com/hackeT/items/fabf6d7f76052f15cb26

# 異常検出

異常検出についての定義や典型的なメソッドについて説明している書籍や記事はたくさんありますので、ここでは割愛いしますが、知っておくべきことは **異常検出** と一概に言っても **何に基づいた異常か** 次第で、たとえ同じデータであっても結果が変わるということです。今回はその一例を紹介したいと思います。

かつてより、異常検出には主に３つのアプローチが唱えられています。
- 外れ値検出
- 変化点検出
- 異常行動検出（異常部位検出）

:::note info
自分が実行しているロジックがこのどれに当たるのかを理解することは重要だと思います。
特に「Outlier」と書いてあるものは外れ値検出であることが多く、今回紹介する
MSTICPyの関数およびSplunkコマンドらも「外れ値検出」に基づいています。
:::

# サンプルデータ

Splunk社が提供しているButtercupgameというWebサーバーへのアクセスログのチュートリアルデータから、いくつかの特徴量フィールドを選択したものをCSVとして保存しました、そのCSVファイルはこちらに配置しています。

https://github.com/Tatsuya-hasegawa/MSTICPy_utils/blob/main/data/buttercupgame_iplocation.csv


# サンプルコード

MSTICPyのライブラリ　timeseries_anomalies_stl　を利用した上記データに対するサンプルコードはこちらからダウンロードできます。

https://github.com/Tatsuya-hasegawa/MSTICPy_utils/blob/main/timeseries_anomalies_stl/msticpy_timeseries_anomalies_stl.ipynb

実行環境
```
Python 3.9.12
msticpy 2.4.0　
statsmodels 0.13.5
```

このNotebookファイルでは`buttercupgame_iplocation.csv`から *["req_time","bytes","Country","status","action"]* のフィールドを抽出し*req_time*を"%Y-%m-%d %H"の日時フォーマットに変更した上で、action="addtocart"に対して１時間ごとのアクセス数をカウントしたものについて異常検出をかけています。

Notebookファイルで今回ポイントとなるのは以下二つのセルです。

### 1. timeseries_anomalies_stl 
```python
from msticpy.analysis.timeseries import timeseries_anomalies_stl
df_count = df_action[['date_clock','count']]
df_count = df_count.set_index('date_clock')
output = timeseries_anomalies_stl(df_count)
```

**msticpy.analysis.timeseries**のライブラリ内の**timeseries_anomalies_stl**関数は
[statsmodels.tsa.seasonal.STL](https://www.statsmodels.org/dev/generated/statsmodels.tsa.seasonal.STL.html#statsmodels.tsa.seasonal.STL)をインポートしている異常検出用の機能です。

STLは*Seasonal and Trend decomposition using Loess*の略であり、データが「トレンド・サイクル」に起因しているのか「季節変動」に起因しているのかを区別するのを手助けします。具体的には、時系列データを「トレンド・サイクル」、「季節変動」、「残余」の3つの成分に分解するためにLoess (Locally estimated scatterplot smoothing)を使用して、3つのコンポーネントの平滑推定値を抽出します。

この例では`output = timeseries_anomalies_stl(df_count)`のようにオプション引数が明示的に付与されておりませんが、実際には、デフォルト値として、以下の値が付与されています。
- seasonal = 7
- period = 24 (Hourly)
- score_threshold = 3

>    Other Parameters
    ----------------
    seasonal: int, optional
        Seasonality period of the input data required for STL.
        Must be an odd integer, and should normally be >= 7 (default).
    period: int, optional
        Periodicity of the the input data. by default 24 (Hourly).
    score_threshold: float, optional
        standard deviation threshold value calculated using Z-score used to flag anomalies,
        by default 3

引用: https://msticpy.readthedocs.io/en/latest/visualization/TimeSeriesAnomalies.html#documentation-of-timeseries-anomalies-stl-function

なお明示的にオプション推定をしなかった理由は、今回利用したButtercupgameのデータが一週間ちょっとのデータであり、seasonalを設定するには短すぎた点とWebサーバーへのアクセスログにおけるperiodは各国からの時差もあるが1日（日次）が妥当だと考えたからです。

これにより、1時間単位のaddtocartアクションのカウント数としては"2022-09-05 16時"のcount 52が異常だと判断されました。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/646242/fa11c50b-93cd-9cce-f2ce-abb994b52333.png)

参考までに追加されたフィールドの説明です。

|フィールド名|説明|
|:---|:---|
|residual|残余|
|trend|トレンド・サイクル|
|seasonal|季節変動|
|weight|外れ値の影響を軽減するために利用される重みづけ|
|baseline|トレンドサイクルと季節変動を加算した値|
|score|残余のZスコア|
|anomaries|score_thresholdオプションに対してscoreが上回ったら`1(異常判定)`になるフラグ|

### 2. display_timeseries_anomolies

```python
from msticpy.nbtools.timeseries import display_timeseries_anomolies

# display_timeseries_anomolies関数の使用上、事前に結果のタイムフィールドにdatetime型を適用し、かつ時系列順に並び替える
output['date_clock'] = pd.to_datetime(output['date_clock'])
output = output.sort_values(by='date_clock')

timeseries_anomalies_plot = display_timeseries_anomolies(
    data=output, 
    y='count',
    time_column='date_clock'
    )
```

このセルでは、obzerved, baseline, anomaliesの３つのデータをプロットするタイムチャートを生成し表示します。ポイントとしてtime_columnに与えるフィールドはdatetime型にしておく必要があります。
また`y='count'`のyは単なる縦軸のラベル名ではなく、Pandasデータフレーム内のobservedで表示する列名を指定します。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/646242/4a46991d-3946-0665-aaf4-0d7ec9e3791d.png)


# Splunkのanomalydetectionコマンドの例

さて、Buttercupgameのデータが長期的なデータではないので季節変動はあまり影響しませんでしたが、異常値が出ました。次にSplunkで同じように異常検出をやってみます。
SplunkのデフォルトコマンドにはSTLを使った異常検出コマンドがないため、
今回はいくつかの異常検出コマンドのラッパーコマンドである `anomalydetection` コマンドを利用しました。
https://docs.splunk.com/Documentation/Splunk/9.0.4/SearchReference/Anomalydetection

```
source="buttercupgame_iplocation.csv" index="msticpy" sourcetype="csv"
| table req_time bytes Country status action
| eval date_clock = strftime(strptime(req_time,"%d/%b/%Y:%H:%M:%S"),"%Y-%m-%d %H")
| eval count = 1
| stats count by date_clock,action
| where action=="addtocart"
| fields date_clock count
| anomalydetection count action=annotate method=histogram
| anomalydetection count action=annotate method=zscore
| anomalydetection count method=iqr action=transform
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/646242/ba92f1cd-0ea1-fc1e-9caf-5e595b6c1ed1.png)

anomalydetectionコマンドのmethod=histgramとmethod=zscoreにて同一の二つの点らがプロットされています。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/646242/27a32c9d-aa55-ceff-84e7-97f833107f68.png)

`method=histgram`は頻度、つまり標準偏差による外れ値判定を行い`method=zscore`も全体から見た発生頻度や平均値からの乖離を用いて異常スコアを計算し、外れ値判定をするものであるため類似しています。そしてこの二つはカウントの最も多い二つを検出しておりピーク検出に近い結果となりました。

最後のmethod=iqrは `outlier`コマンドのラッパーであり、指定したフィールドについて、四分位範囲（IQR）を用いて外れ値を検出するものでしたが、検出はありませんでした。

このようにサイクル変動を気にしない異常検出のアルゴリズムとMSTICPyでSTLを用いたサイクル変動を考慮するものでは、たとえ短いデータでも検出されるものが異なることがわかりました。あとは閾値のパラメータのチューニングや特徴量となるフィールドの選定と加工のトライアンドエラーにより、正しい道を切り開く時間になると思います。この辺をChatGPTをはじめとする生成系AIやDeepLearningなどの特徴量選定およびパラメーターの最適化を自動で行うAIを用いて効率化いきたいものです。

余談ですが、Splunkにてサイクル変動を考慮するコマンドとして [predict](https://docs.splunk.com/Documentation/Splunk/9.0.4/SearchReference/Predict)コマンドというものがありperiodで周期を指定できますのでご興味あれば試してみてください。

今回はMSTICPyの` timeseries_anomalies_stl`関数の紹介でしたので、検出結果への考察はこの辺にします。
SplunkではMLTKを利用することでPythonの `statsmodels`ライブラリを利用できるでしょう。
Python for Scientific ComputingのAppに　numpy, scipy, pandas, scikit-learn と statsmodelsのライブラリが含まれているためです。

# (付録) Splunk MLTKの例

ざっと外れ値検出するMLTKの機能に「Smart Outlier」というものがあります。Density(密度)、つまりカーネル密度推定に基づいて判定します。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/646242/b1953cea-acf3-a841-4d34-719b6064c876.png)

残念ながら、閾値を一番下げても *Total Outliers* が *0* であり、異常検出されたデータはありませんでした。
Warning▲を見るとデータが少なすぎて精度を担保できないとありました。anomalydetectionコマンドが利用するような統計処理ではなく、機械学習アルゴリズムを利用する場合は、ある程度のデータ量と特徴が必要のようです。

# (付録) Splunk DSDLの紹介

Splunk App for Data Science and Deep Learning (DSDL)には、STLを実装している機械学習モデルがありました。こちらはApp内のサンプルデータに対する実行です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/646242/9bc5a211-9985-11a7-97a6-5b076305750d.png)

Outlier（外れ値）の検出には、MTLK同様にカーネル密度推定を用いています
outlier=1.0のところのみ抜粋します
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/646242/122e0746-dc1b-adc0-6553-c215752803be.png)

seasonally_and_trend_decompositionのアルゴリズムに沿って、residual,seasonality,trendが出力されていますね。

本日はこの辺で、ご一読ありがとうございました。
Happy msticpying!


# 引用・参考

[データマイニングによる異常検知](https://www.kyoritsu-pub.co.jp/book/b10010753.html)

[セキュリティエンジニアのための機械学習　第6章　異常検出](https://www.oreilly.co.jp/books/9784873119076/)　

[Splunk anomalydetection command](https://docs.splunk.com/Documentation/Splunk/9.0.4/SearchReference/Anomalydetection)

[Splunk MLTK](https://splunkbase.splunk.com/app/2890)

[Splunk DSDL](https://splunkbase.splunk.com/app/4607)

https://qiita.com/symmr/items/dd07b2e1b515a4023b83