sage/入門機械学習による異常検出勉強ノート6章

FrontPage

2015/09/06からのアクセス回数 3536

ここで紹介したSageワークシートは、以下のURLからダウンロードできます。

http://www15191ue.sakura.ne.jp:8000/home/pub/63/

また、私の公開しているSageのサーバ（http://www15191ue.sakura.ne.jp:8000/）にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。

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

井出本では、Rを使って例題を実行していますが、ここではSage, Pythonをベースに例題を試します。

6章　線形回帰モデルによる異常検知 †

この章でのポイントは、線形回帰モデルを使って、観測量(y', x')の負の対数尤度を異常度とする部分です。

最小二乗法の最尤推定 †

$$ 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のような変数を使って回帰分析します。

sageへの入力:

# RとPandasで必要なユーティリティ
# Rの必要なライブラリ
r('library(ggplot2)')
r('library(jsonlite)')
# RUtilにRとPandasのデータフレームを相互に変換する関数を追加+GgSaveを追加（2015/07/12）
load(DATA + 'RUtil.py')

sageへの入力:

# RのテストデータMASSパッケージのUScribeをPandasのデータフレームに変換
r('library(MASS)')
uscrime = RDf2PandaDf('UScrime')
uscrime.info()

<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

sageへの入力:

# RでのUScrimeの並びを調べるために実行
r('UScrime')

     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
途中省略

sageへの入力:

# UScrimeの説明を表6.1から転記
vars = ['M', 'Ed', 'Po1', 'Po2', 'LF', 'M.F', 'Pop', 'NW', 'U1', 'U2', 'GDP', 'Ineq', 'Prob', 'Time']
abbr = ['14-24歳の男性の割合', '平均就学期間', '1960年における警察予算', '1959年における警察予算', '労働力率', '女性1000人当たりの男性の数', 
'州の人口', '1000人当たりの非白人数', '都市部男性(14-24歳)の失業率', '都市部男性(35-39)の失業率', '州の1人当たりGDP',
 '経済的不平等の度合い', '収監率', '刑務所での平均収監期間']
html.table(zip(vars, abbr), ['記号', '定義'], )

sklearnのRidge回帰を使った最適解 †

リッジ回帰は、sklearnのlinear_modelパッケージのRidgeを使用しました。 sklearnのパッケージはRと比べて使い方が統一されています。

回帰分析の結果求まった予測値(pred)と観測値(y)がどの程度分散しているか プロットしてみます。

sageへの入力:

# データを変数にセット
X = uscrime[vars].values
y = uscrime['y'].values
N = len(y)
M = len(vars)

sageへの入力:

# sklearnのRidge回帰を使って最適解を見つける
from sklearn.linear_model import Ridge
clf = Ridge(alpha=1.0)
clf.fit(X, y)
pred = clf.predict(X)
list_plot(zip(y, pred), figsize=5)

異常度の計算 †

井出本では、リッジ回帰のλについても最適な値を求めていますが、 ここではλ=1.0で計算した結果を使って異常度を計算します。

観測値と回帰結果との残差の自乗は、score関数で求まります。 あとは、係数α（coefs）を使って\(\hat{\sigma}^2_{ridge}\)を求め、 観測値yと変数xの中心値を使って異常度を計算します。

sageへの入力:

# 異常度aを求める
lam = 1.0
coefs = clf.coef_
score = clf.score(X, y)
sig2 = (lam*sum(coefs^2) + score)/N
y_bar = uscrime['y'].mean()
X_bar = uscrime[vars].mean().values
a = (y - y_bar - coefs.dot((X - X_bar).T))^2/sig2
list_plot(a, figsize=5)

計算結果は、井出本とは若干異なりますが、概ねデータの特徴を捉えていると思われます。

sageへの入力:

uscrime['a1'] = a
out = uscrime.sort(['a1'])[vars].head()
#html(out.to_html(classes ="table_form"))
html.table(out.values.tolist(), header= vars)

参考のために、λを0から5まで0.1刻みで計算したscoreの値をプロットしてみました。 1.0でも概ね収束しているようにみえます。

sageへの入力:

# alphaの値を変えてみる => alpha=1.0で十分よい値が求まっていることが分かる
lams = np.arange(0, 5, 5/50)
Res = []
for lam in lams:
    clf = Ridge(alpha=lam)
    clf.fit(X,y)
    Res += [clf.score(X, y)]
list_plot(Res, figsize=5)

コメント †

皆様のご意見、ご希望をお待ちしております。

お名前:

(Input image string)