sage/入門機械学習による異常検出勉強ノート5章
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2015/09/06からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/62/
また、私の公開しているSageのサーバ(http://www15191ue.sak...
アップロードし、実行したり、変更していろいろ動きを試すこ...
井出 剛著の「入門機械学習による異常検出」(以降、井出本と...
井出本では、Rを使って例題を実行していますが、ここではSage...
** 5章 不要な次元を含むデータからの異常検知 [#sd8c8812]
この章でのポイントは、主成分分析で求めた正常な空間からの...
再構成誤差と思われます。
*** 分散最大化による定式化 [#p410c2d5]
井出本のラグランジュ乗数の説明が私が習った方法と違うので...
式を整理します。
PRMLの図12.2から1次元に投影されたD次元ベクトル\(u_1\)を使...
様子をまとめます。
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
\(u_1\)は、単位ベクトルであり、\(u_1^T u_1 = 1\)であり、...
投影されたデータの平均値は、\(u_1^T \bar{x}\)であり、\(\b...
$$
\bar{x} = \frac{1}{N} \sum_{n=1}^N x_n
$$
また、投影されたデータの分散は、以下の様に求まります。
$$
\frac{1}{N} \sum_{n=1}^N \left \{ u_1^T x_n - u_1^T \bar{...
$$
ここで、Sはデータの共分散行列であり、以下の様に定義されま...
$$
S = \frac{1}{N} \sum_{n=1}^N (x_n - \bar{x})(x_n - \bar{x...
$$
投影された分散Sの最大値を求めます。ここで\(u_1^T u_1 = 1\...
以下の様な式を\(u_1\)で微分し、0となる値を求めます。
$$
u_1^T S u_1 + (\lambda_1 - u_1^T u_1)
$$
以下の公式にSが対象行列であることを使うと、第1項の偏微分...
$$
\frac{\partial u_1^T S u_1}{\partial u_1} = (S + S^T)u_1 ...
$$
右辺をまとめると、
$$
\begin{eqnarray}
0 & = & \frac{\partial}{\partial u_1} \left [ u_1^T S u_1...
& = & ( 2S u_1 - 2 \lambda_1 u_1 )
\end{eqnarray}
$$
固有値の定義から、u_1が固有値\(\lambda_1\)の固有ベクトル...
$$
S u_1 = \lambda_1 u_1
$$
上記の式に左から\(u_1^T\)を掛け、\(u_1^T u_1 = 1 \)から、...
$$
u_1^T S u_1 = \lambda_1
$$
*** 準備 [#sef67d4b]
いつものように必要なライブラリを読み込み、テストデータと...
sageへの入力:
#pre{{
# RとPandasで必要なユーティリティ
# Rの必要なライブラリ
r('library(ggplot2)')
r('library(jsonlite)')
# RUtilにRとPandasのデータフレームを相互に変換する関数を...
load(DATA + 'RUtil.py')
}}
sageへの入力:
#pre{{
# RのテストデータCars93をPandasのデータフレームに変換
r('library(MASS)')
cars93 = RDf2PandaDf('Cars93')
}}
*** データ選択 [#a9ace7ad]
Cars93は、93年のアメリカの車の指標と集めたもので、価格...
エンジンサイズ、馬力など、15項目を使用します。
sageへの入力:
#pre{{
# 以下の15変数を選択し、Make(製品名)をインデックスとする
cc = ['Min.Price', 'Price', 'Max.Price', 'MPG.city', 'MPG...
'Length', 'Wheelbase', 'Width', 'Turn.circle', 'Wei...
X = cars93[cc]
X = X.set_index(['Make'])
X.info()
}}
#pre{{
<class 'pandas.core.frame.DataFrame'>
Index: 93 entries, Acura Integra to Volvo 850
Data columns (total 15 columns):
Min.Price 93 non-null float64
Price 93 non-null float64
Max.Price 93 non-null float64
MPG.city 93 non-null int64
MPG.highway 93 non-null int64
EngineSize 93 non-null float64
Horsepower 93 non-null int64
RPM 93 non-null int64
Rev.per.mile 93 non-null int64
Fuel.tank.capacity 93 non-null float64
Length 93 non-null int64
Wheelbase 93 non-null int64
Width 93 non-null int64
Turn.circle 93 non-null int64
Weight 93 non-null int64
dtypes: float64(5), int64(10)
memory usage: 11.6+ KB
}}
*** データのスケーリング [#ld061beb]
データを平均と分散でスケーリングし、データ個数Nをセットし...
sageへの入力:
#pre{{
# 中心にスケーリングしたデータ
Xc = (X - X.mean())/X.std()
N = Xc.shape[0]
}}
Pandasの関数で共分散行列を求めます。これとSageで定義に従...
正しく計算できていることを確認しました。
sageへの入力:
#pre{{
# Pandasで共分散行列Sigを求める
Sig = matrix(Xc.cov().values)
# show(Sig)
}}
sageへの入力:
#pre{{
# 定義に基づきSageで散布行列Sを計算し、正規化して共分散行...
Xcm = matrix(Xc.values).T
S = Xcm*Xcm.T
Sig = S/(N-1)
# show(Sig)
}}
*** 固有値と固有ベクトル [#q4654b44]
Sageのeigenmatrix_leftを使って散布行列Sの固有値と固有ベク...
固有値の値を大きい順にプロットしてみます。2〜3成分程度...
sageへの入力:
#pre{{
# 固有値Lamと固有ベクトルUを求める
(Lam, U) = S.eigenmatrix_left()
# 固有ベクトルのプロット
list_plot(Lam.diagonal(), figsize=5)
}}
&ref(sage0.png);
*** 第1と第2固有ベクトル [#p6162612]
各点を第1と第2固有ベクトルに投影します。
sageへの入力:
#pre{{
# 固有ベクトルの第1成分と第2成分上にプロットする
U12 = U.submatrix(0, 0, 2, 15)
Z = Xcm.T*(U12.T); Z
}}
#pre{{
93 x 2 dense matrix over Real Double Field (use the '.str...
}}
*** 第1と第2固有ベクトル空間での分布 [#j7669960]
観測値を第1と第2固有ベクトル空間での分布をプロットして...
いますが、分布は同じように思われます。
sageへの入力:
#pre{{
list_plot(Z.numpy(), figsize=5)
}}
&ref(sage1.png);
** 異常値の検出 [#j5f4639c]
異常値は、主成分分析を行った固有値ベクトルに投影した値と...
これを再構成残差と呼んでいます。
$$
a_1(x') = (x' - \hat{\mu})^T \left [ I_M -U_m U_m^T \righ...
$$
Sageでの計算では、numpyのarrayにした後、カラムの和(Rのcol...
sageへの入力:
#pre{{
# 要素毎の積を計算するためnumpyデータに変換する
xcm = Xcm.T.numpy()
x2 = (Xcm.T*(U12.T)).numpy()
# 要素の二乗を計算し、カラムの和を取り異常度a1を求める
a1 = (xcm*xcm).sum(1) - (x2*x2).sum(1)
}}
sageへの入力:
#pre{{
list_plot(a1, figsize=5)
}}
&ref(sage2.png);
*** 異常度の高い車種 [#i070ffc7]
異常度a1をXに追加して、異常度の高い順にソートしてトップ6...
井出本の通りの結果となりました。
sageへの入力:
#pre{{
# 異常度をXに追加
X['a1'] = a1
}}
sageへの入力:
#pre{{
# 異常度の高い順に6個出力
cols = ['a1', 'Price', 'MPG.city', 'Horsepower', 'Length']
a6 = X.sort(['a1'], ascending=False).head(6)[cols]
html(a6.to_html(classes ="table_form"))
}}
&ref(table.png);
** カーネル主成分分析 [#zb483c38]
井出本では、カーネル主成分分析の式の導出をしていますが、...
PRMLの説明の方が分かりやすいのでここでは省略します。
井出本ではRのkernlabパッケージを使っていますが、ここではs...
計算してみます。
最初にsklearnのPCAを使って第1主成分と第2主成分にXcをプ...
結果は、先ほど求めた結果と同じになりました。
sageへの入力:
#pre{{
# sklearnのPCAを使ってみる
from sklearn.decomposition import PCA, KernelPCA
pca = PCA()
X_pca = pca.fit_transform(Xc)
list_plot(X_pca[:, 0:2], figsize=5)
}}
&ref(sage3.png);
*** KernelPCAを使ってみる [#ecd0579d]
sklearnのPCAの使い方が分かったので、今度はKernelPCAで同様...
カーネルにはrbfを使いgamma=0.001と0.1の結果を出力してみま...
今度は、井出本と同様の結果が求まりました。KernelPCAは、計...
PCAよりも速く計算することができるので、データ数が多い場合...
sageへの入力:
#pre{{
# sklearnのKernelPCAを使ってみる
# σ=0.001(通常のPCAに近い)、σ=0.1の違いをみる
kpca = KernelPCA(kernel="rbf", gamma=0.001)
X_kpca = kpca.fit_transform(Xc)
g0_001 = list_plot(X_kpca[:, 0:2], figsize=4)
}}
sageへの入力:
#pre{{
kpca = KernelPCA(kernel="rbf", gamma=0.1)
X_kpca = kpca.fit_transform(Xc)
g0_1 = list_plot(X_kpca[:, 0:2], figsize=4)
html.table([['σ=0.001', 'σ=0.1'], [g0_001, g0_1]])
}}
&ref(sage4.png);
** コメント [#d73b071a]
#vote(おもしろかった[2],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2015/09/06からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/62/
また、私の公開しているSageのサーバ(http://www15191ue.sak...
アップロードし、実行したり、変更していろいろ動きを試すこ...
井出 剛著の「入門機械学習による異常検出」(以降、井出本と...
井出本では、Rを使って例題を実行していますが、ここではSage...
** 5章 不要な次元を含むデータからの異常検知 [#sd8c8812]
この章でのポイントは、主成分分析で求めた正常な空間からの...
再構成誤差と思われます。
*** 分散最大化による定式化 [#p410c2d5]
井出本のラグランジュ乗数の説明が私が習った方法と違うので...
式を整理します。
PRMLの図12.2から1次元に投影されたD次元ベクトル\(u_1\)を使...
様子をまとめます。
&ref(http://research.microsoft.com/en-us/um/people/cmbish...
\(u_1\)は、単位ベクトルであり、\(u_1^T u_1 = 1\)であり、...
投影されたデータの平均値は、\(u_1^T \bar{x}\)であり、\(\b...
$$
\bar{x} = \frac{1}{N} \sum_{n=1}^N x_n
$$
また、投影されたデータの分散は、以下の様に求まります。
$$
\frac{1}{N} \sum_{n=1}^N \left \{ u_1^T x_n - u_1^T \bar{...
$$
ここで、Sはデータの共分散行列であり、以下の様に定義されま...
$$
S = \frac{1}{N} \sum_{n=1}^N (x_n - \bar{x})(x_n - \bar{x...
$$
投影された分散Sの最大値を求めます。ここで\(u_1^T u_1 = 1\...
以下の様な式を\(u_1\)で微分し、0となる値を求めます。
$$
u_1^T S u_1 + (\lambda_1 - u_1^T u_1)
$$
以下の公式にSが対象行列であることを使うと、第1項の偏微分...
$$
\frac{\partial u_1^T S u_1}{\partial u_1} = (S + S^T)u_1 ...
$$
右辺をまとめると、
$$
\begin{eqnarray}
0 & = & \frac{\partial}{\partial u_1} \left [ u_1^T S u_1...
& = & ( 2S u_1 - 2 \lambda_1 u_1 )
\end{eqnarray}
$$
固有値の定義から、u_1が固有値\(\lambda_1\)の固有ベクトル...
$$
S u_1 = \lambda_1 u_1
$$
上記の式に左から\(u_1^T\)を掛け、\(u_1^T u_1 = 1 \)から、...
$$
u_1^T S u_1 = \lambda_1
$$
*** 準備 [#sef67d4b]
いつものように必要なライブラリを読み込み、テストデータと...
sageへの入力:
#pre{{
# RとPandasで必要なユーティリティ
# Rの必要なライブラリ
r('library(ggplot2)')
r('library(jsonlite)')
# RUtilにRとPandasのデータフレームを相互に変換する関数を...
load(DATA + 'RUtil.py')
}}
sageへの入力:
#pre{{
# RのテストデータCars93をPandasのデータフレームに変換
r('library(MASS)')
cars93 = RDf2PandaDf('Cars93')
}}
*** データ選択 [#a9ace7ad]
Cars93は、93年のアメリカの車の指標と集めたもので、価格...
エンジンサイズ、馬力など、15項目を使用します。
sageへの入力:
#pre{{
# 以下の15変数を選択し、Make(製品名)をインデックスとする
cc = ['Min.Price', 'Price', 'Max.Price', 'MPG.city', 'MPG...
'Length', 'Wheelbase', 'Width', 'Turn.circle', 'Wei...
X = cars93[cc]
X = X.set_index(['Make'])
X.info()
}}
#pre{{
<class 'pandas.core.frame.DataFrame'>
Index: 93 entries, Acura Integra to Volvo 850
Data columns (total 15 columns):
Min.Price 93 non-null float64
Price 93 non-null float64
Max.Price 93 non-null float64
MPG.city 93 non-null int64
MPG.highway 93 non-null int64
EngineSize 93 non-null float64
Horsepower 93 non-null int64
RPM 93 non-null int64
Rev.per.mile 93 non-null int64
Fuel.tank.capacity 93 non-null float64
Length 93 non-null int64
Wheelbase 93 non-null int64
Width 93 non-null int64
Turn.circle 93 non-null int64
Weight 93 non-null int64
dtypes: float64(5), int64(10)
memory usage: 11.6+ KB
}}
*** データのスケーリング [#ld061beb]
データを平均と分散でスケーリングし、データ個数Nをセットし...
sageへの入力:
#pre{{
# 中心にスケーリングしたデータ
Xc = (X - X.mean())/X.std()
N = Xc.shape[0]
}}
Pandasの関数で共分散行列を求めます。これとSageで定義に従...
正しく計算できていることを確認しました。
sageへの入力:
#pre{{
# Pandasで共分散行列Sigを求める
Sig = matrix(Xc.cov().values)
# show(Sig)
}}
sageへの入力:
#pre{{
# 定義に基づきSageで散布行列Sを計算し、正規化して共分散行...
Xcm = matrix(Xc.values).T
S = Xcm*Xcm.T
Sig = S/(N-1)
# show(Sig)
}}
*** 固有値と固有ベクトル [#q4654b44]
Sageのeigenmatrix_leftを使って散布行列Sの固有値と固有ベク...
固有値の値を大きい順にプロットしてみます。2〜3成分程度...
sageへの入力:
#pre{{
# 固有値Lamと固有ベクトルUを求める
(Lam, U) = S.eigenmatrix_left()
# 固有ベクトルのプロット
list_plot(Lam.diagonal(), figsize=5)
}}
&ref(sage0.png);
*** 第1と第2固有ベクトル [#p6162612]
各点を第1と第2固有ベクトルに投影します。
sageへの入力:
#pre{{
# 固有ベクトルの第1成分と第2成分上にプロットする
U12 = U.submatrix(0, 0, 2, 15)
Z = Xcm.T*(U12.T); Z
}}
#pre{{
93 x 2 dense matrix over Real Double Field (use the '.str...
}}
*** 第1と第2固有ベクトル空間での分布 [#j7669960]
観測値を第1と第2固有ベクトル空間での分布をプロットして...
いますが、分布は同じように思われます。
sageへの入力:
#pre{{
list_plot(Z.numpy(), figsize=5)
}}
&ref(sage1.png);
** 異常値の検出 [#j5f4639c]
異常値は、主成分分析を行った固有値ベクトルに投影した値と...
これを再構成残差と呼んでいます。
$$
a_1(x') = (x' - \hat{\mu})^T \left [ I_M -U_m U_m^T \righ...
$$
Sageでの計算では、numpyのarrayにした後、カラムの和(Rのcol...
sageへの入力:
#pre{{
# 要素毎の積を計算するためnumpyデータに変換する
xcm = Xcm.T.numpy()
x2 = (Xcm.T*(U12.T)).numpy()
# 要素の二乗を計算し、カラムの和を取り異常度a1を求める
a1 = (xcm*xcm).sum(1) - (x2*x2).sum(1)
}}
sageへの入力:
#pre{{
list_plot(a1, figsize=5)
}}
&ref(sage2.png);
*** 異常度の高い車種 [#i070ffc7]
異常度a1をXに追加して、異常度の高い順にソートしてトップ6...
井出本の通りの結果となりました。
sageへの入力:
#pre{{
# 異常度をXに追加
X['a1'] = a1
}}
sageへの入力:
#pre{{
# 異常度の高い順に6個出力
cols = ['a1', 'Price', 'MPG.city', 'Horsepower', 'Length']
a6 = X.sort(['a1'], ascending=False).head(6)[cols]
html(a6.to_html(classes ="table_form"))
}}
&ref(table.png);
** カーネル主成分分析 [#zb483c38]
井出本では、カーネル主成分分析の式の導出をしていますが、...
PRMLの説明の方が分かりやすいのでここでは省略します。
井出本ではRのkernlabパッケージを使っていますが、ここではs...
計算してみます。
最初にsklearnのPCAを使って第1主成分と第2主成分にXcをプ...
結果は、先ほど求めた結果と同じになりました。
sageへの入力:
#pre{{
# sklearnのPCAを使ってみる
from sklearn.decomposition import PCA, KernelPCA
pca = PCA()
X_pca = pca.fit_transform(Xc)
list_plot(X_pca[:, 0:2], figsize=5)
}}
&ref(sage3.png);
*** KernelPCAを使ってみる [#ecd0579d]
sklearnのPCAの使い方が分かったので、今度はKernelPCAで同様...
カーネルにはrbfを使いgamma=0.001と0.1の結果を出力してみま...
今度は、井出本と同様の結果が求まりました。KernelPCAは、計...
PCAよりも速く計算することができるので、データ数が多い場合...
sageへの入力:
#pre{{
# sklearnのKernelPCAを使ってみる
# σ=0.001(通常のPCAに近い)、σ=0.1の違いをみる
kpca = KernelPCA(kernel="rbf", gamma=0.001)
X_kpca = kpca.fit_transform(Xc)
g0_001 = list_plot(X_kpca[:, 0:2], figsize=4)
}}
sageへの入力:
#pre{{
kpca = KernelPCA(kernel="rbf", gamma=0.1)
X_kpca = kpca.fit_transform(Xc)
g0_1 = list_plot(X_kpca[:, 0:2], figsize=4)
html.table([['σ=0.001', 'σ=0.1'], [g0_001, g0_1]])
}}
&ref(sage4.png);
** コメント [#d73b071a]
#vote(おもしろかった[2],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
ページ名:
SmartDoc