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

無限を有限、丸を四角に

多次元球 R 角座標
  • こちらで、\sum_{i=1}^n x_i=1の世界を作っている
    • 有限かつ、辺縁が直線で仕切られている
  • そこで作っている世界はぐるぐる回ったりする
    • それを表すのに、各座標は便利
  • たとえば、掲載図に示すように
  • 限定した世界、角のとがった世界にするために、２つの工夫をする
    • (1)０から無限を有限に納める
    • (2)多次元球面を平面に変える
  • (1)０から無限を有限に納めるのには、１次元の無限に伸びる直線を２次元空間の点(0,1)から眺めることで実現する。
    • 無限遠は角度が\frac{\pi}{2}で、０は角度が０
    • これを、無限の軸の座標→角度→角度を\times k /\frac{\pi}{2}することで、無限遠も含めてすべての座標をその順序変えずに有限の座標に置き換える
  • (2)単位球面上の点は\sum x_i^2=1、これを平面\sum x_i=1に直すのに、xの符号を変えずに2乗して実現する
	# xk=sin(vk)
	# xk-1=cos(vk)*sin(vk-1)
	# xk-2=cos(vk)*cos(vk-1)*sin(vk-2)
	# ...
	# x2=cos(vk)*cos(vk-1)*cos(vk-2)*...*sin(v1)
	# x1=cos(vk)*cos(vk-1)*cos(vk-2)*...*cos(v1)

# デカルト座標から角座標へ
Cartesian2Angular<-function(x){
	r<-sqrt(sum(x^2))
	xst<-x/r
	n<-length(xst)
	t<-rep(0,n-1)
	S<-C<-t
	S[n-1]<-xst[n]
	if(S[n-1]>1)S[n-1]<-1
	if(S[n-1]<(-1))S[n-1]<-(-1)
	t[n-1]<-asin(S[n-1])
	C[n-1]<-cos(t[n-1])
	cumC<-C[n-1]
	for(i in (n-2):1){
		if(cumC!=0){
			S[i]<-xst[i+1]/cumC
			if(S[i]>1)S[i]<-1
			if(S[i]<(-1))S[i]<-(-1)
			t[i]<-asin(S[i])
			C[i]<-cos(t[i])
			cumC<-cumC*C[i]
		}else{
			S[i]<-0
			t[i]<-asin(S[i])
			C[i]<-cos(t[i])
		}
		
	}
	list(r=r,t=t)
}
Cartesian2Angular(c(0,1/sqrt(2),1/sqrt(2)))

# 角座標からデカルト座標へ
Angular2Cartesian<-function(r,t){
	n<-length(t)+1
	C<-cos(t)
	S<-sin(t)
	C2<-cumprod(C[length(C):1])
	x<-c(1,S)
	x[1:(n-1)]<-x[1:(n-1)]*C2[length(C):1]
	x<-x*r
	x
}

# 検算する
x<-runif(4)
a<-Cartesian2Angular(x)

x2<-Angular2Cartesian(a$r,a$t)

x
a
x2


# 無限遠をMaxXに縮める

TangentDist<-function(R,MaxX){
	tmax<-atan(MaxX)
	tan(atan(R)/(pi/2)*tmax)
}
# 無限遠をMinXからMaxXに縮める

TangentDist2<-function(R,MinX,MaxX){
	tmax<-atan(MaxX)
	tmin<-atan(MinX)
	tcurrent<-atan(R)
	tan((tcurrent+pi/2)/(pi)*(tmax-tmin)+tmin)
}

# (0, \infty)への変換には指数関数も使える
TangentDist3<-function(R){
	exp(R)
}

# ある値を符号付きk乗にする
SignSquare<-function(x,k=2){
	sign(x)*(abs(x))^k
}


# ある角座標点を無限遠をMaxXに縮めたうえで、符号付き平方根の座標にする
# 球が単体になる

Sphere2Simplex<-function(r,t,Maxr=1){
	tmpr<-TangentDist(r,Maxr)
	tmpx<-Angular2Cartesian(tmpr,t)
	SignSquare(tmpx)
}
Sphere2Simplex2<-function(r,t,Minr=0,Maxr=1){
	#tmpr<-TangentDist(r,Maxr)
	tmpr<-TangentDist2(r,Minr,Maxr)
	tmpx<-Angular2Cartesian(tmpr,t)
	tmpx^2
}

n<-3
Np<-1000
xs<-matrix(rnorm(n*Np),Np,n)
#xs<-xs/sqrt(apply(xs^2,1,sum))
simplexs<-xs
simplexs2<-xs
simplexs3<-xs
for(i in 1:Np){
	tmp<-Cartesian2Angular(xs[i,])
	#simplexs[i,]<-Angular2Cartesian(tmp$r,tmp$t)
	simplexs[i,]<-Sphere2Simplex(tmp$r,tmp$t)
	#simplexs2[i,]<-TangentDist2(simplexs2[i,],0,1)
	simplexs2[i,]<-TangentDist3(simplexs2[i,])
	tmp2<-Cartesian2Angular(simplexs2[i,])
	simplexs3[i,]<-Sphere2Simplex(tmp2$r,tmp2$t)
}


