入門機械学習による異常検出

井出　剛著の「入門機械学習による異常検出」(以降、井出本と記す)の例題をSageを使ってお復習いします。

2章　正規分布に従うデータからの異常検出

井出本の基本は、データに対するモデルを使って負の対数尤度を求め、それを異常検出関数として使うことです。

多変数正規分布に基づく異常検知

M次元、N個の観測データ$\mathcal{D}$が以下の様な多次元正規分布に従っていると、仮定します。 $$ \mathcal{N}(x|\mu, \Sigma) = \frac{|\Sigma|^{-1/2}}{(2 \pi )^{M/2}} exp \left \{ -\frac{1}{2}(x - \mu)^T \Sigma^{-1} (x - \mu) \right \} $$

対数尤度は、以下の様に計算されます。 $$ L(\mu, \Sigma | \mathcal{D}) = ln \prod_{n=1}^{N} \mathcal{N} (x^{(n)} | \mu, \Sigma) = \sum_{n=1}^N ln \mathcal{N} (x^{(n)} | \mu, \Sigma) $$

これに多次元正規分布$\mathcal{N}$の定義を代入すると、対数尤度Lは以下の様になります。 $$ L(\mu, \Sigma | \mathcal{D}) = -\frac{M N}{2}ln(2 \pi) - \frac{N}{2} ln | \Sigma | -\frac{1}{2} \sum_{n=1}^N (x^{(n)} - \mu)^T \Sigma^{-1} (x^{(n)} - \mu) $$

これを$\mu, \Sigma$で偏微分して、尤度Lを最大にする$\mu, \Sigma$を求めます。

最初にLを$\mu$で偏微分すると、以下の様になります。 $$ \frac{\partial L(\mu, \Sigma | \mathcal{D})}{\partial \mu} = \sum_{n=1}^N \Sigma^{-1} (x^{(n)} - \mu) $$

上記の式が0になる$\hat{\mu}$は以下の様に求まります。 $$ \hat{\mu} = \frac{1}{N} \sum_{n=1}^N x^{(n)} $$

次に、$\Sigma$での偏微分です。これには、以下の関係式を使います。 $$ -ln | \Sigma | = ln | \Sigma^{-1} | $$ と、 $$ (x^{(n)} - \mu)^T \Sigma^{-1} (x^{(n)} - \mu) = Tr \left \{ \Sigma^{-1} (x^{(n)} - \mu) (x^{(n)} - \mu)^T \right \} $$ と、以下の公式を使います。 $$ \frac{\partial}{\partial A} Tr(A B) = \frac{\partial}{\partial A} Tr(B A) = B^T $$ $$ \frac{\partial}{\partial A} ln |A| = (A^{-1})^T $$

偏微分の結果は、以下の様になります。共分散行列$\Sigma$が対象行列であるので、$\Sigma = \Sigma^T$を使っています。 $$ \begin{eqnarray} \frac{\partial L(\mu, \Sigma | \mathcal{D})}{\partial (\Sigma^{-1})} & = & \frac{N}{2}ln | \Sigma^{-1} | - \frac{1}{2} \Sigma^{-1} (x^{(n)} - \mu) (x^{(n)} - \mu)^T \\ & = & \frac{N}{2} \Sigma^T - \frac{1}{2}\sum_{n=1}^N (x^{(n)} - \mu) (x^{(n)} - \mu)^T \\ & = & \frac{N}{2} \Sigma - \frac{1}{2}\sum_{n=1}^N (x^{(n)} - \mu) (x^{(n)} - \mu)^T \end{eqnarray} $$

この値が０になるような$\hat{\Sigma}$は、以下の様に求まります。 $$ \hat{\Sigma} = \frac{1}{N} \sum_{n=1}^N (x^{(n)} - \mu) (x^{(n)} - \mu)^T $$

負の対数尤度を基に、異常度$a(x')$を以下の様に定義します。これは、マハラノビス距離と呼ばれ、「ばらつきが大きい方向の変動は大目にみる」効果があるそうです。 $$ a(x') = ({x'} - \mu)^T \Sigma^{-1} ({x'} - \mu) $$

しかし、実際のRでのプログラムをみると、ちょっと違っています。実行例2.5（Davisに対する異常度の計算部分）は、以下の計算をしています。 $$ a(x') = \sum_{m=1}^M ({x'} - \mu)^2 \Sigma^{-1} $$

