虫が部屋の中を回り続けること - ryamadaのコンピュータ・数学メモ

部屋にハエとかハチとかが入ってくると、飽きもせずにブーンと回り続けることがよくある(というかいつもそう)
どうしてか(理由というよりは、いかにして)
部屋という制約空間(その内側にしか存在できない)
その中で、運動ベクトルがランダムに変わる
ただし、運動ベクトルは勝手に変わることはできなくて、運動ベクトルにランダムな加速度を加えることができるものとする
感覚はあるので、ぶつかる前にランダム加速度を入れる
これは、「制約空間内」の「加速度酔歩」(通常の酔歩は速度ベクトルをランダムに変えるという意味で「速度酔歩」)
加速度酔歩は速度酔歩に比べて「滑らか」であることは、こちらに書いた

こんな感じに
生物の運動がこんな感じだけれど、色々な生物の軌道(複数因子の濃度変化とか)も「因果関係」っぽかったりするけれど、ただの「制約空間内加速度酔歩」であるだけかもしれない。これを「検定」する方法があればよい

# 酔歩データの解釈法

## 格子酔歩関数

```{r}
rw.lattice <- function(n.step,p,x.init = rep(0,p),type=1,m=0,sd=1){
  if(type==1){
    v <- matrix(sample(c(-1,1),n.step*p,replace=TRUE),ncol=p)
  }else if(type==2){
    x. <- matrix(sample(c(-1,1),n.step*p,replace=TRUE),ncol=p)
    s <- sample(1:p,n.step,replace=TRUE)
    v <- matrix(0,n.step,p)
    for(i in 1:p){
      a <- which(s==i)
      v[a,i] <- x.[a,i]
    }
  }else if(type==3){
    v <- matrix(rnorm(n.step*p,m,sd),ncol=p)
    v <- v/sqrt(apply(v^2,1,sum))
  }else if (type==4){
    v <- matrix(rnorm(n.step*p,m,sd),ncol=p)
  }
  v[1,]<-0
  t(t(apply(v,2,cumsum)+x.init))
}
```

## 正規分布による中心束縛のある酔歩関数

```{r}
rw.center.normal <- function(n.step,p,x.init=rep(0,p),m=rep(0,p),sd=rep(1,p),type=1,L.step=1){
  v <- matrix(rnorm(n.step*p),ncol=p)
  v <- t(t(v) * sd + m)
  if(type==1){
    L.step <- rep(L.step[1],n.step)
  }else if(type==2){
    L.step <- abs(rnorm(n.step))
  }else if(type==3){
    # L.step <- L.step
  }
  x <- matrix(0,n.step,p)
  for(i in 2:n.step){
    u <- v[i,] - x[i-1,]
    u. <- L.step[i] * u/sqrt(sum(u^2))
    x[i,] <- x[i-1,] + u.
  }
  t(t(x) + x.init)
}
```

## 球面上酔歩関数

```{r}
rw.sphere.surface <- function(n.step,p,phis=0.1){
  if(length(phis)==1){
    phis <- rep(phis,n.step)
  }
  sphere.brownian(n.step,p,phis)
}
```

## 無限空間の滑らかな酔歩関数

### 正規乱数によるk+1階微分ドライブ
```{r}
rw.norm.kdiff <- function(n.step,p,k,x.init = matrix(0,k+1,p),m=rep(0,p),sd=rep(10^(-p-1),p)){
  ret <- array(0,c(k+1,n.step,p))
  tmp <- matrix(rnorm(n.step*p),ncol=p)
  ret[k+1,,] <- t(t(tmp)*sd+m)
  for(i in k:1){
    tmp <- rbind(x.init[i,],ret[i+1,1:(n.step-1),])
    ret[i,,] <- apply(tmp,2,cumsum)
  }
  ret
}
n.step <- 10^5
p <- 3
k <- 4
par(mfrow=c(2,2))
x <- rw.norm.kdiff(n.step,p,k)
for(i in 1:k){
  matplot(x[i,,],type="l",main=i-1)
}
par(mfcol=c(1,1))
```

## n階差分をとる関数
```{r}
diffwalk.p <- function(x,k){
  n.step <- nrow(x)
  p <- ncol(x)
  ret <- array(0,c(k+1,n.step,p))
  ret[1,,] <- x
  for(i in 1:k){
    ret[i+1,1:(n.step-i),] <- apply(ret[i,1:(n.step-i+1),],2,diff)
  }
  ret
}
```