plot3d(xs[,1],xs[,2],xs[,3],col=rainbow(Np))
open3d()
plot3d(simplexs[,1],simplexs[,2],simplexs[,3],col=rainbow(Np))
open3d()
plot3d(simplexs2[,1],simplexs2[,2],simplexs2[,3],col=rainbow(Np))
open3d()
plot3d(simplexs3[,1],simplexs3[,2],simplexs3[,3],col=rainbow(Np))
  • 正負の無限大を有限に納めた
  • 正のみの世界にするにはどうするか
    • 上記と同じ方法でtangentを使う。そのときに「角度」を-\frac{\pi}{2}から\frac{\pi}{2}にする代わりに0からにする
    • 指数関数を使うこともできる
# 無限遠をMinXからMaxXに縮める

TangentDist2<-function(R,MinX,MaxX){
	tmax<-atan(MaxX)
	tmin<-atan(MinX)
	tcurrent<-atan(R)
	tan((tcurrent+pi/2)/(pi)*(tmax-tmin)+tmin)
}
  • この処理を加えることで、以下が可能
    • 正の世界へ
    • 正の無限大を有限に抑え
    • \sum x_i=1の平面内に納める
      • tangentを使う

http://www.genome.med.kyoto-u.ac.jp/StatGenet/lectures/2010/positiveWorld2.png

      • 指数関数を使う

http://www.genome.med.kyoto-u.ac.jp/StatGenet/lectures/2010/positiveWorld3.png

  • 掲載図のソースは以下の通り
# トーラスの前に多次元球を描こう
# k次元球の座標はオイラーの多次元角座標を用いて、k-1個の角パラメタと
# １個の半径パラメタとの合わせてk個のパラメタで表せる
# 単位k次元球の場合は、半径が１なので、
# 角座標(k-1個の角)ｄ表される
# このk次元球面の直交座標系の座標を出す関数は以下の通り
# vはk-1個の角度のベクトル

sphereCoords<-function(v){
	# すべての角度のコサインとサインをまず出す
	C<-cos(v)
	S<-sin(v)
	# xk=sin(vk)
	# xk-1=cos(vk)*sin(vk-1)
	# xk-2=cos(vk)*cos(vk-1)*sin(vk-2)
	# ...
	# x2=cos(vk)*cos(vk-1)*cos(vk-2)*...*sin(v1)
	# x1=cos(vk)*cos(vk-1)*cos(vk-2)*...*cos(v1)
	# なので以下の通り
	ret<-c(1,S)

	for(i in 1:length(v)){
		for(j in i:length(v)){
			ret[i]<-ret[i]*C[j]
		}
	}
	return(ret)
}

# これを用いて、(k次元)トーラスの座標を角座標で作る
# まず、ドーナツの主要な輪について、回転角を与える
# 次に、主要な輪の点を中心として、その位置ベクトルを含む平面であって
# 輪を乗せている面と垂直な平面を考え
# そこにドーナツの実の部分の輪を描く
# 第２の角度はこの実の部分の輪における角度である
# ただし、トーラスの場合には、輪の階層ごとに異なる半径があるのに対し
# 多次元球の場合は、角座標を上げて行っても半径は同じであることに注意
# 普通のドーナツの場合はここまで
# 次元を上げるときは、作業を繰り返す

# トーラスk次

# パラメタは以下の通り
# いわゆるトーラス
# 周期の重層性の度数
k<-2 # kを増やしてもよいようにパラメタ化してあるが、ドーナツはk=2

# トーラスでは、主要な輪について処理した後、
# その位置を中心として次の輪の分の座標を足しこむが、
# そうしなければいけないというわけでもないので、
# 輪の強さをだんだん小さくするための係数を incrRatioとする
incrRatio<-1
# トーラス上の直線は、ドーナツの主要な輪をある周波数で回りつつ、ドーナツの実の部分の小さな輪もある周波数で回る
# その周波数を与えるためのパラメタがPs
Ps<-c(7,11)

# ドットを打つ時点を作る
n<-1000
t<-seq(from=0,to=1,length.out=n)*2*pi



MultDimTorus3<-function(v,k){
	C<-cos(v)
	S<-sin(v)
	ret<-rep(0,length(v)+1)
	preret<-ret
	for(i in 1:length(v)){
		if(i==1){
			ret[1]<-C[1]
			ret[2]<-S[1]
		}else{
			# 次元を上げながら座標を作り上げていくときに
			# 前段階に増えて分の座標を確認する
			# 新たな輪の分は、次のルールで足しこむ
			# 前段階までに使われていない軸には新な輪のコサインを与える
			# 残りの足しこみ分はサインだが、それは、前段階までに使われた軸に
			# 分配する必要がある
			# その分配が「前段階に増えた分の座標」に応じて行われる
			incr<-ret-preret 
			preret<-ret
			incr<-incr/sqrt(sum(incr^2))
			ret<-k*ret+C[i]*incr
			ret[i+1]<-S[i]
			
		}
	}
	return(ret)
}

	# xk=sin(vk)
	# xk-1=cos(vk)*sin(vk-1)
	# xk-2=cos(vk)*cos(vk-1)*sin(vk-2)
	# ...
	# x2=cos(vk)*cos(vk-1)*cos(vk-2)*...*sin(v1)
	# x1=cos(vk)*cos(vk-1)*cos(vk-2)*...*cos(v1)
