異常検知と変化検知： 10章メモ（疎構造学習による変化検知）

以下の赤い本の10章を読みます。キャラクターは適当です。誤りがありましたらご指摘いただけますと幸いです。

異常検知と変化検知 (機械学習プロフェッショナルシリーズ)
井手剛杉山将

講談社 2015-08-08
売り上げランキング : 56236

Amazonで詳しく見る by G-Tools

f:id:cookie-box:20180513082851p:plain:w60

先の9章では「『時系列が変化した』とはどのように判定するのか」という決してトリビアルではない問題について、基本的な発想の例として「正常というより異常らしい度合い（対数尤度）のランプ関数の累積和が閾値に達したら変化したと判定する（累積和法）」という解法と、「部分時系列の、予め定めた部分時系列群のうち最近傍標本までの対数距離が大きい場合に異常部位と判定する（近傍法）」という解法を与えました。また、「少し前の状況と今の状況が違うか否か」という形で変化検知問題を抽象的に定式化し、その最も重要な実装例の1つである特異スペクトル変換法を扱いましたね。

f:id:cookie-box:20180513082908p:plain:w60

特異スペクトル変換法は式 (9.8) の過去側の分布／現在側の分布どちらの $p$ にもフォンミーゼス・フィッシャー分布（vMF分布）を仮定したものだけど、vMF分布は長さ1のベクトルの分布であるのに対して、部分時系列はそれ自体は長さ1のベクトルじゃない。部分時系列を並べた行列を特異値分解して左特異ベクトルをとる、という操作をすることで長さ1のベクトルにするんだったね。時系列の変化度とは抽象的には「少し前の状況（分布）と今の状況（分布）のKL情報量」と定義されたけど、特異スペクトル変換法は「少し前の状況」及び「今の状況」を「1本（ないしは数本）の単位ベクトル」で表現するという点で、仮定する分布というよりこの変数変換自体が独特であるように感じるね。でもこの「元データ → 元データを特異値分解して出てきた特異値の大きい数本の左特異ベクトル」という変換がノイズ除去に相当するみたいだね。

続く10章は、多変量時系列の異常検知において各変数の寄与度を知るにはどうすればいいのか、という話のようですね…察するに、ブラスバンドのハーモニーが乱れたときに、トランペットが悪いのかクラリネットが悪いのかチューバが悪いのか、各パートの異常への寄与度を知りたいといった感じですね。

なんでそんな吹奏楽部がぎすぎすしそうな例を考えたの…。

出力値が常に変動するような時系列から異常を検知するには、変数間の依存関係に着目するとよいみたいですね。トランペットとクラリネットはどちらも高音のメロディ楽器なので、両者の振幅の相関は高いだろうとか、相関が乱れたらどちらかが演奏を誤っているだろうとかそんな感じでしょうか。というかトランペットとクラリネットって実際同じようなタイミングで鳴るんでしょうか。

知らないよ…。

ただ変数間の依存関係に着目するには「直接相関」と「間接相関」の区別が重要ということで…130～131ページにある「都市ごとの教会の数と殺人件数の相関（どちらも都市の人口に比例しているだけ）」のような例って岩波DS3でも見ましたね。

小学生の握力と算数ドリルの点数の疑似相関だったっけ。学年が交絡因子で。

ああそういうのでしたね。でも、算数ってノートにねばり強く書き続ける力だと思うんですよね。つまり、握力が高まるほど長時間ノートに力強く数式を書き続けられるようになって、算数力も向上するのではないでしょうか。

話ややこしくしなくていいよ…。

132ページからもしばらく時系列からは離れて（？）、確率ベクトルの対マルコフグラフをデータからどう描けばいいのか、という話が続きますね。それでもうこの本では多変量正規分布を仮定したガウス型グラフィカルモデルのみ取り扱うようですね。132ページの「多変量正規分布というめがねをかけて」という表現がわかるようでわかりませんが置いておきますね。…あれ、データに多変量正規分布を仮定したときの共分散行列の推定値ってこの (2.4) 式でいいんでしたっけ。最近別件でアンサンブルカルマンフィルタを実装したんですけど、そちらではフィルタを実行するときアンサンブル粒子から以下の下側の式のように分母が N-1 の共分散行列を計算して、これを用いてカルマンゲインを計算したと記憶しているんですが…。

$\displaystyle \hat{\Sigma} \equiv \frac{1}{N} \sum_{n=1}^N (x^{(n)} - \hat{\mu})(x^{(n)} - \hat{\mu})^{\top} \tag{2.4}$ $\displaystyle \hat{\Sigma} \equiv \frac{1}{N - 1} \sum_{n=1}^N (x^{(n)} - \hat{\mu})(x^{(n)} - \hat{\mu})^{\top}$

前者が最尤推定量で後者が不偏推定量だね。これらの違いは、雑に言うと以下のような感じかな…。

