![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
2015/09/06からのアクセス回数 5390
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
http://www15191ue.sakura.ne.jp:8000/home/pub/62/
また、私の公開しているSageのサーバ（http://www15191ue.sakura.ne.jp:8000/）にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。
井出 剛著の「入門機械学習による異常検出」(以降、井出本と記す)の例題をSageを使ってお復習いします。
井出本では、Rを使って例題を実行していますが、ここではSage, Pythonをベースに例題を試します。
この章でのポイントは、主成分分析で求めた正常な空間からのずれの距離を異常度とする 再構成誤差と思われます。
井出本のラグランジュ乗数の説明が私が習った方法と違うので、ここでは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を使用します。
sageへの入力:
# RとPandasで必要なユーティリティ
# Rの必要なライブラリ
r('library(ggplot2)')
r('library(jsonlite)')
# RUtilにRとPandasのデータフレームを相互に変換する関数を追加+GgSaveを追加（2015/07/12）
load(DATA + 'RUtil.py')
sageへの入力:
# RのテストデータCars93をPandasのデータフレームに変換
r('library(MASS)')
cars93 = RDf2PandaDf('Cars93')
Cars93は、９３年のアメリカの車の指標と集めたもので、価格（Min.Price, Price, Max.Price）と燃費（MPG.city, MPG.highway）、 エンジンサイズ、馬力など、１５項目を使用します。
sageへの入力:
# 以下の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()
<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
データを平均と分散でスケーリングし、データ個数Nをセットします。
sageへの入力:
# 中心にスケーリングしたデータ Xc = (X - X.mean())/X.std() N = Xc.shape[0]
Pandasの関数で共分散行列を求めます。これとSageで定義に従って計算した共分散行列Sigの値を比較し、 正しく計算できていることを確認しました。
sageへの入力:
# Pandasで共分散行列Sigを求める Sig = matrix(Xc.cov().values) # show(Sig)
sageへの入力:
# 定義に基づきSageで散布行列Sを計算し、正規化して共分散行列Sigを表示して上記と一致することを確認 Xcm = matrix(Xc.values).T S = Xcm*Xcm.T Sig = S/(N-1) # show(Sig)
Sageのeigenmatrix_leftを使って散布行列Sの固有値と固有ベクトルを求めます。
固有値の値を大きい順にプロットしてみます。２〜３成分程度で十分表現できることが分かります。
sageへの入力:
# 固有値Lamと固有ベクトルUを求める (Lam, U) = S.eigenmatrix_left() # 固有ベクトルのプロット list_plot(Lam.diagonal(), figsize=5)
各点を第１と第２固有ベクトルに投影します。
sageへの入力:
# 固有ベクトルの第1成分と第2成分上にプロットする U12 = U.submatrix(0, 0, 2, 15) Z = Xcm.T*(U12.T); Z
93 x 2 dense matrix over Real Double Field (use the '.str()' method to see the entries)
観測値を第１と第２固有ベクトル空間での分布をプロットしてみます。井出本の分布と比べると第１成分が正負逆になって いますが、分布は同じように思われます。
sageへの入力:
list_plot(Z.numpy(), figsize=5)
異常値は、主成分分析を行った固有値ベクトルに投影した値と観測値の残差と井出本では定義しています。 これを再構成残差と呼んでいます。 $$ 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)を使っています。
sageへの入力:
# 要素毎の積を計算するためnumpyデータに変換する xcm = Xcm.T.numpy() x2 = (Xcm.T*(U12.T)).numpy() # 要素の二乗を計算し、カラムの和を取り異常度a1を求める a1 = (xcm*xcm).sum(1) - (x2*x2).sum(1)
sageへの入力:
list_plot(a1, figsize=5)
異常度a1をXに追加して、異常度の高い順にソートしてトップ６を取ってみると、 井出本の通りの結果となりました。
sageへの入力:
# 異常度をXに追加 X['a1'] = a1
sageへの入力:
# 異常度の高い順に6個出力 cols = ['a1', 'Price', 'MPG.city', 'Horsepower', 'Length'] a6 = X.sort(['a1'], ascending=False).head(6)[cols] html(a6.to_html(classes ="table_form"))
井出本では、カーネル主成分分析の式の導出をしていますが、潜在変数の導入理由や考え方は、 PRMLの説明の方が分かりやすいのでここでは省略します。
井出本ではRのkernlabパッケージを使っていますが、ここではsklearnのKernelPCAを使って 計算してみます。
最初にsklearnのPCAを使って第１主成分と第２主成分にXcをプロットしてみます。 結果は、先ほど求めた結果と同じになりました。
sageへの入力:
# sklearnのPCAを使ってみる from sklearn.decomposition import PCA, KernelPCA pca = PCA() X_pca = pca.fit_transform(Xc) list_plot(X_pca[:, 0:2], figsize=5)
sklearnのPCAの使い方が分かったので、今度はKernelPCAで同様の計算をしてみます。 カーネルにはrbfを使いgamma=0.001と0.1の結果を出力してみます。
今度は、井出本と同様の結果が求まりました。KernelPCAは、計算に使用する固有値が少ない場合、 PCAよりも速く計算することができるので、データ数が多い場合には有効です。
sageへの入力:
# 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への入力:
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]])
皆様のご意見、ご希望をお待ちしております。
