曲線の次元

• 多変量の時系列データは、多次元空間に曲線を描く
• 観測データ(が十分に精度がよい、と、とりあえず仮定して)から
• 速度ベクトルをだし
• そのあとは、観測各時点における動標構Moving frameベクトルを差分から算出し
• また、その時点における1〜n-1次曲率を算出する
• いくつかのパターンで行う
  • 回転ベクトルによるらせん
  • トポロジーを作る曲線
  • ロトカ-ヴォルテラ(的な)曲線
• らせんもトポロジーを描く曲線も、曲率は一定
• ロトカ-ヴォルテラは一定ではないが、あるなめらかな動きをしている上位少数の曲率と、ごちゃごちゃの下位の曲率(乱雑なだけ)とに分けられる
• 曲線の中の規則性のある部分は、曲率として微分可能なものであるとすると、微分可能様に引き出された曲率の数が、この曲線の「有意味な次元」と言える(？？)

# 回転行列を作る
NormalBase<-function(n){ # n次元
	I<-X<-diag(rep(1,n))
# ペアワイズに適当な角度を定める
	thetas<-runif(n*(n-1)/2)*2*pi
	T<-matrix(0,n,n)
	T[lower.tri(T)]<-thetas
# ペアワイズに適当に回転させることを全ペアについて実施
	for(i in 1:(n-1)){
		for(j in (i+1):n){
			R<-I
			R[i,i]<-R[j,j]<-cos(T[j,i])
			R[i,j]<-sin(T[j,i])
			R[j,i]<--R[i,j]
			X<-R%*%X
		}
	}
	X
}

# Moving frameを算出するための関数


# 直交単位ベクトルを取り出す
PerpendicVector<-function(v1,v2){
	r1<-sqrt(sum(v1^2))
	r2<-sqrt(sum(v2^2))
	ret<-rep(0,length(v2))
	ct<-0
	if(r1*r2!=0){
		ip<-sum(v1*v2)
		ct<-ip/(r1*r2)
		ret<-v2/r2-v1/r1*ct
		ret<-ret/sqrt(sum(ret^2))
	}
	list(vector=ret,cos=ct)
}
# Moving frameを作る、曲率も出す
MovingFrame<-function(xs,n=NULL){
	if(is.null(n)){
		n<-length(xs[1,])
	}
	Es<-array(0,c(n,length(xs[,1]),length(xs[1,])))
	Cs<-matrix(0,n,length(xs[,1]))
	Es[1,,]<-(xs)/sqrt(apply(xs^2,1,sum))
	current<-xs
	for(i in 2:n){
		tmp<-apply(current,2,diff)
		current<-tmp
		for(j in 1:length(tmp[,1])){
			tmpout<-PerpendicVector(Es[i-1,j,],tmp[j,])
			Es[i,j,]<-tmpout$vector
			Cs[i-1,j]<-tmpout$cos
		}
	}
	list(Es=Es,Cs=Cs)
}

# トーラスのRot
n<-4
Rot<-matrix(
c(0.003640657,0.02308859,0.05751062,-0.99807124,
0.542523942,-0.63809206,-0.54456963,-0.04416122,
 -0.062159929,0.61777600,-0.78327751,-0.03106943,
 0.837729468,0.45897523,0.29430051,0.03063146),4,4,byrow=TRUE)


# シミュレーション

E.out<-eigen(Rot)

V<-E.out[[2]]
V2<-solve(V)
S<-diag(E.out[[1]])

Niter<-10000
xs<-matrix(0,Niter,n)
xs[1,]<-runif(n)

dt<-0.005

T<-runif(Niter-1)
T<-sort(T)

dhk<-NormalBase(n)
v<-dhk[1,]
for(i in 2:Niter){
	#xs[i,]<-xs[i-1,]+v*dt*T[i-1]
	xs[i,]<-xs[i-1,]+v*dt
	
	#dhk<-(V%*%S^(dt)%*%V2) %*% dhk
	dhk<-(V%*%diag(E.out[[1]]^(dt))%*%V2) %*% dhk
	# 単位ベクトル補正
	#ll<-apply(dhk^2,1,sum)
	#dhk<-((dhk)/sqrt(ll))
	v<-dhk[1,]
}



library(rgl)
plot3d(xs[,1],xs[,2],xs[,3],col=rainbow(Niter))
library(gtools)
co<-combinations(n,3)
for(i in 1:length(co[,1])){
	open3d()
	plot3d(xs[,co[i,1]],xs[,co[i,2]],xs[,co[i,3]],col=rainbow(Niter))

}

matplot(xs,type="l")


xs<-Re(xs)


vs<-apply(xs,2,diff)

Es<-MovingFrame(vs)
matplot(t(Es[[2]]),type="l")

plot(as.data.frame(t(Es[[2]])))