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

体積計算

単体 行列式 正単体 特異値分解 
SimplexVolume<-function(x,Factorial=FALSE){
	n<-length(x[,1])
	#d<-t(x[2:n,])-x[1,]
	d<-apply(x,2,FUN="diff")
	if(Factorial){
		ret<-log(abs(det(d))) - lfactorial(n-1)
	}else{
		ret<-log(abs(det(d)))
	}
	return(ret)
}

x1<-runif(2)
x2<-runif(2)
xx<-matrix(c(0,0,x1,x2),ncol=2,byrow=TRUE)

exp(SimplexVolume(xx,Factorial=TRUE))

d<-sqrt(sum(x1-x2)^2)
d<-sqrt(sum(x1^2))
h<-sqrt(sum(x2^2))
t<-acos(sum(x1*x2)/(d*h))
h2<-h*sin(t)
d*h2/2
det(matrix(c(x1,x2-x1),2,2))/2
  • 正単体の体積は\frac{\sqrt{n+1}}{n!}に比例する
    • たとえば、k次元ユークリッド空間の正規直交基底の各軸の上の点で、原点から距離１の点k個をつなぐとk-1次正単体ができるが、この正単体の１辺の長さは\sqrt{2}であるのだが、このような正単体の体積は\frac{\sqrt{k}{(k-1)!}そのものである
    • こちら
  • 高次元空間にあるk個の点は(k-1)次元空間に納まるが、そのk個の点を頂点とする(k-1)次元空間凸包の体積を求めるときには、次元を下げる必要もあって、そのために特異値分解を以下のように使えて(これについてはこちら)、それによって、高次元体積の計算ができる(ただし、以下のソースではFactorial=FALSEをデフォルトにすることで、単体の体積ではなく、単体の１点を基準にk-1個のベクトルが作る「ひし形のお化け」の体積(の対数)を計算している。Factorial=TRUEにすると、単体の体積になる。
  • このソースを使うと、単体の体積
SimplexVolume<-function(x,Factorial=FALSE){
	n<-length(x[,1])
	#d<-t(x[2:n,])-x[1,]
	d<-apply(x,2,FUN="diff")
	if(Factorial){
		ret<-log(abs(det(d))) - lfactorial(n-1)
	}else{
		ret<-log(abs(det(d)))
	}
	return(ret)
}


SimplexVolume2<-function(x,Factorial=FALSE){
	n<-length(x[,1])
	tmpX<-t(x[2:n,])-x[1,]
	if(n==2){
		ret<-log(sqrt(sum(tmpX^2)))
	}else{
		svd.out<-svd(t(tmpX))
		#print(svd.out)
		#print(tmpX)
		#print(svd.out[[3]])
		tmpX<-t(tmpX)%*%svd.out[[3]]
		#print(tmpX)

		ret<-SimplexVolume(rbind(rep(0,length(tmpX[1,])),tmpX),Factorial=Factorial)
	}
	return(ret)
}



k<-3
df<-5
Npt<-5

X<-matrix(runif(Npt*df),ncol=df)
X<-diag(rep(1,Npt))
library(gtools)

cmb<-combinations(Npt,k)
cmb
Ncmb<-length(cmb[,1])
w<-rep(0,Ncmb)

for(i in 1:Ncmb){
	# k x df matrix
	tmpX<-X[cmb[i,],]
	tmpn<-length(tmpX[,1])
	# (k-1) x df matrix
	tmpX<-t(tmpX[2:tmpn,])-tmpX[1,]

	if(tmpn==2){
		w[i]<-log(sqrt(sum(tmpX^2)))
	}else{
		svd.out<-svd(t(tmpX))
		print(svd.out)
		print(tmpX)
		print(svd.out[[3]])
		tmpX<-t(tmpX)%*%svd.out[[3]]
		print(tmpX)

		w[i]<-SimplexVolume(rbind(rep(0,length(tmpX[1,])),tmpX),Factorial=TRUE)
	}
	
}

exp(w)

k<-3
df<-2
Npt<-10

X<-matrix(runif(Npt*df),ncol=df)
#X<-diag(rep(1,Npt))
library(gtools)

cmb<-combinations(Npt,k)
cmb
Ncmb<-length(cmb[,1])
w<-rep(0,Ncmb)

for(i in 1:Ncmb){
	# k x df matrix
	tmpX<-X[cmb[i,],]
	w[i]<-SimplexVolume2(tmpX,Factorial=TRUE)
}

exp(w)
cmb2<-cmb
print(cmb2)
selected<-rep(0,Npt)
h<-cmb[which(w==min(w)),]
selected[h]<-1
ss<-which(selected==1)
for(i in 1:length(ss)){
	cmb2[which(cmb2==ss[i])]<-0
	cmb2<-t(apply(cmb2,1,sort))
}

while(sum(selected)<Npt){
	candidates<-which((cmb2[,k-1]==0) & (cmb2[,k]!=0))
	wc<-w[candidates]
	minw<-candidates[which(wc==min(wc))]
	h<-rbind(h,cmb[minw,])
	selected[cmb[minw,]]<-1
	newselect<-cmb2[minw,k]
	#Hx<-c(Hx,newselect)
	cmb2[which(cmb2==newselect)]<-0
	cmb2<-t(apply(cmb2,1,sort))
	print(cmb2)
}

print(h)


dM<-dist(X)
stree<-spantree(dM)

stree