2015-02-07

酔歩＠多次元球面

単位球面上の点から、ランダムに方向を選んで、球面上の位置を変える
方向に関して、単位球が存在している空間に関するランダムな方向ベクトルを定め、その方向に指定角だけ移動することにする
かなり制約がきつい酔歩だけれど、次元が少し上がると、どこにどういう『強い制約』があるのか、パッと見はわからない。このパッと見はわからない(けれど、本当はある)制約をどうしたら見破ることができるか、どこまでは見破ることができるか、その信じるべき根拠とその強さはどう測るか…

angular.move <- function(x,y,theta){
	x. <- x/sqrt(sum(x^2))
	y. <- y/sqrt(sum(y^2))
	ctheta <- (sum(x. * y.))
	# 方向ベクトルのうち位置ベクトルと垂直な成分
	tmp <- y. - x.*ctheta
	# 垂直ベクトルを単位ベクトル化
	tmp <- tmp/sqrt(sum(tmp^2))
	# 一歩の角から、一歩先の単位球面点を通るベクトルを位置ベクトルと方向ベクトルの線形和に表すための係数を算出
	# 移動空間である単位球と、位置ベクトルを中心とする単位球とに関して、位置ベクトルと、方向ベクトルとが作る２次元平面を切り出して、そこで２次元幾何計算をする
	Q <- tan(2*theta)
	a <- 1/(Q+1)
	b <- Q/(Q+1)
	# 位置ベクトルと方向単位ベクトルとの線形和
	z  <- a*x. + b*tmp
	# それを位置ベクトルを中心とする単位球面上に伸ばす
	z <- z/sqrt(sum(z^2))
	# 移動先の点と移動空間中心とを結ぶ直線上の点
	ret <- x + z
	# それを移動空間上に射影
	ret/sqrt(sum(ret^2))*sqrt(sum(x^2))

}

sphere.brownian <- function(n.step,p=3,x.init=c(1,rep(0,p-1)),phis=rep(0.01,n.step)){
	R.sp <- function(n,p){
		m <- matrix(rnorm(n*p),ncol=p)
		d <- sqrt(apply(m^2,1,sum))
		m/d
	}
	# 方向ベクトルをn.step発生させる
	S <- R.sp(n.step,p)
	# 位置座標軌道を格納するX
	X <- S
	for(i in 2:n.step){
		X[i,] <- angular.move(X[i-1,],S[i,],phis[i])
	}
	X
}

X <- sphere.brownian(n.step=10000,p=3,phis=rep(0.01,10000))

plot3d(X)

以下は旧版ソース

# p次元単位球面上のn個の乱点座標
R.sp <- function(n,p){
	m <- matrix(rnorm(n*p),ncol=p)
	d <- sqrt(apply(m^2,1,sum))
	m/d
}
p <- 5 # 次元

n.step <- 10^4 # 点の数・歩数-1

# 方向ベクトルをn.step発生させる
S <- R.sp(n.step,p)
# １歩の「角度」を定角にしたり、ばらつかせたり…
phis <- rep(0.01,n.step)
phis <- abs(rnorm(n.step,0,0.01))
# 位置座標軌道を格納するX
X <- S
for(i in 2:n.step){
	# 位置ベクトルと方向ベクトルのコサイン
	ctheta <- (sum(X[i-1,] * S[i,]))
	# 方向ベクトルのうち位置ベクトルと垂直な成分
	tmp <- S[i,] - X[i-1,]*ctheta
	# 垂直ベクトルを単位ベクトル化
	tmp <- tmp/sqrt(sum(tmp^2))
	# 一歩の角から、一歩先の単位球面点を通るベクトルを位置ベクトルと方向ベクトルの線形和に表すための係数を算出
	# 移動空間である単位球と、位置ベクトルを中心とする単位球とに関して、位置ベクトルと、方向ベクトルとが作る２次元平面を切り出して、そこで２次元幾何計算をする
	Q <- tan(2*phis[i])
	a <- 1/(Q+1)
	b <- Q/(Q+1)
	# 位置ベクトルと方向単位ベクトルとの線形和
	z  <- a*X[i-1,] + b*tmp
	# それを位置ベクトルを中心とする単位球面上に伸ばす
	z <- z/sqrt(sum(z^2))
	# 移動先の点と移動空間中心とを結ぶ直線上の点
	X[i,] <- X[i-1,] + z
	# それを移動空間上に射影
	X[i,] <- X[i,]/sqrt(sum(X[i,]^2))
}
library(rgl)
if(p>=3){
	plot3d(X[,1:3])
}else{
	plot(X)
}
plot(as.data.frame(X))

matplot(X,type="l")

angs <- rep(0,n.step-1)
X1 <- X[1:(n.step-1),]
X2 <- X[2:n.step,]
X12 <- X1*X2
angs <- apply(X12,1,sum)
range(acos(angs)) # 一歩の角の範囲