ryamadaの遺伝学・遺伝統計学メモ このページをアンテナに追加 RSSフィード

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

2011-02-10 偏微分方程式を作る さらに

ryamada222011-02-10

[][]ローレンツアトラクタ風

  • 3次元
  • 回転は2つの中心の2次元球(円)への収束
  • 2つの円の半径は同じにしてやると、中心を結ぶベクトルと中心をめぐる回転面のなす角との関係が適当だと、いわゆるローレンツアトラクタ風の絵が描ける

http://www.genome.med.kyoto-u.ac.jp/StatGenet/testRY20110208/%E3%83%AD%E3%83%BC%E3%83%AC%E3%83%B3%E3%83%84%E3%82%A2%E3%83%88%E3%83%A9%E3%82%AF%E3%82%BF%E9%A2%A8.png

> p
           [,1]       [,2]     [,3]
[1,]  0.1840166 -0.2398052 0.684940
[2,] -0.3749743  1.4440024 1.824822
> Rs
[1] 0.9714508 0.9714508
> dimRots
[1] 2 2
> Rots
[[1]]
          [,1]       [,2] [,3]
[1,] 0.1328534 -0.9911357    0
[2,] 0.9911357  0.1328534    0
[3,] 0.0000000  0.0000000    1

[[2]]
          [,1]       [,2] [,3]
[1,] 0.1328534 -0.9911357    0
[2,] 0.9911357  0.1328534    0
[3,] 0.0000000  0.0000000    1

> RotsSub
[[1]]
           [,1]       [,2]       [,3]
[1,] -0.4272632  0.7044868 -0.5666961
[2,] -0.2547384 -0.6951955 -0.6721692
[3,] -0.8674990 -0.1428339  0.4764914

[[2]]
           [,1]       [,2]       [,3]
[1,] 0.46354856 -0.3160336  0.8277956
[2,] 0.02531692 -0.9291258 -0.3688962
[3,] 0.88570976  0.1919585 -0.4226939

> 
library(rgl)
 NormalBase<-function(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
}

n<-3
np<-2

p<-matrix(rnorm(np*n),np,n)

Niter<-1000
dt<-0.01

# 回転球の半径が、点間距離の1/2よりみじかくなるようにする
# そうすると、球がつながった形に収束する
dists<-c()
for(i in 1:(np-1)){
	for(j in (i+1):np){
		dists<-c(dists,sqrt((p[i,]-p[j,])^2))
	}
}


Rs<-runif(np,min=0,max=max(dists)/2)
Rs<-runif(np,min=max(dists)/2,max=max(dists)*1)
Rs<-runif(np)
Rs[2]<-Rs[1]

# 回転は最大でn次元
dimRots<-sample(2:n,np,replace=TRUE)
#dimRots<-rep(2,np)
# 亜次元回転の方向を決める行列
RotsSub<-NULL
RotsSubInv<-NULL
Rots<-NULL
e.outs<-NULL
Rdts<-NULL

for(i in 1:np){
	RotsSub[[i]]<-NormalBase(n)
	RotsSubInv[[i]]<-solve(RotsSub[[i]])
	Rots[[i]]<-diag(rep(1,n))
	Rots[[i]][1:dimRots[i],1:dimRots[i]]<-NormalBase(dimRots[i])
	if(i==2){
		Rots[[i]]<-Rots[[1]]
	}
	e.outs[[i]]<-eigen(Rots[[i]])
	Rdts[[i]]<-(e.outs[[i]][[2]])%*%diag((e.outs[[i]][[1]])^dt)%*%solve(e.outs[[i]][[2]])
}
Nrep<-20
xssum<-NULL
col<-c()
xssum<-p 
# 集中点も描かせる
# 回転球の中心は「赤(col=2)」
col<-rep(2,np)
#k<--runif(1)*10
k2<-3
k3<-0.1 # 亜球への収束とそれ以外の次元の収束の比
fracattraction<-1
for(rep in 1:Nrep){
	xs<-matrix(0,Niter,n)
	xs[1,]<-runif(n)*10
	#xs[1,]<-apply(p,2,sum)/np + rnorm(n)*0.001
	#xs[1,]<-p[sample(1:np,1),]+rnorm(n)*0.001
	#xs[1,]<-xs[1,]/sqrt(sum(xs[1,]^2))*0.1
	for(i in 2:Niter){
		v<-rep(0,n)
		vs<-matrix(0,np,n)
		vsRot<-vs
		for(j in 1:np){
			#vs[j,]<--(xs[i-1,]-p[j,])
			# 回転亜球の回転軸に関して回して
			# そののちに回転亜球の回転をし
			# 軸に関して戻し回転をした後
			# オリジナルの点からの移動分を求める
			tmp<-xs[i-1,]-p[j,]
			tmp<-RotsSub[[j]]%*%tmp
			tmpvs<--tmp
			tmpL<-sqrt(sum(tmpvs[1:dimRots[j]]^2))
			tmpvs[1:dimRots[j]]<-(tmpL-Rs[j])/tmpL*tmpvs[1:dimRots[j]]*k3
			vs[j,]<-RotsSubInv[[j]]%*%tmpvs
			vsRot[j,]<-Re((RotsSubInv[[j]])%*%(Rdts[[j]]%*%(tmp)))-(xs[i-1,]-p[j,])			
		}
		# 各中心までの距離
		tmpl<-sqrt(apply(vs^2,1,sum))
		# 各中心までの距離が「与えられた半径」とどれくらい違うか
		#tmpl<-tmpl-Rs
		# 各球の表面からへと近づくが
		# 遠いと速く、近いとゆっくりと近づき、離れない
		stvs<-matrix(0,np,n)
		if(prod(tmpl)!=0){
			stvs<-sign(tmpl)*vs/tmpl^k2
		}	
		
		tmpv<-apply(stvs,2,sum)
		tmpv<-tmpv/sqrt(sum(tmpv^2))

		v<-fracattraction*tmpv*(cumprod(tmpl)[length(tmpl)])*dt
		# 球表面に近いとその球面の回転成分が大きくなる
		vsRot<-vsRot/tmpl
		tmpvRot<-apply(vsRot,2,sum)
		v<-v+tmpvRot
#v<-tmpv*(cumprod(tmpl)[length(tmpl)])^k
		xs[i,]<-xs[i-1,]+v
	}
	
	xssum<-rbind(xssum,xs)
	col<-c(col,rep(rep,Niter))
}
cex<-rep(0.1,length(col))
cex[1:Nrep]<-3
plot3d(xssum[,1],xssum[,2],xssum[,3],col=col)

filename<-"attractor"
M <- par3d("userMatrix")
play3d( par3dinterp( userMatrix=list(M,
                                     rotate3d(M, pi/2, 1, 0, 0),
                                     rotate3d(M, pi/2, 0, 1, 0) ) ), 
        duration=4 )