値としては、Wikipediaにある共分散行列が対角行列であるときの以下の式に近いと思われます。 $$ d(\vec{x}, \vec{y}) = \sqrt{\sum_{i=1}^P \frac{(x_i - y_i)^2}{\sigma_i^2}} $$

データを使って異常度を計算する

Rのcarパッケージに入っているデータDavisを使って、多次元データの異常度を求めてみましょう。

準備

いつものように必要なパッケージを読込ます。

データの性格を知る

Rのデータをpandaのデータフレームに変換して、データの性格を調べます。 最初にinfoとdescribeで個数や欠損値、平均、分散などの統計情報を求めます。

次にDavisの体重のヒストグラムと体重（weight）と身長（height）の散布図をプロットし、おおざっぱな性質を調べます。

Sageとpandasの機能を使って共分散行列を求める

Sageとpandas(numpy)の機能を使って共分散を求めみましょう。

平均mxを求めてデータの平均からの差Xcを求めます。pandasではこれらをmeanとマイナス記号で処理することができます。

数式では共分散行列$\hat{\Sigma}$は、以下の様に定義されていますが、 $$ \hat{\mu} = \frac{1}{N}\sum_{n=1}^N x^{(n)}, \hat{\Sigma} = \frac{1}{N} \sum_{n=1}^N (x^{(n)} - \hat{\mu})(x^{(n)} - \hat{\mu})^T $$ Xcが列ベクトルではなく、行ベクトルなので、以下の様に計算しています。 $$ \hat{\Sigma} = \frac{1}{N} Xc^T Xc $$

numpyの固有の理由だと思いますが、マトリックスの積に.dot()関数を使います。この辺はちょっと違和感があります。

Pnadasの機能を使うと共分散行列は一発で計算できるのですが欠損値の影響なのか値が一致しません。

Sageのマトリックスを使った異常度の計算

Sageのmatrixを使った異常度の計算方法を以下に示します。

理論をベースとしているSageと実験値を処理するためのPandas, numpyをどうやって融合するのかがSage使いの腕の見せ所です。

と言ってもやっていることは簡単です。マトリックスの要素単位の処理をするには、 Sageのmatrixにnumpy()メソッドで、numpyのマトリックスに変換するだけです。

Rやnumpyでは多くのデータを一括して処理することが多く、要素単位の処理に四則演算子が 使われますが、Sageではマトリックスやベクトルの定義に沿って計算が行われます。

Pandaを使った異常度の計算

pandaとnumpyの機能を使って異常値を計算します。ここでは、$\Sigma^{-1} Xc$がsolveの解であることを 使っています。また、次数の和には、pythonのsumを使っています。

出力結果では、１つだけ突出した異常データが見つかります。

全体的な処理から以降ではPandaを使った分散行列の方式を使います。

マハラノビス=タグチ法

異常値の検出では、全体に対する個々の観測値の異常を検出するのに使用しましたが、 個々の観測値の変数値の異常を検出するには、マハラノビス=タグチ法を使用します。

M変数のなかからいくつかの変数を選び、1変数当たりの異常度を以下の様に検出します。 $$ SN_q = -10 log_{10} \left \{ \frac{1}{N'} \sum_{n=1}^{N'} \frac{1}{a_q(x'^{(n)})/M_q} \right \} $$

MASSパッケージのroadデータで実験

マハラノビス=タグチ法の実行例2.6をSageを使って試してみましょう。

MASSパッケージのroadデータには、アメリカ２６州の交通事故死亡者数deaths、運転者数drivers、 人口密度popden、郊外地区の道路長rural、1月における1日の最高気温の平均値temp、1年度との燃料消費量 fuelの６変数を記録しています。

データの前段階でroadをdriversで割って対数化する必要があるのですが、pandaで処理する方法が見つからず、 Rで前処理することにしました（残念）。

取り込んだデータは、infoで個数と欠損値の有無を調べ、内容を調べます。 $rowという変な名前が存在するので、stateに置き替えます。

1変数当たりの異常度とCalifのSN比を求める

異常度の計算と同様に異常度を求め、次元数で割って1変数当たりの異常度を求めます。

最も異常度の大きい、CalifのデータをXcから取り出し、SN比を計算し、棒グラフで表示します。

表示順はdeaths fuel popden rural tempで、図2.7(b)と並びが異なりますが、 値は同じ結果になっています。

この結果から、異常に寄与している変数は、通事故死亡者数deaths、1年度との燃料消費量 fuel、人口密度popdenであることが分かります。 この結果は、車が多く使用される過密な大都市で交通事故死亡者が増加するという一般的な 常識と一致します。