手元にデータがある。色々な多変量正規分布の中から、手元のデータを生成する確率が最も高いものを選ぶ。これは対数尤度を最大化すればよい（計算略）。
手元にデータがある。このデータがある多変量正規分布から生成されたということを踏まえる＝つまり、手元のデータはたまたま手元にあるように生成されただけで、別の世界線では同じ多変量正規分布から生成された別のデータになっていたかもしれない。これを踏まえた上で、色々な多変量正規分布の中から、手元のデータを生成する確率が最も高いものを選ぶ。手元のデータを母集団と見做すならば共分散行列は式 (2.4) の標本共分散行列でいいのだが、「標本平均がたまたま観測された量である」という仮定の下では、この「標本平均（データに対して最も都合がよい場所に仮定した平均）の周りのばらつき」でばらつきを見積もるのはやや過小評価になる（証明略）。

最尤推定量の方は、データが fix された上でそれを多変量正規分布にぼかして考えたらどうなる、って感じで、不偏推定量の方は、データは真の分布からたまたまこのようにサンプリングされただけって感じがするんだよね。

異常検知におけるガウス型グラフィカルモデルの推定では、手元のデータを生成した最も尤もらしい共分散行列を知りたいので最尤推定量を求める。それを通して手元のデータの変数間の相関構造を知りたいだけなので、サンプリングが違ったらどうだっただろう、という考慮はそぐわない。
アンサンブルカルマンフィルタでは状態分布を代替的にアンサンブル粒子で表現していて、フィルタのステップではアンサンブル粒子を共分散行列に翻訳したい。このとき、現在のアンサンブル粒子はたまたま今回の計算のために生成しただけなので、現在のアンサンブル粒子を生成した最も尤もらしい共分散行列を知りたいのではない。色々なアンサンブル粒子の選び方があることを踏まえて、不偏推定量を求める。

長々と書いてみたけど N が大きかったら N と N-1 の違いなんて大したことないし、そもそも10章通して N-1 を N に読み換えても話変わらない気がするんだけどね…。

お疲れさまです。ところで、133ページを読むと、結局 $x_i$ と $x_j$ の間に直接相関がないことと、精度行列の $(i,j)$ 成分がゼロであることが同値なんですね？ということは、「教会の数」「殺人件数」「人口」のデータの場合、「教会の数」と「殺人件数」の組合せに対応する精度行列の要素がゼロになるということですが、本当にちゃんとそうなるんでしょうか…。疑似相関とはいえ相関はあるんですよね？あとこれって一意に定まらなさそうな気もするんですよね。だって「教会の数」「殺人件数」「人口」のどれがニワトリでタマゴかなんてデータは語らないですよね？

今回のケースでは、132ページのグラフみたいに、「人口」軸で輪切りにしたときに「教会の数」と「殺人件数」は無相関だ、というのはわかるよね。「教会の数」軸での輪切りや「殺人件数」軸での輪切りではそうならないんじゃないかな。雑な例だけど実験してみるね。まず、0番目（ Python だから便宜上 0 から始めるね）の変数から1 番目の変数と2番目の変数それぞれにグラフの辺があって、1番目の変数と2番目の変数の間にはグラフの辺がない例。

# -*- coding: utf-8 -*- 
import numpy as np

# 標本共分散行列
def calc_cov_mat(z_en):
    z_hat = z_en - np.average(z_en, axis=0)
    v = np.outer(np.zeros(z_en.shape[1]), np.zeros(z_en.shape[1]))
    for i_en in range(z_en.shape[0]):
        v = v + np.outer(z_hat[i_en], z_hat[i_en])
    v = v / z_en.shape[0]
    return v

n_en = 100
x0_en = np.random.randn(n_en, 1)
x1_en = x0_en + 0.01 * np.random.randn(n_en, 1)
x2_en = -1.0 * x0_en + 0.01 * np.random.randn(n_en, 1)
x_en = np.concatenate([x0_en, x1_en, x2_en], axis=1)

Sigma = calc_cov_mat(x_en)
Lambda = np.linalg.inv(Sigma)

print(Sigma / Sigma[0, 0])
print(Lambda / Lambda[0, 0])

[[ 1.          0.99916607 -1.00078985]
 [ 0.99916607  0.99843717 -0.9999501 ]
 [-1.00078985 -0.9999501   1.00167383]]
[[ 1.         -0.4735883   0.52634417]
 [-0.4735883   0.44918519 -0.02475815]
 [ 0.52634417 -0.02475815  0.50121079]]

精度行列 Lambda の (1, 2) 要素＝(2, 1) 要素は小さいね。試しに変数の役割を入れ替えて、x2―x0、x2―x1 にグラフの辺がある場合もやってみるね。

x2_en = np.random.randn(n_en, 1)
x0_en = x2_en + 0.01 * np.random.randn(n_en, 1)
x1_en = -1.0 * x2_en + 0.01 * np.random.randn(n_en, 1)
x_en = np.concatenate([x0_en, x1_en, x2_en], axis=1)

Sigma = calc_cov_mat(x_en)
Lambda = np.linalg.inv(Sigma)

print(Sigma / Sigma[0, 0])
print(Lambda / Lambda[0, 0])

