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

データの生成〜時系列データ解析をやってみる

時系列
  • 参考書は、前記事でリンクさせていただいたブログ記事と

カオス時系列解析の基礎と応用

カオス時系列解析の基礎と応用

ランダムウォークとくりこみ群―確率論から数理物理学へ (新しい解析学の流れ)

ランダムウォークとくりこみ群―確率論から数理物理学へ (新しい解析学の流れ)

  • イントロ
    • ここで言う時系列解析というのは、p種類の変数が作るp次元空間中の点の移動軌跡に関する解析のことである
    • 時間軸に関して必ず１つの点を取り、１個未満でも２個以上でもない。また、軌跡は連続である
    • この「連続」制約は時系列解析の特徴であって、軌跡を1次元多様体とし、また、時間軸に関する１時点１点という制約が、１次元多様体の広がり方に制約を与えている
    • これを、m次元空間中の[tex:n
  • データを作ろう
    • 乱数によってのみ作られる軌道データを作成する
    • 乱数だけを使っているので、その部分については、「何のパターンやルールがなく」それ以外のところに生じていて、解析によって見出されて欲しいのは、「乱数生成時に採用した決まり」である
      • いくつかの関数を用意する
## 格子酔歩関数

```{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))
}
rw.weighted.sample.2 <- function(X,u,fx,x.current){
  tmp <- fx(X,u,x.current)
  if(sum(tmp)==0){
		return(NULL)
	}
	tmp <- tmp/sum(tmp)
	return(sample(1:length(X[u[1],,1]),1,prob=tmp))
}
```

```{r}
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))
}
```
```{r}

w.fx.gravity.2 <- function(X,u=3,x.current){
  m <- x.current[1,]
  if(sum(m^2)==0){
    return(rep(1,length(X[1,,1])))
  }
  tmp <- (X[u,,] %*% (-m))/(sqrt(apply(X[u,,]^2,1,sum))*sqrt(sum(m^2)))
  #tmp2 <- apply(tmp^2,2,sum)
  exp((tmp-min(tmp)))
}

w.fx.gravity.2.2 <- function(X,u=c(3,2),x.current){
  m <- x.current[1,]
  if(sum(m^2)==0){
    return(rep(1,length(X[1,,1])))
  }
  tmp <- (X[u[1],,] %*% (-m))/(sqrt(apply(X[u[1],,]^2,1,sum))*sqrt(sum(m^2)))
  tmp2 <- apply(X[u[2],,]^2,1,sum)
  #tmp2 <- apply(tmp^2,2,sum)
  exp((tmp-min(tmp))) * exp(-(tmp2-min(tmp2)))
}

```


```{r}
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)))
}
my.rw.2 <- 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.2(tmp.X,u,fx,X[,i-1,])
			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 <- c(3,2) # 加速度は原点方向に向きやすい、速度の絶対値は0になりやすい
fx <- w.fx.gravity.2.2

X5 <- my.rw.2(n.step,p,k,n.neighbor,u,scale=0.2,fx=fx)
plot3d(X5$X[1,,],type="l")
matplot(X5$t,X5$X[1,,],type="l")
```