verum ipsum factum

sudillap's blog

「金塊か、キノコ料理か」(外れ値検出問題)を解く[1クラスサポートベクターマシン]

サポートベクターマシンの一種である1クラスサポートベクターマシンで外れ値を見つけてみます。

1クラスサポートベクターマシンでデータを入力空間から特徴空間に写像すると、入力空間で孤立しているデータは特徴空間内の原点近くに写像されます。この性質から支持超平面より原点側にあるデータを外れ値とみなせばいいことになります。

ついでに述べておきますと、クラス分類問題に用いられるサポートベクターマシンは教師あり学習ですが、1クラスサポートベクターマシンは教師なし学習です。したがって、外れ値検出のための教師データは不要です。


ところで、1クラスサポートベクターマシンには次に述べる2つのパラメータがあり、それらをユーザーが指定する必要があります。

$\sigma$
カーネルにRBF(radial basis function)を用いているため。
$\nu$
1クラスサポートベクターマシンは$\nu$サポートベクターマシンにもとづいています。$\nu$によりデータに占める外れ値の割合の上限を指定できます。

Rスクリプトとその結果を次に示します。SVMのパラメータのうち、$\sigma$は自動推定機能で求め、$\nu$は0.01~0.2まで0.01刻みで指定しています。
実行結果から1番、87番、100番が外れ値の候補となります。ただし、このスクリプトとは異なったパラメータ($\sigma$、$\nu$)を指定するとまた異なった結果が得られると思います。

library(kernlab)
for(i in 1:20){
  nu <- 0.01*i
  mdl <- ksvm(as.matrix(gold),type="one-svc",kernel="rbfdot",nu=nu)
  outliers <- SVindex(mdl)[which(coef(mdl)==1.0)]
  if(length(outliers)>0){
    print(paste("nu:",nu," -> ",paste(outliers,collapse=" / ")))
  }
}
[1] "nu: 0.15  ->  87 / 100"
[1] "nu: 0.16  ->  1 / 87 / 100"
[1] "nu: 0.18  ->  1 / 87 / 100"
[1] "nu: 0.19  ->  1 / 7 / 16 / 18 / 21 / 22 / 63 / 87 / 91 / 94 / 99 / 100"
[1] "nu: 0.2  ->  1 / 7 / 16 / 18 / 21 / 22 / 63 / 87 / 91 / 94 / 99 / 100"

「金塊か、キノコ料理か」(外れ値検出問題)を解く[クラスター分析]

外れ値とは他のデータから離れているという意味なので、次のようにクラスター分析を用いれは見つけられそうです。

  • データを一つしか含まない孤立したクラスターに分類されたデータ
  • クラスターの中心から離れた場所にあるデータ

ここではクラスター分析の代表的な手法であるk平均法(k-means)階層的方法を用います。

結果を先に述べると、87番目と100番目のデータのどちらかが外れ値と推測されます。

階層的クラスタリング

標準化後のデータをクラスタリングした結果は下図(デンドログラム)のとおりです。クラスタ間距離の測り方によって結果は異なりますが、全般的に100番目と87番目のデータが外れ値のようです。

gold.clust <- hclust(dist(gold)^2, method="average") #平均法の場合
plot(gold.clust,cex=0.8)

f:id:sudillap:20130325105552p:plainf:id:sudillap:20130325105600p:plainf:id:sudillap:20130325105606p:plainf:id:sudillap:20130325105609p:plainf:id:sudillap:20130325105615p:plainf:id:sudillap:20130325105616p:plainf:id:sudillap:20130325105619p:plain

k平均法1

3次元散布図から、このデータは大きく2つのクラスターに分けることができます。
k平均法でデータを2つに分けたとき、各クラスターの中心から最も離れているデータを外れ値とします。

次の図の青色のデータが外れ値候補上位5つ(1,2,4,87,100番)となります。

n.cluster <- 2
km <- kmeans(gold,n.cluster,nstart=30)
centers <- km$centers[km$cluster,]
outliers <- order(rowSums((gold-centers)^2),decreasing=TRUE)[1:5]
outliers
[1]   1   2 100  87   4

f:id:sudillap:20130325110312p:plain

k平均法2

ここでもk平均法を用いますが、別のアイデアで外れ値を探してみます。

クラスター数を増やしていけば、ある時点で外れ値のみを含むクラスターが生成されると考えられます。
データを一つだけ含むクラスターが複数個できた時点のクラスター情報をもとに、外れ値を探します。

Rスクリプトとその実行結果は次のとおりです。73番目と100番目を外れ値とみなせそうです。

max.outlier <- 3
outliers <- c()
for(iter in 1:20){
  for(i in 2:40){
    gold.kmeans <- kmeans(gold,i,nstart=30)
    found <- which(gold.kmeans$size==1)
    if(length(found)>=max.outlier){
      break
    }
  }
  outliers <- c(outliers, which(gold.kmeans$cluster %in% found))
}
table(outliers)
outliers
 22  73  79  87 100 
 12  19   3  11  20 

「金塊か、キノコ料理か」(外れ値検出問題)を解く[主成分分析]

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

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

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

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

普通の主成分分析

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

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

f:id:sudillap:20130324232204p:plain

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

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

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

  • 最小体積楕円体(minimum volume ellipsoid)[MASS]
  • 最小共分散行列式(minimum covariance determinant)[MASS]
  • ROBPCA(ROBust method for Principal Component Analysis)[rrcov]

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

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

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

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

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

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

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

