sage/入門機械学習による異常検出勉強ノート6章
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2015/09/06からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/63/
また、私の公開しているSageのサーバ(http://www15191ue.sak...
アップロードし、実行したり、変更していろいろ動きを試すこ...
井出 剛著の「入門機械学習による異常検出」(以降、井出本と...
井出本では、Rを使って例題を実行していますが、ここではSage...
** 6章 線形回帰モデルによる異常検知 [#he82f37c]
この章でのポイントは、線形回帰モデルを使って、観測量(y', ...
*** 最小二乗法の最尤推定 [#rbdf96f2]
$$
f(x) = \alpha_0 + \alpha^T x = \alpha_0 + \alpha_1 x_1 + ...
$$
係数αは、M+1個をもち、これをデータ\(\mathcal{D}\)から以下...
$$
\begin{eqnarray}
p(y | x) & = & \mathcal{N} (y | \alpha_0 + \alpha^T x, \s...
& = & \frac{1}{\sqrt{2 \pi \sigma^2}} exp \left [ - \frac...
\end{eqnarray}
$$
未知のパラメータ\(\alpha_0, \alpha, \sigma^2\)に対する尤...
$$
p(\mathcal{D} | \alpha_0, \alpha, \sigma^2) = \prod_{n=1}...
$$
これの対数尤度は、以下の様になります。
$$
L(\alpha_0, \alpha, \sigma^2 | \mathcal{D}) = - \frac{N}{...
$$
最初に\(\alpha_0\)で微分して、その値が0となる\(\hat{\alph...
$$
\hat{\alpha_0} = \frac{1}{N} \left [ y^{(n)} - \alpha^T x...
$$
ここで、\(\bar{x}, \bar{y}\)は、変数xと観測値yの平均です。
$$
\bar{y} = \frac{1}{N} \sum_{n=1}^N y^{(n)}, \bar{x} = \fr...
$$
求まった\(\alpha_0\)を定数尤度の式に入れ、\(\alpha\)で微...
$$
0 = \frac{1}{2} \frac{\partial}{\partial \alpha} \sum_{n=...
$$
ここで、変数xと観測値yと平均の差(中心化)を\(\tilde{X}, ...
$$
0 = \frac{1}{2} \frac{\partial}{\partial \alpha} \sum_{n=...
$$
整理すると、
$$
0 = \sum_{n=1}^N \tilde{y} \tilde{x}^T - \alpha^T \left (...
$$
これを\(\alpha\)について解くと、以下の様に求まります。
$$
\hat{\alpha} = \left ( \tilde{X}^T \tilde{X} \right )^{-1...
$$
同様に\(\sigma^{-2}\)で微分すると、
$$
\sigma^2 = \frac{1}{N} \sum_{n=1}^N \left \{ \tilde{y} - ...
$$
*** 異常度の定義 [#n9d29588]
観測量\((y', x')\)についての異常度は、この観測量に対する...
$$
a(y', x') = \frac{1}{\hat{\sigma^2}} \left [ y' - \bar{y}...
$$
*** リッジ回帰を使って異常検出 [#b84cb52f]
実際の観測量では、すべて変数xが独立ではなく互いに関連があ...
そこで、井出本ではリッジ回帰を使って\(\alpha, \sigma^2\)...
以下の様に求まります。
$$
\hat{\alpha}_{ridge} = \left [ \tilde{X}^T \tilde{X} + \l...
$$
分散\(\hat{\sigma}^2_{ridge}\)は、以下の様になります。
$$
\hat{\sigma}^2_{ridge} = \frac{1}{N} \left \{ \hat{\lambd...
$$
*** 準備 [#b6b0468e]
いつものように必要なライブラリを読み込み、テストデータと...
UScribeは、都市の犯罪率yで、以下の表6.1のような変数を使っ...
sageへの入力:
#pre{{
# RとPandasで必要なユーティリティ
# Rの必要なライブラリ
r('library(ggplot2)')
r('library(jsonlite)')
# RUtilにRとPandasのデータフレームを相互に変換する関数を...
load(DATA + 'RUtil.py')
}}
sageへの入力:
#pre{{
# RのテストデータMASSパッケージのUScribeをPandasのデータ...
r('library(MASS)')
uscrime = RDf2PandaDf('UScrime')
uscrime.info()
}}
#pre{{
<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への入力:
#pre{{
# RでのUScrimeの並びを調べるために実行
r('UScrime')
}}
#pre{{
M So Ed Po1 Po2 LF M.F Pop NW U1 U2 GDP Ineq ...
1 151 1 91 58 56 510 950 33 301 108 41 394 261 0....
2 143 0 113 103 95 583 1012 13 102 96 36 557 194 0....
途中省略
}}
sageへの入力:
#pre{{
# UScrimeの説明を表6.1から転記
vars = ['M', 'Ed', 'Po1', 'Po2', 'LF', 'M.F', 'Pop', 'NW'...
abbr = ['14-24歳の男性の割合', '平均就学期間', '1960年に...
'州の人口', '1000人当たりの非白人数', '都市部男性(14-24歳...
'経済的不平等の度合い', '収監率', '刑務所での平均収監期...
html.table(zip(vars, abbr), ['記号', '定義'], )
}}
&ref(th_tbl-6.1.jpg);
*** sklearnのRidge回帰を使った最適解 [#ecee9962]
リッジ回帰は、sklearnのlinear_modelパッケージのRidgeを使...
sklearnのパッケージはRと比べて使い方が統一されています。
回帰分析の結果求まった予測値(pred)と観測値(y)がどの程度分...
プロットしてみます。
sageへの入力:
#pre{{
# データを変数にセット
X = uscrime[vars].values
y = uscrime['y'].values
N = len(y)
M = len(vars)
}}
sageへの入力:
#pre{{
# 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)
}}
&ref(sage0.png);
*** 異常度の計算 [#le810f42]
井出本では、リッジ回帰のλについても最適な値を求めています...
ここではλ=1.0で計算した結果を使って異常度を計算します。
観測値と回帰結果との残差の自乗は、score関数で求まります。
あとは、係数α(coefs)を使って\(\hat{\sigma}^2_{ridge}\)...
観測値yと変数xの中心値を使って異常度を計算します。
sageへの入力:
#pre{{
# 異常度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)
}}
&ref(sage1.png);
計算結果は、井出本とは若干異なりますが、概ねデータの特徴...
sageへの入力:
#pre{{
uscrime['a1'] = a
out = uscrime.sort(['a1'])[vars].head()
#html(out.to_html(classes ="table_form"))
html.table(out.values.tolist(), header= vars)
}}
&ref(th_tbl-a1.jpg);
参考のために、λを0から5まで0.1刻みで計算したscoreの値をプ...
1.0でも概ね収束しているようにみえます。
sageへの入力:
#pre{{
# 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)
}}
&ref(sage2.png);
** コメント [#q37a201f]
#vote(おもしろかった[4],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2015/09/06からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/63/
また、私の公開しているSageのサーバ(http://www15191ue.sak...
アップロードし、実行したり、変更していろいろ動きを試すこ...
井出 剛著の「入門機械学習による異常検出」(以降、井出本と...
井出本では、Rを使って例題を実行していますが、ここではSage...
** 6章 線形回帰モデルによる異常検知 [#he82f37c]
この章でのポイントは、線形回帰モデルを使って、観測量(y', ...
*** 最小二乗法の最尤推定 [#rbdf96f2]
$$
f(x) = \alpha_0 + \alpha^T x = \alpha_0 + \alpha_1 x_1 + ...
$$
係数αは、M+1個をもち、これをデータ\(\mathcal{D}\)から以下...
$$
\begin{eqnarray}
p(y | x) & = & \mathcal{N} (y | \alpha_0 + \alpha^T x, \s...
& = & \frac{1}{\sqrt{2 \pi \sigma^2}} exp \left [ - \frac...
\end{eqnarray}
$$
未知のパラメータ\(\alpha_0, \alpha, \sigma^2\)に対する尤...
$$
p(\mathcal{D} | \alpha_0, \alpha, \sigma^2) = \prod_{n=1}...
$$
これの対数尤度は、以下の様になります。
$$
L(\alpha_0, \alpha, \sigma^2 | \mathcal{D}) = - \frac{N}{...
$$
最初に\(\alpha_0\)で微分して、その値が0となる\(\hat{\alph...
$$
\hat{\alpha_0} = \frac{1}{N} \left [ y^{(n)} - \alpha^T x...
$$
ここで、\(\bar{x}, \bar{y}\)は、変数xと観測値yの平均です。
$$
\bar{y} = \frac{1}{N} \sum_{n=1}^N y^{(n)}, \bar{x} = \fr...
$$
求まった\(\alpha_0\)を定数尤度の式に入れ、\(\alpha\)で微...
$$
0 = \frac{1}{2} \frac{\partial}{\partial \alpha} \sum_{n=...
$$
ここで、変数xと観測値yと平均の差(中心化)を\(\tilde{X}, ...
$$
0 = \frac{1}{2} \frac{\partial}{\partial \alpha} \sum_{n=...
$$
整理すると、
$$
0 = \sum_{n=1}^N \tilde{y} \tilde{x}^T - \alpha^T \left (...
$$
これを\(\alpha\)について解くと、以下の様に求まります。
$$
\hat{\alpha} = \left ( \tilde{X}^T \tilde{X} \right )^{-1...
$$
同様に\(\sigma^{-2}\)で微分すると、
$$
\sigma^2 = \frac{1}{N} \sum_{n=1}^N \left \{ \tilde{y} - ...
$$
*** 異常度の定義 [#n9d29588]
観測量\((y', x')\)についての異常度は、この観測量に対する...
$$
a(y', x') = \frac{1}{\hat{\sigma^2}} \left [ y' - \bar{y}...
$$
*** リッジ回帰を使って異常検出 [#b84cb52f]
実際の観測量では、すべて変数xが独立ではなく互いに関連があ...
そこで、井出本ではリッジ回帰を使って\(\alpha, \sigma^2\)...
以下の様に求まります。
$$
\hat{\alpha}_{ridge} = \left [ \tilde{X}^T \tilde{X} + \l...
$$
分散\(\hat{\sigma}^2_{ridge}\)は、以下の様になります。
$$
\hat{\sigma}^2_{ridge} = \frac{1}{N} \left \{ \hat{\lambd...
$$
*** 準備 [#b6b0468e]
いつものように必要なライブラリを読み込み、テストデータと...
UScribeは、都市の犯罪率yで、以下の表6.1のような変数を使っ...
sageへの入力:
#pre{{
# RとPandasで必要なユーティリティ
# Rの必要なライブラリ
r('library(ggplot2)')
r('library(jsonlite)')
# RUtilにRとPandasのデータフレームを相互に変換する関数を...
load(DATA + 'RUtil.py')
}}
sageへの入力:
#pre{{
# RのテストデータMASSパッケージのUScribeをPandasのデータ...
r('library(MASS)')
uscrime = RDf2PandaDf('UScrime')
uscrime.info()
}}
#pre{{
<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への入力:
#pre{{
# RでのUScrimeの並びを調べるために実行
r('UScrime')
}}
#pre{{
M So Ed Po1 Po2 LF M.F Pop NW U1 U2 GDP Ineq ...
1 151 1 91 58 56 510 950 33 301 108 41 394 261 0....
2 143 0 113 103 95 583 1012 13 102 96 36 557 194 0....
途中省略
}}
sageへの入力:
#pre{{
# UScrimeの説明を表6.1から転記
vars = ['M', 'Ed', 'Po1', 'Po2', 'LF', 'M.F', 'Pop', 'NW'...
abbr = ['14-24歳の男性の割合', '平均就学期間', '1960年に...
'州の人口', '1000人当たりの非白人数', '都市部男性(14-24歳...
'経済的不平等の度合い', '収監率', '刑務所での平均収監期...
html.table(zip(vars, abbr), ['記号', '定義'], )
}}
&ref(th_tbl-6.1.jpg);
*** sklearnのRidge回帰を使った最適解 [#ecee9962]
リッジ回帰は、sklearnのlinear_modelパッケージのRidgeを使...
sklearnのパッケージはRと比べて使い方が統一されています。
回帰分析の結果求まった予測値(pred)と観測値(y)がどの程度分...
プロットしてみます。
sageへの入力:
#pre{{
# データを変数にセット
X = uscrime[vars].values
y = uscrime['y'].values
N = len(y)
M = len(vars)
}}
sageへの入力:
#pre{{
# 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)
}}
&ref(sage0.png);
*** 異常度の計算 [#le810f42]
井出本では、リッジ回帰のλについても最適な値を求めています...
ここではλ=1.0で計算した結果を使って異常度を計算します。
観測値と回帰結果との残差の自乗は、score関数で求まります。
あとは、係数α(coefs)を使って\(\hat{\sigma}^2_{ridge}\)...
観測値yと変数xの中心値を使って異常度を計算します。
sageへの入力:
#pre{{
# 異常度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)
}}
&ref(sage1.png);
計算結果は、井出本とは若干異なりますが、概ねデータの特徴...
sageへの入力:
#pre{{
uscrime['a1'] = a
out = uscrime.sort(['a1'])[vars].head()
#html(out.to_html(classes ="table_form"))
html.table(out.values.tolist(), header= vars)
}}
&ref(th_tbl-a1.jpg);
参考のために、λを0から5まで0.1刻みで計算したscoreの値をプ...
1.0でも概ね収束しているようにみえます。
sageへの入力:
#pre{{
# 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)
}}
&ref(sage2.png);
** コメント [#q37a201f]
#vote(おもしろかった[4],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
ページ名:
SmartDoc