トポロジカルデータアナリシス：TDAパッケージを使ってみる

導入

とあることがきっかけで、とっても久しぶりにRでTDA（トポロジカルデータアナリシス）をしてみました。だいぶいろんなことを忘れていたので、単純な例を使ったメモを残しておきます。

トポロジカルデータアナリシスとは

とてもざっくりいうと、位相幾何学という数学の知見をつかって、データから「形」の情報を抽出するような手法になっています。

導入は、過去の記事にも載っていますので、こちらも参照していただけるとうれしいです。
tekenuko.hatenablog.com

「形」の情報の抽出のしかたですが、データ点のまわりにある半径の球をとり、その球たちの重なりによって点の間に線を引く、といった方法でデータ点から図形を見立てます。ただし、データに付随するノイズなどにロバストな形で図形的な情報を抜き出すために、半径を連続的に変化させた場合の構造の変化を見るといったテクニックを用います。

これらの簡単な解説は、一応過去の記事でも紹介しています（わかりづらかったらすいません）。
tekenuko.hatenablog.com

半径を連続的に変化させた場合、図形が生まれて、亡くなる…という動きが起こります。その一連を動きをプロットしたものが、パーシステント図と呼ばれます。今回、具体例を用いてパーシステント図がどうなるかを説明したいと思います。

使用パッケージ

ここでは、Rでトポロジカルデータアナリシスを行うためのパッケージである{TDA}の関数を使います。以下のようにしてインストールします。Rcppなど、いくつか依存パッケージがあり、それらのバージョンなどでエラーが出ることがあります。その場合は、依存パッケージを入れ直すなどするなどの対処をする必要があります（うろ覚え）。

# TDAパッケージ
if(!require(TDA)){
  install.packages("TDA", quiet = TRUE)
}
require(TDA)

使い方（マニュアル）はこちら：
https://cran.r-project.org/web/packages/TDA/vignettes/article.pdf

日本語の解説だと、ALBERTさんの記事が参考になります。
blog.albert2005.co.jp

サンプルデータ

大きな円と小さな円が繋がっている＋微小なノイズがのっているような人工データを生成します。ここでは、{TDA}パッケージにある、半径1の円上に沿ってランダムサンプリングを行う関数circleUnif()をうまく組み合わせます。

X1 <- circleUnif(100) * 2 + rnorm(100, 0, 0.03)
X2 <- circleUnif(50) * 0.3 + c(1.65,1.65) + + rnorm(50, 0, 0.02)
X <- rbind(X1, X2)
X <- data.frame(x = X[, 1], y = X[, 2])

可視化すると以下のようになります。

ggplot(X) + 
  geom_point(aes(x = x, y = y))

f:id:tekenuko:20171101001723p:plain

パーシステント図

パーシステント図を描きます。図形の構築法は、ヴィートリス・リップス複体を利用します。ヴィートリス・リップス複体を含めた構築法の概要は、過去記事を参照いただけると幸いです。
tekenuko.hatenablog.com

{TDA}パッケージでは、以下のような関数を使って実現します。

# パーシステント図の大元（パーシステントホモロジー）を計算
Diag_Cir <- ripsDiag(X = X, maxdimension = 1, maxscale = 5)

Xはデータ点、maxdimensionは何次元の三角形の移り変わりを見るかを指定するオプションです。今回は1としており、穴が空いているかの移り変わりを見ています。maxscaleは図を描くときの軸の最大値の値を指定するオプションです。パーシステント図を作るために動かす半径のパラメータは、本来無限大まで取れますが、変化は実質どこかで止まるはずで、その止まるパラメータ点をmaxscaleとしてスケーリングします。

可視化すると以下のようになります。{TDA}に内蔵されているplot関数を使うのが基本ですが、ここでは{ggplot2}で可視化を試みました。

# ggplot2で描くためにripsDiagの結果を整形
dimension_Cir <- Diag_Cir$diagram[,1]
birth_Cir <- Diag_Cir$diagram[,2]
death_Cir <- Diag_Cir$diagram[,3]
Diag_Cir_DF <- data.frame(cbind(Dimension = dimension_Cir, 
                              Birth = birth_Cir,
                              Death = death_Cir))

# plot
ggplot(Diag_Cir_DF) + 
  geom_point(aes(x = Birth, y = Death, colour = as.factor(Dimension)), size = 4) + 
  xlim(0,5) + ylim(0,5) + 
  geom_abline(slope = 1) + 
  labs(color = "Dimension") + 
  scale_color_manual(name="Dimension",
                     values=c("0"="black","1"="#F8766D"))

f:id:tekenuko:20171101002745p:plain

赤い点が穴に対応しており、対角線から距離が離れているほど半径の変化で安定的な穴であることを表しています。今回、赤い点は相対的に安定的なものが2点、不安定なものが1点出現しています。前者はデータ点を大きな円（対角線から大きく離れている点）と小さな円（対角線からちょっと離れている点）に対応しており、後者はノイズ起源でたまたま穴の空いた図形が生まれてすぐ消えていったものだと考えられます。このように、パーシステント図を利用すると、データ点の形に付随した量を抽出することができます。