Cartesian2Angular(c(0,1/sqrt(2),1/sqrt(2)))

Cartesian2Angular(c(0,1/sqrt(2),1/sqrt(2)))

# 各座標からデカルト座標へ
Angular2Cartesian<-function(r,t){
	n<-length(t)+1
	#tmpt<-(t+pi)%%(2*pi)-pi
	C<-abs(cos(t))
	S<-sin(t)
	C2<-cumprod(C[length(C):1])
	x<-c(1,S)
	x[1:(n-1)]<-x[1:(n-1)]*C2[length(C):1]
	x[1]<-x[1]*sign(cos(t[1]))
	x<-x*r
	x
}


# 検算する
x<-runif(4)
a<-Cartesian2Angular(x)

x2<-Angular2Cartesian(a$r,a$t)

x
a
x2


# 球面の動き
library(rgl)
n<-5
Niter<-10000
xs<-matrix(0,Niter,n)
ts<-matrix(0,Niter,n-1)
ts[1,]<-runif(n-1)
r<-1
xs[1,]<-Angular2Cartesian(r,ts[1,])
ps<-rnorm(n-1)
dt<-0.01
for(i in 2:Niter){
	ts[i,]<-ts[i-1,]+ps*dt
	xs[i,]<-Angular2Cartesian(r,ts[i,])

}
plot3d(xs[,1],xs[,2],xs[,3],cex=0.1,col=gray((1:Niter/Niter)))
matplot(xs,type="l")


# 球面の動きができたので、これを使う

n<-3
m<-3
Ctr<-runif(n)
R<-runif(1)*min(Ctr,1-Ctr)

incrRatio<-0
Ps<-sqrt(sample(2:100,n-1,replace=TRUE))
Ps<-Ps%%(2*pi)


Nrep<-5
insideoutRatio<-0.5
xssum<-NULL
col<-c()
Niter<-1000
dt<-pi/100
k<-1
dt2<-R*0.1
dt2<-0.01
for(rep in 1:Nrep){
	col<-c(col,1,rep(rep,Niter-1))
	xs<-matrix(0,Niter,n)
	rs<-rep(Niter)
	ts<-matrix(0,Niter,n-1)
	xs[1,]<-runif(n)*100
	tmpr<-runif(1)
	if(tmpr<insideoutRatio){
		tmpr<-1/tmpr
	}
	#tmpr<-1000
	xs[1,]<-Ctr+xs[1,]/sqrt(sum(xs[1,]^2))*R*tmpr
	
	tmp<-Cartesian2Angular(xs[1,]-Ctr)
	rs[1]<-tmp$r
	ts[1,]<-tmp$t
	for(i in 2:Niter){
		DifVec<-xs[i-1,]-Ctr
		DifVecLen<-sqrt(sum(DifVec^2))
		ts[i,]<-ts[i-1,]+Ps*dt
		rs[i]<-rs[i-1]-dt2*(DifVecLen-R)
		
		xs[i,]<-Angular2Cartesian(rs[i],ts[i,])+Ctr
		Dif<-xs[i,]-xs[i-1,]
		Dif<-Dif*xs[i-1,]
		#xs[i,]<-xs[i-1,]+Dif
		tmp<-Cartesian2Angular(xs[i,]-Ctr)
		#rs[i]<-tmp$r
		#ts[i,]<-tmp$t
	}
	xssum<-rbind(xssum,xs)
}
xlim<-ylim<-zlim<-c(min(xssum),max(xssum))
plot3d(xssum[,1],xssum[,2],xssum[,3],col=col,xlim=xlim,ylim=ylim,zlim=zlim)


# xssumをピラミッドに納めよう
xssum2<-xssum
for(i in 1:length(xssum[,1])){
	#tmpxssum2<-TangentDist2(xssum[i,],0,1)
	tmpxssum2<-TangentDist3(xssum[i,])
	tmp2<-Cartesian2Angular(tmpxssum2)
	xssum2[i,]<-Sphere2Simplex(tmp2$r,tmp2$t)
	
}

xlim<-ylim<-zlim<-c(min(xssum),max(xssum))
plot3d(xssum[,1],xssum[,2],xssum[,3],col=col,xlim=xlim,ylim=ylim,zlim=zlim)
open3d()
xlim<-ylim<-zlim<-c(min(xssum2),max(xssum2))
plot3d(xssum2[,1],xssum2[,2],xssum2[,3],col=col,xlim=xlim,ylim=ylim,zlim=zlim)