ROBPCA

標準化したデータにROBPCAを適用した結果を示します(非標準化データの結果もほぼ同じなので省略)。縦軸の値が最も大きい87番目のデータが最も外れ値の可能性が高いですね。実のところ1,16,63番目のデータも外れ値となりますが、外れ値は一つしかないことがわかっていますので、87番目を外れ値としました。
この図の見方や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は次の3つのステップからなります。

1:元のデータを高々$n-1$次元(=rank($X$))の部分空間に変換します。
2:共分散行列を求めた後、データによく適合する部分空間の次元$k$を選択します。
3:データを$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$は要素がすべて1の$n$次元縦ベクトル、$\hat{\mu}$は位置推定値(頑健中心)で$p$次元縦ベクトルです。

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


まず、データを次の4種類(正常なデータと3種類の外れ値)に分類します(下図参照)。

正常観測値(regular observations)
主成分部分空間(図の2次元平面)に近接している均質な群
好てこ点(good leverage points)
主成分空間に近接していますが、正常観測値からは離れている点(図の点1と点4)
直交外れ値(orthogonal outliers)
主成分空間の直交方向にかなり離れていますが、主成分空間に射影すると見分けられない点(図の点5)
悪てこ点(bad leverage points)
主成分空間の直交方向にかなり離れていることに加え、主成分空間に射影すると他の点からも離れている点(図の点2,点3)

f:id:sudillap:20130325103001p:plain

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

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

  • 左下:正常観測値(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.

「金塊か、キノコ料理か」(外れ値検出問題)を解く[はじめに]

ITエンジニアのための実務スキル評価サービスCodeIQ機械学習分野の問題を眺めていたら、「金塊か、キノコ料理か」(by naoya_tさん)という問題がありました。

おもしろそうなので、解答締め切りは過ぎていますが挑戦してみました。


この問題の挑戦受付は終了しているため、CodeIQのページでは詳細を知ることはもうできませんが、Tech総研のページ「「機械学習基礎」簡単な問題を解いて理解しよう!後篇」では問題の詳細(+データ)と出題者のnaoya_tさんによる解答を読むことができます。

CodeIQの「回答評価のポイント」欄で「異常(外れ値)検出で、選んではいけない玉を見つけ出してもらいます。」と書かれているようにこれは外れ値を探す問題です。次の問題文(抜粋)にあるように100個のデータの中から外れ値を一つ見つけられたら正解となります。

「ここに100個の玉がある。
 この中に1つだけ貴重な石でできた玉がある。
 大切な客人ならば、それを選び出すことができるはずだ。
 大きさも重さも調べたいだけ調べてもらって構わない。」
 
  というわけで、N君への今回の依頼はこれだ。

100個の玉のうち、貴重な石でできた玉がどれなのか教えてほしい。
判定が微妙なものがあれば、怪しい玉も含めて3つぐらい教えてほしい。

データがこちらからダウンロードできるので次にあげた方法を使って外れ値検出にトライしてみました。


ほかの方もこの課題に挑戦されているようです。たとえば、

データ概要

ファイル名はhundred.txtで、100個の玉のデータが含まれています。最初の玉5つ分のデータは次のとおりです。左からそれぞれ、重さ[g]、比熱[J/kg・K]、光の反射率[%](玉はすべて半径1.0cmの球体)だそうです。

69.613 129.070 52.111
70.670 128.161 52.446
72.303 128.450 52.853
73.759 127.522 51.786
74.085 129.067 53.352

早速、散布図行列を表示してデータの分布を見ます。図の座標軸V1はデータの1列目に相当(V2、V3も同様)します。
V2とV3の散布図を見ると、なんとなく左下隅のデータが外れ値っぽいですが、これだけではなんとも言えません。
f:id:sudillap:20130324204740p:plain

つぎに3次元散布図を描いてみます(以下は2方向から見た図)。散布図を回転させてみるとわかるのですが、外れ値の候補がいくつか見つかります。ちなみに正解を先に言うと、87番目のデータが外れ値です。
f:id:sudillap:20130324205035p:plainf:id:sudillap:20130324205030p:plain

散布図からでも外れ値の候補をいくつか見つけることはできますが、これだけでは面白くありません。そこで、本記事では代表的な機械学習を使って外れ値を見つけます。

なお、散布図を見ると分かりますが、次元ごとにデータが取る範囲が異なります。このままでは色々面倒なので、基本的に標準化後のデータ分析をしています。
Tech総研の記事でnaoya_tさんは

【註】3要素の分布がそれぞれ異なるのであらかじめ正規化するのが望ましいですがそのまま突っ込んでも結果が変わらないようなデータにしてあります。

と書かれていますが、(問題にも依りますが)一般的には標準化した方がいいと思われます。
このデータではあまり問題は起こりませんが、石の重さをグラム[g]ではなくKgで測定されていた場合を考えると、いろいろ面倒なことが起こりそうですね。

まずは、データ読み込みのRスクリプトを書いておきます。読み込んだデータを標準化し、データフレームgoldに格納しています。

gold <- read.table("hundred.txt", header=FALSE)
gold <- as.data.frame(scale(as.matrix(gold), center = TRUE, scale = TRUE))