[[ 1.         -0.99804501  0.99911337]
 [-0.99804501  0.99624298 -0.99724056]
 [ 0.99911337 -0.99724056  0.99831405]]
[[ 1.          0.08005913 -0.92082764]
 [ 0.08005913  1.15942445  1.07805449]
 [-0.92082764  1.07805449  1.99854651]]

確かに今度は (0, 1) 要素＝(1, 0) 要素が小さいですね。偏相関係数をきちんととれるものなんですね。それで、134ページからは、この非対角要素に偏相関係数（に比例する量）をもつ精度行列をなるべく疎にしたいという話になるんですね。精度行列において閾値未満の非対角要素をゼロに切り捨てるというアイデアは駄目で、そもそも共分散行列が正則でなく精度行列が求まらなかったり、精度行列が求まっても閾値に敏感だったり切り捨てた結果の精度行列が正定値でなくなったりして多変量正規分布の精度行列ではなくなってしまうので駄目だとありますね。じゃあなんでそのアイデアを出してみたんですかね。

いや、素朴な発想が上手くいかないと断っておくことも大事だからね？何にせよ、そういう要望なら目的関数に正則化項を入れて多変量正規分布をフィッティングするんだろうね。

なぜネタバレするんです…。おっしゃる通り、尤度関数に精度行列の要素の和のラプラス分布を乗じた上でこれを最大化するんですね。対数をとると単純に対数尤度関数から $\rho || \Lambda||_1$ をマイナスすることになります（ $\rho$ はラプラス分布の尺度の逆数で、 $|| \Lambda||_1$ は精度行列の各要素の絶対値の和ですね）。この項のことをしばしばL1正則化項と呼称すると。しかし、このL1正則化項の導入によって、精度行列を推定するのに単に標本共分散行列の逆行列を取ればよいというわけにはいかなくなってしまいました。どうやって精度行列を求めればいいんでしょう…。

精度行列に関する（L1正則化項付き）対数尤度の偏微分をとればいいんじゃない。

しかし137ページに $\rho>0$ の場合は「偏微分＝0」が解析的に解けないとありますね。 $\rho<0$ だったら解けるんでしょうか…そもそも今回は精度行列を疎にしたいという目的より必ず $\rho>0$ なので、 $\rho<0$ の場合に興味はありませんが…。それで、解析的に解けないのでどうするかというと… $j$ 次元目以外固定してしまう？ $j$ 次元目を一番後ろに寄せて、精度行列、精度行列の逆行列、標本共分散行列を137ページのようにブロック表記して、そうすると各行列が4つの部分に分解されて、うち2つは転置の関係なので、元の「偏微分＝0」は (10.15) (10.16) (10.17) の3式になるわけですね…それで138～139ページの手続きにしたがえば元の精度行列の $j$ 次元目の行・列以外を固定したもとで目的関数を極大値にできますが…それを全ての次元について、ぐるぐると繰り返せば目的の精度行列が求まると…このように、元データを多変量正規分布のめがねでみてかつ疎な精度行列を得ようとする手法を「グラフィカルラッソ」というんですね。うーん、手順はわかるんですが、これ最適解に収束するんですか？ある次元以外固定してしまおう、というのがなんか心配なんですが…。

その辺の議論はこの本にはなさそうだね。参考文献を参照かな。

…まあいいです。何らかの参照データから疎な精度行列が求まり、変数間の依存関係がわかったとしましょう。これを利用する異常検知手法が 10.5 節に2例紹介されています。1例目は、ある多変量の観測値が得られたときに、その $j$ 次元目の要素の異常度は、 $j$ 次元目以外の要素を given とした条件付き確率密度の対数のマイナス1倍とする（外れ値解析）。

つまり、時系列の $j$ 次元目だけに着目してよく観測される量だったとしても、 $j$ 次元目以外がこの状況のときには観測されない量だ、という場合は異常と判定されるようなイメージだね。

2つ目の例は1点の観測値だけでなくまとまった観測値に対する異常判定ですね。このときは、参照データと観測データのそれぞれの分布を求めて、両者間のKL情報量を $j$ 次元目以外について積分したものが異常度ということです。どちらの例も特にグラフィカルラッソや、何なら多変量正規分布でなくてもよさそうですね。ただ相関構造をシンプルにしておかなければ有効な解析ができないということなのですかね。多変量正規分布の場合の (10.28) → (10.29)、(10.30) → (10.31) の式変形は後で追っておきましょう。…143ページに「 $\mathcal{D}$ と $\mathcal{D'}$ の立場を入れ替えた定義も可能です」とありますが、確かにKL情報量は非対称ですが、実際過去側の分布と現在側の分布のどちらをどちらにすべきなんでしょう？

9章を読んだときにも似た話をしたね。KL情報量は対数密度比の期待値だから、過去側の分布で期待値をとるのか現在側の分布で期待値をとるのかって話だと思うけど…過去側の分布が基準と考えれば (10.30) 式の向きが自然なんだろうけど。異常時に分布がどう変化するかにもよるだろうね。

（他の章のメモがあれば）つづく