BigQueryで異常検知(複数系列に対応するよう改修)

概要

前回まで、異常検知系では指標の実装をメインに書いてきました(2記事だけですが)。一方で、前回のブログの終わりでも書きましたが、別パターンのsampleデータの生成と複数系列への対応も行いたいです。
BigQueryで時系列データから異常部位検出 (EMA 指数移動平均) - アクセス解析担当者のBigQuery日記

指標系を一通り作成し、改めて別パターンサンプルデータの生成・複数系列への対応をしようと考えていました。が、先にこちらをやってしまい、今後異常検知系指標をどんどん追加していき、それぞれ比較をした方が各指標の比較や活用がしやすくなると思いましたので、今回はこれをやっていきます。

今回のキーワード

1. sampleデータの生成
2. 以前作成したクエリの改修

方針

1. 複数系列のsampleデータ生成を行います。
   1. 各系列には、その系列を表すid列、時系列を表す_time列、データとして_data列をもたせます。
   2. また各系列はスケールが変わるよう、rand()で生成した乱数を系列ごとに定数倍します。
2. クエリを、上記データに適用できるよう改修します。
3. 最後に各系列と、直近の異常度の推移を確認します。

クエリ

1.複数系列のsampleデータ生成

WITH
  sample AS (
    WITH
    sample_step1 AS (
      SELECT
        id,
        generate_array(1,100) AS _time_arr #乱数生成用の箱兼_time列用の数値生成。第2引数が系列の長さ
      FROM
        UNNEST(generate_array(1,10)) AS id #第2引数が系列の数
      )
      SELECT
        id,
        _time,
        RAND() * id AS _data #系列ごとに定数倍。今回は簡単のためid倍する。
      FROM
        sample_step1,sample_step1._time_arr　AS _time
    )
    SELECT
      *
    FROM
      sample

2.クエリを、上記データに適用できるよう改修

k近傍法

BigQueryで時系列データから異常部位検出 (k近傍法) - アクセス解析担当者のBigQuery日記

EMA

BigQueryで時系列データから異常部位検出 (EMA 指数移動平均) - アクセス解析担当者のBigQuery日記

CREATE TEMPORARY FUNCTION calc_min_distance(part_arr ARRAY<FLOAT64>,part_arr_arr ARRAY<STRUCT< inner_arr ARRAY<FLOAT64>>>) RETURNS FLOAT64
LANGUAGE js AS """
  function calc_distance(x,y) {
    if(x.length != y.length) {
      return null
    } else {
      var d = 0
      for(var i = 0;i<x.length;i++){
        d += Math.pow(x[i] - y[i],2)
      }
      return Math.sqrt(d)
    }
  }
  
  var ret = null
  if(part_arr_arr == null){
    return ret
  }
  for(var i = 0 ; i < part_arr_arr.length;i++){
    var distance = calc_distance(part_arr,part_arr_arr[i].inner_arr)
    if(distance!=null){
      if(ret==null || ret > distance){
        ret = distance
      }
    }
  }
  
  return ret
""";
WITH
  sample AS (
    WITH
    sample_step1 AS (
      SELECT
        id,
        generate_array(1,100) AS _time_arr #乱数生成用の箱兼_time列用の数値生成。第2引数が系列の長さ
      FROM
        UNNEST(generate_array(1,10)) AS id #第2引数が系列の数
      )
      SELECT
        id,
        _time,
        RAND() * id AS _data #系列ごとに定数倍。今回は簡単のためid倍する。
      FROM
        sample_step1,sample_step1._time_arr　AS _time
    ),
  #####################
  #k近傍法による異常検知の実装
  #####################
  km1 AS (
     #部分系列の生成
  SELECT
     *,
     ARRAY_AGG(_data) OVER(PARTITION BY id ORDER BY _time asc ROWS BETWEEN 6 PRECEDING AND CURRENT ROW) AS part_arr # 6 PRECEDINGを、使用したい系列の長さにより変更してください
  FROM sample
  ),
  km2 AS (
    #比較したいウィンドウで、部分系列の系列を生成します。
    SELECT
      *,
      ARRAY_AGG(STRUCT(part_arr)) OVER(PARTITION BY id ORDER BY _time asc ROWS BETWEEN 30 PRECEDING AND 1 PRECEDING)　AS part_arr_arr # 30 PRECEDINGを、比較したいウィンドウサイズにより変更してください
    FROM
      km1
  ),
  km_result AS (
  #k近傍法結果をkm_min_distanceとして出力
    SELECT
       id,
       _time,
       _data,
      calc_min_distance(part_arr,part_arr_arr) AS km_min_distance
    FROM
      km2
  ),
  #####################
  #EMA(指数移動平均) n=7
  #####################
  ema1 AS (
      SELECT
        *,
        ARRAY_REVERSE(ARRAY_AGG(_data) OVER(partition by id ORDER BY _time asc ROWS BETWEEN 6 PRECEDING AND CURRENT ROW)) AS part_arr #自身とn-1行前までをARRAYにする。また自身が初めに来るようにソートを逆順にする。
      FROM
        km_result #km_resultサブクエリを引き継ぐ
  ),
  ema_result AS (
    SELECT
      *,
      (#サブクエリでEMAを計算する
        SELECT
          SUM((2/(7+1))*POW((1-2/(7+1)),k) * f_i_k) #数式部分
         FROM
          UNNEST(part_arr) AS f_i_k WITH OFFSET AS k #ARRAYの位置も一緒に出力。数式で言うkになる。
       ) AS ema
    FROM 
      ema1
  )

SELECT
  * except(part_arr)
FROM
  ema_result
WHERE id IN(1,5,10) #可視化用にデータを絞る

解説

sampleデータ生成については特に解説しない。可視化のパートでグラフ出します。

k近傍法、EMAについて

主だった変更点としては、どちらも部分系列生成のウィンドウ関数で、

partition by id
||<l
を追加し、
>|sql|
 #以前はorder by id asc
order by _time asc

に変更しました。
以前書いたsqlでは、1系列でidを時間の区切りを表すのに使っていましたが、今回のsqlで言うとそれが_time列になるからです。紛らわしくて申し訳ありません。

可視化と所感

各系列のデータ

id 1,5,10 のデータのみ出力しています。
sampleデータ生成のsql内で

RAND() * id

というしているので、それぞれ分布が異なっています。

k近傍法

1系列で見ていた時との違いとして、idが大きい系列ほど、異常度が高くなってしまっていることがわかります。

EMA

k近傍法同様、idが大きい系列ほど、高くなっています。

結論

よってどちらも何かしら正規化した指標が別途必要であると考えており、今後は別指標の作成と併せてその辺もやりたいと思います。先は長い。

以上です。