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

無限を有限に　丸を四角に

R 球 多次元球 ロトカ-ヴォルテラ 正単体
  • こちらで、多次元の無限空間を有限な角張った空間に納めることについて書いた
  • それは、そもそも、「この世の偏微分方程式」が動く世界をそのような空間にしたかったから(昨日の記事がその関係)
  • このことを復習する
    • n個の変数が\sum_{i=1}^n x_i=1を満足する(保存則)
    • n個の変数はx_i > 0もしくはx_i \ge 0(存在則)
    • このとき、n個の変数をn個のデカルト座標軸に取れば、保存則と存在則を満足する部分空間は、n-1次元多様体であるn-1正単体(こちら)
    • ここで、n-1次元の無限空間に「アトラクタ」や「軌道」「運動」をシミュレーションする(それは、無限空間に取り扱いやすい性質を持つ図形(球とか、(常)らせんとか)を用いると便利)
    • その上で、無限空間をn-1正単体に圧縮してやれば、「アトラクタ」「軌道」「運動」が『この世』としての部分空間に納められる
  • さて。
  • まず、n次元空間をn-1正単体に縮める方法を考える
    • ２段階ある
      • 無限空間を有限にする
      • 球を単体にする
    • 無限を有限に
      • 無限空間は中心からの距離が一定な球殻が無限距離まで広がっているような空間であるとする
      • 距離を有限にする
        • tangentを用いて、無限大をMaxXという有限正値に、0を0にマッピングする(以下のソース)
TangentDist<-function(R,MaxX){
	tmax<-atan(MaxX)
	tan(atan(R)/(pi/2)*tmax)
}
    • 球を単体に
      • 球とは、\sum_{i=1}^n x_i^2=R^2
      • 正単体とは、\sum_{i=1}^n |x_i|^1 = R^2
      • この球にこの正単体は、n個の点で内接する
      • この球に内接し、正単体よりも大きく、『球』のように丸味を帯びていて、正単体のように、n方向への突出傾向をもつような図形は\sum_{i=1}^n |x_i|^k = R^2; 1 < k < 2で表せる
      • Rが大きいとき(R=Rmax)に球であって、Rが小さいとき(R=Rmin)に正単体になるように空間圧縮をするとすれば、\sum_{i=1}^n |x_i|^ k = R^2; k = 1 + \frac{R-Rmin}{Rmax-Rmin}のように線形性を持ってkを決めてもよい。
      • n-1正単体では、 -\frac{1}{n-1} \le x_i \le 1となっているので、以下のようにして縮めることができる
# nc-1 正単体(nc個の頂点)の頂点ベクトルを作る
CategoryVector<-
function (nc = 3) 
{
    df <- nc - 1
    d <- df + 1
    diagval <- 1:d
    diagval <- sqrt((df + 1)/df) * sqrt((df - diagval + 1)/(df - 
        diagval + 2))
    others <- -diagval/(df - (0:(d - 1)))
    m <- matrix(rep(others, df + 1), nrow = df + 1, byrow = TRUE)
    diag(m) <- diagval
    m[upper.tri(m)] <- 0
    as.matrix(m[, 1:df])
}
# 無限遠をMaxXに縮めるときに距離Rをいくつに縮めるか
TangentDist<-function(R,MaxX){
	tmax<-atan(MaxX)
	tan(atan(R)/(pi/2)*tmax)
}
# n+1個の値のうち、一番小さい値とそのベクトルＩＤを出す
minind<-function(x){
	which(x==min(x))[1]
}
minindval<-function(x){
	list(ind=minind(x),v=min(x))
}

Sphere2Simplex<-function(x){
	if(!is.matrix(x)){
		x<-matrix(x,nr=1)
	}
	n<-length(x[1,])
	# n単体のn+1 個の頂点座標
	cv<-CategoryVector(n+1)
	# n+1個の頂点ベクトルへの射影の長さ
	X<-x%*%t(cv)
	# ベクトルxがどの頂点ベクトルに対してその射影を最短とするかを調べる
	# その頂点ベクトル以外の頂点が作る「n-1次元面」をxベクトルは通過する
	MinIndVal<-matrix(unlist(apply(X,1,minindval)),ncol=2,byrow=TRUE)
	# xのノルム
	r<-sqrt(apply(x^2,1,sum))
	# 無限遠を1/nに縮めるときにノルムrの点はどのくらいのノルムに縮めるか
	tmpr<-TangentDist(r,1)
	# 最長の縮み位置に対するtmprの比率
	tmpRatio<-tmpr/(1/(n))
	# 1乗から2乗の間の何乗にするか
	tmpk<-1+tmpRatio
	# 最終的に何倍に縮めるか
	ratio<-(tmpr^2/apply(abs(x)^tmpk,1,sum))^(1/tmpk)
	x*ratio
}

xs<-matrix(0,100,3)
xs[1,]<-c(0,0,1)
for(i in 2:100){
	xs[i,]<-xs[i-1,]*1.1
}
xsc<-Sphere2Simplex(xs)

