2011-10-02

フェノタイプの整理メモ１６　確率的に考える１２

ネットワークの因子たち(推移行列編)
$\frac{d}{dt}\mathbf{x} = \mathbf{M} \mathbf{x}$ と表されるような要素数n( $\mathbf{x}$ の長さ)の微分方程式の関係を考える
ここで、 $n \times n$ 行列 $\mathbf{M}$ を正規直交基底が作るそれとすれば、このような微分方程式において $\mathbf{x}$ は原点を中心としたn次元球面上を運動する
また、 $\mathbf{M}$ の行列式が1である(0ではない)から、 $\mathbf{Mx}=\mathbf{0}$ を満足するような $\mathbf{x}$ は存在する。すなわち、そのような $\mathbf{x}$ を不動点として持つ。不動点は $\mathbf{x},-\mathbf{x}$ の２点である
任意の点は、この微分方程式に従って、この不動点に収束していく。その様子が以下の図であり、それを実行するソースである

library(GPArotation)

# n次元空間で
# 原点を中心として、点Pが
# ベクトルVを接ベクトルとするような回転を表す行列

# デカルト座標を各座標に直す関数を使う
# デカルト座標から角座標へ
Cartesian2AngularX<-function(x){
	n<-length(x)
	if(n==1){
		return(list(r=abs(x),t=acos(sign(x))))
	}else if(n==2){
		tmpt<-0
		if(x[1]!=0){
			tmpt<-acos(x[1]/sqrt(sum(x^2)))*sign(x[2])
			if(sign(x[2])==0){
				if(x[1]>0){
					tmpt<-0
				}else if(x[1]<0){
					tmpt<-pi
				}
			}
			#tmpt<-atan(x[2]/x[1])
		}else{
			if(x[2]>0){
				tmpt<-pi/2
			}else if(x[2]<0){
				tmpt<--pi/2
			}
		}
		return(list(r=sqrt(sum(x^2)),t=tmpt))
	}
	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)
}

# (1,0,...,0)ベクトルを与えたベクトルに移す回転行列(もしくはその逆行列)を返す

VectorRotation<-function(P,inv=FALSE){ # inv=TRUEは逆行列
	n<-length(P)
	#print(n)
	if(n==1){
		return(matrix(c(1),1,1))
	}
	Q<-Cartesian2AngularX(P)

	M<-diag(n)
	M2<-M

	for(i in 1:length(Q[[2]])){
		tmp<-M
		tmp[i,i]<--sin(Q[[2]][i])
		tmp[i,i+1]<-cos(Q[[2]][i])
		tmp[i+1,i]<-cos(Q[[2]][i])
		tmp[i+1,i+1]<-sin(Q[[2]][i])
		M2<-tmp%*%M2
	}
	M2<-t(M2)[,n:1]
	
	S<-sign(M2[,1]*P)
	M2<-M2*S
	if(inv){
		M2<-solve(M2)
	}
	M2
}

# ２つのベクトルを取って
# その第１ベクトルを(1,0,...,0)に写し
# 第２ベクトルを(1,0,...,0)と(0,1,0,...,0)とが張る面に移す
# P,Vを第１・２軸平面に移す行列はinv=TRUE
# その逆はinv=FALSE
XYRotation<-function(P,V,inv=TRUE){
	# Pを(1,0,...,0)軸に回転する行列
	R1<-VectorRotation(P,inv=TRUE)
	# VそR1で回転して
	V2<-R1%*%V
	# その第１軸以外の成分を取り出して、この部分の回転を求める
	tmpV<-V2[2:length(V)]
	tmpR2<-VectorRotation(tmpV,inv=TRUE)
	R2<-diag(length(V))
	R2[2:n,2:n]<-tmpR2
	ret<-R2%*%R1
	if(!inv)ret<-solve(ret)
	ret
}

# 次元
n<-3
# 正規直交基底行列を作る
M<-Random.Start(n)
# 繰返し数(始点を変えて何度か、実施)
Nrep<-100
# 一つの始点から、追跡する回数
Niter<-1000

# プロットに揺らぎを入れるための条件
commondrift<-0.00
raredrift<-0.0
rarefreq<-0.005

# 揺らぎ部分のプロットを省くための変数
plotstart<-1


# プロットの角速度幅
dt<-0.01
# シミュレーションを行うために、x1=x2平面での回転をするが、そのための行列
S<-diag(n)
S1<-matrix(c(cos(dt),sin(dt),-sin(dt),cos(dt)),2,2)
S[1:2,1:2]<-S1


# シミュレーション記録
# 座標の記録とプロット用の色の記録
Xsum<-NULL
tmpcol<-c()


for(ir in 1:Nrep){# 始点ごとのループ

	# この始点での座標の記録
	X<-matrix(0,Niter,n)

	# 初期座標
	X[1,]<-runif(n)
	X[1,]<-X[1,]/sqrt(sum(X[1,]^2))
	# シミュレーション
	for(i in 2:Niter){
		P<-X[i-1,]
		Q<-M%*%P
		PQ<-Q-P
		# PとQを通る平面をx1=x2平面に移す行列
		R<-XYRotation(P,PQ)
		tmpP<-R%*%P
		tmpQ<-S%*%tmpP
		newQ<-solve(R)%*%tmpQ
		
		currentdrift<-commondrift
		if(runif(1)<rarefreq)currentdrift<-raredrift
		X[i,]<-newQ+rnorm(n)*currentdrift
		X[i,]<-X[i,]/sqrt(sum(X[i,]^2))
	}
	Xsum<-rbind(Xsum,X[plotstart:length(X[,1]),])
	tmpcol<-c(tmpcol,rep(ir,(Niter-plotstart+1)))
}

normX<-sqrt(sum(Xsum[1,]^2))

NrandPt<-10000
XX<-matrix(rnorm(NrandPt*n),NrandPt,n)
XX<-XX/sqrt(apply(XX^2,1,sum))*normX

XXX<-rbind(XX,Xsum)
col<-c(rep("gray",length(XX[,1])),tmpcol)
plot3d(XXX,col=col)
#open3d()
lims<-range(Xsum)
plot3d(Xsum,col=tmpcol,xlim=lims,ylim=lims,zlim=lims)