Brownian excursion

• Brownian excursionというのがある
• どういうことかというと、１次元数直線の原点から正の方向に出発し、酔歩する。ただし、原点に戻ってしまったら終わり、というルール
• このブラウン運動的な小旅行(excursion)、ブラウン散歩は、「木」のルートから『木の外側』をたどって歩くのと同じである、という話
• ブラウン散歩の位置座標を表す関数をgとすれば、gの値が増えているときというのは、木の枝が成長しているとき
• gが小さくなっているときは、ルートに向かって歩いているとき
• gが減少から増加に転ずるときに、分岐することにする
• そうすると、この関数gによって木が一意に決まる
• gの最大値は、木のルートからの最遠点になる
• gが単調増加して、ピークになって単調減少するような場合には、ルートから枝が一本伸びているだけの構造で、それを往復することになる
• これを用いて離散ブラウン散歩データから木構造を作ろう…とがんばってみて、「できた！」と思ったのだが、そしてできていたのだが
• できた瞬間に、馬鹿みたいに簡単な方法があることに気づくことは、ままあって…
• ブラウン散歩の各点の間に距離を計算する。ただし、その距離は、時刻s,tの位置、g(s),g(t)と、時刻s~tの間のもっとも値の小さいgの値、min(g,s,t)とを使って、g(s)+g(t)-2min(g,s,t)で定まる
• こうすると、すべての点ペアの木構造上の距離になっている
• すべての点ペアの距離が行列で作れるので、neighbor joiningを使えば、節ノードも適切に発生しつつ、木構造が得られる(今の場合、もとの距離データが正確に木になっているので、必ず正解が返ってくる)
• それだけ！

library(ape)

my.tree.dist <- function(g){
	n <- length(g)
	ret <- matrix(0,n,n)
	for(i in 1:(n-1)){
		for(j in (i+1):n){
			ret[i,j] <- ret[j,i] <- g[i]+g[j]-2*min(g[i:j])
		}
	}
	ret
}
# 自由に酔歩し、始点終点を0にしつつ、途中で負になったら正に反転することでブラウン散歩の離散データを作る
n <- 200
t <- c(0,sort(runif(n)),1)

s <- rnorm(length(t))
S <- cumsum(s)
S <- S-S[1]
S <- S-t*S[length(S)]
Splus <- abs(S)

d.mat.2 <- my.tree.dist(Splus)
nj.out2 <- nj(d.mat.2)
plot(nj.out2,"f",show.tip.label=FALSE)

• 延々と書いたコードもいつか役に立つかもしれないので貼っておく

# 増加・減少を交互に起こさせるために、上記のブラウン散歩で増加が連続したり
# 減少が連続したりしたら、途中に長さ0のステップを入れる
t2 <- t[1:2]
Splus2 <- Splus[1:2]
Splus.diff <- diff(Splus)
for(i in 2:length(Splus.diff)){
	if(sign(Splus.diff[i-1])==sign(Splus.diff[i])){
		t2 <- c(t2,t2[length(t2)],t[i+1])
		Splus2 <- c(Splus2,Splus2[length(Splus2)],Splus[i+1])
	}else{
		t2 <- c(t2,t[i+1])
		Splus2 <- c(Splus2,Splus[i+1])
	}
}
matplot(t,cbind(S,Splus),type="l")
points(t2,Splus2,col=3)


Splus.diff2 <- diff(Splus2)
# 増減が必ず交互なので、増加〜エッジ相当と
# 減少、後戻りして、次の分岐点を決める
# この２通りに分けられる
#Elen <- matrix(Splus.diff2,nrow=2)[1,]
N <- length(Splus2)
# エッジの長さなど
Edges <- t(matrix(Splus2[-N],nrow=2))
Elen <- Edges[,2]-Edges[,1]
# 個々のエッジの分岐点がどのエッジの上にあるかを調べる
Branch <- rep(0,length(Elen))
Branch <- (1:length(Elen))-1

for(i in 2:length(Branch)){
	Inside <- which((Edges[,1]-Edges[i,1])*(Edges[,2]-Edges[i,1]) < 0)
	if(length(which(Inside < i))>0){
		Branch[i] <- max(Inside[which(Inside < i)])
	}
	
}
plot(t2,Splus2,type="b")
abline(h=Edges[,1],col=2)
Branch
# 一つずつエッジを調べる
# そのエッジから生えるエッジを確認しつつ
# そのエッジを分割して、節ノードを発生する
# エッジ分割をしたら、元のエッジは消去する
vs <- c(1,2)
es <- matrix(c(1,2),nrow=1)
# 節を考慮しながら、個々のエッジの長さを格納
elen <- c(Elen[1])
# その差表を進めるために、追加エッジの始点終点ノードのIDがいくつになったかを格納する
stend.E <- matrix(0,nrow(Edges),2)
stend.E[1,] <- c(1,2)

for(i in 1:length(Elen)){
	selected <- which(Branch==i)
	if(length(selected)==0){
		next
	}
	selected.st <- Edges[selected,1]
	vs.on.this.edge <- sort(unique(c(Edges[i,],selected.st)))
	tmp.n <- length(vs.on.this.edge)-2
	new.vs <- c()
	if(tmp.n>0){
		new.vs <- length(vs)+(1:tmp.n)
	}
	vid.on.this.edge <- c(stend.E[i,1],new.vs,stend.E[i,2])
	if(tmp.n>0){
		vs <- c(vs,new.vs)
		remove.edge <- which(es[,1]==stend.E[i,1] & es[,2]==stend.E[i,2])
		new.edge <- cbind(c(es[remove.edge,1],new.vs),c(new.vs,es[remove.edge,2]))
		new.edge.len <- diff(sort(unique(c(Edges[i,],selected.st))))
		es <- rbind(es,new.edge)
		es <- es[-remove.edge,]
		elen <- c(elen,new.edge.len)
		elen <- elen[-remove.edge]
	}
	for(j in 1:length(selected)){
		#if(Elen[selected[j]]!=0){
			tmp.v <- which(vs.on.this.edge == selected.st[j])
			tmp.new.v <- length(vs)+1
			vs <- c(vs,tmp.new.v)
			es <- rbind(es,c(tmp.v,tmp.new.v))
			stend.E[selected[j],] <- c(tmp.v,tmp.new.v)
			elen <- c(elen,Elen[selected[j]])
		#}

	}
	tmp.g <- graph.edgelist(es,directed=FALSE)

}

g3 <- graph.edgelist(es,directed=FALSE)
E(g3)$weight <- elen
#plot(g3,vertex.size = 0.1,vertex.label="")
plot(g3)
# グラフ上の距離をノードペアワイズで取り出せば
d.mat <- shortest.paths(g3,weights=elen)
# nj法で木に戻せる
nj.out <- nj(d.mat)
# 同じ木の絵が描ける
plot(nj.out)