入門機械学習による異常検出
井出 剛著の「入門機械学習による異常検出」(以降、井出本と記す)の例題をSageを使ってお復習いします。
6章 線形回帰モデルによる異常検知
この章でのポイントは、線形回帰モデルを使って、観測量(y', x')の負の対数尤度を異常度とする部分です。
最小二乗法の最尤推定
線形回帰では、関数fを以下のような1次関数で表現します。 $$ f(x) = \alpha_0 + \alpha^T x = \alpha_0 + \alpha_1 x_1 + \cdots + \alpha_M x_M $$
係数αは、M+1個をもち、これをデータ$\mathcal{D}$から以下の確立モデルを使って求めます。 $$ \begin{eqnarray} p(y | x) & = & \mathcal{N} (y | \alpha_0 + \alpha^T x, \sigma^2) \\ & = & \frac{1}{\sqrt{2 \pi \sigma^2}} exp \left [ - \frac{1}{2\sigma^2}(y - \alpha_0 - \alpha^T x)^2 \right ] \end{eqnarray} $$
未知のパラメータ$\alpha_0, \alpha, \sigma^2$に対する尤度は、以下の様に定義されます。 $$ p(\mathcal{D} | \alpha_0, \alpha, \sigma^2) = \prod_{n=1}^N \mathcal{N}(y^{(n)} | \alpha_0, \alpha x^{(n)}, \sigma^2) $$ これの対数尤度は、以下の様になります。 $$ L(\alpha_0, \alpha, \sigma^2 | \mathcal{D}) = - \frac{N}{2} ln(2\pi\sigma^2) - \frac{1}{2\sigma^2} \sum_{n-=1}^N \left [ y^{(n)} - \alpha_0 - \alpha^T x^{(n)} \right ]^2 $$
最初に$\alpha_0$で微分して、その値が0となる$\hat{\alpha_0}$を求めると、以下の様になります。 $$ \hat{\alpha_0} = \frac{1}{N} \left [ y^{(n)} - \alpha^T x^{(n)} \right ] = \bar{y} - \alpha^T \bar{x} $$ ここで、$\bar{x}, \bar{y}$は、変数xと観測値yの平均です。 $$ \bar{y} = \frac{1}{N} \sum_{n=1}^N y^{(n)}, \bar{x} = \frac{1}{N} \sum_{n=1}^N x^{(n)} $$
求まった$\alpha_0$を定数尤度の式に入れ、$\alpha$で微分し0（勾配0）とすると、 $$ 0 = \frac{1}{2} \frac{\partial}{\partial \alpha} \sum_{n=1}^N \left [ y^{(n)} - \bar{y} - \alpha (x^{(n)} - \bar{x}) \right]^2 $$ ここで、変数xと観測値yと平均の差（中心化）を$\tilde{X}, \tilde{y}$で表すと、 $$ 0 = \frac{1}{2} \frac{\partial}{\partial \alpha} \sum_{n=1}^N \left [ \tilde{y} - \alpha \tilde{X} \right]^2 = \sum_{n=1}^N \left \{ \tilde{y} - \alpha^T \tilde{X} \right \} \tilde{X}^T $$ 整理すると、 $$ 0 = \sum_{n=1}^N \tilde{y} \tilde{x}^T - \alpha^T \left ( \sum_{n=1}^N \tilde{X} \tilde{X}^T \right ) $$ これを$\alpha$について解くと 、以下の様に求まります。 $$ \hat{\alpha} = \left ( \tilde{X}^T \tilde{X} \right )^{-1} \tilde{X}^T \tilde{y} $$
同様に$\sigma^{-2}$で微分すると、 $$ \sigma^2 = \frac{1}{N} \sum_{n=1}^N \left \{ \tilde{y} - \hat{\alpha}^T \tilde{X} \right \} ^2 $$
異常度の定義
観測量$(y', x')$についての異常度は、この観測量に対する負の対数尤度から、以下の様になります。 $$ a(y', x') = \frac{1}{\hat{\sigma^2}} \left [ y' - \bar{y} - \hat{\alpha}^T (x' - \bar{x}) \right ] ^2 $$
リッジ回帰を使って異常検出
実際の観測量では、すべて変数xが独立ではなく互いに関連があるため、行列のランクが不足して逆行列が求まりません。
そこで、井出本ではリッジ回帰を使って$\alpha, \sigma^2$を求めています。リッジ回帰での$\hat{\alpha}_{ridge}$は、 以下の様に求まります。 $$ \hat{\alpha}_{ridge} = \left [ \tilde{X}^T \tilde{X} + \lambda I_M \right ]^{-1} \tilde{X}^T \tilde{y} $$
分散$\hat{\sigma}^2_{ridge}$は、以下の様になります。 $$ \hat{\sigma}^2_{ridge} = \frac{1}{N} \left \{ \hat{\lambda} \hat{\alpha}^T_{ridge} \hat{\alpha}_{ridge} + \sum_{n=1}^N \left [ \tilde{y}^{(n)} - \hat{\alpha}^T_{ridge} \tilde{x}^{(n)} \right ] \right \} $$
準備
いつものように必要なライブラリを読み込み、テストデータとしてRのMASSパッケージに含まれているUScribeを使用します。
UScribeは、都市の犯罪率yで、以下の表6.1のような変数を使って回帰分析します。
|
|
<class 'pandas.core.frame.DataFrame'> Int64Index: 47 entries, 0 to 46 Data columns (total 16 columns): Ed 47 non-null int64 GDP 47 non-null int64 Ineq 47 non-null int64 LF 47 non-null int64 M 47 non-null int64 M.F 47 non-null int64 NW 47 non-null int64 Po1 47 non-null int64 Po2 47 non-null int64 Pop 47 non-null int64 Prob 47 non-null float64 So 47 non-null int64 Time 47 non-null float64 U1 47 non-null int64 U2 47 non-null int64 y 47 non-null int64 dtypes: float64(2), int64(14) memory usage: 6.2 KB <class 'pandas.core.frame.DataFrame'> Int64Index: 47 entries, 0 to 46 Data columns (total 16 columns): Ed 47 non-null int64 GDP 47 non-null int64 Ineq 47 non-null int64 LF 47 non-null int64 M 47 non-null int64 M.F 47 non-null int64 NW 47 non-null int64 Po1 47 non-null int64 Po2 47 non-null int64 Pop 47 non-null int64 Prob 47 non-null float64 So 47 non-null int64 Time 47 non-null float64 U1 47 non-null int64 U2 47 non-null int64 y 47 non-null int64 dtypes: float64(2), int64(14) memory usage: 6.2 KB |
M So Ed Po1 Po2 LF M.F Pop NW U1 U2 GDP Ineq Prob Time y 1 151 1 91 58 56 510 950 33 301 108 41 394 261 0.084602 26.2011 791 2 143 0 113 103 95 583 1012 13 102 96 36 557 194 0.029599 25.2999 1635 途中省略 M So Ed Po1 Po2 LF M.F Pop NW U1 U2 GDP Ineq Prob Time y 1 151 1 91 58 56 510 950 33 301 108 41 394 261 0.084602 26.2011 791 2 143 0 113 103 95 583 1012 13 102 96 36 557 194 0.029599 25.2999 1635 途中省略 |
|
sklearnのRidge回帰を使った最適解
リッジ回帰は、sklearnのlinear_modelパッケージのRidgeを使用しました。 sklearnのパッケージはRと比べて使い方が統一されています。
回帰分析の結果求まった予測値(pred)と観測値(y)がどの程度分散しているか プロットしてみます。
|
|
|
異常度の計算
井出本では、リッジ回帰のλについても最適な値を求めていますが、 ここではλ=1.0で計算した結果を使って異常度を計算します。
観測値と回帰結果との残差の自乗は、score関数で求まります。 あとは、係数α（coefs）を使って$\hat{\sigma}^2_{ridge}$を求め、 観測値yと変数xの中心値を使って異常度を計算します。
|
計算結果は、井出本とは若干異なりますが、概ねデータの特徴を捉えていると思われます。
|
参考のために、λを0から5まで0.1刻みで計算したscoreの値をプロットしてみました。 1.0でも概ね収束しているようにみえます。
|
|
|