## 無限空間酔歩とn階差分
```{r}
n.step <- 10^4
p <- 2
x.list <- list()
for(i in 1:4){
  x.list[[i]] <- rw.lattice(n.step,p,type=i)
}

xlim <- ylim <- range(unlist(x.list))
plot(x.list[[1]],type="l",col=1,xlim=xlim,ylim=ylim)
for(i in 2:4){
  points(x.list[[i]],type="l",col=i)
}
```
無限空間酔歩(格子酔歩、正規乱数酔歩)の場合、
1階差分と2階差分には関係があり
0階差分と1階差分には関係がない
```{r}
k <- 4
diffs <- lapply(x.list,diffwalk.p,k)
```
```{r}
par(mfrow=c(2,2))
for(i in 1:4){
  plot(diffs[[i]][2,,1],diffs[[i]][3,,1],main="1 vs. 2")
}
par(mfrow=c(1,1))
```
```{r}
par(mfrow=c(2,2))
for(i in 1:4){
  plot(diffs[[i]][1,,1],diffs[[i]][2,,1],main="0 vs. 1")
}
par(mfrow=c(1,1))
```

```{r}
par(mfrow=c(2,2))
for(i in 1:4){
  plot(diffs[[i]][1,,1],diffs[[i]][3,,1],main="0 vs. 2")
}
par(mfrow=c(1,1))
```
```{r}
par(mfrow=c(2,2))
for(i in 1:k){
  plot(diffs[[2]][1,,1],diffs[[2]][i+1,,1],main=paste("0 vs. ",i))
}
par(mfrow=c(1,1))
```
## 正規分布による中心束縛のある酔歩とn階差分

```{r}
n.step <- 10^4
p <- 2
x.list <- list()
for(i in 1:2){
  x.list[[i]] <- rw.center.normal(n.step,p,type=i)
}

xlim <- ylim <- range(unlist(x.list))
plot(x.list[[1]],type="l",col=1,xlim=xlim,ylim=ylim)
points(x.list[[2]],type="l",col=2)
```
```{r}
k <- 4
diffs <- lapply(x.list,diffwalk.p,k)
```
```{r}
par(mfrow=c(2,2))
for(i in 1:2){
  plot(diffs[[i]][2,,1],diffs[[i]][3,,1],main="1 vs. 2")
}
par(mfrow=c(1,1))
```
```{r}
par(mfrow=c(2,2))
for(i in 1:2){
  plot(diffs[[i]][1,,1],diffs[[i]][2,,1],main="0 vs. 1")
}
par(mfrow=c(1,1))
```

```{r}
par(mfrow=c(2,2))
for(i in 1:2){
  plot(diffs[[i]][1,,1],diffs[[i]][3,,1],main="0 vs. 2")
}
par(mfrow=c(1,1))
```
```{r}
par(mfrow=c(2,2))
for(i in 1:k){
  plot(diffs[[2]][1,,1],diffs[[2]][i+1,,1],main=paste("0 vs. ",i))
}
par(mfrow=c(1,1))
```

## 球面酔歩のn階差分
```{r}
n.step <- 1000
p <- 3
#x <- rw.sphere.surface(n.step,p,phis=abs(rnorm(n.step)))
x <- rw.sphere.surface(n.step,p,phis=0.1)
x.list <- list(x)
xlim <- ylim <- range(unlist(x.list))
matplot(x.list[[1]],type="l")
```
```{r,rgl=TRUE}
plot3d(x.list[[1]],type="l")
```
```{r}
k <- 4
diffs <- lapply(x.list,diffwalk.p,k)
```
```{r}
par(mfrow=c(2,2))
for(i in 1:k){
  plot(diffs[[1]][1,,1],diffs[[1]][i+1,,1],main=paste("0 vs. ",i))
}
par(mfrow=c(1,1))
```

```{r}
p <- 3
n.step <- 10^4
np <- n.step*p
n.neighbor <- 10^2

k <- 2

X <- array(0,c(k+1,n.step,p))
# X[,1,] should be given
R <- array(rnorm(np*n.neighbor),c(n.step,n.neighbor,p))*0.01

for(i in 2:n.step){
  tmp.X <- array(0,c(k+1,n.neighbor,p))
	tmp.X[k+1,,] <- R[i,,]
	for(j in k:1){
		tmp.X[j,,] <- t(t(tmp.X[j+1,,]) + X[j,i-1,])
	}
	tmp.val <- apply(tmp.X[1,,]^2,1,sum)
	tmp.pr <- exp(-(tmp.val-min(tmp.val)))
	tmp.pr <- tmp.pr/sum(tmp.pr)
	s <- sample(1:n.neighbor,1,prob=tmp.pr)
	X[1:k,i,] <- tmp.X[1:k,s,]
}

plot(X[1,,1:2],type="l")
```
```{r,rgl=TRUE}
plot3d(X[1,,],type="l")
```

