遺伝学・遺伝統計学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

phomパッケージ in R

R トポロジー トポロジー統計 Persistent topology phom
  • CRANのphomページ
  • phomパッケージが準拠しているのは、こちらのToplogy and Data by G. Carlsson
  • phomパッケージがすること
    • 空間の点の集合として得られたデータを用いる
    • (単体的)複体を構成
    • 複体のホモロジーを計算する
  • phomパッケージの複体構成法
    • Vietoris-Rips
    • Lazy-Witness
  • Vietoris-Rips
    • ２点間の「距離」が指定値以下であるような２点間のみを連結する
    • できた森から(単体的)複体の条件に合うように、２次元以上の複体の構成を決める
  • Lazy-Witness
    • Vietoris-Ripsよりノイズの影響に対して強い
    • 比較的少数の代表点がトポロジーを決めるのでは、という考えに基づいた構成法
    • ある閾値によって、「入っている・入っていない」の判断をしながら構成する
    • 構成しながら、「すでに構成体に入っている点」は相手にせず、作業行程でつながっていない点を入れていく
  • Homologyを計算する
  • Homology計算では
    • dim 0は位相空間の数、dim 1は、穴、dim 2はさらに高次の穴…
    • interval値が大きいほど「本物らしさ」が大きい
  • ２次元円の中空の度合いを変化させたときの、次元別のintervalの値の変化を見よう
    • 右に行くにつれ、円周になっていく
    • 最上段はその様子
    • 中段は、dim 0 のinterval
    • 下段は、dim 1 のinterval
    • 一つながりの位相空間であることは、dim 0 のintervalが比較的大きいものが１つだけある(しかも、中空の程度に寄らず)ことが示されている
    • dim 1(円)が一つあることが、右に行くほど確実になっていることがdim 1のintervalの大きくなることで示されている


  • ３次元球とその中空化。点の数が少ないと、３次元球の中空というよりも、２次元円がたくさんある(３次元球面上の円がたくさん見つかる)ような推定結果になる(Npt=100とかだと)
  • 点の数を増やすとdim 2 が１つ見つかるようになる

library(phom)

# 球、中空球
k<-2 # 空間の次元
Npt<-100 # 点の数

# 外側の円の半径
r1<-1
# パラメタ化する内側の円の半径
r2s<-seq(from=0,to=1,length=5)
# pHom()関数での閾値
max.dim<-2
max.f<-0.7
intervals<-list()
par(ask=FALSE)
par(mfcol=c(3,5))
for(i in 1:length(r2s)){
	r2<-r2s[i]
	rs<-runif(Npt,min=r2,max=r1)
	d<-matrix(rnorm(k*Npt),ncol=k)
	d<-d/sqrt(apply(d^2,1,sum))*rs
	intervals[[i]]<-pHom(d,max.dim,max.f,mode="lw",metric="euclidean")
	plot(d)
	for(j in 0:max.dim){
		title.st<-paste("dim",j)
		plotBarcodeDiagram(intervals[[i]],j,max.f,title=title.st)
	}

}
par(mfcol=c(1,1))

# 外側の円の半径
r1<-1
# パラメタ化する内側の円の半径
r2s<-seq(from=0,to=1,length=5)
# pHom()関数での閾値
max.dim<-2
max.f<-0.8
intervals<-list()
par(ask=FALSE)
par(mfcol=c(2,5))
for(i in 1:length(r2s)){
	r2<-r2s[i]
	rs<-runif(Npt,min=r2,max=r1)
	d<-matrix(rnorm(k*Npt),ncol=k)
	d<-d/sqrt(apply(d^2,1,sum))*rs
	intervals[[i]]<-pHom(d,max.dim,max.f,mode="lw",metric="euclidean")
	plot(d)
	for(j in 0:max.dim){
		title.st<-paste("dim",j)
		#plotBarcodeDiagram(intervals[[i]],j,max.f,title=title.st)
		

	}
	plotPersistenceDiagram(intervals[[i]],max.dim,max.f,title=title.st)
}
par(mfcol=c(1,1))

# ３次元にしてみる
library(phom)

# 球、中空球
k<-3 # 空間の次元
Npt<-500 # 点の数

# 外側の円の半径
r1<-1
# パラメタ化する内側の円の半径
r2s<-seq(from=0,to=1,length=5)
# pHom()関数での閾値
max.dim<-k
max.f<-0.6
intervals<-list()
par(ask=FALSE)
par(mfcol=c(2,5))
for(i in 1:length(r2s)){
	r2<-r2s[i]
	rs<-runif(Npt,min=r2,max=r1)
	d<-matrix(rnorm(k*Npt),ncol=k)
	d<-d/sqrt(apply(d^2,1,sum))*rs
	intervals[[i]]<-pHom(d,max.dim,max.f,mode="lw",metric="euclidean")
	plot(d)
	for(j in 0:max.dim){
		title.st<-paste("dim",j)
		#plotBarcodeDiagram(intervals[[i]],j,max.f,title=title.st)
		

	}
	plotPersistenceDiagram(intervals[[i]],max.dim,max.f,title=title.st)
}
par(mfcol=c(1,1))