「金塊か、キノコ料理か」（外れ値検出問題）を解く［主成分分析］

主成分分析を用いて外れ値を見つけましょう。

ここでの方針は、主成分分析を使って３次元から２次元に縮約したデータをもとに外れ値で見つようということです。
スクリープロットを見ると、主成分２つで元データの９０％以上の情報を保持していますので２次元に縮約しても問題ありません。

結論を先に言うと、頑健主成分分析の一種であるROBPCA（後述）の結果から判断して、８７番目のデータ（正解と同じ！）が外れ値となります。通常の主成分分析では検出は無理でした。

さて、主成分分析と一口に言っても実は色々な変種が存在します。ここでは、一般の本に載っている通常の主成分分析だけではなく、これとは異なる３種類の頑健主成分分析(robust principal component analysis)という一般にはあまり馴染みがない手法も用いて外れ値を探します。

普通の主成分分析

３次元の玉データを２次元に縮約すると次の図が得られます。１００番のデータが外れ値っぽいですね。

gold.pc <- princomp(gold)
summary(gold.pc,loadings=TRUE)

f:id:sudillap:20130324232204p:plain

ちなみにデータを標準化しない場合は次のようになり、今度は７３番、１００番のデータが怪しいですね。
f:id:sudillap:20130324232607p:plain

実は普通の主成分分析は外れ値の影響を受けやすいことが知られています（分散、共分散が外れ値の影響を受けやすいため）。明らかに外れ値を含んでいるデータに対して主成分分析を用いること自体、良くなかったのかもしれません。

外れ値の影響を受けやすいことが問題なのであれば、その影響を受けにくい（頑健な）主成分分析を使えば良いかもしれません。
そこで次の頑健主成分分析を用いてみます。カギ括弧内は対応するRのライブラリ名です。

最小体積楕円体（minimum volume ellipsoid）［MASS］
最小共分散行列式（minimum covariance determinant）［MASS］
ROBPCA（ROBust method for Principal Component Analysis）［rrcov］

最小体積楕円体による主成分分析

最小体積楕円体により求めた頑健な分散・共分散行列を用いて主成分分析を行います。標準化したデータと非標準化データによる結果は次のとおりです。外れ値らしい点はいくつかありますが１つだけ選ぶことは困難です。

library(MASS)
gold.pc.rob <- princomp(gold,covmat=cov.mve(gold))
summary(gold.pc.rob,loadings=TRUE)

f:id:sudillap:20130324234337p:plain f:id:sudillap:20130324234340p:plain

最小共分散行列式による主成分分析

最小体積楕円体による主成分分析の結果は最小体積楕円体の場合とほとんど同じなので省略します。

library(MASS)
gold.pc.rob <- princomp(gold,covmat=cov.mcd(gold))
summary(gold.pc.rob,loadings=TRUE)

ROBPCA

標準化したデータにROBPCAを適用した結果を示します（非標準化データの結果もほぼ同じなので省略）。縦軸の値が最も大きい８７番目のデータが最も外れ値の可能性が高いですね。実のところ１，１６，６３番目のデータも外れ値となりますが、外れ値は一つしかないことがわかっていますので、８７番目を外れ値としました。
この図の見方やROBPCAの概要に興味のある方は後述するアルゴリズム概要を参考にしてください。

library(rrcov)
gold.pc.h <- PcaHubert(gold,k=2)
summary(gold.pc.h)
plot(gold.pc.h)

f:id:sudillap:20130324234421p:plain

ROBPCAアルゴリズムの概要

概要のみを書きます。詳細をお知りになりたい方は参考文献を参照してください。

データ数を$n$、次元数を$p$とし、元のデータは$n \times p$行列$X$に格納されているとします。
このとき、ROBPCAは次の３つのステップからなります。

１：元のデータを高々$n-1$次元（＝rank($X$)）の部分空間に変換します。
２：共分散行列を求めた後、データによく適合する部分空間の次元$k$を選択します。
３：データを$k$次元部分空間に射影し、この空間で頑健な位置行列(location matrix)と散乱行列(scatter matrix)を推定します。さらにこの行列から$k$個の非ゼロ固有値$l_1, \ldots, l_k$と固有ベクトルを求めます（固有ベクトルはまとめて$p \times k$行列$P$とします）。
このとき、主成分得点（$n \times k$行列$T$）は
$$
T = (X - 1 \hat{\mu}^T)P
$$
と計算されます。ここで$1$は要素がすべて１の$n$次元縦ベクトル、$\hat{\mu}$は位置推定値（頑健中心）で$p$次元縦ベクトルです。

ROBPCA主成分プロットの読み方

まず、データを次の４種類（正常なデータと３種類の外れ値）に分類します（下図参照）。

正常観測値（regular observations）: 主成分部分空間（図の２次元平面）に近接している均質な群

好てこ点(good leverage points): 主成分空間に近接していますが、正常観測値からは離れている点（図の点１と点４）

直交外れ値（orthogonal outliers）: 主成分空間の直交方向にかなり離れていますが、主成分空間に射影すると見分けられない点（図の点５）

悪てこ点(bad leverage points): 主成分空間の直交方向にかなり離れていることに加え、主成分空間に射影すると他の点からも離れている点（図の点２，点３）

f:id:sudillap:20130325103001p:plain

ROBPCAの結果をもとにRの関数を使うと、これらの４種類の点を区別できる診断プロットを表示することができます。このページの最初の方で示したROBPCAの結果がこの診断プロットです（再掲しておきます）。
f:id:sudillap:20130324234421p:plain

診断プロットは４つの区画に分割され、それぞれ上述した４種類のデータ種別に対応しています。すなわち、区画とデータは次にように対応付けられます。

左下：正常観測値（regular observations）
左上：直交外れ値（orthogonal outliers）
右下：好てこ点(good leverage points)
右上：悪てこ点(bad leverage points)

診断プロットの横軸、縦軸はそれぞれ頑健得点距離（robust score distance）、直交距離（orthogonal distance）で次の式で定義されます。
$$
\begin{align}
\text{score distance} &= \sqrt{\sum_{j=1}^{k}\frac{t_{ij}^2}{l_j}}\\
\text{orthogonal distance} &= ||x_i-\hat{\mu}-Pt_i^T||
\end{align}
$$

参考文献：

Hubert,M., Rousseeuw,P.J. and Vanden Branden,K. (2004), ROBPCA: a new approach to robust principal component analysis. Technometrics, 45, 301-320.