Brownian excursionとBrownian contunuum random tree

• 0からスタートして0に戻るブラウン運動はBrownian bridge
• 正の値だけをとるBrownian bridgeはBrownian excursion(ブラウン散歩)
• Rでシミュレーションするには、e1071パッケージに実装されたBrownian bridgeを何度も回して、excursion条件を満足するのができたらそれを返す、というやり方(参考：こちら)
• その上で、ブラウン散歩の値についてという距離定義を入れることによって、ブラウン散歩から木が得られるそれがBrownian continuum random tree

Bt <- rnorm(500, 0, 1)
Bt <- cumsum(Bt)
Bt <- rnorm(500, 0, 1)
Bt <- Bt[Bt > 0 & Bt < 1]
Bt <- cumsum(Bt)
Bt <- rnorm(500, 0, 1/sqrt(500))
  Bt <- cumsum(Bt)
  Bt <- ts(Bt, start = 1/500, frequency = 500)
  ts(abs(Bt - time(Bt) * as.vector(Bt)[500]), start = 1/500, 
     frequency = 500)
     
library(e1071)
succeeded <- FALSE
while(!succeeded)
{
  bridge <- rbridge(end = 1, frequency = 500)
  succeeded=all(bridge>=0)
}
plot(bridge)


while(!succeeded)
{
  bridge <- rbridge(end = 1, frequency = 500)
  succeeded=all(bridge>=0)||all(bridge<=0)
}
bridge = abs(bridge)
plot(bridge)

# library(e1071)
my.rexcursion <- function(frequency = 1000, end = 1){
	succeeded <- FALSE
	while(!succeeded)
	{
		bridge <- rbridge(end = end, frequency = frequency)
		succeeded=all(bridge>=0)||all(bridge<=0)
	}
	return(c(0,c(abs(bridge))))
}
my.rwiener <- function(frequency = 1000, end = 1){
	c(0,cumsum(rnorm(end * frequency)/sqrt(frequency)))
}
br.ex <- my.rexcursion(frequency=100)
plot(br.ex)
length(br.ex)


fr <- 10
t <- (1:fr)/fr
br.ex <- my.rexcursion(frequency=fr)
br.z <- my.rwiener(frequency=fr)

library(rgl)
library(ape)
plot3d(t,br.ex,br.z,type="l")

my.Rtree.dist <- function(g,s,t){
	g[s] + g[t] - 2 * min(g[s:t])
}

my.Rtree.dist.mat <- function(g){
	L <- length(g)
	D.mat <- matrix(0,L,L)
	
	for(i in 1:(L-1)){
		for(j in (i+1):L){
			D.mat[i,j] <- D.mat[j,i] <- my.Rtree.dist(g,i,j)
		}
	}
	return(D.mat)
}
D.mat <- my.Rtree.dist.mat(br.ex)

tr <- nj(D.mat)

plot(nj(my.Rtree.dist.mat(br.ex)),type="u",show.tip.label=FALSE)

my.exUpDown <- function(br.ex){
	L <- length(br.ex)
	d <- diff(br.ex)
	#print(d)
	d.sign <- sign(d)
	#print(d.sign)
	d.d.sign <- diff(d.sign)
	dd.positive <- which(d.d.sign>0) + 1
	dd.negative <- which(d.d.sign<0) + 1
	#print(dd.positive)
	#print(dd.negative)
	col <- rep(1,L)
	col[dd.positive] <- 2
	col[dd.negative] <- 3
	plot(br.ex,col=col,pch=20,type="b")

}
my.exUpDown(br.ex)


my.ex2tr(br.ex)