2012-01-25

複数要因の近さの関係と立体角

R 立体角相関

関係するのはこちらの話しやこちらの話し
- 三角不等式の一般次元版や球面版に関してはこちら
２次元での角度について
- ２次元平面に原点を中心とする単位円を置き、単位円周上に２点を取る
- ２点を結ぶ弧の長さは、円の中心と円周上の点を結んでできる２つの半径がなす角 $t$ を用いて、 $t$ と表される
- 円周の１周分の「弧」は $t_{\Omega}=2 \pi$ とされる
- 円周全体に対応する角度を１としたときの「角度の大きさ」は $s=\frac{t}{t_{\Omega}}$
- のときは直角に対応する
  - これを $s_{\perp}$ と書くことにする
一般次元での角度について
- k次元空間に原点を中心とする単位球を置き、単位球面上にk点を取る
  - 話を簡単にするために、k個の点が作るk個の位置ベクトルは線形独立であるとする
- k点が単位球面上に切り取る部分球面の「(長さ〜面積〜)体積」を $t^{(k)}$ とする
- 球全体の表面積は $t_{\Omega}^{(k)}=\frac{2\pi^{\frac{k}{2}}}{\Gamma(\frac{k}{2})}$
- $k=2$ 次元のときのように、球面全体に対応する相対値 $s^{(k)}=\frac{t^{(k)}}{t_{\Omega}^{(k)}}$ を定める
- このとき、 $t^{(k)}$ をk本の単位ベクトルの「角度」とし $s^{(k)}$ をその相対値と考える
- k本の単位ベクトルがすべて相互に直交しているとき、 $s^{(k)}_{\perp}=\frac{1}{2^k}$ となる
相互に独立とは限らないk個の要素の相互関係を定量することにしよう
- ペアワイズな「相関」を見るとすれば $\frac{k(k-1)}{2}$ ペアについて「相関」の程度を量で表すことになる
- ペアワイズに限らないとすれば要素数k個の集合のべき集合の要素数 $2^k$ について考慮することになる
- べき集合のうち、空集合と要素数１個の集合は「関係」に関する情報を持っているわけではないとして除外することにすれば、 $2^k-(1+k)$ 個の「関係の強さの量」を問題にすればよいことになる
- $k$ 次元空間の $k$ 個の単位ベクトルによって、要素を表すことができれば、上記の考え方によって、べき集合の要素ごとに、「多次元弧」の「長さ〜面積〜体積」を知ればよいことになる
- この計算方法が知られているのかいないのか、知らないけれど(球面三角形(こちらとか、こちらとか、の面積の一般化のはず(こちらのHyperspherical volume elementとか？？？))、モンテカルロ的に知ることはできそうだ

library(sets)
library(GPArotation)
library(MCMCpack)

# ベクトルが張る亜空間への垂線の足とそれをベクトルの線形結合で表したときの係数を返す
PerpFoot<-function(V,A,P=NULL){
	if(is.null(P))P<-rep(0,length(A))
	
	x<-solve(V%*%t(V))%*%V%*%(A-P)
	B<-P+t(V)%*%x
	return(list(X=B,k=x))
}


SolidAngle<-function(M,Niter=10000){
	k<-length(M[1,])
	# 要素を集合扱いにして、べき集合の要素を表す行列を作る
	s<-as.set(1:k)
	# べき集合
	S<-set_power(s)
	# べき集合の要素を表す行列
	Smat<-matrix(0,2^k,k)
	cnt<-1
	for(i in S){
		tmp<-as.integer(i)
		Smat[cnt,tmp]<-1
		cnt<-cnt+1
	}
	# べき集合の要素である部分集合の要素数を表すベクトル
	sumSmat<-apply(Smat,1,sum)

	# 乱点発生
	R.sp<-RandomSphere(df=k,r=1,n=Niter)
	# InsidesはInsidesの割合を
	# Insides2は「直角」に照らした比を格納
	#Minv<-solve(t(M))
	#ts<-Minv%*%t(R.sp)
	#ts.sign<-sign(ts)
	Insides<-rep(0,2^k)
	# べき集合の要素ごとに乱点の垂線の足を求め
	# 垂線の足が、要素ベクトルの正の線形和であるかどうかを確認
	#for(i in 1:(2^k)){
	for(i in (1+k+1):(2^k)){
		foot<-PerpFoot(M[which(Smat[i,]==1),],t(R.sp))
		#print(foot)
		tmp<-(foot$k)>=0
		#tmp<-(ts.sign*Smat[i,])>=0
		tmp2<-apply(tmp,2,prod)
		Insides[i]<-sum(tmp2)/length(tmp2)
	}
	Insides2<-Insides*(2^sumSmat)
	return(list(Insides=Insides,InsidesRatio=Insides2))

}

# 次元
k<-5

# k個のベクトルをランダムに作る
# 相互に直交したkベクトルを作る場合
M<-Random.Start(k)
# それなりに偏ったkベクトルを作る場合(行ベクトルが要素に対応)
#M<-as.matrix(rdirichlet(k,rep(1,k)))
#M<-M/sqrt(apply(M^2,1,sum))
#M<-matrix(c(1,0,0,-1/2,sqrt(3)/2,0,0,1/sqrt(2),1/sqrt(2)),byrow=TRUE,ncol=3)

out<-SolidAngle(M,Niter=10000)
par(mfcol=c(1,2))
plot(out$Insides)

plot(out$InsidesRatio)
par(mfcol=c(1,1))

out$InsidesRatio[(k+2):(k+1+k*(k-1)/2)]
# ペアワイズに限って言えば、相関係数と
# InsidesRatioをpi/2倍したもののコサインとが一致する
M%*%t(M)
cos(out$InsidesRatio[(k+2):(k+1+k*(k-1)/2)]*pi/2)

とはいえ、k次元空間のk本の単位ベクトルはペアワイズな角度の情報 $\frac{k(k-1)}{2}$ だけを使って、再構成できるから、自由度は $\frac{k(k-1)}{2}$ のはずで、ペアワイズな角度のみから、すべての、より高次な「面積〜体積」は求められるはずなのだが…

#次元を決めて
k<-4
# 適当に単位ベクトルの束である正方行列を作る
R<-matrix(rnorm(k^2),k,k)
R<-R/sqrt(apply(R^2,1,sum))
# その内積
T<-R%*%t(R)
# 次の手順で、この内積を守らせたベクトルセットを作ろう
# 一意に作れる
X<-matrix(0,k,k)
X[1,1]<-1
for(i in 2:k){
	tmp<-X[1:(i-1),1:(i-1)]
	tmp2<-solve(tmp)
	y<-tmp2%*%T[i,1:(i-1)]
	y2<-sqrt(1-sum(y^2))
	X[i,1:i]<-c(y,y2)
}

X