```{r}

rw.weighted.sample <- function(X,k,fx){
  tmp <- fx(X[k,,])
	tmp <- tmp/sum(tmp)
	sample(1:length(X[k,,1]),1,prob=tmp)
}

w.fx.norm <- function(X,m=rep(0,ncol(X)),sd=rep(1,ncol(X))){
	tmp <- (t(X)-m)/sd
	tmp2 <- apply(tmp^2,2,sum)
	exp(-(tmp2-min(tmp2)))
}
```
```{r}
p <- 2
n.step <- 10^3
np <- n.step*p
n.neighbor <- 10^2

k <- 2

X <- array(0,c(k+1,n.step,p))
# X[,1,] should be given
R <- array(rnorm(np*n.neighbor),c(n.step,n.neighbor,p))*0.01

for(i in 2:n.step){
  tmp.X <- array(0,c(k+1,n.neighbor,p))
	tmp.X[k+1,,] <- R[i,,]
	for(j in k:1){
		tmp.X[j,,] <- t(t(tmp.X[j+1,,]) + X[j,i-1,])
	}
	#tmp.val <- apply(tmp.X[1,,]^2,1,sum)
	#tmp.pr <- exp(-(tmp.val-min(tmp.val)))
	#tmp.pr <- tmp.pr/sum(tmp.pr)
	#s <- sample(1:n.neighbor,1,prob=tmp.pr)
	s <- rw.weighted.sample(tmp.X,1,w.fx.norm)
	X[1:k,i,] <- tmp.X[1:k,s,]
}

plot(X[1,,],type="l")
```
```{r}
rw.weighted.sample <- function(X,fx){
  tmp <- fx(X)
	if(sum(tmp)==0){
		return(NULL)
	}
	tmp <- tmp/sum(tmp)
	return(sample(1:length(X[,1]),1,prob=tmp))
}


w.fx.norm <- function(X,m=rep(0,ncol(X)),sd=rep(1,ncol(X))){
	tmp <- (t(X)-m)/sd
	tmp2 <- apply(tmp^2,2,sum)
	exp(-(tmp2-min(tmp2)))
}
w.fx.sphere <- function(X,m=rep(0,ncol(X)),r=1){
	tmp <- t(X)-m
	tmp2 <- apply(tmp^2,2,sum)
	as.numeric(tmp2 < r^2)
}
w.fx.simplex <- function(X,m=rep(0,ncol(X)),r=1){
  tmp <- t(X)-m
  tmp0 <- tmp >=0
  #tmp1 <- tmp <=r
  tmpsum <- apply(tmp,2,sum) <= r
  #tmp01 <- rbind(tmp0,tmp1,tmpsum)
  tmp01 <- rbind(tmp0,tmpsum)
	as.numeric(apply(tmp01,2,prod))
}

my.rw <- function(n.step,p,k,n.neighbor,u,scale=0.01,fx){
  X <- array(0,c(k+1,n.step,p))
	y.hx <- rep(1,n.step)
  np <- n.step*p
  #X[,1,] <- X.init
	R <- array(rnorm(np*n.neighbor),c(n.step,n.neighbor,p))*scale

	for(i in 2:n.step){
		tmp.X <- array(0,c(k+1,n.neighbor,p))
		loop <- TRUE
		y <- 1
		while(loop){
			tmp.X[k+1,,] <- R[i,,]*y
			for(j in k:1){
  			#tmp.X[j,,] <- t(t(tmp.X[j+1,,])*y + X[j,i-1,])
  			tmp.X[j,,] <- t(t(tmp.X[j+1,,]) + X[j,i-1,])
			}
			
			s <- rw.weighted.sample(tmp.X[u,,],fx)
			if(is.null(s)){
				y <- y*2
			}else{
				X[,i,] <- tmp.X[,s,]
  			y.hx[i] <- y
  			#y.hx[i] <- 1
				loop <- FALSE
			}
		}
	
	}
	return(list(X=X,t=cumsum(y.hx)))
}
```



```{r}
n.step <- 10^3
n.neighbor <- 10^2
p <- 3
k <- 2
u <- 1
fx <- w.fx.norm
X <- my.rw(n.step,p,k,n.neighbor,u,fx=fx)
fx <- w.fx.sphere
X.init <- matrix(0,k+1,p)
X.init[1,] <- rep(1/p,p)
X2 <- my.rw(n.step,p,k,n.neighbor,u,scale=0.01,fx=fx)

plot3d(X2$X[1,,],type="l")
matplot(X2$t,X2$X[1,,],type="l")
```
```{r}
n.step <- 10^4
n.neighbor <- 10^2
p <- 2
k <- 2
u <- 1
fx <- w.fx.simplex
X3 <- my.rw(n.step,p,k,n.neighbor,u,fx=fx)
X3.1 <- X3$X[1,,]
X3.1 <- cbind(X3.1,1-apply(X3.1,1,sum))
cv <- CategoryVector(p+1)
plot(X3.1 %*% cv,type="l",asp=1,xlim=range(cv[,1]),ylim=range(cv[,2]))
points(cv[c(1,2,2,3,3,1),],type="l",col=2)
```