2014-09-06

曲線状のデータ

こちらの記事とそれに至る１か月ほどの記事群で、しばらく、曲線や曲面と言った、「幾何」的な扱いをいじってみた
さて、それを「観測データの解釈」に持ってくるとはどういうことか、を考えることにする
曲線を考えるためには「微分可能」な「つながったもの」がd>=2次元で得られる必要がある。このdには時間の１次元が含まれることもある
曲線がd次元で得られたとしよう
- 一つの連続量観測項目を経時的に測定すれば、空間と時間とを合わせてd=2次元の曲線であると考えることとする
- 二つの連続量観測項目を経時的に測定すればd=3次元の曲線であると考えることとする
観測とは、点の測定であるという立場に立つ(というかそういうことなのだろうけれど)と、d次元空間データ点に順序をつけて離散点を得ることが観測である
全部でk回の観察を、順序を持ってとったとする(時系列データの場合は、kの順序を時刻が決めている
取られたデータは、 $k \times d$ 個の座標の値がある。そのほかの情報はk個の点に順序があること、それらを結ぶと『曲線になるだろう』とみなすこと(→(滑らかに)結ぶ、というのは、そのように予想するというモデルを当てはめることになる)
以下の話は、kが小さいと差分幅の大きさを無視できなくなって面倒くさいので、ひとまずkは十分大きく、観測点が作る折れ線は曲線らしく見えるくらいであるとする(そうでないときは、そうでないことに対する対応策をとればよい)
曲線を考えるときは、「曲線上の点」の曲率・捩率を考えればよい。その考え方で、どのような情報が得られていることになるのかを考えてみる( $k \times d$ の座標値とkの順序、という情報が、曲線扱いをするとどうなるか、を考えてみるということ)
まず、「観測点」自体はd個の値でできているから、 $k \times d$ 個の値は情報としてある
それ以外の情報はすべて、k個の点に順序があって、曲線になっている、ということから生じる情報である
曲線を「曲線よりも高い次元空間の中に収められてオブジェクト」と考えるのが「座標」の考え方であり、「曲線のみを考えて、それが収められた空間と無関係にとらえる」のが「曲線論」なので、その「曲線だけを用いて得られる情報」が、「それ以外の情報」である
「曲線」には、すべての「点」にMoving Frame(動標構)がある。これは、進行方向ベクトル、進行方向ベクトルの微分のうち、進行方向ベクトルによらない成分の方向のベクトル、さらに２階の微分について、直交成分を…とやって、d-1階微分をしてd-1個の単位ベクトルをとり、最後にそのd-1ベクトルと併せて正規直交基底をなすように１つの単位ベクトルを決める。このようにして曲線上の点を原点とする正規直交基底が出来上がるが、この基底は「この点における、曲線の特徴」を定めている。これがMoving Frame。ちなみに、この $d^2$ 行列は最初の単位ベクトルを決めるときにd-1の自由度、次のを決めるときは長さが１という制約と先に決めたベクトルと直交という制約があるので自由度d-2、、、、ということで結局 $\frac{d(d-1)}{2}$ の自由度がある
さらに、それぞれの点の「Moving Frame」だけだと、どれくらい素早く曲がっているのかがわからない。その方向変化の速さを表すのが、スカラーで与えられる。それが「曲率・捩率・さらに高階の曲率」。これはd-1個とれる。ちなみに(一般化した)曲率は、曲線がどれだけ伸びる間にどれだけ曲がるか、という大きさを表す
結局、曲線の内在的な情報は $\frac{d(d-1)}{2} + d-1$ の自由度を観測点ごとに持ち、観測点は、曲線の収められた空間において $d$ の自由度を持つから、すべての点の持つ情報は $k\times (\frac{d(d-1)}{2} + d-1 + d) = k \times \frac{(d+1)(d+2)-4}{2}$

ds < -1:10
(d^2+3*d-2)/2

> (d^2+3*d-2)/2
 [1]  1  4  8 13 19 26 34 43 53 64

次元１の観測曲線の自由度が $k*1$ というのは、どんな観測がなされようとも、すべて１次元空間にある直線にその、どこに「観測点」が乗っているかしか問題ではなく、k個の観測座標を決めることが曲線を規定するすべてである、ということに相当している
関連してぐちゃぐちゃ書いたソース

# X 
my.Moving.Frame <- function(X,k=length(X[1,])){
	n <- length(X[,1]) # No. points
	d <- length(X[1,]) # Dimension
	# inter-point distance
	L <- rep(0,n-1)
	for(i in 1:(n-1)){
		L[i] <- sqrt(sum((X[i,]-X[i+1,])^2))
	}
	diff.X <- list()
	diff.X[[1]] <- X
	diff.L <- list()
	diff.L[[1]] <- L
	for(i in 2:k){
		diff.X[[i]] <- apply(diff.X[[i-1]],2,diff)
		tmp.n.pt <- length(diff.X[[i]][,1])
		tmp.L <- rep(0,tmp.n.pt-1)
		for(j in 1:(tmp.n.pt-1)){
			tmp.L[j] <- sqrt(sum((diff.X[[i]][j,]-diff.X[[i]][j+1,])^2))
		}
		diff.L[[i]] <- tmp.L
	}
	mf <- list() # Moving Frame
	mf[[2]] <- diff.X[[2]]/sqrt(apply(diff.X[[2]]^2,1,sum))
	if(k>2){
		for(i in 3:k){
			n.pt <- length(diff.X[[i]][,1])
			ip <- apply(diff.X[[i]] * mf[[i-1]][1:n.pt,],1,sum)
			mf[[i]] <- diff.X[[i]] - ip * mf[[i-1]][1:n.pt,]
			mf[[i]] <- mf[[i]]/sqrt(apply(mf[[i]]^2,1,sum))
		}
	}
	n.pt <- length(diff.X[[k]][,1])
	mf[[1]] <- matrix(0,n.pt,d)
	for(i in 1:n.pt){
		tmp.mat <- matrix(0,d-1,d-1)
		tmp.vec <- rep(0,d-1)
		for(j in 2:d){
			tmp <- mf[[j]][i,]
			tmp.mat[j-1,] <- tmp[1:(d-1)]
			tmp.vec[j-1] <- tmp[d]
		}
		mf[[1]][i,] <- c(-solve(tmp.mat) %*% tmp.vec,1)
	}
	mf[[1]] <- mf[[1]]/sqrt(apply(mf[[1]]^2,1,sum))
	
	Cs <- matrix(0,d-1,length(diff.L[[1]]))
	for(i in 1:(d-1)){
		n.pt <- length(diff.L[[i+1]])
		Cs[i,1:n.pt] <- diff.L[[i+1]][1:n.pt]/L[1:n.pt]
	}
	return(list(L=L,mf=mf,diff.X=diff.X,diff.L=diff.L,Cs=Cs))
}

library(rgl)

# １本の対数らせんをもう一度作る
t <- seq(from=0,to=1,length=100)*3
krA <- -3
kzA <- 0.25
thetaA <- 3.5
X <- exp(krA*t) * cos(thetaA*t)
Y <- exp(krA*t) * sin(thetaA*t)
Z <- exp(kzA*t)
#Z <- sin(t*10)
XYZ <- cbind(X,Y,Z)




out.mmf <- my.Moving.Frame(XYZ)
out.mmf2 <- my.Moving.Frame(XYZ[,1:2])
head(out.mmf2$Cs)
plot(out.mmf2$Cs)

XYZ.st <- XYZ
XYZ.st <- rbind(XYZ,rep(min(XYZ),3))
XYZ.st <- rbind(XYZ.st,rep(max(XYZ),3))

plot3d(XYZ.st)

s <- sample(1:length(out.mmf$mf[[length(XYZ[1,])]][,1]),length(out.mmf$mf[[length(XYZ[1,])]][,1])*0.5)
for(i in 1:length(out.mmf$mf)){
	#for(j in 1:length(out.mmf$mf[[i]][,1])){
	for(j in s){
#if(runif(1)<0.1){
lines3d(rbind(XYZ[j,],XYZ[j,] + out.mmf$mf[[i]][j,]*10),col=i+1)
#lines3d(rbind(XYZ[s,],XYZ[s,] + E2[s,]*10),col=3)
#lines3d(rbind(XYZ[s,],XYZ[s,] + E3[s,]*10),col=4)
#}
	}
}




matplot(t(out.mmf$Cs),type="l")

plot(atan(XYZ[,1]/XYZ[,2])[-1],out.mmf$Cs[1,],type="l")