#!/usr/bin/env python # coding: utf-8 # # Hiroshi TAKEMOTO # (take.pwave@gmail.com) # #

入門機械学習による異常検出

#

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

# # #

5章　不要な次元を含むデータからの異常検知

#

# この章でのポイントは、主成分分析で求めた正常な空間からのずれの距離を異常度とする # 再構成誤差と思われます。 #

# # # #

分散最大化による定式化

#

# 井出本のラグランジュ乗数の説明が私が習った方法と違うので、ここではPRMLの12章の定式化を元に # 式を整理します。 #

#

# PRMLの図12.2から1次元に投影されたD次元ベクトル$u_1$を使って、分散最大化による定式化の # 様子をまとめます。 #

# #

# $u_1$は、単位ベクトルであり、$u_1^T u_1 = 1$であり、各データ点は$x_n$は、$u_1^T x_n$として投影されます。（$\tilde{x}_n$の点） # 投影されたデータの平均値は、$u_1^T \bar{x}$であり、$\bar{x}$はサンプルデータの平均です。 # $$ # \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{x} \right \} = u_1^T S u_1 # $$ # ここで、Sはデータの共分散行列であり、以下の様に定義されます。 # $$ # S = \frac{1}{N} \sum_{n=1}^N (x_n - \bar{x})(x_n - \bar{x})^T # $$ # 投影された分散Sの最大値を求めます。ここで$u_1^T u_1 = 1$の制約が加わり、これをラグランジュ乗数を導入して、 # 以下の様な式を$u_1$で微分し、０となる値を求めます。 # $$ # 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 = 2 S u_1 # $$ # 右辺をまとめると、 # $$ # \begin{eqnarray} # # 0 & = & \frac{\partial}{\partial u_1} \left [ u_1^T S u_1 - \lambda_1(1 - u_1^T u_1) \right] \\ # & = & ( 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 $から、分散$\lambda_1$（固有値）は以下の様に求まります。 # $$ # u_1^T S u_1 = \lambda_1 # $$ #

# # #

準備

#

# いつものように必要なライブラリを読み込み、テストデータとしてRのMASSパッケージに含まれているCars93を使用します。 #

# # In[1]: get_ipython().run_cell_magic('HTML', '', '\n') # In[2]: # 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 get_ipython().run_line_magic('matplotlib', 'inline') # jupyter用のdisplayメソッド from IPython.display import display, Latex, HTML, Math, JSON # sageユーティリティ load('script/sage_util.py') # Rユーティリティ load('script/RUtil.py') # In[3]: # RのテストデータCars93をPandasのデータフレームに変換 r('library(MASS)') cars93 = RDf2PandaDf('Cars93') # #

データ選択

#

# Cars93は、９３年のアメリカの車の指標と集めたもので、価格（Min.Price, Price, Max.Price）と燃費（MPG.city, MPG.highway）、 # エンジンサイズ、馬力など、１５項目を使用します。 #

# # In[4]: # 以下の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をセットします。 #

# # In[5]: # 中心にスケーリングしたデータ Xc = (X - X.mean())/X.std() N = Xc.shape[0] # #

# Pandasの関数で共分散行列を求めます。これとSageで定義に従って計算した共分散行列Sigの値を比較し、 # 正しく計算できていることを確認しました。 #

# # In[6]: # Pandasで共分散行列Sigを求める Sig = matrix(Xc.cov().values) # show(Sig) # In[7]: # 定義に基づきSageで散布行列Sを計算し、正規化して共分散行列Sigを表示して上記と一致することを確認 Xcm = matrix(Xc.values).T S = Xcm*Xcm.T Sig = S/(N-1) # show(Sig) # #

固有値と固有ベクトル

#

# Sageのeigenmatrix_leftを使って散布行列Sの固有値と固有ベクトルを求めます。 #

#

# 固有値の値を大きい順にプロットしてみます。２〜３成分程度で十分表現できることが分かります。 #

# # In[8]: # 固有値Lamと固有ベクトルUを求める (Lam, U) = S.eigenmatrix_left() # 固有ベクトルのプロット list_plot(Lam.diagonal(), figsize=4, zorder=2) # #

第１と第２固有ベクトル

#

# 各点を第１と第２固有ベクトルに投影します。 #

# # In[9]: # 固有ベクトルの第1成分と第2成分上にプロットする U12 = U.submatrix(0, 0, 2, 15) Z = Xcm.T*(U12.T); Z # #

第１と第２固有ベクトル空間での分布

#

# 観測値を第１と第２固有ベクトル空間での分布をプロットしてみます。井出本の分布と比べると第１成分が正負逆になって # いますが、分布は同じように思われます。 #

# # In[10]: 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)を使っています。 #

# # In[11]: # 要素毎の積を計算するためnumpyデータに変換する xcm = Xcm.T.numpy() x2 = (Xcm.T*(U12.T)).numpy() # 要素の二乗を計算し、カラムの和を取り異常度a1を求める a1 = (xcm*xcm).sum(1) - (x2*x2).sum(1) # In[12]: list_plot(a1, figsize=4, zorder=2) # #

異常度の高い車種

#

# 異常度a1をXに追加して、異常度の高い順にソートしてトップ5を取ってみると、 # 井出本の通りの結果となりました。 #

# # In[13]: # 異常度をXに追加 X['a1'] = a1 # In[14]: # 異常度の高い順に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を使って第１主成分と第２主成分にXcをプロットしてみます。 # 結果は、先ほど求めた結果と同じになりました。 #

# # In[15]: # 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) # #

KernelPCAを使ってみる

#

# sklearnのPCAの使い方が分かったので、今度はKernelPCAで同様の計算をしてみます。 # カーネルにはrbfを使いgamma=0.001と0.1の結果を出力してみます。 #

#

# 今度は、井出本と同様の結果が求まりました。KernelPCAは、計算に使用する固有値が少ない場合、 # PCAよりも速く計算することができるので、データ数が多い場合には有効です。 #

# # In[16]: # 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) # In[17]: 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']) # In[ ]: