井出 剛著の「入門機械学習による異常検出」(以降、井出本と記す)の例題をSageを使ってお復習いします。
この章でのポイントは、主成分分析で求めた正常な空間からのずれの距離を異常度とする 再構成誤差と思われます。
井出本のラグランジュ乗数の説明が私が習った方法と違うので、ここではPRMLの12章の定式化を元に 式を整理します。
PRMLの図12.2から1次元に投影されたD次元ベクトル$u_1$を使って、分散最大化による定式化の 様子をまとめます。
いつものように必要なライブラリを読み込み、テストデータとしてRのMASSパッケージに含まれているCars93を使用します。
%%HTML
<link rel="stylesheet" type="text/css" href="css/sage_table_form.css">
# RとPandasのデータフレームを相互に変換する関数を読み込む
# Rの必要なライブラリ
r('library(ggplot2)')
r('library(jsonlite)')
# python用のパッケージ
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
# jupyter用のdisplayメソッド
from IPython.display import display, Latex, HTML, Math, JSON
# sageユーティリティ
load('script/sage_util.py')
# Rユーティリティ
load('script/RUtil.py')
# RのテストデータCars93をPandasのデータフレームに変換
r('library(MASS)')
cars93 = RDf2PandaDf('Cars93')
Cars93は、93年のアメリカの車の指標と集めたもので、価格(Min.Price, Price, Max.Price)と燃費(MPG.city, MPG.highway)、 エンジンサイズ、馬力など、15項目を使用します。
# 以下の15変数を選択し、Make(製品名)をインデックスとする
cc = ['Min.Price', 'Price', 'Max.Price', 'MPG.city', 'MPG.highway', 'EngineSize', 'Horsepower', 'RPM', 'Rev.per.mile', 'Fuel.tank.capacity', 'Length', 'Wheelbase', 'Width', 'Turn.circle', 'Weight', 'Make']
X = cars93[cc]
X = X.set_index(['Make'])
X.info()
データを平均と分散でスケーリングし、データ個数Nをセットします。
# 中心にスケーリングしたデータ
Xc = (X - X.mean())/X.std()
N = Xc.shape[0]
Pandasの関数で共分散行列を求めます。これとSageで定義に従って計算した共分散行列Sigの値を比較し、 正しく計算できていることを確認しました。
# Pandasで共分散行列Sigを求める
Sig = matrix(Xc.cov().values)
# show(Sig)
# 定義に基づきSageで散布行列Sを計算し、正規化して共分散行列Sigを表示して上記と一致することを確認
Xcm = matrix(Xc.values).T
S = Xcm*Xcm.T
Sig = S/(N-1)
# show(Sig)
Sageのeigenmatrix_leftを使って散布行列Sの固有値と固有ベクトルを求めます。
固有値の値を大きい順にプロットしてみます。2〜3成分程度で十分表現できることが分かります。
# 固有値Lamと固有ベクトルUを求める
(Lam, U) = S.eigenmatrix_left()
# 固有ベクトルのプロット
list_plot(Lam.diagonal(), figsize=4, zorder=2)
各点を第1と第2固有ベクトルに投影します。
# 固有ベクトルの第1成分と第2成分上にプロットする
U12 = U.submatrix(0, 0, 2, 15)
Z = Xcm.T*(U12.T); Z
観測値を第1と第2固有ベクトル空間での分布をプロットしてみます。井出本の分布と比べると第1成分が正負逆になって いますが、分布は同じように思われます。
list_plot(Z.numpy(), figsize=4, zorder=2)
異常値は、主成分分析を行った固有値ベクトルに投影した値と観測値の残差と井出本では定義しています。 これを再構成残差と呼んでいます。 $$ a_1(x') = (x' - \hat{\mu})^T \left [ I_M -U_m U_m^T \right ] (x' -\hat{\mu}) $$
Sageでの計算では、numpyのarrayにした後、カラムの和(Rのcolsumsの代わり)を求めるため、sum(1)を使っています。
# 要素毎の積を計算するためnumpyデータに変換する
xcm = Xcm.T.numpy()
x2 = (Xcm.T*(U12.T)).numpy()
# 要素の二乗を計算し、カラムの和を取り異常度a1を求める
a1 = (xcm*xcm).sum(1) - (x2*x2).sum(1)
list_plot(a1, figsize=4, zorder=2)
異常度a1をXに追加して、異常度の高い順にソートしてトップ5を取ってみると、 井出本の通りの結果となりました。
# 異常度をXに追加
X['a1'] = a1
# 異常度の高い順に6個出力
cols = ['a1', 'Price', 'MPG.city', 'Horsepower', 'Length']
a5 = X.sort_values(by=['a1'], ascending=False)[cols]
a5.head()
井出本では、カーネル主成分分析の式の導出をしていますが、潜在変数の導入理由や考え方は、 PRMLの説明の方が分かりやすいのでここでは省略します。
井出本ではRのkernlabパッケージを使っていますが、ここではsklearnのKernelPCAを使って 計算してみます。
最初にsklearnのPCAを使って第1主成分と第2主成分にXcをプロットしてみます。 結果は、先ほど求めた結果と同じになりました。
# sklearnのPCAを使ってみる
from sklearn.decomposition import PCA, KernelPCA
pca = PCA()
X_pca = pca.fit_transform(Xc)
list_plot(X_pca[:, 0:2], figsize=4, zorder=2)
sklearnのPCAの使い方が分かったので、今度はKernelPCAで同様の計算をしてみます。 カーネルにはrbfを使いgamma=0.001と0.1の結果を出力してみます。
今度は、井出本と同様の結果が求まりました。KernelPCAは、計算に使用する固有値が少ない場合、 PCAよりも速く計算することができるので、データ数が多い場合には有効です。
# 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, zorder=2)
kpca = KernelPCA(kernel="rbf", gamma=0.1)
X_kpca = kpca.fit_transform(Xc)
g0_1 = list_plot(X_kpca[:, 0:2], figsize=4, zorder=2)
Table2Html([[_to_png(g0_001), _to_png(g0_1)]], header=['σ=0.001', 'σ=0.1'])