plot(xsc[,1])
  • これを用いて、n次元正単体空間を動くn+1変数の軌道をn次元空間に描かせよう

http://www.genome.med.kyoto-u.ac.jp/StatGenet/testRY20110208/%E6%AD%A3%E5%8D%98%E4%BD%93%E7%A9%BA%E9%96%93%E8%BB%8C%E9%81%93.png

# nc-1 正単体(nc個の頂点)の頂点ベクトルを作る
CategoryVector<-
function (nc = 3) 
{
    df <- nc - 1
    d <- df + 1
    diagval <- 1:d
    diagval <- sqrt((df + 1)/df) * sqrt((df - diagval + 1)/(df - 
        diagval + 2))
    others <- -diagval/(df - (0:(d - 1)))
    m <- matrix(rep(others, df + 1), nrow = df + 1, byrow = TRUE)
    diag(m) <- diagval
    m[upper.tri(m)] <- 0
    as.matrix(m[, 1:df])
}
# 無限遠をMaxXに縮めるときに距離Rをいくつに縮めるか
TangentDist<-function(R,MaxX){
	tmax<-atan(MaxX)
	tan(atan(R)/(pi/2)*tmax)
}
# n+1個の値のうち、一番小さい値とそのベクトルＩＤを出す
minind<-function(x){
	which(x==min(x))[1]
}
minindval<-function(x){
	list(ind=minind(x),v=min(x))
}

Sphere2Simplex<-function(x){
	if(!is.matrix(x)){
		x<-matrix(x,nr=1)
	}
	n<-length(x[1,])
	# n単体のn+1 個の頂点座標
	cv<-CategoryVector(n+1)
	# n+1個の頂点ベクトルへの射影の長さ
	X<-x%*%t(cv)
	# ベクトルxがどの頂点ベクトルに対してその射影を最短とするかを調べる
	# その頂点ベクトル以外の頂点が作る「n-1次元面」をxベクトルは通過する
	MinIndVal<-matrix(unlist(apply(X,1,minindval)),ncol=2,byrow=TRUE)
	# xのノルム
	r<-sqrt(apply(x^2,1,sum))
	# 無限遠を1/nに縮めるときにノルムrの点はどのくらいのノルムに縮めるか
	tmpr<-TangentDist(r,1)
	# 最長の縮み位置に対するtmprの比率
	tmpRatio<-tmpr/(1/(n))
	# 1乗から2乗の間の何乗にするか
	tmpk<-1+tmpRatio
	# 最終的に何倍に縮めるか
	ratio<-(tmpr^2/apply(abs(x)^tmpk,1,sum))^(1/tmpk)
	x*ratio
}

xs<-matrix(0,100,3)
xs[1,]<-c(0,0,1)
for(i in 2:100){
	xs[i,]<-xs[i-1,]*1.1
}
xsc<-Sphere2Simplex(xs)

plot(xsc[,1])


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)*5)

# 回転は最大で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])
	
	e.outs[[i]]<-eigen(Rots[[i]])
	Rdts[[i]]<-(e.outs[[i]][[2]])%*%diag((e.outs[[i]][[1]])^dt)%*%solve(e.outs[[i]][[2]])
}
Nrep<-10
xssum<-NULL
col<-c()
xssum<-p 
# 集中点も描かせる
# 回転球の中心は「赤(col=2)」
col<-rep(2,np)
#k<--runif(1)*10
k2<-3
fracattraction<-0
for(rep in 1:Nrep){
	xs<-matrix(0,Niter,n)
	xs[1,]<-rnorm(n)
	#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]]
			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<-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 )
        
        
        
        
        
        
SimplexXs<-Sphere2Simplex(xssum)
# 正単体は中心が原点なので、全体に(1/(n-1),1/(n-1),...,1/(n-1))だけ平行移動する
SimplexXs<-SimplexXs+1/(n-1)

plot3d(SimplexXs[,1],SimplexXs[,2],SimplexXs[,3],col=col)

logo

図書室

日記全体の目次
(記事一覧)

疾患遺伝子解析を学ぶためのカタログ(H18 3月改訂)の目次

 教科書

京都大学大学院医学研究科集中講義 全８限Part I(2007)
京都大学大学院医学研究科集中講義 全８限Part II(2008)
京都大学大学院医学研究科統計遺伝学講義 全８限(2009)

以下の関連サイトへのリンク道案内ページ
遺伝子多型解析Wiki
東大医科研HGCゲノム機能解析分野(旧)
理研リウマチチーム解析資料ページ


疾患遺伝子解析を学ぶためのカタログ(H18 3月改訂)
遺伝統計学のための統計学基礎 全９限
Haploviewを用いた連鎖不平衡マッピング解析実習 全９限
Arlequinを用いた集団遺伝学実習 全７限
一から始める遺伝子多型・形質関連解析
多数のGWAS/GWASメタアナリシスとその周辺解析用のアプリケーションへのリンク

複合遺伝性疾患とは

----------News----------
20万アクセス
2008/07/20
40万アクセス
2010/09/12


-------------------------