井出 剛著の「入門機械学習による異常検出」(以降、井出本と記す)の例題をSageを使ってお復習いします。
この章でのポイントは、主成分分析で求めた正常な空間からのずれの距離を異常度とする 再構成誤差と思われます。
井出本のラグランジュ乗数の説明が私が習った方法と違うので、ここでは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$で微分し、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 = 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を使用します。
|
|
Cars93は、93年のアメリカの車の指標と集めたもので、価格(Min.Price, Price, Max.Price)と燃費(MPG.city, MPG.highway)、 エンジンサイズ、馬力など、15項目を使用します。
<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 <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をセットします。
|
Pandasの関数で共分散行列を求めます。これとSageで定義に従って計算した共分散行列Sigの値を比較し、 正しく計算できていることを確認しました。
|
|
Sageのeigenmatrix_leftを使って散布行列Sの固有値と固有ベクトルを求めます。
固有値の値を大きい順にプロットしてみます。2〜3成分程度で十分表現できることが分かります。
|
各点を第1と第2固有ベクトルに投影します。
93 x 2 dense matrix over Real Double Field (use the '.str()' method to see the entries) 93 x 2 dense matrix over Real Double Field (use the '.str()' method to see the entries) |
観測値を第1と第2固有ベクトル空間での分布をプロットしてみます。井出本の分布と比べると第1成分が正負逆になって いますが、分布は同じように思われます。
|
異常値は、主成分分析を行った固有値ベクトルに投影した値と観測値の残差と井出本では定義しています。 これを再構成残差と呼んでいます。 $$ 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)を使っています。
|
|
異常度a1をXに追加して、異常度の高い順にソートしてトップ6を取ってみると、 井出本の通りの結果となりました。
|
|
井出本では、カーネル主成分分析の式の導出をしていますが、潜在変数の導入理由や考え方は、 PRMLの説明の方が分かりやすいのでここでは省略します。
井出本ではRのkernlabパッケージを使っていますが、ここではsklearnのKernelPCAを使って 計算してみます。
最初にsklearnのPCAを使って第1主成分と第2主成分にXcをプロットしてみます。 結果は、先ほど求めた結果と同じになりました。
|
sklearnのPCAの使い方が分かったので、今度はKernelPCAで同様の計算をしてみます。 カーネルにはrbfを使いgamma=0.001と0.1の結果を出力してみます。
今度は、井出本と同様の結果が求まりました。KernelPCAは、計算に使用する固有値が少ない場合、 PCAよりも速く計算することができるので、データ数が多い場合には有効です。
|
|
|