2014-06-21

Dirichlet diffusion trees

先日の記事の続き
作り方
原点から第一の要素がランダムウォークする。単位時間後とすれば、ランダムウォークするので、平均0、分散がある値(1としてもよいだろう)の点(t=1,x)に到達する
第二以降の要素は、原点から出発し、「分離」イベントが起きるまでは、すでにある道をたどる。複数の道があるときはどれを選ぶかは等確率。したがって、i番目の要素は、i-1個の先行要素のどれかを「親」としてついて歩く
「分離イベント」がおきると、そこを始点にランダムウォークする。時刻t (0<= t <= 1)で分離したとすると、その始点からは、平均0、分散が残時間(1-t)に比例した値の正規分布に従う分だけ、ずれて(t=1,x')に到達する
「分離」イベントは、「ついて歩いて」いるときに「一緒に歩いている先行要素」の数が多ければ多いほど「離れにくい」。先ほど「親」を一人、選んだが、実際には、「一緒に歩いている先行要素」はみんな親で、「すべての親と別れないといけない」から
このようなルールを入れると、道の先に行けば行くほど、分離しやすくなる
一方、分離イベントを「道連れ人数」によらず、時刻tの関数としてやることで、分離のパターンを作ることもできる。たとえば、時刻0で必ず分離する(ただの独立ランダムウォーク)とか、時刻0で定確率で分離するがそれ以外は分離しない(複数カテゴリに分かれる)とか
実際の分離は、この時刻の関数と、道連れ人数のしがらみで決まる
以下は、この「関数」と「しがらみ」を無視して(その部分が面倒なので)、分離木を順次成長させる部分だけをアルゴリズムにしたもの

library(igraph)
# 要素数
n.pt <- 17
# 最初の要素のグラフは原点とその要素をつないだだけ
el <- matrix(c(1,2),nrow=1)
# 葉ノードが要素なので、そのIDを格納
Leaves <- c(2)
# ルートノード
Root <- c(1)
# グラフオブジェクト化
g <- graph.edgelist(el)
# エッジの重みは、その始点・終点ノードが対応する時刻の差。最初のエッジのそれは１
E(g)$weight <- c(1)
plot(g)
g.hx <- list()
g.hx[[1]] <- g
par(mfcol=c(2,2))
# 第２要素から順にグラフを成長させる
for(i in 2:n.pt){
# 「親」を選ぶ
	tmp.leave <- Leaves[sample(1:length(Leaves),1)]
	#tmp.leave <- sample(Leaves,1)
# 原点から親へのパスをエッジリストでチェック
	tmp.path <- get.shortest.paths(g, Root, to=tmp.leave, mode = c("out"), weights = NULL, output=c("epath"))
# 現時点のエッジリスト
	el.now <- get.edgelist(g)
# ノード数
	vnum <- vcount(g)
# エッジの重み
	e.weight <- E(g)$weight
# 親へのパスのエッジを一つ等確率で選ぶ
	tmp.edge <- tmp.path$epath[[1]][sample(1:length(tmp.path$epath[[1]]),1)]
# 選ばれたエッジの始点終点ノードの確認
	tmp.edge.stend <- el.now[tmp.edge,]
# 選ばれたエッジの途中に新規ノードを作り、そこから新しい葉ノードへのエッジをはやす
	el.next <- el.now[-tmp.edge,]
	el.next <- rbind(el.next,matrix(c(tmp.edge.stend[1],vnum+1,vnum+1,tmp.edge.stend[2],vnum+1,vnum+2),byrow=TRUE,ncol=2))
	Leaves <- c(Leaves,vnum+2)
	g <- graph.edgelist(el.next)
	weight.now <- e.weight[tmp.edge]
	e.weight.next <- e.weight[-tmp.edge]
# 新しい３エッジの重み、分割してできたエッジの重みの和は元のエッジの重みであり、葉へのパスの重みの和は１であることを反映
	r <- runif(1)
	e.weight.next <- c(e.weight.next, weight.now * c(r,1-r),1)
	E(g)$weight <- e.weight.next
	tmp.path <- get.shortest.paths(g, Root, to=vnum+2, mode = c("out"), weights = NULL, output=c("epath"))
	tmp.path.weight <- e.weight.next[tmp.path$epath[[1]]]
	last.weight <- 1-sum(tmp.path.weight[-length(tmp.path.weight)])
	e.weight.next[length(e.weight.next)] <- last.weight
	E(g)$weight <- e.weight.next
	g.hx[[i]] <- g
	
	col <- rep(5,vnum+2)
	col[Leaves] <- 3
	col[1] <- 6
# ルートからの時刻を計算
	dist.from.Root <- shortest.paths(g,1)
# エッジの長さに応じて、ドリフト距離を正規乱数で作成	
	tmp.dist <- rep(0,length(e.weight.next))
	for(j in 1:length(tmp.dist)){
		tmp.dist[j] <- rnorm(1,0,sqrt(e.weight.next[j]))
	}
	tmp.g <- graph.edgelist(el.next)
	E(tmp.g)$weight <- tmp.dist-min(tmp.dist)
	dist.random <- shortest.paths(tmp.g,1)
	mean.dist.random <- mean(dist.random)
	dist.random <- dist.random -mean.dist.random
	dist.random[1] <- 0
	coords <- matrix(c(dist.random,dist.from.Root),ncol=2)
	plot(g.hx[[i]],layout=coords,vertex.color = col)
}
par(mfcol=c(1,1))